Actual source code: dtds.c

petsc-3.7.5 2017-01-01
Report Typos and Errors
  1: #include <petsc/private/petscdsimpl.h> /*I "petscds.h" I*/

  3: PetscClassId PETSCDS_CLASSID = 0;

  5: PetscFunctionList PetscDSList              = NULL;
  6: PetscBool         PetscDSRegisterAllCalled = PETSC_FALSE;

 10: /*@C
 11:   PetscDSRegister - Adds a new PetscDS implementation

 13:   Not Collective

 15:   Input Parameters:
 16: + name        - The name of a new user-defined creation routine
 17: - create_func - The creation routine itself

 19:   Notes:
 20:   PetscDSRegister() may be called multiple times to add several user-defined PetscDSs

 22:   Sample usage:
 23: .vb
 24:     PetscDSRegister("my_ds", MyPetscDSCreate);
 25: .ve

 27:   Then, your PetscDS type can be chosen with the procedural interface via
 28: .vb
 29:     PetscDSCreate(MPI_Comm, PetscDS *);
 30:     PetscDSSetType(PetscDS, "my_ds");
 31: .ve
 32:    or at runtime via the option
 33: .vb
 34:     -petscds_type my_ds
 35: .ve

 37:   Level: advanced

 39: .keywords: PetscDS, register
 40: .seealso: PetscDSRegisterAll(), PetscDSRegisterDestroy()

 42: @*/
 43: PetscErrorCode PetscDSRegister(const char sname[], PetscErrorCode (*function)(PetscDS))
 44: {

 48:   PetscFunctionListAdd(&PetscDSList, sname, function);
 49:   return(0);
 50: }

 54: /*@C
 55:   PetscDSSetType - Builds a particular PetscDS

 57:   Collective on PetscDS

 59:   Input Parameters:
 60: + prob - The PetscDS object
 61: - name - The kind of system

 63:   Options Database Key:
 64: . -petscds_type <type> - Sets the PetscDS type; use -help for a list of available types

 66:   Level: intermediate

 68: .keywords: PetscDS, set, type
 69: .seealso: PetscDSGetType(), PetscDSCreate()
 70: @*/
 71: PetscErrorCode PetscDSSetType(PetscDS prob, PetscDSType name)
 72: {
 73:   PetscErrorCode (*r)(PetscDS);
 74:   PetscBool      match;

 79:   PetscObjectTypeCompare((PetscObject) prob, name, &match);
 80:   if (match) return(0);

 82:   PetscDSRegisterAll();
 83:   PetscFunctionListFind(PetscDSList, name, &r);
 84:   if (!r) SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscDS type: %s", name);

 86:   if (prob->ops->destroy) {
 87:     (*prob->ops->destroy)(prob);
 88:     prob->ops->destroy = NULL;
 89:   }
 90:   (*r)(prob);
 91:   PetscObjectChangeTypeName((PetscObject) prob, name);
 92:   return(0);
 93: }

 97: /*@C
 98:   PetscDSGetType - Gets the PetscDS type name (as a string) from the object.

100:   Not Collective

102:   Input Parameter:
103: . prob  - The PetscDS

105:   Output Parameter:
106: . name - The PetscDS type name

108:   Level: intermediate

110: .keywords: PetscDS, get, type, name
111: .seealso: PetscDSSetType(), PetscDSCreate()
112: @*/
113: PetscErrorCode PetscDSGetType(PetscDS prob, PetscDSType *name)
114: {

120:   PetscDSRegisterAll();
121:   *name = ((PetscObject) prob)->type_name;
122:   return(0);
123: }

127: static PetscErrorCode PetscDSView_Ascii(PetscDS prob, PetscViewer viewer)
128: {
129:   PetscViewerFormat format;
130:   PetscInt          f;
131:   PetscErrorCode    ierr;

134:   PetscViewerGetFormat(viewer, &format);
135:   PetscViewerASCIIPrintf(viewer, "Discrete System with %d fields\n", prob->Nf);
136:   PetscViewerASCIIPushTab(viewer);
137:   for (f = 0; f < prob->Nf; ++f) {
138:     PetscObject  obj;
139:     PetscClassId id;
140:     const char  *name;
141:     PetscInt     Nc;

143:     PetscDSGetDiscretization(prob, f, &obj);
144:     PetscObjectGetClassId(obj, &id);
145:     PetscObjectGetName(obj, &name);
146:     PetscViewerASCIIPrintf(viewer, "Field %s", name ? name : "<unknown>");
147:     if (id == PETSCFE_CLASSID)      {
148:       PetscFEGetNumComponents((PetscFE) obj, &Nc);
149:       PetscViewerASCIIPrintf(viewer, " FEM");
150:     } else if (id == PETSCFV_CLASSID) {
151:       PetscFVGetNumComponents((PetscFV) obj, &Nc);
152:       PetscViewerASCIIPrintf(viewer, " FVM");
153:     }
154:     else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
155:     if (Nc > 1) {PetscViewerASCIIPrintf(viewer, "%d components", Nc);}
156:     else        {PetscViewerASCIIPrintf(viewer, "%d component ", Nc);}
157:     if (prob->implicit[f]) {PetscViewerASCIIPrintf(viewer, " (implicit)");}
158:     else                   {PetscViewerASCIIPrintf(viewer, " (explicit)");}
159:     if (prob->adjacency[f*2+0]) {
160:       if (prob->adjacency[f*2+1]) {PetscViewerASCIIPrintf(viewer, " (adj FVM++)");}
161:       else                        {PetscViewerASCIIPrintf(viewer, " (adj FVM)");}
162:     } else {
163:       if (prob->adjacency[f*2+1]) {PetscViewerASCIIPrintf(viewer, " (adj FEM)");}
164:       else                        {PetscViewerASCIIPrintf(viewer, " (adj FUNKY)");}
165:     }
166:     PetscViewerASCIIPrintf(viewer, "\n");
167:     if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
168:       if (id == PETSCFE_CLASSID)      {PetscFEView((PetscFE) obj, viewer);}
169:       else if (id == PETSCFV_CLASSID) {PetscFVView((PetscFV) obj, viewer);}
170:     }
171:   }
172:   PetscViewerASCIIPopTab(viewer);
173:   return(0);
174: }

178: /*@C
179:   PetscDSView - Views a PetscDS

181:   Collective on PetscDS

183:   Input Parameter:
184: + prob - the PetscDS object to view
185: - v  - the viewer

187:   Level: developer

189: .seealso PetscDSDestroy()
190: @*/
191: PetscErrorCode PetscDSView(PetscDS prob, PetscViewer v)
192: {
193:   PetscBool      iascii;

198:   if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) prob), &v);}
200:   PetscObjectTypeCompare((PetscObject) v, PETSCVIEWERASCII, &iascii);
201:   if (iascii) {PetscDSView_Ascii(prob, v);}
202:   if (prob->ops->view) {(*prob->ops->view)(prob, v);}
203:   return(0);
204: }

208: /*@
209:   PetscDSSetFromOptions - sets parameters in a PetscDS from the options database

211:   Collective on PetscDS

213:   Input Parameter:
214: . prob - the PetscDS object to set options for

216:   Options Database:

218:   Level: developer

220: .seealso PetscDSView()
221: @*/
222: PetscErrorCode PetscDSSetFromOptions(PetscDS prob)
223: {
224:   const char    *defaultType;
225:   char           name[256];
226:   PetscBool      flg;

231:   if (!((PetscObject) prob)->type_name) {
232:     defaultType = PETSCDSBASIC;
233:   } else {
234:     defaultType = ((PetscObject) prob)->type_name;
235:   }
236:   PetscDSRegisterAll();

238:   PetscObjectOptionsBegin((PetscObject) prob);
239:   PetscOptionsFList("-petscds_type", "Discrete System", "PetscDSSetType", PetscDSList, defaultType, name, 256, &flg);
240:   if (flg) {
241:     PetscDSSetType(prob, name);
242:   } else if (!((PetscObject) prob)->type_name) {
243:     PetscDSSetType(prob, defaultType);
244:   }
245:   if (prob->ops->setfromoptions) {(*prob->ops->setfromoptions)(prob);}
246:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
247:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) prob);
248:   PetscOptionsEnd();
249:   PetscDSViewFromOptions(prob, NULL, "-petscds_view");
250:   return(0);
251: }

255: /*@C
256:   PetscDSSetUp - Construct data structures for the PetscDS

258:   Collective on PetscDS

260:   Input Parameter:
261: . prob - the PetscDS object to setup

263:   Level: developer

265: .seealso PetscDSView(), PetscDSDestroy()
266: @*/
267: PetscErrorCode PetscDSSetUp(PetscDS prob)
268: {
269:   const PetscInt Nf = prob->Nf;
270:   PetscInt       dim, work, NcMax = 0, NqMax = 0, f;

275:   if (prob->setup) return(0);
276:   /* Calculate sizes */
277:   PetscDSGetSpatialDimension(prob, &dim);
278:   prob->totDim = prob->totDimBd = prob->totComp = 0;
279:   PetscCalloc4(Nf+1,&prob->off,Nf+1,&prob->offDer,Nf+1,&prob->offBd,Nf+1,&prob->offDerBd);
280:   PetscMalloc4(Nf,&prob->basis,Nf,&prob->basisDer,Nf,&prob->basisBd,Nf,&prob->basisDerBd);
281:   for (f = 0; f < Nf; ++f) {
282:     PetscFE         feBd = (PetscFE) prob->discBd[f];
283:     PetscObject     obj;
284:     PetscClassId    id;
285:     PetscQuadrature q;
286:     PetscInt        Nq = 0, Nb, Nc;

288:     PetscDSGetDiscretization(prob, f, &obj);
289:     PetscObjectGetClassId(obj, &id);
290:     if (id == PETSCFE_CLASSID)      {
291:       PetscFE fe = (PetscFE) obj;

293:       PetscFEGetQuadrature(fe, &q);
294:       PetscFEGetDimension(fe, &Nb);
295:       PetscFEGetNumComponents(fe, &Nc);
296:       PetscFEGetDefaultTabulation(fe, &prob->basis[f], &prob->basisDer[f], NULL);
297:     } else if (id == PETSCFV_CLASSID) {
298:       PetscFV fv = (PetscFV) obj;

300:       PetscFVGetQuadrature(fv, &q);
301:       Nb   = 1;
302:       PetscFVGetNumComponents(fv, &Nc);
303:       PetscFVGetDefaultTabulation(fv, &prob->basis[f], &prob->basisDer[f], NULL);
304:     } else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
305:     prob->off[f+1]    = Nc     + prob->off[f];
306:     prob->offDer[f+1] = Nc*dim + prob->offDer[f];
307:     if (q) {PetscQuadratureGetData(q, NULL, &Nq, NULL, NULL);}
308:     NqMax          = PetscMax(NqMax, Nq);
309:     NcMax          = PetscMax(NcMax, Nc);
310:     prob->totDim  += Nb*Nc;
311:     prob->totComp += Nc;
312:     if (feBd) {
313:       PetscFEGetDimension(feBd, &Nb);
314:       PetscFEGetNumComponents(feBd, &Nc);
315:       PetscFEGetDefaultTabulation(feBd, &prob->basisBd[f], &prob->basisDerBd[f], NULL);
316:       prob->totDimBd += Nb*Nc;
317:       prob->offBd[f+1]    = Nc     + prob->offBd[f];
318:       prob->offDerBd[f+1] = Nc*dim + prob->offDerBd[f];
319:     }
320:   }
321:   work = PetscMax(prob->totComp*dim, PetscSqr(NcMax*dim));
322:   /* Allocate works space */
323:   PetscMalloc5(prob->totComp,&prob->u,prob->totComp,&prob->u_t,prob->totComp*dim,&prob->u_x,dim,&prob->x,work,&prob->refSpaceDer);
324:   PetscMalloc6(NqMax*NcMax,&prob->f0,NqMax*NcMax*dim,&prob->f1,NqMax*NcMax*NcMax,&prob->g0,NqMax*NcMax*NcMax*dim,&prob->g1,NqMax*NcMax*NcMax*dim,&prob->g2,NqMax*NcMax*NcMax*dim*dim,&prob->g3);
325:   if (prob->ops->setup) {(*prob->ops->setup)(prob);}
326:   prob->setup = PETSC_TRUE;
327:   return(0);
328: }

332: static PetscErrorCode PetscDSDestroyStructs_Static(PetscDS prob)
333: {

337:   PetscFree4(prob->off,prob->offDer,prob->offBd,prob->offDerBd);
338:   PetscFree4(prob->basis,prob->basisDer,prob->basisBd,prob->basisDerBd);
339:   PetscFree5(prob->u,prob->u_t,prob->u_x,prob->x,prob->refSpaceDer);
340:   PetscFree6(prob->f0,prob->f1,prob->g0,prob->g1,prob->g2,prob->g3);
341:   return(0);
342: }

346: static PetscErrorCode PetscDSEnlarge_Static(PetscDS prob, PetscInt NfNew)
347: {
348:   PetscObject      *tmpd, *tmpdbd;
349:   PetscBool        *tmpi, *tmpa;
350:   PetscPointFunc   *tmpobj, *tmpf;
351:   PetscPointJac    *tmpg, *tmpgp, *tmpgt;
352:   PetscBdPointFunc *tmpfbd;
353:   PetscBdPointJac  *tmpgbd;
354:   PetscRiemannFunc *tmpr;
355:   void            **tmpctx;
356:   PetscInt          Nf = prob->Nf, f, i;
357:   PetscErrorCode    ierr;

360:   if (Nf >= NfNew) return(0);
361:   prob->setup = PETSC_FALSE;
362:   PetscDSDestroyStructs_Static(prob);
363:   PetscMalloc4(NfNew, &tmpd, NfNew, &tmpdbd, NfNew, &tmpi, NfNew*2, &tmpa);
364:   for (f = 0; f < Nf; ++f) {tmpd[f] = prob->disc[f]; tmpdbd[f] = prob->discBd[f]; tmpi[f] = prob->implicit[f]; for (i = 0; i < 2; ++i) tmpa[f*2+i] = prob->adjacency[f*2+i];}
365:   for (f = Nf; f < NfNew; ++f) {PetscContainerCreate(PetscObjectComm((PetscObject) prob), (PetscContainer *) &tmpd[f]);
366:     tmpdbd[f] = NULL; tmpi[f] = PETSC_TRUE; tmpa[f*2+0] = PETSC_FALSE; tmpa[f*2+1] = PETSC_TRUE;}
367:   PetscFree4(prob->disc, prob->discBd, prob->implicit, prob->adjacency);
368:   prob->Nf        = NfNew;
369:   prob->disc      = tmpd;
370:   prob->discBd    = tmpdbd;
371:   prob->implicit  = tmpi;
372:   prob->adjacency = tmpa;
373:   PetscCalloc7(NfNew, &tmpobj, NfNew*2, &tmpf, NfNew*NfNew*4, &tmpg, NfNew*NfNew*4, &tmpgp, NfNew*NfNew*4, &tmpgt, NfNew, &tmpr, NfNew, &tmpctx);
374:   for (f = 0; f < Nf; ++f) tmpobj[f] = prob->obj[f];
375:   for (f = 0; f < Nf*2; ++f) tmpf[f] = prob->f[f];
376:   for (f = 0; f < Nf*Nf*4; ++f) tmpg[f] = prob->g[f];
377:   for (f = 0; f < Nf*Nf*4; ++f) tmpgp[f] = prob->gp[f];
378:   for (f = 0; f < Nf; ++f) tmpr[f] = prob->r[f];
379:   for (f = 0; f < Nf; ++f) tmpctx[f] = prob->ctx[f];
380:   for (f = Nf; f < NfNew; ++f) tmpobj[f] = NULL;
381:   for (f = Nf*2; f < NfNew*2; ++f) tmpf[f] = NULL;
382:   for (f = Nf*Nf*4; f < NfNew*NfNew*4; ++f) tmpg[f] = NULL;
383:   for (f = Nf*Nf*4; f < NfNew*NfNew*4; ++f) tmpgp[f] = NULL;
384:   for (f = Nf*Nf*4; f < NfNew*NfNew*4; ++f) tmpgt[f] = NULL;
385:   for (f = Nf; f < NfNew; ++f) tmpr[f] = NULL;
386:   for (f = Nf; f < NfNew; ++f) tmpctx[f] = NULL;
387:   PetscFree7(prob->obj, prob->f, prob->g, prob->gp, prob->gt, prob->r, prob->ctx);
388:   prob->obj = tmpobj;
389:   prob->f   = tmpf;
390:   prob->g   = tmpg;
391:   prob->gp  = tmpgp;
392:   prob->gt  = tmpgt;
393:   prob->r   = tmpr;
394:   prob->ctx = tmpctx;
395:   PetscCalloc2(NfNew*2, &tmpfbd, NfNew*NfNew*4, &tmpgbd);
396:   for (f = 0; f < Nf*2; ++f) tmpfbd[f] = prob->fBd[f];
397:   for (f = 0; f < Nf*Nf*4; ++f) tmpgbd[f] = prob->gBd[f];
398:   for (f = Nf*2; f < NfNew*2; ++f) tmpfbd[f] = NULL;
399:   for (f = Nf*Nf*4; f < NfNew*NfNew*4; ++f) tmpgbd[f] = NULL;
400:   PetscFree2(prob->fBd, prob->gBd);
401:   prob->fBd = tmpfbd;
402:   prob->gBd = tmpgbd;
403:   return(0);
404: }

408: /*@
409:   PetscDSDestroy - Destroys a PetscDS object

411:   Collective on PetscDS

413:   Input Parameter:
414: . prob - the PetscDS object to destroy

416:   Level: developer

418: .seealso PetscDSView()
419: @*/
420: PetscErrorCode PetscDSDestroy(PetscDS *prob)
421: {
422:   PetscInt       f;

426:   if (!*prob) return(0);

429:   if (--((PetscObject)(*prob))->refct > 0) {*prob = 0; return(0);}
430:   ((PetscObject) (*prob))->refct = 0;
431:   PetscDSDestroyStructs_Static(*prob);
432:   for (f = 0; f < (*prob)->Nf; ++f) {
433:     PetscObjectDereference((*prob)->disc[f]);
434:     PetscObjectDereference((*prob)->discBd[f]);
435:   }
436:   PetscFree4((*prob)->disc, (*prob)->discBd, (*prob)->implicit, (*prob)->adjacency);
437:   PetscFree7((*prob)->obj,(*prob)->f,(*prob)->g,(*prob)->gp,(*prob)->gt,(*prob)->r,(*prob)->ctx);
438:   PetscFree2((*prob)->fBd,(*prob)->gBd);
439:   if ((*prob)->ops->destroy) {(*(*prob)->ops->destroy)(*prob);}
440:   PetscHeaderDestroy(prob);
441:   return(0);
442: }

446: /*@
447:   PetscDSCreate - Creates an empty PetscDS object. The type can then be set with PetscDSSetType().

449:   Collective on MPI_Comm

451:   Input Parameter:
452: . comm - The communicator for the PetscDS object

454:   Output Parameter:
455: . prob - The PetscDS object

457:   Level: beginner

459: .seealso: PetscDSSetType(), PETSCDSBASIC
460: @*/
461: PetscErrorCode PetscDSCreate(MPI_Comm comm, PetscDS *prob)
462: {
463:   PetscDS   p;

468:   *prob  = NULL;
469:   PetscDSInitializePackage();

471:   PetscHeaderCreate(p, PETSCDS_CLASSID, "PetscDS", "Discrete System", "PetscDS", comm, PetscDSDestroy, PetscDSView);

473:   p->Nf    = 0;
474:   p->setup = PETSC_FALSE;

476:   *prob = p;
477:   return(0);
478: }

482: /*@
483:   PetscDSGetNumFields - Returns the number of fields in the DS

485:   Not collective

487:   Input Parameter:
488: . prob - The PetscDS object

490:   Output Parameter:
491: . Nf - The number of fields

493:   Level: beginner

495: .seealso: PetscDSGetSpatialDimension(), PetscDSCreate()
496: @*/
497: PetscErrorCode PetscDSGetNumFields(PetscDS prob, PetscInt *Nf)
498: {
502:   *Nf = prob->Nf;
503:   return(0);
504: }

508: /*@
509:   PetscDSGetSpatialDimension - Returns the spatial dimension of the DS

511:   Not collective

513:   Input Parameter:
514: . prob - The PetscDS object

516:   Output Parameter:
517: . dim - The spatial dimension

519:   Level: beginner

521: .seealso: PetscDSGetNumFields(), PetscDSCreate()
522: @*/
523: PetscErrorCode PetscDSGetSpatialDimension(PetscDS prob, PetscInt *dim)
524: {

530:   *dim = 0;
531:   if (prob->Nf) {
532:     PetscObject  obj;
533:     PetscClassId id;

535:     PetscDSGetDiscretization(prob, 0, &obj);
536:     PetscObjectGetClassId(obj, &id);
537:     if (id == PETSCFE_CLASSID)      {PetscFEGetSpatialDimension((PetscFE) obj, dim);}
538:     else if (id == PETSCFV_CLASSID) {PetscFVGetSpatialDimension((PetscFV) obj, dim);}
539:     else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", 0);
540:   }
541:   return(0);
542: }

546: /*@
547:   PetscDSGetTotalDimension - Returns the total size of the approximation space for this system

549:   Not collective

551:   Input Parameter:
552: . prob - The PetscDS object

554:   Output Parameter:
555: . dim - The total problem dimension

557:   Level: beginner

559: .seealso: PetscDSGetNumFields(), PetscDSCreate()
560: @*/
561: PetscErrorCode PetscDSGetTotalDimension(PetscDS prob, PetscInt *dim)
562: {

567:   PetscDSSetUp(prob);
569:   *dim = prob->totDim;
570:   return(0);
571: }

575: /*@
576:   PetscDSGetTotalBdDimension - Returns the total size of the boundary approximation space for this system

578:   Not collective

580:   Input Parameter:
581: . prob - The PetscDS object

583:   Output Parameter:
584: . dim - The total boundary problem dimension

586:   Level: beginner

588: .seealso: PetscDSGetNumFields(), PetscDSCreate()
589: @*/
590: PetscErrorCode PetscDSGetTotalBdDimension(PetscDS prob, PetscInt *dim)
591: {

596:   PetscDSSetUp(prob);
598:   *dim = prob->totDimBd;
599:   return(0);
600: }

604: /*@
605:   PetscDSGetTotalComponents - Returns the total number of components in this system

607:   Not collective

609:   Input Parameter:
610: . prob - The PetscDS object

612:   Output Parameter:
613: . dim - The total number of components

615:   Level: beginner

617: .seealso: PetscDSGetNumFields(), PetscDSCreate()
618: @*/
619: PetscErrorCode PetscDSGetTotalComponents(PetscDS prob, PetscInt *Nc)
620: {

625:   PetscDSSetUp(prob);
627:   *Nc = prob->totComp;
628:   return(0);
629: }

633: /*@
634:   PetscDSGetDiscretization - Returns the discretization object for the given field

636:   Not collective

638:   Input Parameters:
639: + prob - The PetscDS object
640: - f - The field number

642:   Output Parameter:
643: . disc - The discretization object

645:   Level: beginner

647: .seealso: PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetBdDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
648: @*/
649: PetscErrorCode PetscDSGetDiscretization(PetscDS prob, PetscInt f, PetscObject *disc)
650: {
654:   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
655:   *disc = prob->disc[f];
656:   return(0);
657: }

661: /*@
662:   PetscDSGetBdDiscretization - Returns the boundary discretization object for the given field

664:   Not collective

666:   Input Parameters:
667: + prob - The PetscDS object
668: - f - The field number

670:   Output Parameter:
671: . disc - The boundary discretization object

673:   Level: beginner

675: .seealso: PetscDSSetBdDiscretization(), PetscDSAddBdDiscretization(), PetscDSGetDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
676: @*/
677: PetscErrorCode PetscDSGetBdDiscretization(PetscDS prob, PetscInt f, PetscObject *disc)
678: {
682:   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
683:   *disc = prob->discBd[f];
684:   return(0);
685: }

689: /*@
690:   PetscDSSetDiscretization - Sets the discretization object for the given field

692:   Not collective

694:   Input Parameters:
695: + prob - The PetscDS object
696: . f - The field number
697: - disc - The discretization object

699:   Level: beginner

701: .seealso: PetscDSGetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
702: @*/
703: PetscErrorCode PetscDSSetDiscretization(PetscDS prob, PetscInt f, PetscObject disc)
704: {

710:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
711:   PetscDSEnlarge_Static(prob, f+1);
712:   if (prob->disc[f]) {PetscObjectDereference(prob->disc[f]);}
713:   prob->disc[f] = disc;
714:   PetscObjectReference(disc);
715:   {
716:     PetscClassId id;

718:     PetscObjectGetClassId(disc, &id);
719:     if (id == PETSCFV_CLASSID) {
720:       prob->implicit[f]      = PETSC_FALSE;
721:       prob->adjacency[f*2+0] = PETSC_TRUE;
722:       prob->adjacency[f*2+1] = PETSC_FALSE;
723:     }
724:   }
725:   return(0);
726: }

730: /*@
731:   PetscDSSetBdDiscretization - Sets the boundary discretization object for the given field

733:   Not collective

735:   Input Parameters:
736: + prob - The PetscDS object
737: . f - The field number
738: - disc - The boundary discretization object

740:   Level: beginner

742: .seealso: PetscDSGetBdDiscretization(), PetscDSAddBdDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
743: @*/
744: PetscErrorCode PetscDSSetBdDiscretization(PetscDS prob, PetscInt f, PetscObject disc)
745: {

751:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
752:   PetscDSEnlarge_Static(prob, f+1);
753:   if (prob->discBd[f]) {PetscObjectDereference(prob->discBd[f]);}
754:   prob->discBd[f] = disc;
755:   PetscObjectReference(disc);
756:   return(0);
757: }

761: /*@
762:   PetscDSAddDiscretization - Adds a discretization object

764:   Not collective

766:   Input Parameters:
767: + prob - The PetscDS object
768: - disc - The boundary discretization object

770:   Level: beginner

772: .seealso: PetscDSGetDiscretization(), PetscDSSetDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
773: @*/
774: PetscErrorCode PetscDSAddDiscretization(PetscDS prob, PetscObject disc)
775: {

779:   PetscDSSetDiscretization(prob, prob->Nf, disc);
780:   return(0);
781: }

785: /*@
786:   PetscDSAddBdDiscretization - Adds a boundary discretization object

788:   Not collective

790:   Input Parameters:
791: + prob - The PetscDS object
792: - disc - The boundary discretization object

794:   Level: beginner

796: .seealso: PetscDSGetBdDiscretization(), PetscDSSetBdDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
797: @*/
798: PetscErrorCode PetscDSAddBdDiscretization(PetscDS prob, PetscObject disc)
799: {

803:   PetscDSSetBdDiscretization(prob, prob->Nf, disc);
804:   return(0);
805: }

809: /*@
810:   PetscDSGetImplicit - Returns the flag for implicit solve for this field. This is just a guide for IMEX

812:   Not collective

814:   Input Parameters:
815: + prob - The PetscDS object
816: - f - The field number

818:   Output Parameter:
819: . implicit - The flag indicating what kind of solve to use for this field

821:   Level: developer

823: .seealso: PetscDSSetImplicit(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetBdDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
824: @*/
825: PetscErrorCode PetscDSGetImplicit(PetscDS prob, PetscInt f, PetscBool *implicit)
826: {
830:   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
831:   *implicit = prob->implicit[f];
832:   return(0);
833: }

837: /*@
838:   PetscDSSetImplicit - Set the flag for implicit solve for this field. This is just a guide for IMEX

840:   Not collective

842:   Input Parameters:
843: + prob - The PetscDS object
844: . f - The field number
845: - implicit - The flag indicating what kind of solve to use for this field

847:   Level: developer

849: .seealso: PetscDSGetImplicit(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetBdDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
850: @*/
851: PetscErrorCode PetscDSSetImplicit(PetscDS prob, PetscInt f, PetscBool implicit)
852: {
855:   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
856:   prob->implicit[f] = implicit;
857:   return(0);
858: }

862: /*@
863:   PetscDSGetAdjacency - Returns the flags for determining variable influence

865:   Not collective

867:   Input Parameters:
868: + prob - The PetscDS object
869: - f - The field number

871:   Output Parameter:
872: + useCone    - Flag for variable influence starting with the cone operation
873: - useClosure - Flag for variable influence using transitive closure

875:   Note: See the discussion in DMPlexGetAdjacencyUseCone() and DMPlexGetAdjacencyUseClosure()

877:   Level: developer

879: .seealso: PetscDSSetAdjacency(), DMPlexGetAdjacencyUseCone(), DMPlexGetAdjacencyUseClosure(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetBdDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
880: @*/
881: PetscErrorCode PetscDSGetAdjacency(PetscDS prob, PetscInt f, PetscBool *useCone, PetscBool *useClosure)
882: {
887:   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
888:   *useCone    = prob->adjacency[f*2+0];
889:   *useClosure = prob->adjacency[f*2+1];
890:   return(0);
891: }

895: /*@
896:   PetscDSSetAdjacency - Set the flags for determining variable influence

898:   Not collective

900:   Input Parameters:
901: + prob - The PetscDS object
902: . f - The field number
903: . useCone    - Flag for variable influence starting with the cone operation
904: - useClosure - Flag for variable influence using transitive closure

906:   Note: See the discussion in DMPlexGetAdjacencyUseCone() and DMPlexGetAdjacencyUseClosure()

908:   Level: developer

910: .seealso: PetscDSGetAdjacency(), DMPlexGetAdjacencyUseCone(), DMPlexGetAdjacencyUseClosure(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetBdDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
911: @*/
912: PetscErrorCode PetscDSSetAdjacency(PetscDS prob, PetscInt f, PetscBool useCone, PetscBool useClosure)
913: {
916:   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
917:   prob->adjacency[f*2+0] = useCone;
918:   prob->adjacency[f*2+1] = useClosure;
919:   return(0);
920: }

924: PetscErrorCode PetscDSGetObjective(PetscDS prob, PetscInt f,
925:                                    void (**obj)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
926:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
927:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
928:                                                 PetscReal t, const PetscReal x[], PetscScalar obj[]))
929: {
933:   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
934:   *obj = prob->obj[f];
935:   return(0);
936: }

940: PetscErrorCode PetscDSSetObjective(PetscDS prob, PetscInt f,
941:                                    void (*obj)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
942:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
943:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
944:                                                PetscReal t, const PetscReal x[], PetscScalar obj[]))
945: {

951:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
952:   PetscDSEnlarge_Static(prob, f+1);
953:   prob->obj[f] = obj;
954:   return(0);
955: }

959: /*@C
960:   PetscDSGetResidual - Get the pointwise residual function for a given test field

962:   Not collective

964:   Input Parameters:
965: + prob - The PetscDS
966: - f    - The test field number

968:   Output Parameters:
969: + f0 - integrand for the test function term
970: - f1 - integrand for the test function gradient term

972:   Note: We are using a first order FEM model for the weak form:

974:   \int_\Omega \phi f_0(u, u_t, \nabla u, x, t) + \nabla\phi \cdot {\vec f}_1(u, u_t, \nabla u, x, t)

976: The calling sequence for the callbacks f0 and f1 is given by:

978: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
979: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
980: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
981: $    PetscReal t, const PetscReal x[], PetscScalar f0[])

983: + dim - the spatial dimension
984: . Nf - the number of fields
985: . uOff - the offset into u[] and u_t[] for each field
986: . uOff_x - the offset into u_x[] for each field
987: . u - each field evaluated at the current point
988: . u_t - the time derivative of each field evaluated at the current point
989: . u_x - the gradient of each field evaluated at the current point
990: . aOff - the offset into a[] and a_t[] for each auxiliary field
991: . aOff_x - the offset into a_x[] for each auxiliary field
992: . a - each auxiliary field evaluated at the current point
993: . a_t - the time derivative of each auxiliary field evaluated at the current point
994: . a_x - the gradient of auxiliary each field evaluated at the current point
995: . t - current time
996: . x - coordinates of the current point
997: - f0 - output values at the current point

999:   Level: intermediate

1001: .seealso: PetscDSSetResidual()
1002: @*/
1003: PetscErrorCode PetscDSGetResidual(PetscDS prob, PetscInt f,
1004:                                   void (**f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1005:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1006:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1007:                                               PetscReal t, const PetscReal x[], PetscScalar f0[]),
1008:                                   void (**f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1009:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1010:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1011:                                               PetscReal t, const PetscReal x[], PetscScalar f1[]))
1012: {
1015:   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
1018:   return(0);
1019: }

1023: /*@C
1024:   PetscDSSetResidual - Set the pointwise residual function for a given test field

1026:   Not collective

1028:   Input Parameters:
1029: + prob - The PetscDS
1030: . f    - The test field number
1031: . f0 - integrand for the test function term
1032: - f1 - integrand for the test function gradient term

1034:   Note: We are using a first order FEM model for the weak form:

1036:   \int_\Omega \phi f_0(u, u_t, \nabla u, x, t) + \nabla\phi \cdot {\vec f}_1(u, u_t, \nabla u, x, t)

1038: The calling sequence for the callbacks f0 and f1 is given by:

1040: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1041: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1042: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1043: $    PetscReal t, const PetscReal x[], PetscScalar f0[])

1045: + dim - the spatial dimension
1046: . Nf - the number of fields
1047: . uOff - the offset into u[] and u_t[] for each field
1048: . uOff_x - the offset into u_x[] for each field
1049: . u - each field evaluated at the current point
1050: . u_t - the time derivative of each field evaluated at the current point
1051: . u_x - the gradient of each field evaluated at the current point
1052: . aOff - the offset into a[] and a_t[] for each auxiliary field
1053: . aOff_x - the offset into a_x[] for each auxiliary field
1054: . a - each auxiliary field evaluated at the current point
1055: . a_t - the time derivative of each auxiliary field evaluated at the current point
1056: . a_x - the gradient of auxiliary each field evaluated at the current point
1057: . t - current time
1058: . x - coordinates of the current point
1059: - f0 - output values at the current point

1061:   Level: intermediate

1063: .seealso: PetscDSGetResidual()
1064: @*/
1065: PetscErrorCode PetscDSSetResidual(PetscDS prob, PetscInt f,
1066:                                   void (*f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1067:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1068:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1069:                                              PetscReal t, const PetscReal x[], PetscScalar f0[]),
1070:                                   void (*f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1071:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1072:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1073:                                              PetscReal t, const PetscReal x[], PetscScalar f1[]))
1074: {

1081:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1082:   PetscDSEnlarge_Static(prob, f+1);
1083:   prob->f[f*2+0] = f0;
1084:   prob->f[f*2+1] = f1;
1085:   return(0);
1086: }

1090: /*@C
1091:   PetscDSGetJacobian - Get the pointwise Jacobian function for given test and basis field

1093:   Not collective

1095:   Input Parameters:
1096: + prob - The PetscDS
1097: . f    - The test field number
1098: - g    - The field number

1100:   Output Parameters:
1101: + g0 - integrand for the test and basis function term
1102: . g1 - integrand for the test function and basis function gradient term
1103: . g2 - integrand for the test function gradient and basis function term
1104: - g3 - integrand for the test function gradient and basis function gradient term

1106:   Note: We are using a first order FEM model for the weak form:

1108:   \int_\Omega \phi g_0(u, u_t, \nabla u, x, t) \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi

1110: The calling sequence for the callbacks g0, g1, g2 and g3 is given by:

1112: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1113: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1114: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1115: $    PetscReal t, const PetscReal u_tShift, const PetscReal x[], PetscScalar g0[])

1117: + dim - the spatial dimension
1118: . Nf - the number of fields
1119: . uOff - the offset into u[] and u_t[] for each field
1120: . uOff_x - the offset into u_x[] for each field
1121: . u - each field evaluated at the current point
1122: . u_t - the time derivative of each field evaluated at the current point
1123: . u_x - the gradient of each field evaluated at the current point
1124: . aOff - the offset into a[] and a_t[] for each auxiliary field
1125: . aOff_x - the offset into a_x[] for each auxiliary field
1126: . a - each auxiliary field evaluated at the current point
1127: . a_t - the time derivative of each auxiliary field evaluated at the current point
1128: . a_x - the gradient of auxiliary each field evaluated at the current point
1129: . t - current time
1130: . u_tShift - the multiplier a for dF/dU_t
1131: . x - coordinates of the current point
1132: - g0 - output values at the current point

1134:   Level: intermediate

1136: .seealso: PetscDSSetJacobian()
1137: @*/
1138: PetscErrorCode PetscDSGetJacobian(PetscDS prob, PetscInt f, PetscInt g,
1139:                                   void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1140:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1141:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1142:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g0[]),
1143:                                   void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1144:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1145:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1146:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g1[]),
1147:                                   void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1148:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1149:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1150:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g2[]),
1151:                                   void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1152:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1153:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1154:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g3[]))
1155: {
1158:   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
1159:   if ((g < 0) || (g >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, prob->Nf);
1164:   return(0);
1165: }

1169: /*@C
1170:   PetscDSSetJacobian - Set the pointwise Jacobian function for given test and basis fields

1172:   Not collective

1174:   Input Parameters:
1175: + prob - The PetscDS
1176: . f    - The test field number
1177: . g    - The field number
1178: . g0 - integrand for the test and basis function term
1179: . g1 - integrand for the test function and basis function gradient term
1180: . g2 - integrand for the test function gradient and basis function term
1181: - g3 - integrand for the test function gradient and basis function gradient term

1183:   Note: We are using a first order FEM model for the weak form:

1185:   \int_\Omega \phi g_0(u, u_t, \nabla u, x, t) \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi

1187: The calling sequence for the callbacks g0, g1, g2 and g3 is given by:

1189: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1190: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1191: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1192: $    PetscReal t, const PetscReal x[], PetscScalar g0[])

1194: + dim - the spatial dimension
1195: . Nf - the number of fields
1196: . uOff - the offset into u[] and u_t[] for each field
1197: . uOff_x - the offset into u_x[] for each field
1198: . u - each field evaluated at the current point
1199: . u_t - the time derivative of each field evaluated at the current point
1200: . u_x - the gradient of each field evaluated at the current point
1201: . aOff - the offset into a[] and a_t[] for each auxiliary field
1202: . aOff_x - the offset into a_x[] for each auxiliary field
1203: . a - each auxiliary field evaluated at the current point
1204: . a_t - the time derivative of each auxiliary field evaluated at the current point
1205: . a_x - the gradient of auxiliary each field evaluated at the current point
1206: . t - current time
1207: . u_tShift - the multiplier a for dF/dU_t
1208: . x - coordinates of the current point
1209: - g0 - output values at the current point

1211:   Level: intermediate

1213: .seealso: PetscDSGetJacobian()
1214: @*/
1215: PetscErrorCode PetscDSSetJacobian(PetscDS prob, PetscInt f, PetscInt g,
1216:                                   void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1217:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1218:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1219:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g0[]),
1220:                                   void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1221:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1222:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1223:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g1[]),
1224:                                   void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1225:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1226:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1227:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g2[]),
1228:                                   void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1229:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1230:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1231:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g3[]))
1232: {

1241:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1242:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
1243:   PetscDSEnlarge_Static(prob, PetscMax(f, g)+1);
1244:   prob->g[(f*prob->Nf + g)*4+0] = g0;
1245:   prob->g[(f*prob->Nf + g)*4+1] = g1;
1246:   prob->g[(f*prob->Nf + g)*4+2] = g2;
1247:   prob->g[(f*prob->Nf + g)*4+3] = g3;
1248:   return(0);
1249: }

1253: /*@C
1254:   PetscDSHasJacobianPreconditioner - Signals that a Jacobian preconditioner matrix has been set

1256:   Not collective

1258:   Input Parameter:
1259: . prob - The PetscDS

1261:   Output Parameter:
1262: . hasJacPre - flag that pointwise function for Jacobian preconditioner matrix has been set

1264:   Level: intermediate

1266: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1267: @*/
1268: PetscErrorCode PetscDSHasJacobianPreconditioner(PetscDS prob, PetscBool *hasJacPre)
1269: {
1270:   PetscInt f, g, h;

1274:   *hasJacPre = PETSC_FALSE;
1275:   for (f = 0; f < prob->Nf; ++f) {
1276:     for (g = 0; g < prob->Nf; ++g) {
1277:       for (h = 0; h < 4; ++h) {
1278:         if (prob->gp[(f*prob->Nf + g)*4+h]) *hasJacPre = PETSC_TRUE;
1279:       }
1280:     }
1281:   }
1282:   return(0);
1283: }

1287: /*@C
1288:   PetscDSGetJacobianPreconditioner - Get the pointwise Jacobian preconditioner function for given test and basis field. If this is missing, the system matrix is used to build the preconditioner.

1290:   Not collective

1292:   Input Parameters:
1293: + prob - The PetscDS
1294: . f    - The test field number
1295: - g    - The field number

1297:   Output Parameters:
1298: + g0 - integrand for the test and basis function term
1299: . g1 - integrand for the test function and basis function gradient term
1300: . g2 - integrand for the test function gradient and basis function term
1301: - g3 - integrand for the test function gradient and basis function gradient term

1303:   Note: We are using a first order FEM model for the weak form:

1305:   \int_\Omega \phi g_0(u, u_t, \nabla u, x, t) \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi

1307: The calling sequence for the callbacks g0, g1, g2 and g3 is given by:

1309: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1310: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1311: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1312: $    PetscReal t, const PetscReal u_tShift, const PetscReal x[], PetscScalar g0[])

1314: + dim - the spatial dimension
1315: . Nf - the number of fields
1316: . uOff - the offset into u[] and u_t[] for each field
1317: . uOff_x - the offset into u_x[] for each field
1318: . u - each field evaluated at the current point
1319: . u_t - the time derivative of each field evaluated at the current point
1320: . u_x - the gradient of each field evaluated at the current point
1321: . aOff - the offset into a[] and a_t[] for each auxiliary field
1322: . aOff_x - the offset into a_x[] for each auxiliary field
1323: . a - each auxiliary field evaluated at the current point
1324: . a_t - the time derivative of each auxiliary field evaluated at the current point
1325: . a_x - the gradient of auxiliary each field evaluated at the current point
1326: . t - current time
1327: . u_tShift - the multiplier a for dF/dU_t
1328: . x - coordinates of the current point
1329: - g0 - output values at the current point

1331:   Level: intermediate

1333: .seealso: PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1334: @*/
1335: PetscErrorCode PetscDSGetJacobianPreconditioner(PetscDS prob, PetscInt f, PetscInt g,
1336:                                   void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1337:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1338:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1339:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g0[]),
1340:                                   void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1341:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1342:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1343:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g1[]),
1344:                                   void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1345:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1346:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1347:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g2[]),
1348:                                   void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1349:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1350:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1351:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g3[]))
1352: {
1355:   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
1356:   if ((g < 0) || (g >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, prob->Nf);
1361:   return(0);
1362: }

1366: /*@C
1367:   PetscDSSetJacobianPreconditioner - Set the pointwise Jacobian preconditioner function for given test and basis fields. If this is missing, the system matrix is used to build the preconditioner.

1369:   Not collective

1371:   Input Parameters:
1372: + prob - The PetscDS
1373: . f    - The test field number
1374: . g    - The field number
1375: . g0 - integrand for the test and basis function term
1376: . g1 - integrand for the test function and basis function gradient term
1377: . g2 - integrand for the test function gradient and basis function term
1378: - g3 - integrand for the test function gradient and basis function gradient term

1380:   Note: We are using a first order FEM model for the weak form:

1382:   \int_\Omega \phi g_0(u, u_t, \nabla u, x, t) \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi

1384: The calling sequence for the callbacks g0, g1, g2 and g3 is given by:

1386: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1387: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1388: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1389: $    PetscReal t, const PetscReal x[], PetscScalar g0[])

1391: + dim - the spatial dimension
1392: . Nf - the number of fields
1393: . uOff - the offset into u[] and u_t[] for each field
1394: . uOff_x - the offset into u_x[] for each field
1395: . u - each field evaluated at the current point
1396: . u_t - the time derivative of each field evaluated at the current point
1397: . u_x - the gradient of each field evaluated at the current point
1398: . aOff - the offset into a[] and a_t[] for each auxiliary field
1399: . aOff_x - the offset into a_x[] for each auxiliary field
1400: . a - each auxiliary field evaluated at the current point
1401: . a_t - the time derivative of each auxiliary field evaluated at the current point
1402: . a_x - the gradient of auxiliary each field evaluated at the current point
1403: . t - current time
1404: . u_tShift - the multiplier a for dF/dU_t
1405: . x - coordinates of the current point
1406: - g0 - output values at the current point

1408:   Level: intermediate

1410: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobian()
1411: @*/
1412: PetscErrorCode PetscDSSetJacobianPreconditioner(PetscDS prob, PetscInt f, PetscInt g,
1413:                                   void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1414:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1415:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1416:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g0[]),
1417:                                   void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1418:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1419:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1420:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g1[]),
1421:                                   void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1422:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1423:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1424:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g2[]),
1425:                                   void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1426:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1427:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1428:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g3[]))
1429: {

1438:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1439:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
1440:   PetscDSEnlarge_Static(prob, PetscMax(f, g)+1);
1441:   prob->gp[(f*prob->Nf + g)*4+0] = g0;
1442:   prob->gp[(f*prob->Nf + g)*4+1] = g1;
1443:   prob->gp[(f*prob->Nf + g)*4+2] = g2;
1444:   prob->gp[(f*prob->Nf + g)*4+3] = g3;
1445:   return(0);
1446: }

1450: /*@C
1451:   PetscDSHasDynamicJacobian - Signals that a dynamic Jacobian, dF/du_t, has been set

1453:   Not collective

1455:   Input Parameter:
1456: . prob - The PetscDS

1458:   Output Parameter:
1459: . hasDynJac - flag that pointwise function for dynamic Jacobian has been set

1461:   Level: intermediate

1463: .seealso: PetscDSGetDynamicJacobian(), PetscDSSetDynamicJacobian(), PetscDSGetJacobian()
1464: @*/
1465: PetscErrorCode PetscDSHasDynamicJacobian(PetscDS prob, PetscBool *hasDynJac)
1466: {
1467:   PetscInt f, g, h;

1471:   *hasDynJac = PETSC_FALSE;
1472:   for (f = 0; f < prob->Nf; ++f) {
1473:     for (g = 0; g < prob->Nf; ++g) {
1474:       for (h = 0; h < 4; ++h) {
1475:         if (prob->gt[(f*prob->Nf + g)*4+h]) *hasDynJac = PETSC_TRUE;
1476:       }
1477:     }
1478:   }
1479:   return(0);
1480: }

1484: /*@C
1485:   PetscDSGetDynamicJacobian - Get the pointwise dynamic Jacobian, dF/du_t, function for given test and basis field

1487:   Not collective

1489:   Input Parameters:
1490: + prob - The PetscDS
1491: . f    - The test field number
1492: - g    - The field number

1494:   Output Parameters:
1495: + g0 - integrand for the test and basis function term
1496: . g1 - integrand for the test function and basis function gradient term
1497: . g2 - integrand for the test function gradient and basis function term
1498: - g3 - integrand for the test function gradient and basis function gradient term

1500:   Note: We are using a first order FEM model for the weak form:

1502:   \int_\Omega \phi g_0(u, u_t, \nabla u, x, t) \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi

1504: The calling sequence for the callbacks g0, g1, g2 and g3 is given by:

1506: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1507: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1508: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1509: $    PetscReal t, const PetscReal u_tShift, const PetscReal x[], PetscScalar g0[])

1511: + dim - the spatial dimension
1512: . Nf - the number of fields
1513: . uOff - the offset into u[] and u_t[] for each field
1514: . uOff_x - the offset into u_x[] for each field
1515: . u - each field evaluated at the current point
1516: . u_t - the time derivative of each field evaluated at the current point
1517: . u_x - the gradient of each field evaluated at the current point
1518: . aOff - the offset into a[] and a_t[] for each auxiliary field
1519: . aOff_x - the offset into a_x[] for each auxiliary field
1520: . a - each auxiliary field evaluated at the current point
1521: . a_t - the time derivative of each auxiliary field evaluated at the current point
1522: . a_x - the gradient of auxiliary each field evaluated at the current point
1523: . t - current time
1524: . u_tShift - the multiplier a for dF/dU_t
1525: . x - coordinates of the current point
1526: - g0 - output values at the current point

1528:   Level: intermediate

1530: .seealso: PetscDSSetJacobian()
1531: @*/
1532: PetscErrorCode PetscDSGetDynamicJacobian(PetscDS prob, PetscInt f, PetscInt g,
1533:                                          void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1534:                                                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1535:                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1536:                                                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g0[]),
1537:                                          void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1538:                                                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1539:                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1540:                                                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g1[]),
1541:                                          void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1542:                                                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1543:                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1544:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g2[]),
1545:                                          void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1546:                                                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1547:                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1548:                                                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g3[]))
1549: {
1552:   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
1553:   if ((g < 0) || (g >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, prob->Nf);
1558:   return(0);
1559: }

1563: /*@C
1564:   PetscDSSetDynamicJacobian - Set the pointwise dynamic Jacobian, dF/du_t, function for given test and basis fields

1566:   Not collective

1568:   Input Parameters:
1569: + prob - The PetscDS
1570: . f    - The test field number
1571: . g    - The field number
1572: . g0 - integrand for the test and basis function term
1573: . g1 - integrand for the test function and basis function gradient term
1574: . g2 - integrand for the test function gradient and basis function term
1575: - g3 - integrand for the test function gradient and basis function gradient term

1577:   Note: We are using a first order FEM model for the weak form:

1579:   \int_\Omega \phi g_0(u, u_t, \nabla u, x, t) \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi

1581: The calling sequence for the callbacks g0, g1, g2 and g3 is given by:

1583: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1584: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1585: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1586: $    PetscReal t, const PetscReal x[], PetscScalar g0[])

1588: + dim - the spatial dimension
1589: . Nf - the number of fields
1590: . uOff - the offset into u[] and u_t[] for each field
1591: . uOff_x - the offset into u_x[] for each field
1592: . u - each field evaluated at the current point
1593: . u_t - the time derivative of each field evaluated at the current point
1594: . u_x - the gradient of each field evaluated at the current point
1595: . aOff - the offset into a[] and a_t[] for each auxiliary field
1596: . aOff_x - the offset into a_x[] for each auxiliary field
1597: . a - each auxiliary field evaluated at the current point
1598: . a_t - the time derivative of each auxiliary field evaluated at the current point
1599: . a_x - the gradient of auxiliary each field evaluated at the current point
1600: . t - current time
1601: . u_tShift - the multiplier a for dF/dU_t
1602: . x - coordinates of the current point
1603: - g0 - output values at the current point

1605:   Level: intermediate

1607: .seealso: PetscDSGetJacobian()
1608: @*/
1609: PetscErrorCode PetscDSSetDynamicJacobian(PetscDS prob, PetscInt f, PetscInt g,
1610:                                          void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1611:                                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1612:                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1613:                                                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g0[]),
1614:                                          void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1615:                                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1616:                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1617:                                                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g1[]),
1618:                                          void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1619:                                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1620:                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1621:                                                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g2[]),
1622:                                          void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1623:                                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1624:                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1625:                                                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g3[]))
1626: {

1635:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1636:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
1637:   PetscDSEnlarge_Static(prob, PetscMax(f, g)+1);
1638:   prob->gt[(f*prob->Nf + g)*4+0] = g0;
1639:   prob->gt[(f*prob->Nf + g)*4+1] = g1;
1640:   prob->gt[(f*prob->Nf + g)*4+2] = g2;
1641:   prob->gt[(f*prob->Nf + g)*4+3] = g3;
1642:   return(0);
1643: }

1647: /*@C
1648:   PetscDSGetRiemannSolver - Returns the Riemann solver for the given field

1650:   Not collective

1652:   Input Arguments:
1653: + prob - The PetscDS object
1654: - f    - The field number

1656:   Output Argument:
1657: . r    - Riemann solver

1659:   Calling sequence for r:

1661: $ r(PetscInt dim, PetscInt Nf, const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscScalar flux[], void *ctx)

1663: + dim  - The spatial dimension
1664: . Nf   - The number of fields
1665: . x    - The coordinates at a point on the interface
1666: . n    - The normal vector to the interface
1667: . uL   - The state vector to the left of the interface
1668: . uR   - The state vector to the right of the interface
1669: . flux - output array of flux through the interface
1670: - ctx  - optional user context

1672:   Level: intermediate

1674: .seealso: PetscDSSetRiemannSolver()
1675: @*/
1676: PetscErrorCode PetscDSGetRiemannSolver(PetscDS prob, PetscInt f,
1677:                                        void (**r)(PetscInt dim, PetscInt Nf, const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscScalar flux[], void *ctx))
1678: {
1681:   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
1683:   *r = prob->r[f];
1684:   return(0);
1685: }

1689: /*@C
1690:   PetscDSSetRiemannSolver - Sets the Riemann solver for the given field

1692:   Not collective

1694:   Input Arguments:
1695: + prob - The PetscDS object
1696: . f    - The field number
1697: - r    - Riemann solver

1699:   Calling sequence for r:

1701: $ r(PetscInt dim, PetscInt Nf, const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscScalar flux[], void *ctx)

1703: + dim  - The spatial dimension
1704: . Nf   - The number of fields
1705: . x    - The coordinates at a point on the interface
1706: . n    - The normal vector to the interface
1707: . uL   - The state vector to the left of the interface
1708: . uR   - The state vector to the right of the interface
1709: . flux - output array of flux through the interface
1710: - ctx  - optional user context

1712:   Level: intermediate

1714: .seealso: PetscDSGetRiemannSolver()
1715: @*/
1716: PetscErrorCode PetscDSSetRiemannSolver(PetscDS prob, PetscInt f,
1717:                                        void (*r)(PetscInt dim, PetscInt Nf, const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscScalar flux[], void *ctx))
1718: {

1724:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1725:   PetscDSEnlarge_Static(prob, f+1);
1726:   prob->r[f] = r;
1727:   return(0);
1728: }

1732: PetscErrorCode PetscDSGetContext(PetscDS prob, PetscInt f, void **ctx)
1733: {
1736:   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
1738:   *ctx = prob->ctx[f];
1739:   return(0);
1740: }

1744: PetscErrorCode PetscDSSetContext(PetscDS prob, PetscInt f, void *ctx)
1745: {

1750:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1751:   PetscDSEnlarge_Static(prob, f+1);
1752:   prob->ctx[f] = ctx;
1753:   return(0);
1754: }

1758: /*@C
1759:   PetscDSGetBdResidual - Get the pointwise boundary residual function for a given test field

1761:   Not collective

1763:   Input Parameters:
1764: + prob - The PetscDS
1765: - f    - The test field number

1767:   Output Parameters:
1768: + f0 - boundary integrand for the test function term
1769: - f1 - boundary integrand for the test function gradient term

1771:   Note: We are using a first order FEM model for the weak form:

1773:   \int_\Gamma \phi {\vec f}_0(u, u_t, \nabla u, x, t) \cdot \hat n + \nabla\phi \cdot {\overleftrightarrow f}_1(u, u_t, \nabla u, x, t) \cdot \hat n

1775: The calling sequence for the callbacks f0 and f1 is given by:

1777: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1778: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1779: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1780: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[])

1782: + dim - the spatial dimension
1783: . Nf - the number of fields
1784: . uOff - the offset into u[] and u_t[] for each field
1785: . uOff_x - the offset into u_x[] for each field
1786: . u - each field evaluated at the current point
1787: . u_t - the time derivative of each field evaluated at the current point
1788: . u_x - the gradient of each field evaluated at the current point
1789: . aOff - the offset into a[] and a_t[] for each auxiliary field
1790: . aOff_x - the offset into a_x[] for each auxiliary field
1791: . a - each auxiliary field evaluated at the current point
1792: . a_t - the time derivative of each auxiliary field evaluated at the current point
1793: . a_x - the gradient of auxiliary each field evaluated at the current point
1794: . t - current time
1795: . x - coordinates of the current point
1796: . n - unit normal at the current point
1797: - f0 - output values at the current point

1799:   Level: intermediate

1801: .seealso: PetscDSSetBdResidual()
1802: @*/
1803: PetscErrorCode PetscDSGetBdResidual(PetscDS prob, PetscInt f,
1804:                                     void (**f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1805:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1806:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1807:                                                 PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[]),
1808:                                     void (**f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1809:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1810:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1811:                                                 PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f1[]))
1812: {
1815:   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
1818:   return(0);
1819: }

1823: /*@C
1824:   PetscDSSetBdResidual - Get the pointwise boundary residual function for a given test field

1826:   Not collective

1828:   Input Parameters:
1829: + prob - The PetscDS
1830: . f    - The test field number
1831: . f0 - boundary integrand for the test function term
1832: - f1 - boundary integrand for the test function gradient term

1834:   Note: We are using a first order FEM model for the weak form:

1836:   \int_\Gamma \phi {\vec f}_0(u, u_t, \nabla u, x, t) \cdot \hat n + \nabla\phi \cdot {\overleftrightarrow f}_1(u, u_t, \nabla u, x, t) \cdot \hat n

1838: The calling sequence for the callbacks f0 and f1 is given by:

1840: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1841: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1842: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1843: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[])

1845: + dim - the spatial dimension
1846: . Nf - the number of fields
1847: . uOff - the offset into u[] and u_t[] for each field
1848: . uOff_x - the offset into u_x[] for each field
1849: . u - each field evaluated at the current point
1850: . u_t - the time derivative of each field evaluated at the current point
1851: . u_x - the gradient of each field evaluated at the current point
1852: . aOff - the offset into a[] and a_t[] for each auxiliary field
1853: . aOff_x - the offset into a_x[] for each auxiliary field
1854: . a - each auxiliary field evaluated at the current point
1855: . a_t - the time derivative of each auxiliary field evaluated at the current point
1856: . a_x - the gradient of auxiliary each field evaluated at the current point
1857: . t - current time
1858: . x - coordinates of the current point
1859: . n - unit normal at the current point
1860: - f0 - output values at the current point

1862:   Level: intermediate

1864: .seealso: PetscDSGetBdResidual()
1865: @*/
1866: PetscErrorCode PetscDSSetBdResidual(PetscDS prob, PetscInt f,
1867:                                     void (*f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1868:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1869:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1870:                                                PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[]),
1871:                                     void (*f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1872:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1873:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1874:                                                PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f1[]))
1875: {

1880:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1881:   PetscDSEnlarge_Static(prob, f+1);
1884:   return(0);
1885: }

1889: /*@C
1890:   PetscDSGetBdJacobian - Get the pointwise boundary Jacobian function for given test and basis field

1892:   Not collective

1894:   Input Parameters:
1895: + prob - The PetscDS
1896: . f    - The test field number
1897: - g    - The field number

1899:   Output Parameters:
1900: + g0 - integrand for the test and basis function term
1901: . g1 - integrand for the test function and basis function gradient term
1902: . g2 - integrand for the test function gradient and basis function term
1903: - g3 - integrand for the test function gradient and basis function gradient term

1905:   Note: We are using a first order FEM model for the weak form:

1907:   \int_\Gamma \phi {\vec g}_0(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \cdot \hat n \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \hat n \cdot \nabla \psi

1909: The calling sequence for the callbacks g0, g1, g2 and g3 is given by:

1911: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1912: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1913: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1914: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar g0[])

1916: + dim - the spatial dimension
1917: . Nf - the number of fields
1918: . uOff - the offset into u[] and u_t[] for each field
1919: . uOff_x - the offset into u_x[] for each field
1920: . u - each field evaluated at the current point
1921: . u_t - the time derivative of each field evaluated at the current point
1922: . u_x - the gradient of each field evaluated at the current point
1923: . aOff - the offset into a[] and a_t[] for each auxiliary field
1924: . aOff_x - the offset into a_x[] for each auxiliary field
1925: . a - each auxiliary field evaluated at the current point
1926: . a_t - the time derivative of each auxiliary field evaluated at the current point
1927: . a_x - the gradient of auxiliary each field evaluated at the current point
1928: . t - current time
1929: . u_tShift - the multiplier a for dF/dU_t
1930: . x - coordinates of the current point
1931: . n - normal at the current point
1932: - g0 - output values at the current point

1934:   Level: intermediate

1936: .seealso: PetscDSSetBdJacobian()
1937: @*/
1938: PetscErrorCode PetscDSGetBdJacobian(PetscDS prob, PetscInt f, PetscInt g,
1939:                                     void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1940:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1941:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1942:                                                 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscScalar g0[]),
1943:                                     void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1944:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1945:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1946:                                                 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscScalar g1[]),
1947:                                     void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1948:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1949:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1950:                                                 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscScalar g2[]),
1951:                                     void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1952:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1953:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1954:                                                 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscScalar g3[]))
1955: {
1958:   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
1959:   if ((g < 0) || (g >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, prob->Nf);
1964:   return(0);
1965: }

1969: /*@C
1970:   PetscDSSetBdJacobian - Set the pointwise boundary Jacobian function for given test and basis field

1972:   Not collective

1974:   Input Parameters:
1975: + prob - The PetscDS
1976: . f    - The test field number
1977: . g    - The field number
1978: . g0 - integrand for the test and basis function term
1979: . g1 - integrand for the test function and basis function gradient term
1980: . g2 - integrand for the test function gradient and basis function term
1981: - g3 - integrand for the test function gradient and basis function gradient term

1983:   Note: We are using a first order FEM model for the weak form:

1985:   \int_\Gamma \phi {\vec g}_0(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \cdot \hat n \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \hat n \cdot \nabla \psi

1987: The calling sequence for the callbacks g0, g1, g2 and g3 is given by:

1989: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1990: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1991: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1992: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar g0[])

1994: + dim - the spatial dimension
1995: . Nf - the number of fields
1996: . uOff - the offset into u[] and u_t[] for each field
1997: . uOff_x - the offset into u_x[] for each field
1998: . u - each field evaluated at the current point
1999: . u_t - the time derivative of each field evaluated at the current point
2000: . u_x - the gradient of each field evaluated at the current point
2001: . aOff - the offset into a[] and a_t[] for each auxiliary field
2002: . aOff_x - the offset into a_x[] for each auxiliary field
2003: . a - each auxiliary field evaluated at the current point
2004: . a_t - the time derivative of each auxiliary field evaluated at the current point
2005: . a_x - the gradient of auxiliary each field evaluated at the current point
2006: . t - current time
2007: . u_tShift - the multiplier a for dF/dU_t
2008: . x - coordinates of the current point
2009: . n - normal at the current point
2010: - g0 - output values at the current point

2012:   Level: intermediate

2014: .seealso: PetscDSGetBdJacobian()
2015: @*/
2016: PetscErrorCode PetscDSSetBdJacobian(PetscDS prob, PetscInt f, PetscInt g,
2017:                                     void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2018:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2019:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2020:                                                PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscScalar g0[]),
2021:                                     void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2022:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2023:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2024:                                                PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscScalar g1[]),
2025:                                     void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2026:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2027:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2028:                                                PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscScalar g2[]),
2029:                                     void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2030:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2031:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2032:                                                PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscScalar g3[]))
2033: {

2042:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2043:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
2044:   PetscDSEnlarge_Static(prob, PetscMax(f, g)+1);
2045:   prob->gBd[(f*prob->Nf + g)*4+0] = g0;
2046:   prob->gBd[(f*prob->Nf + g)*4+1] = g1;
2047:   prob->gBd[(f*prob->Nf + g)*4+2] = g2;
2048:   prob->gBd[(f*prob->Nf + g)*4+3] = g3;
2049:   return(0);
2050: }

2054: /*@
2055:   PetscDSGetFieldOffset - Returns the offset of the given field in the full space basis

2057:   Not collective

2059:   Input Parameters:
2060: + prob - The PetscDS object
2061: - f - The field number

2063:   Output Parameter:
2064: . off - The offset

2066:   Level: beginner

2068: .seealso: PetscDSGetBdFieldOffset(), PetscDSGetNumFields(), PetscDSCreate()
2069: @*/
2070: PetscErrorCode PetscDSGetFieldOffset(PetscDS prob, PetscInt f, PetscInt *off)
2071: {
2072:   PetscInt       g;

2078:   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
2079:   *off = 0;
2080:   for (g = 0; g < f; ++g) {
2081:     PetscFE  fe = (PetscFE) prob->disc[g];
2082:     PetscInt Nb, Nc;

2084:     PetscFEGetDimension(fe, &Nb);
2085:     PetscFEGetNumComponents(fe, &Nc);
2086:     *off += Nb*Nc;
2087:   }
2088:   return(0);
2089: }

2093: /*@
2094:   PetscDSGetBdFieldOffset - Returns the offset of the given field in the full space boundary basis

2096:   Not collective

2098:   Input Parameters:
2099: + prob - The PetscDS object
2100: - f - The field number

2102:   Output Parameter:
2103: . off - The boundary offset

2105:   Level: beginner

2107: .seealso: PetscDSGetFieldOffset(), PetscDSGetNumFields(), PetscDSCreate()
2108: @*/
2109: PetscErrorCode PetscDSGetBdFieldOffset(PetscDS prob, PetscInt f, PetscInt *off)
2110: {
2111:   PetscInt       g;

2117:   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
2118:   *off = 0;
2119:   for (g = 0; g < f; ++g) {
2120:     PetscFE  fe = (PetscFE) prob->discBd[g];
2121:     PetscInt Nb, Nc;

2123:     PetscFEGetDimension(fe, &Nb);
2124:     PetscFEGetNumComponents(fe, &Nc);
2125:     *off += Nb*Nc;
2126:   }
2127:   return(0);
2128: }

2132: /*@
2133:   PetscDSGetComponentOffset - Returns the offset of the given field on an evaluation point

2135:   Not collective

2137:   Input Parameters:
2138: + prob - The PetscDS object
2139: - f - The field number

2141:   Output Parameter:
2142: . off - The offset

2144:   Level: beginner

2146: .seealso: PetscDSGetBdFieldOffset(), PetscDSGetNumFields(), PetscDSCreate()
2147: @*/
2148: PetscErrorCode PetscDSGetComponentOffset(PetscDS prob, PetscInt f, PetscInt *off)
2149: {
2150:   PetscInt       g;

2156:   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
2157:   *off = 0;
2158:   for (g = 0; g < f; ++g) {
2159:     PetscFE  fe = (PetscFE) prob->disc[g];
2160:     PetscInt Nc;

2162:     PetscFEGetNumComponents(fe, &Nc);
2163:     *off += Nc;
2164:   }
2165:   return(0);
2166: }

2170: /*@
2171:   PetscDSGetComponentOffsets - Returns the offset of each field on an evaluation point

2173:   Not collective

2175:   Input Parameter:
2176: . prob - The PetscDS object

2178:   Output Parameter:
2179: . offsets - The offsets

2181:   Level: beginner

2183: .seealso: PetscDSGetBdFieldOffset(), PetscDSGetNumFields(), PetscDSCreate()
2184: @*/
2185: PetscErrorCode PetscDSGetComponentOffsets(PetscDS prob, PetscInt *offsets[])
2186: {
2190:   *offsets = prob->off;
2191:   return(0);
2192: }

2196: /*@
2197:   PetscDSGetComponentDerivativeOffsets - Returns the offset of each field derivative on an evaluation point

2199:   Not collective

2201:   Input Parameter:
2202: . prob - The PetscDS object

2204:   Output Parameter:
2205: . offsets - The offsets

2207:   Level: beginner

2209: .seealso: PetscDSGetBdFieldOffset(), PetscDSGetNumFields(), PetscDSCreate()
2210: @*/
2211: PetscErrorCode PetscDSGetComponentDerivativeOffsets(PetscDS prob, PetscInt *offsets[])
2212: {
2216:   *offsets = prob->offDer;
2217:   return(0);
2218: }

2222: /*@
2223:   PetscDSGetComponentBdOffsets - Returns the offset of each field on a boundary evaluation point

2225:   Not collective

2227:   Input Parameter:
2228: . prob - The PetscDS object

2230:   Output Parameter:
2231: . offsets - The offsets

2233:   Level: beginner

2235: .seealso: PetscDSGetBdFieldOffset(), PetscDSGetNumFields(), PetscDSCreate()
2236: @*/
2237: PetscErrorCode PetscDSGetComponentBdOffsets(PetscDS prob, PetscInt *offsets[])
2238: {
2242:   *offsets = prob->offBd;
2243:   return(0);
2244: }

2248: /*@
2249:   PetscDSGetComponentBdDerivativeOffsets - Returns the offset of each field derivative on a boundary evaluation point

2251:   Not collective

2253:   Input Parameter:
2254: . prob - The PetscDS object

2256:   Output Parameter:
2257: . offsets - The offsets

2259:   Level: beginner

2261: .seealso: PetscDSGetBdFieldOffset(), PetscDSGetNumFields(), PetscDSCreate()
2262: @*/
2263: PetscErrorCode PetscDSGetComponentBdDerivativeOffsets(PetscDS prob, PetscInt *offsets[])
2264: {
2268:   *offsets = prob->offDerBd;
2269:   return(0);
2270: }

2274: /*@C
2275:   PetscDSGetTabulation - Return the basis tabulation at quadrature points for the volume discretization

2277:   Not collective

2279:   Input Parameter:
2280: . prob - The PetscDS object

2282:   Output Parameters:
2283: + basis - The basis function tabulation at quadrature points
2284: - basisDer - The basis function derivative tabulation at quadrature points

2286:   Level: intermediate

2288: .seealso: PetscDSGetBdTabulation(), PetscDSCreate()
2289: @*/
2290: PetscErrorCode PetscDSGetTabulation(PetscDS prob, PetscReal ***basis, PetscReal ***basisDer)
2291: {

2296:   PetscDSSetUp(prob);
2299:   return(0);
2300: }

2304: /*@C
2305:   PetscDSGetBdTabulation - Return the basis tabulation at quadrature points for the boundary discretization

2307:   Not collective

2309:   Input Parameter:
2310: . prob - The PetscDS object

2312:   Output Parameters:
2313: + basis - The basis function tabulation at quadrature points
2314: - basisDer - The basis function derivative tabulation at quadrature points

2316:   Level: intermediate

2318: .seealso: PetscDSGetTabulation(), PetscDSCreate()
2319: @*/
2320: PetscErrorCode PetscDSGetBdTabulation(PetscDS prob, PetscReal ***basis, PetscReal ***basisDer)
2321: {

2326:   PetscDSSetUp(prob);
2329:   return(0);
2330: }

2334: PetscErrorCode PetscDSGetEvaluationArrays(PetscDS prob, PetscScalar **u, PetscScalar **u_t, PetscScalar **u_x)
2335: {

2340:   PetscDSSetUp(prob);
2344:   return(0);
2345: }

2349: PetscErrorCode PetscDSGetWeakFormArrays(PetscDS prob, PetscScalar **f0, PetscScalar **f1, PetscScalar **g0, PetscScalar **g1, PetscScalar **g2, PetscScalar **g3)
2350: {

2355:   PetscDSSetUp(prob);
2362:   return(0);
2363: }

2367: PetscErrorCode PetscDSGetRefCoordArrays(PetscDS prob, PetscReal **x, PetscScalar **refSpaceDer)
2368: {

2373:   PetscDSSetUp(prob);
2376:   return(0);
2377: }

2381: /*@
2382:   PetscDSCopyEquations - Copy all pointwise function pointers to the new problem

2384:   Not collective

2386:   Input Parameter:
2387: . prob - The PetscDS object

2389:   Output Parameter:
2390: . newprob - The PetscDS copy

2392:   Level: intermediate

2394: .seealso: PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
2395: @*/
2396: PetscErrorCode PetscDSCopyEquations(PetscDS prob, PetscDS newprob)
2397: {
2398:   PetscInt       Nf, Ng, f, g;

2404:   PetscDSGetNumFields(prob, &Nf);
2405:   PetscDSGetNumFields(newprob, &Ng);
2406:   if (Nf != Ng) SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_SIZ, "Number of fields must match %D != %D", Nf, Ng);
2407:   for (f = 0; f < Nf; ++f) {
2408:     PetscPointFunc   obj;
2409:     PetscPointFunc   f0, f1;
2410:     PetscPointJac    g0, g1, g2, g3;
2411:     PetscBdPointFunc f0Bd, f1Bd;
2412:     PetscBdPointJac  g0Bd, g1Bd, g2Bd, g3Bd;
2413:     PetscRiemannFunc r;

2415:     PetscDSGetObjective(prob, f, &obj);
2416:     PetscDSGetResidual(prob, f, &f0, &f1);
2417:     PetscDSGetBdResidual(prob, f, &f0Bd, &f1Bd);
2418:     PetscDSGetRiemannSolver(prob, f, &r);
2419:     PetscDSSetObjective(newprob, f, obj);
2420:     PetscDSSetResidual(newprob, f, f0, f1);
2421:     PetscDSSetBdResidual(newprob, f, f0Bd, f1Bd);
2422:     PetscDSSetRiemannSolver(newprob, f, r);
2423:     for (g = 0; g < Nf; ++g) {
2424:       PetscDSGetJacobian(prob, f, g, &g0, &g1, &g2, &g3);
2425:       PetscDSGetBdJacobian(prob, f, g, &g0Bd, &g1Bd, &g2Bd, &g3Bd);
2426:       PetscDSSetJacobian(newprob, f, g, g0, g1, g2, g3);
2427:       PetscDSSetBdJacobian(newprob, f, g, g0Bd, g1Bd, g2Bd, g3Bd);
2428:     }
2429:   }
2430:   return(0);
2431: }

2435: static PetscErrorCode PetscDSDestroy_Basic(PetscDS prob)
2436: {
2438:   return(0);
2439: }

2443: static PetscErrorCode PetscDSInitialize_Basic(PetscDS prob)
2444: {
2446:   prob->ops->setfromoptions = NULL;
2447:   prob->ops->setup          = NULL;
2448:   prob->ops->view           = NULL;
2449:   prob->ops->destroy        = PetscDSDestroy_Basic;
2450:   return(0);
2451: }

2453: /*MC
2454:   PETSCDSBASIC = "basic" - A discrete system with pointwise residual and boundary residual functions

2456:   Level: intermediate

2458: .seealso: PetscDSType, PetscDSCreate(), PetscDSSetType()
2459: M*/

2463: PETSC_EXTERN PetscErrorCode PetscDSCreate_Basic(PetscDS prob)
2464: {
2465:   PetscDS_Basic *b;
2466:   PetscErrorCode      ierr;

2470:   PetscNewLog(prob, &b);
2471:   prob->data = b;

2473:   PetscDSInitialize_Basic(prob);
2474:   return(0);
2475: }