Actual source code: dtds.c

petsc-3.15.0 2021-04-05
Report Typos and Errors
  1: #include <petsc/private/petscdsimpl.h>

  3: PetscClassId PETSCDS_CLASSID = 0;

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

  8: /* A PetscDS (Discrete System) encodes a set of equations posed in a discrete space, which represents a set of
  9:    nonlinear continuum equations. The equations can have multiple fields, each field having a different
 10:    discretization. In addition, different pieces of the domain can have different field combinations and equations.

 12:    The DS provides the user a description of the approximation space on any given cell. It also gives pointwise
 13:    functions representing the equations.

 15:    Each field is associated with a label, marking the cells on which it is supported. Note that a field can be
 16:    supported on the closure of a cell not in the label due to overlap of the boundary of neighboring cells. The DM
 17:    then creates a DS for each set of cells with identical approximation spaces. When assembling, the user asks for
 18:    the space associated with a given cell. DMPlex uses the labels associated with each DS in the default integration loop.
 19: */

 21: /*@C
 22:   PetscDSRegister - Adds a new PetscDS implementation

 24:   Not Collective

 26:   Input Parameters:
 27: + name        - The name of a new user-defined creation routine
 28: - create_func - The creation routine itself

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

 33:   Sample usage:
 34: .vb
 35:     PetscDSRegister("my_ds", MyPetscDSCreate);
 36: .ve

 38:   Then, your PetscDS type can be chosen with the procedural interface via
 39: .vb
 40:     PetscDSCreate(MPI_Comm, PetscDS *);
 41:     PetscDSSetType(PetscDS, "my_ds");
 42: .ve
 43:    or at runtime via the option
 44: .vb
 45:     -petscds_type my_ds
 46: .ve

 48:   Level: advanced

 50:    Not available from Fortran

 52: .seealso: PetscDSRegisterAll(), PetscDSRegisterDestroy()

 54: @*/
 55: PetscErrorCode PetscDSRegister(const char sname[], PetscErrorCode (*function)(PetscDS))
 56: {

 60:   PetscFunctionListAdd(&PetscDSList, sname, function);
 61:   return(0);
 62: }

 64: /*@C
 65:   PetscDSSetType - Builds a particular PetscDS

 67:   Collective on prob

 69:   Input Parameters:
 70: + prob - The PetscDS object
 71: - name - The kind of system

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

 76:   Level: intermediate

 78:    Not available from Fortran

 80: .seealso: PetscDSGetType(), PetscDSCreate()
 81: @*/
 82: PetscErrorCode PetscDSSetType(PetscDS prob, PetscDSType name)
 83: {
 84:   PetscErrorCode (*r)(PetscDS);
 85:   PetscBool      match;

 90:   PetscObjectTypeCompare((PetscObject) prob, name, &match);
 91:   if (match) return(0);

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

 97:   if (prob->ops->destroy) {
 98:     (*prob->ops->destroy)(prob);
 99:     prob->ops->destroy = NULL;
100:   }
101:   (*r)(prob);
102:   PetscObjectChangeTypeName((PetscObject) prob, name);
103:   return(0);
104: }

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

109:   Not Collective

111:   Input Parameter:
112: . prob  - The PetscDS

114:   Output Parameter:
115: . name - The PetscDS type name

117:   Level: intermediate

119:    Not available from Fortran

121: .seealso: PetscDSSetType(), PetscDSCreate()
122: @*/
123: PetscErrorCode PetscDSGetType(PetscDS prob, PetscDSType *name)
124: {

130:   PetscDSRegisterAll();
131:   *name = ((PetscObject) prob)->type_name;
132:   return(0);
133: }

135: static PetscErrorCode PetscDSView_Ascii(PetscDS prob, PetscViewer viewer)
136: {
137:   PetscViewerFormat  format;
138:   const PetscScalar *constants;
139:   PetscInt           numConstants, f;
140:   PetscErrorCode     ierr;

143:   PetscViewerGetFormat(viewer, &format);
144:   PetscViewerASCIIPrintf(viewer, "Discrete System with %d fields\n", prob->Nf);
145:   PetscViewerASCIIPushTab(viewer);
146:   PetscViewerASCIIPrintf(viewer, "  cell total dim %D total comp %D\n", prob->totDim, prob->totComp);
147:   if (prob->isHybrid) {PetscViewerASCIIPrintf(viewer, "  hybrid cell\n");}
148:   for (f = 0; f < prob->Nf; ++f) {
149:     DSBoundary      b;
150:     PetscObject     obj;
151:     PetscClassId    id;
152:     PetscQuadrature q;
153:     const char     *name;
154:     PetscInt        Nc, Nq, Nqc;

156:     PetscDSGetDiscretization(prob, f, &obj);
157:     PetscObjectGetClassId(obj, &id);
158:     PetscObjectGetName(obj, &name);
159:     PetscViewerASCIIPrintf(viewer, "Field %s", name ? name : "<unknown>");
160:     PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
161:     if (id == PETSCFE_CLASSID)      {
162:       PetscFEGetNumComponents((PetscFE) obj, &Nc);
163:       PetscFEGetQuadrature((PetscFE) obj, &q);
164:       PetscViewerASCIIPrintf(viewer, " FEM");
165:     } else if (id == PETSCFV_CLASSID) {
166:       PetscFVGetNumComponents((PetscFV) obj, &Nc);
167:       PetscFVGetQuadrature((PetscFV) obj, &q);
168:       PetscViewerASCIIPrintf(viewer, " FVM");
169:     }
170:     else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", f);
171:     if (Nc > 1) {PetscViewerASCIIPrintf(viewer, "%D components", Nc);}
172:     else        {PetscViewerASCIIPrintf(viewer, "%D component ", Nc);}
173:     if (prob->implicit[f]) {PetscViewerASCIIPrintf(viewer, " (implicit)");}
174:     else                   {PetscViewerASCIIPrintf(viewer, " (explicit)");}
175:     if (q) {
176:       PetscQuadratureGetData(q, NULL, &Nqc, &Nq, NULL, NULL);
177:       PetscViewerASCIIPrintf(viewer, " (Nq %D Nqc %D)", Nq, Nqc);
178:     }
179:     PetscViewerASCIIPrintf(viewer, " %D-jet", prob->jetDegree[f]);
180:     PetscViewerASCIIPrintf(viewer, "\n");
181:     PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
182:     PetscViewerASCIIPushTab(viewer);
183:     if (id == PETSCFE_CLASSID)      {PetscFEView((PetscFE) obj, viewer);}
184:     else if (id == PETSCFV_CLASSID) {PetscFVView((PetscFV) obj, viewer);}
185:     PetscViewerASCIIPopTab(viewer);

187:     for (b = prob->boundary; b; b = b->next) {
188:       PetscInt c, i;

190:       if (b->field != f) continue;
191:       PetscViewerASCIIPushTab(viewer);
192:       PetscViewerASCIIPrintf(viewer, "Boundary %s (%s) %s\n", b->name, b->labelname, DMBoundaryConditionTypes[b->type]);
193:       if (!b->numcomps) {
194:         PetscViewerASCIIPrintf(viewer, "  all components\n");
195:       } else {
196:         PetscViewerASCIIPrintf(viewer, "  components: ");
197:         PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
198:         for (c = 0; c < b->numcomps; ++c) {
199:           if (c > 0) {PetscViewerASCIIPrintf(viewer, ", ");}
200:           PetscViewerASCIIPrintf(viewer, "%D", b->comps[c]);
201:         }
202:         PetscViewerASCIIPrintf(viewer, "\n");
203:         PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
204:       }
205:       PetscViewerASCIIPrintf(viewer, "  ids: ");
206:       PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
207:       for (i = 0; i < b->numids; ++i) {
208:         if (i > 0) {PetscViewerASCIIPrintf(viewer, ", ");}
209:         PetscViewerASCIIPrintf(viewer, "%D", b->ids[i]);
210:       }
211:       PetscViewerASCIIPrintf(viewer, "\n");
212:       PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
213:       PetscViewerASCIIPopTab(viewer);
214:     }
215:   }
216:   PetscDSGetConstants(prob, &numConstants, &constants);
217:   if (numConstants) {
218:     PetscViewerASCIIPrintf(viewer, "%D constants\n", numConstants);
219:     PetscViewerASCIIPushTab(viewer);
220:     for (f = 0; f < numConstants; ++f) {PetscViewerASCIIPrintf(viewer, "%g\n", (double) PetscRealPart(constants[f]));}
221:     PetscViewerASCIIPopTab(viewer);
222:   }
223:   PetscWeakFormView(prob->wf, viewer);
224:   PetscViewerASCIIPopTab(viewer);
225:   return(0);
226: }

228: /*@C
229:    PetscDSViewFromOptions - View from Options

231:    Collective on PetscDS

233:    Input Parameters:
234: +  A - the PetscDS object
235: .  obj - Optional object
236: -  name - command line option

238:    Level: intermediate
239: .seealso:  PetscDS, PetscDSView, PetscObjectViewFromOptions(), PetscDSCreate()
240: @*/
241: PetscErrorCode  PetscDSViewFromOptions(PetscDS A,PetscObject obj,const char name[])
242: {

247:   PetscObjectViewFromOptions((PetscObject)A,obj,name);
248:   return(0);
249: }

251: /*@C
252:   PetscDSView - Views a PetscDS

254:   Collective on prob

256:   Input Parameter:
257: + prob - the PetscDS object to view
258: - v  - the viewer

260:   Level: developer

262: .seealso PetscDSDestroy()
263: @*/
264: PetscErrorCode PetscDSView(PetscDS prob, PetscViewer v)
265: {
266:   PetscBool      iascii;

271:   if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) prob), &v);}
273:   PetscObjectTypeCompare((PetscObject) v, PETSCVIEWERASCII, &iascii);
274:   if (iascii) {PetscDSView_Ascii(prob, v);}
275:   if (prob->ops->view) {(*prob->ops->view)(prob, v);}
276:   return(0);
277: }

279: /*@
280:   PetscDSSetFromOptions - sets parameters in a PetscDS from the options database

282:   Collective on prob

284:   Input Parameter:
285: . prob - the PetscDS object to set options for

287:   Options Database:
288: + -petscds_type <type>     : Set the DS type
289: . -petscds_view <view opt> : View the DS
290: . -petscds_jac_pre         : Turn formation of a separate Jacobian preconditioner on and off
291: . -bc_<name> <ids>         : Specify a list of label ids for a boundary condition
292: - -bc_<name>_comp <comps>  : Specify a list of field components to constrain for a boundary condition

294:   Level: developer

296: .seealso PetscDSView()
297: @*/
298: PetscErrorCode PetscDSSetFromOptions(PetscDS prob)
299: {
300:   DSBoundary     b;
301:   const char    *defaultType;
302:   char           name[256];
303:   PetscBool      flg;

308:   if (!((PetscObject) prob)->type_name) {
309:     defaultType = PETSCDSBASIC;
310:   } else {
311:     defaultType = ((PetscObject) prob)->type_name;
312:   }
313:   PetscDSRegisterAll();

315:   PetscObjectOptionsBegin((PetscObject) prob);
316:   for (b = prob->boundary; b; b = b->next) {
317:     char       optname[1024];
318:     PetscInt   ids[1024], len = 1024;
319:     PetscBool  flg;

321:     PetscSNPrintf(optname, sizeof(optname), "-bc_%s", b->name);
322:     PetscMemzero(ids, sizeof(ids));
323:     PetscOptionsIntArray(optname, "List of boundary IDs", "", ids, &len, &flg);
324:     if (flg) {
325:       b->numids = len;
326:       PetscFree(b->ids);
327:       PetscMalloc1(len, &b->ids);
328:       PetscArraycpy(b->ids, ids, len);
329:     }
330:     len = 1024;
331:     PetscSNPrintf(optname, sizeof(optname), "-bc_%s_comp", b->name);
332:     PetscMemzero(ids, sizeof(ids));
333:     PetscOptionsIntArray(optname, "List of boundary field components", "", ids, &len, &flg);
334:     if (flg) {
335:       b->numcomps = len;
336:       PetscFree(b->comps);
337:       PetscMalloc1(len, &b->comps);
338:       PetscArraycpy(b->comps, ids, len);
339:     }
340:   }
341:   PetscOptionsFList("-petscds_type", "Discrete System", "PetscDSSetType", PetscDSList, defaultType, name, 256, &flg);
342:   if (flg) {
343:     PetscDSSetType(prob, name);
344:   } else if (!((PetscObject) prob)->type_name) {
345:     PetscDSSetType(prob, defaultType);
346:   }
347:   PetscOptionsBool("-petscds_jac_pre", "Discrete System", "PetscDSUseJacobianPreconditioner", prob->useJacPre, &prob->useJacPre, &flg);
348:   if (prob->ops->setfromoptions) {(*prob->ops->setfromoptions)(prob);}
349:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
350:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) prob);
351:   PetscOptionsEnd();
352:   if (prob->Nf) {PetscDSViewFromOptions(prob, NULL, "-petscds_view");}
353:   return(0);
354: }

356: /*@C
357:   PetscDSSetUp - Construct data structures for the PetscDS

359:   Collective on prob

361:   Input Parameter:
362: . prob - the PetscDS object to setup

364:   Level: developer

366: .seealso PetscDSView(), PetscDSDestroy()
367: @*/
368: PetscErrorCode PetscDSSetUp(PetscDS prob)
369: {
370:   const PetscInt Nf   = prob->Nf;
371:   PetscBool      hasH = PETSC_FALSE;
372:   PetscInt       dim, dimEmbed, NbMax = 0, NcMax = 0, NqMax = 0, NsMax = 1, f;

377:   if (prob->setup) return(0);
378:   /* Calculate sizes */
379:   PetscDSGetSpatialDimension(prob, &dim);
380:   PetscDSGetCoordinateDimension(prob, &dimEmbed);
381:   prob->totDim = prob->totComp = 0;
382:   PetscMalloc2(Nf,&prob->Nc,Nf,&prob->Nb);
383:   PetscCalloc2(Nf+1,&prob->off,Nf+1,&prob->offDer);
384:   PetscMalloc2(Nf,&prob->T,Nf,&prob->Tf);
385:   for (f = 0; f < Nf; ++f) {
386:     PetscObject     obj;
387:     PetscClassId    id;
388:     PetscQuadrature q = NULL;
389:     PetscInt        Nq = 0, Nb, Nc;

391:     PetscDSGetDiscretization(prob, f, &obj);
392:     if (prob->jetDegree[f] > 1) hasH = PETSC_TRUE;
393:     if (!obj) {
394:       /* Empty mesh */
395:       Nb = Nc = 0;
396:       prob->T[f] = prob->Tf[f] = NULL;
397:     } else {
398:       PetscObjectGetClassId(obj, &id);
399:       if (id == PETSCFE_CLASSID)      {
400:         PetscFE fe = (PetscFE) obj;

402:         PetscFEGetQuadrature(fe, &q);
403:         PetscFEGetDimension(fe, &Nb);
404:         PetscFEGetNumComponents(fe, &Nc);
405:         PetscFEGetCellTabulation(fe, prob->jetDegree[f], &prob->T[f]);
406:         PetscFEGetFaceTabulation(fe, prob->jetDegree[f], &prob->Tf[f]);
407:       } else if (id == PETSCFV_CLASSID) {
408:         PetscFV fv = (PetscFV) obj;

410:         PetscFVGetQuadrature(fv, &q);
411:         PetscFVGetNumComponents(fv, &Nc);
412:         Nb   = Nc;
413:         PetscFVGetCellTabulation(fv, &prob->T[f]);
414:         /* TODO: should PetscFV also have face tabulation? Otherwise there will be a null pointer in prob->basisFace */
415:       } else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
416:     }
417:     prob->Nc[f]       = Nc;
418:     prob->Nb[f]       = Nb;
419:     prob->off[f+1]    = Nc     + prob->off[f];
420:     prob->offDer[f+1] = Nc*dim + prob->offDer[f];
421:     if (q) {PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);}
422:     NqMax          = PetscMax(NqMax, Nq);
423:     NbMax          = PetscMax(NbMax, Nb);
424:     NcMax          = PetscMax(NcMax, Nc);
425:     prob->totDim  += Nb;
426:     prob->totComp += Nc;
427:     /* There are two faces for all fields but the cohesive field on a hybrid cell */
428:     if (prob->isHybrid && (f < Nf-1)) prob->totDim += Nb;
429:   }
430:   /* Allocate works space */
431:   if (prob->isHybrid) NsMax = 2;
432:   PetscMalloc3(NsMax*prob->totComp,&prob->u,NsMax*prob->totComp,&prob->u_t,NsMax*prob->totComp*dimEmbed + (hasH ? NsMax*prob->totComp*dimEmbed*dimEmbed : 0),&prob->u_x);
433:   PetscMalloc5(dimEmbed,&prob->x,NbMax*NcMax,&prob->basisReal,NbMax*NcMax*dimEmbed,&prob->basisDerReal,NbMax*NcMax,&prob->testReal,NbMax*NcMax*dimEmbed,&prob->testDerReal);
434:   PetscMalloc6(NsMax*NqMax*NcMax,&prob->f0,NsMax*NqMax*NcMax*dimEmbed,&prob->f1,
435:                       NsMax*NsMax*NqMax*NcMax*NcMax,&prob->g0,NsMax*NsMax*NqMax*NcMax*NcMax*dimEmbed,&prob->g1,
436:                       NsMax*NsMax*NqMax*NcMax*NcMax*dimEmbed,&prob->g2,NsMax*NsMax*NqMax*NcMax*NcMax*dimEmbed*dimEmbed,&prob->g3);
437:   if (prob->ops->setup) {(*prob->ops->setup)(prob);}
438:   prob->setup = PETSC_TRUE;
439:   return(0);
440: }

442: static PetscErrorCode PetscDSDestroyStructs_Static(PetscDS prob)
443: {

447:   PetscFree2(prob->Nc,prob->Nb);
448:   PetscFree2(prob->off,prob->offDer);
449:   PetscFree2(prob->T,prob->Tf);
450:   PetscFree3(prob->u,prob->u_t,prob->u_x);
451:   PetscFree5(prob->x,prob->basisReal, prob->basisDerReal,prob->testReal,prob->testDerReal);
452:   PetscFree6(prob->f0,prob->f1,prob->g0,prob->g1,prob->g2,prob->g3);
453:   return(0);
454: }

456: static PetscErrorCode PetscDSEnlarge_Static(PetscDS prob, PetscInt NfNew)
457: {
458:   PetscObject      *tmpd;
459:   PetscBool        *tmpi;
460:   PetscInt         *tmpk;
461:   PetscPointFunc   *tmpup;
462:   PetscSimplePointFunc *tmpexactSol,  *tmpexactSol_t;
463:   void                **tmpexactCtx, **tmpexactCtx_t;
464:   void            **tmpctx;
465:   PetscInt          Nf = prob->Nf, f;
466:   PetscErrorCode    ierr;

469:   if (Nf >= NfNew) return(0);
470:   prob->setup = PETSC_FALSE;
471:   PetscDSDestroyStructs_Static(prob);
472:   PetscMalloc3(NfNew, &tmpd, NfNew, &tmpi, NfNew, &tmpk);
473:   for (f = 0; f < Nf; ++f) {tmpd[f] = prob->disc[f]; tmpi[f] = prob->implicit[f]; tmpk[f] = prob->jetDegree[f];}
474:   for (f = Nf; f < NfNew; ++f) {tmpd[f] = NULL; tmpi[f] = PETSC_TRUE; tmpk[f] = 1;}
475:   PetscFree3(prob->disc, prob->implicit, prob->jetDegree);
476:   PetscWeakFormSetNumFields(prob->wf, NfNew);
477:   prob->Nf        = NfNew;
478:   prob->disc      = tmpd;
479:   prob->implicit  = tmpi;
480:   prob->jetDegree = tmpk;
481:   PetscCalloc2(NfNew, &tmpup, NfNew, &tmpctx);
482:   for (f = 0; f < Nf; ++f) tmpup[f] = prob->update[f];
483:   for (f = 0; f < Nf; ++f) tmpctx[f] = prob->ctx[f];
484:   for (f = Nf; f < NfNew; ++f) tmpup[f] = NULL;
485:   for (f = Nf; f < NfNew; ++f) tmpctx[f] = NULL;
486:   PetscFree2(prob->update, prob->ctx);
487:   prob->update = tmpup;
488:   prob->ctx = tmpctx;
489:   PetscCalloc4(NfNew, &tmpexactSol, NfNew, &tmpexactCtx, NfNew, &tmpexactSol_t, NfNew, &tmpexactCtx_t);
490:   for (f = 0; f < Nf; ++f) tmpexactSol[f] = prob->exactSol[f];
491:   for (f = 0; f < Nf; ++f) tmpexactCtx[f] = prob->exactCtx[f];
492:   for (f = 0; f < Nf; ++f) tmpexactSol_t[f] = prob->exactSol_t[f];
493:   for (f = 0; f < Nf; ++f) tmpexactCtx_t[f] = prob->exactCtx_t[f];
494:   for (f = Nf; f < NfNew; ++f) tmpexactSol[f] = NULL;
495:   for (f = Nf; f < NfNew; ++f) tmpexactCtx[f] = NULL;
496:   for (f = Nf; f < NfNew; ++f) tmpexactSol_t[f] = NULL;
497:   for (f = Nf; f < NfNew; ++f) tmpexactCtx_t[f] = NULL;
498:   PetscFree4(prob->exactSol, prob->exactCtx, prob->exactSol_t, prob->exactCtx_t);
499:   prob->exactSol = tmpexactSol;
500:   prob->exactCtx = tmpexactCtx;
501:   prob->exactSol_t = tmpexactSol_t;
502:   prob->exactCtx_t = tmpexactCtx_t;
503:   return(0);
504: }

506: /*@
507:   PetscDSDestroy - Destroys a PetscDS object

509:   Collective on prob

511:   Input Parameter:
512: . prob - the PetscDS object to destroy

514:   Level: developer

516: .seealso PetscDSView()
517: @*/
518: PetscErrorCode PetscDSDestroy(PetscDS *ds)
519: {
520:   PetscInt       f;
521:   DSBoundary     next;

525:   if (!*ds) return(0);

528:   if (--((PetscObject)(*ds))->refct > 0) {*ds = NULL; return(0);}
529:   ((PetscObject) (*ds))->refct = 0;
530:   if ((*ds)->subprobs) {
531:     PetscInt dim, d;

533:     PetscDSGetSpatialDimension(*ds, &dim);
534:     for (d = 0; d < dim; ++d) {PetscDSDestroy(&(*ds)->subprobs[d]);}
535:   }
536:   PetscFree((*ds)->subprobs);
537:   PetscDSDestroyStructs_Static(*ds);
538:   for (f = 0; f < (*ds)->Nf; ++f) {
539:     PetscObjectDereference((*ds)->disc[f]);
540:   }
541:   PetscFree3((*ds)->disc, (*ds)->implicit, (*ds)->jetDegree);
542:   PetscWeakFormDestroy(&(*ds)->wf);
543:   PetscFree2((*ds)->update,(*ds)->ctx);
544:   PetscFree4((*ds)->exactSol,(*ds)->exactCtx,(*ds)->exactSol_t,(*ds)->exactCtx_t);
545:   if ((*ds)->ops->destroy) {(*(*ds)->ops->destroy)(*ds);}
546:   next = (*ds)->boundary;
547:   while (next) {
548:     DSBoundary b = next;

550:     next = b->next;
551:     PetscFree(b->comps);
552:     PetscFree(b->ids);
553:     PetscFree(b->name);
554:     PetscFree(b->labelname);
555:     PetscFree(b);
556:   }
557:   PetscFree((*ds)->constants);
558:   PetscHeaderDestroy(ds);
559:   return(0);
560: }

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

565:   Collective

567:   Input Parameter:
568: . comm - The communicator for the PetscDS object

570:   Output Parameter:
571: . ds   - The PetscDS object

573:   Level: beginner

575: .seealso: PetscDSSetType(), PETSCDSBASIC
576: @*/
577: PetscErrorCode PetscDSCreate(MPI_Comm comm, PetscDS *ds)
578: {
579:   PetscDS        p;

584:   *ds  = NULL;
585:   PetscDSInitializePackage();

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

589:   p->Nf           = 0;
590:   p->setup        = PETSC_FALSE;
591:   p->numConstants = 0;
592:   p->constants    = NULL;
593:   p->dimEmbed     = -1;
594:   p->useJacPre    = PETSC_TRUE;
595:   PetscWeakFormCreate(comm, &p->wf);

597:   *ds = p;
598:   return(0);
599: }

601: /*@
602:   PetscDSGetNumFields - Returns the number of fields in the DS

604:   Not collective

606:   Input Parameter:
607: . prob - The PetscDS object

609:   Output Parameter:
610: . Nf - The number of fields

612:   Level: beginner

614: .seealso: PetscDSGetSpatialDimension(), PetscDSCreate()
615: @*/
616: PetscErrorCode PetscDSGetNumFields(PetscDS prob, PetscInt *Nf)
617: {
621:   *Nf = prob->Nf;
622:   return(0);
623: }

625: /*@
626:   PetscDSGetSpatialDimension - Returns the spatial dimension of the DS, meaning the topological dimension of the discretizations

628:   Not collective

630:   Input Parameter:
631: . prob - The PetscDS object

633:   Output Parameter:
634: . dim - The spatial dimension

636:   Level: beginner

638: .seealso: PetscDSGetCoordinateDimension(), PetscDSGetNumFields(), PetscDSCreate()
639: @*/
640: PetscErrorCode PetscDSGetSpatialDimension(PetscDS prob, PetscInt *dim)
641: {

647:   *dim = 0;
648:   if (prob->Nf) {
649:     PetscObject  obj;
650:     PetscClassId id;

652:     PetscDSGetDiscretization(prob, 0, &obj);
653:     if (obj) {
654:       PetscObjectGetClassId(obj, &id);
655:       if (id == PETSCFE_CLASSID)      {PetscFEGetSpatialDimension((PetscFE) obj, dim);}
656:       else if (id == PETSCFV_CLASSID) {PetscFVGetSpatialDimension((PetscFV) obj, dim);}
657:       else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", 0);
658:     }
659:   }
660:   return(0);
661: }

663: /*@
664:   PetscDSGetCoordinateDimension - Returns the coordinate dimension of the DS, meaning the dimension of the space into which the discretiaztions are embedded

666:   Not collective

668:   Input Parameter:
669: . prob - The PetscDS object

671:   Output Parameter:
672: . dimEmbed - The coordinate dimension

674:   Level: beginner

676: .seealso: PetscDSSetCoordinateDimension(), PetscDSGetSpatialDimension(), PetscDSGetNumFields(), PetscDSCreate()
677: @*/
678: PetscErrorCode PetscDSGetCoordinateDimension(PetscDS prob, PetscInt *dimEmbed)
679: {
683:   if (prob->dimEmbed < 0) SETERRQ(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONGSTATE, "No coordinate dimension set for this DS");
684:   *dimEmbed = prob->dimEmbed;
685:   return(0);
686: }

688: /*@
689:   PetscDSSetCoordinateDimension - Set the coordinate dimension of the DS, meaning the dimension of the space into which the discretiaztions are embedded

691:   Logically collective on prob

693:   Input Parameters:
694: + prob - The PetscDS object
695: - dimEmbed - The coordinate dimension

697:   Level: beginner

699: .seealso: PetscDSGetCoordinateDimension(), PetscDSGetSpatialDimension(), PetscDSGetNumFields(), PetscDSCreate()
700: @*/
701: PetscErrorCode PetscDSSetCoordinateDimension(PetscDS prob, PetscInt dimEmbed)
702: {
705:   if (dimEmbed < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Coordinate dimension must be non-negative, not %D", dimEmbed);
706:   prob->dimEmbed = dimEmbed;
707:   return(0);
708: }

710: /*@
711:   PetscDSGetHybrid - Returns the flag for a hybrid (cohesive) cell

713:   Not collective

715:   Input Parameter:
716: . prob - The PetscDS object

718:   Output Parameter:
719: . isHybrid - The flag

721:   Level: developer

723: .seealso: PetscDSSetHybrid(), PetscDSCreate()
724: @*/
725: PetscErrorCode PetscDSGetHybrid(PetscDS prob, PetscBool *isHybrid)
726: {
730:   *isHybrid = prob->isHybrid;
731:   return(0);
732: }

734: /*@
735:   PetscDSSetHybrid - Set the flag for a hybrid (cohesive) cell

737:   Not collective

739:   Input Parameters:
740: + prob - The PetscDS object
741: - isHybrid - The flag

743:   Level: developer

745: .seealso: PetscDSGetHybrid(), PetscDSCreate()
746: @*/
747: PetscErrorCode PetscDSSetHybrid(PetscDS prob, PetscBool isHybrid)
748: {
751:   prob->isHybrid = isHybrid;
752:   return(0);
753: }

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

758:   Not collective

760:   Input Parameter:
761: . prob - The PetscDS object

763:   Output Parameter:
764: . dim - The total problem dimension

766:   Level: beginner

768: .seealso: PetscDSGetNumFields(), PetscDSCreate()
769: @*/
770: PetscErrorCode PetscDSGetTotalDimension(PetscDS prob, PetscInt *dim)
771: {

776:   PetscDSSetUp(prob);
778:   *dim = prob->totDim;
779:   return(0);
780: }

782: /*@
783:   PetscDSGetTotalComponents - Returns the total number of components in this system

785:   Not collective

787:   Input Parameter:
788: . prob - The PetscDS object

790:   Output Parameter:
791: . dim - The total number of components

793:   Level: beginner

795: .seealso: PetscDSGetNumFields(), PetscDSCreate()
796: @*/
797: PetscErrorCode PetscDSGetTotalComponents(PetscDS prob, PetscInt *Nc)
798: {

803:   PetscDSSetUp(prob);
805:   *Nc = prob->totComp;
806:   return(0);
807: }

809: /*@
810:   PetscDSGetDiscretization - Returns the discretization object for the given field

812:   Not collective

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

818:   Output Parameter:
819: . disc - The discretization object

821:   Level: beginner

823: .seealso: PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
824: @*/
825: PetscErrorCode PetscDSGetDiscretization(PetscDS prob, PetscInt f, PetscObject *disc)
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:   *disc = prob->disc[f];
832:   return(0);
833: }

835: /*@
836:   PetscDSSetDiscretization - Sets the discretization object for the given field

838:   Not collective

840:   Input Parameters:
841: + prob - The PetscDS object
842: . f - The field number
843: - disc - The discretization object

845:   Level: beginner

847: .seealso: PetscDSGetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
848: @*/
849: PetscErrorCode PetscDSSetDiscretization(PetscDS prob, PetscInt f, PetscObject disc)
850: {

856:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
857:   PetscDSEnlarge_Static(prob, f+1);
858:   PetscObjectDereference(prob->disc[f]);
859:   prob->disc[f] = disc;
860:   PetscObjectReference(disc);
861:   if (disc) {
862:     PetscClassId id;

864:     PetscObjectGetClassId(disc, &id);
865:     if (id == PETSCFE_CLASSID) {
866:       PetscDSSetImplicit(prob, f, PETSC_TRUE);
867:     } else if (id == PETSCFV_CLASSID) {
868:       PetscDSSetImplicit(prob, f, PETSC_FALSE);
869:     }
870:     PetscDSSetJetDegree(prob, f, 1);
871:   }
872:   return(0);
873: }

875: /*@
876:   PetscDSGetWeakForm - Returns the weak form object

878:   Not collective

880:   Input Parameter:
881: . ds - The PetscDS object

883:   Output Parameter:
884: . wf - The weak form object

886:   Level: beginner

888: .seealso: PetscDSSetWeakForm(), PetscDSGetNumFields(), PetscDSCreate()
889: @*/
890: PetscErrorCode PetscDSGetWeakForm(PetscDS ds, PetscWeakForm *wf)
891: {
895:   *wf = ds->wf;
896:   return(0);
897: }

899: /*@
900:   PetscDSSetWeakForm - Sets the weak form object

902:   Not collective

904:   Input Parameters:
905: + ds - The PetscDS object
906: - wf - The weak form object

908:   Level: beginner

910: .seealso: PetscDSGetWeakForm(), PetscDSGetNumFields(), PetscDSCreate()
911: @*/
912: PetscErrorCode PetscDSSetWeakForm(PetscDS ds, PetscWeakForm wf)
913: {

919:   PetscObjectDereference((PetscObject) ds->wf);
920:   ds->wf = wf;
921:   PetscObjectReference((PetscObject) wf);
922:   PetscWeakFormSetNumFields(wf, ds->Nf);
923:   return(0);
924: }

926: /*@
927:   PetscDSAddDiscretization - Adds a discretization object

929:   Not collective

931:   Input Parameters:
932: + prob - The PetscDS object
933: - disc - The boundary discretization object

935:   Level: beginner

937: .seealso: PetscDSGetDiscretization(), PetscDSSetDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
938: @*/
939: PetscErrorCode PetscDSAddDiscretization(PetscDS prob, PetscObject disc)
940: {

944:   PetscDSSetDiscretization(prob, prob->Nf, disc);
945:   return(0);
946: }

948: /*@
949:   PetscDSGetQuadrature - Returns the quadrature, which must agree for all fields in the DS

951:   Not collective

953:   Input Parameter:
954: . prob - The PetscDS object

956:   Output Parameter:
957: . q - The quadrature object

959: Level: intermediate

961: .seealso: PetscDSSetImplicit(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
962: @*/
963: PetscErrorCode PetscDSGetQuadrature(PetscDS prob, PetscQuadrature *q)
964: {
965:   PetscObject    obj;
966:   PetscClassId   id;

970:   *q = NULL;
971:   if (!prob->Nf) return(0);
972:   PetscDSGetDiscretization(prob, 0, &obj);
973:   PetscObjectGetClassId(obj, &id);
974:   if      (id == PETSCFE_CLASSID) {PetscFEGetQuadrature((PetscFE) obj, q);}
975:   else if (id == PETSCFV_CLASSID) {PetscFVGetQuadrature((PetscFV) obj, q);}
976:   else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", 0);
977:   return(0);
978: }

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

983:   Not collective

985:   Input Parameters:
986: + prob - The PetscDS object
987: - f - The field number

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

992:   Level: developer

994: .seealso: PetscDSSetImplicit(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
995: @*/
996: PetscErrorCode PetscDSGetImplicit(PetscDS prob, PetscInt f, PetscBool *implicit)
997: {
1001:   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);
1002:   *implicit = prob->implicit[f];
1003:   return(0);
1004: }

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

1009:   Not collective

1011:   Input Parameters:
1012: + prob - The PetscDS object
1013: . f - The field number
1014: - implicit - The flag indicating what kind of solve to use for this field

1016:   Level: developer

1018: .seealso: PetscDSGetImplicit(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
1019: @*/
1020: PetscErrorCode PetscDSSetImplicit(PetscDS prob, PetscInt f, PetscBool implicit)
1021: {
1024:   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);
1025:   prob->implicit[f] = implicit;
1026:   return(0);
1027: }

1029: /*@
1030:   PetscDSGetJetDegree - Returns the highest derivative for this field equation, or the k-jet that the discretization needs to tabulate.

1032:   Not collective

1034:   Input Parameters:
1035: + ds - The PetscDS object
1036: - f  - The field number

1038:   Output Parameter:
1039: . k  - The highest derivative we need to tabulate

1041:   Level: developer

1043: .seealso: PetscDSSetJetDegree(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
1044: @*/
1045: PetscErrorCode PetscDSGetJetDegree(PetscDS ds, PetscInt f, PetscInt *k)
1046: {
1050:   if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
1051:   *k = ds->jetDegree[f];
1052:   return(0);
1053: }

1055: /*@
1056:   PetscDSSetJetDegree - Set the highest derivative for this field equation, or the k-jet that the discretization needs to tabulate.

1058:   Not collective

1060:   Input Parameters:
1061: + ds - The PetscDS object
1062: . f  - The field number
1063: - k  - The highest derivative we need to tabulate

1065:   Level: developer

1067: .seealso: PetscDSGetJetDegree(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
1068: @*/
1069: PetscErrorCode PetscDSSetJetDegree(PetscDS ds, PetscInt f, PetscInt k)
1070: {
1073:   if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
1074:   ds->jetDegree[f] = k;
1075:   return(0);
1076: }

1078: PetscErrorCode PetscDSGetObjective(PetscDS ds, PetscInt f,
1079:                                    void (**obj)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1080:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1081:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1082:                                                 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[]))
1083: {
1084:   PetscPointFunc *tmp;
1085:   PetscInt        n;
1086:   PetscErrorCode  ierr;

1091:   if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
1092:   PetscWeakFormGetObjective(ds->wf, NULL, 0, f, &n, &tmp);
1093:   *obj = tmp ? tmp[0] : NULL;
1094:   return(0);
1095: }

1097: PetscErrorCode PetscDSSetObjective(PetscDS ds, PetscInt f,
1098:                                    void (*obj)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1099:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1100:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1101:                                                PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[]))
1102: {

1108:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1109:   PetscWeakFormSetIndexObjective(ds->wf, NULL, 0, f, 0, obj);
1110:   return(0);
1111: }

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

1116:   Not collective

1118:   Input Parameters:
1119: + ds - The PetscDS
1120: - f  - The test field number

1122:   Output Parameters:
1123: + f0 - integrand for the test function term
1124: - f1 - integrand for the test function gradient term

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

1128:   \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)

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

1132: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1133: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1134: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1135: $    PetscReal t, const PetscReal x[], PetscScalar f0[])

1137: + dim - the spatial dimension
1138: . Nf - the number of fields
1139: . uOff - the offset into u[] and u_t[] for each field
1140: . uOff_x - the offset into u_x[] for each field
1141: . u - each field evaluated at the current point
1142: . u_t - the time derivative of each field evaluated at the current point
1143: . u_x - the gradient of each field evaluated at the current point
1144: . aOff - the offset into a[] and a_t[] for each auxiliary field
1145: . aOff_x - the offset into a_x[] for each auxiliary field
1146: . a - each auxiliary field evaluated at the current point
1147: . a_t - the time derivative of each auxiliary field evaluated at the current point
1148: . a_x - the gradient of auxiliary each field evaluated at the current point
1149: . t - current time
1150: . x - coordinates of the current point
1151: . numConstants - number of constant parameters
1152: . constants - constant parameters
1153: - f0 - output values at the current point

1155:   Level: intermediate

1157: .seealso: PetscDSSetResidual()
1158: @*/
1159: PetscErrorCode PetscDSGetResidual(PetscDS ds, PetscInt f,
1160:                                   void (**f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1161:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1162:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1163:                                               PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]),
1164:                                   void (**f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1165:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1166:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1167:                                               PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
1168: {
1169:   PetscPointFunc *tmp0, *tmp1;
1170:   PetscInt        n0, n1;
1171:   PetscErrorCode  ierr;

1175:   if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
1176:   PetscWeakFormGetResidual(ds->wf, NULL, 0, f, &n0, &tmp0, &n1, &tmp1);
1177:   *f0  = tmp0 ? tmp0[0] : NULL;
1178:   *f1  = tmp1 ? tmp1[0] : NULL;
1179:   return(0);
1180: }

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

1185:   Not collective

1187:   Input Parameters:
1188: + ds - The PetscDS
1189: . f  - The test field number
1190: . f0 - integrand for the test function term
1191: - f1 - integrand for the test function gradient term

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

1195:   \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)

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

1199: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1200: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1201: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1202: $    PetscReal t, const PetscReal x[], PetscScalar f0[])

1204: + dim - the spatial dimension
1205: . Nf - the number of fields
1206: . uOff - the offset into u[] and u_t[] for each field
1207: . uOff_x - the offset into u_x[] for each field
1208: . u - each field evaluated at the current point
1209: . u_t - the time derivative of each field evaluated at the current point
1210: . u_x - the gradient of each field evaluated at the current point
1211: . aOff - the offset into a[] and a_t[] for each auxiliary field
1212: . aOff_x - the offset into a_x[] for each auxiliary field
1213: . a - each auxiliary field evaluated at the current point
1214: . a_t - the time derivative of each auxiliary field evaluated at the current point
1215: . a_x - the gradient of auxiliary each field evaluated at the current point
1216: . t - current time
1217: . x - coordinates of the current point
1218: . numConstants - number of constant parameters
1219: . constants - constant parameters
1220: - f0 - output values at the current point

1222:   Level: intermediate

1224: .seealso: PetscDSGetResidual()
1225: @*/
1226: PetscErrorCode PetscDSSetResidual(PetscDS ds, PetscInt f,
1227:                                   void (*f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1228:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1229:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1230:                                              PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]),
1231:                                   void (*f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1232:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1233:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1234:                                              PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
1235: {

1242:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1243:   PetscWeakFormSetIndexResidual(ds->wf, NULL, 0, f, 0, f0, 0, f1);
1244:   return(0);
1245: }

1247: /*@C
1248:   PetscDSHasJacobian - Signals that Jacobian functions have been set

1250:   Not collective

1252:   Input Parameter:
1253: . prob - The PetscDS

1255:   Output Parameter:
1256: . hasJac - flag that pointwise function for the Jacobian has been set

1258:   Level: intermediate

1260: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1261: @*/
1262: PetscErrorCode PetscDSHasJacobian(PetscDS ds, PetscBool *hasJac)
1263: {

1268:   PetscWeakFormHasJacobian(ds->wf, hasJac);
1269:   return(0);
1270: }

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

1275:   Not collective

1277:   Input Parameters:
1278: + ds - The PetscDS
1279: . f  - The test field number
1280: - g  - The field number

1282:   Output Parameters:
1283: + g0 - integrand for the test and basis function term
1284: . g1 - integrand for the test function and basis function gradient term
1285: . g2 - integrand for the test function gradient and basis function term
1286: - g3 - integrand for the test function gradient and basis function gradient term

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

1290:   \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

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

1294: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1295: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1296: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1297: $    PetscReal t, const PetscReal u_tShift, const PetscReal x[], PetscScalar g0[])

1299: + dim - the spatial dimension
1300: . Nf - the number of fields
1301: . uOff - the offset into u[] and u_t[] for each field
1302: . uOff_x - the offset into u_x[] for each field
1303: . u - each field evaluated at the current point
1304: . u_t - the time derivative of each field evaluated at the current point
1305: . u_x - the gradient of each field evaluated at the current point
1306: . aOff - the offset into a[] and a_t[] for each auxiliary field
1307: . aOff_x - the offset into a_x[] for each auxiliary field
1308: . a - each auxiliary field evaluated at the current point
1309: . a_t - the time derivative of each auxiliary field evaluated at the current point
1310: . a_x - the gradient of auxiliary each field evaluated at the current point
1311: . t - current time
1312: . u_tShift - the multiplier a for dF/dU_t
1313: . x - coordinates of the current point
1314: . numConstants - number of constant parameters
1315: . constants - constant parameters
1316: - g0 - output values at the current point

1318:   Level: intermediate

1320: .seealso: PetscDSSetJacobian()
1321: @*/
1322: PetscErrorCode PetscDSGetJacobian(PetscDS ds, PetscInt f, PetscInt g,
1323:                                   void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1324:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1325:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1326:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1327:                                   void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1328:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1329:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1330:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1331:                                   void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1332:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1333:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1334:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1335:                                   void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1336:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1337:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1338:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1339: {
1340:   PetscPointJac *tmp0, *tmp1, *tmp2, *tmp3;
1341:   PetscInt       n0, n1, n2, n3;

1346:   if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
1347:   if ((g < 0) || (g >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, ds->Nf);
1348:   PetscWeakFormGetJacobian(ds->wf, NULL, 0, f, g, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3);
1349:   *g0  = tmp0 ? tmp0[0] : NULL;
1350:   *g1  = tmp1 ? tmp1[0] : NULL;
1351:   *g2  = tmp2 ? tmp2[0] : NULL;
1352:   *g3  = tmp3 ? tmp3[0] : NULL;
1353:   return(0);
1354: }

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

1359:   Not collective

1361:   Input Parameters:
1362: + ds - The PetscDS
1363: . f  - The test field number
1364: . g  - The field number
1365: . g0 - integrand for the test and basis function term
1366: . g1 - integrand for the test function and basis function gradient term
1367: . g2 - integrand for the test function gradient and basis function term
1368: - g3 - integrand for the test function gradient and basis function gradient term

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

1372:   \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

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

1376: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1377: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1378: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1379: $    PetscReal t, const PetscReal x[], PetscScalar g0[])

1381: + dim - the spatial dimension
1382: . Nf - the number of fields
1383: . uOff - the offset into u[] and u_t[] for each field
1384: . uOff_x - the offset into u_x[] for each field
1385: . u - each field evaluated at the current point
1386: . u_t - the time derivative of each field evaluated at the current point
1387: . u_x - the gradient of each field evaluated at the current point
1388: . aOff - the offset into a[] and a_t[] for each auxiliary field
1389: . aOff_x - the offset into a_x[] for each auxiliary field
1390: . a - each auxiliary field evaluated at the current point
1391: . a_t - the time derivative of each auxiliary field evaluated at the current point
1392: . a_x - the gradient of auxiliary each field evaluated at the current point
1393: . t - current time
1394: . u_tShift - the multiplier a for dF/dU_t
1395: . x - coordinates of the current point
1396: . numConstants - number of constant parameters
1397: . constants - constant parameters
1398: - g0 - output values at the current point

1400:   Level: intermediate

1402: .seealso: PetscDSGetJacobian()
1403: @*/
1404: PetscErrorCode PetscDSSetJacobian(PetscDS ds, PetscInt f, PetscInt g,
1405:                                   void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1406:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1407:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1408:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1409:                                   void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1410:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1411:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1412:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1413:                                   void (*g2)(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[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1417:                                   void (*g3)(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[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1421: {

1430:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1431:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
1432:   PetscWeakFormSetIndexJacobian(ds->wf, NULL, 0, f, g, 0, g0, 0, g1, 0, g2, 0, g3);
1433:   return(0);
1434: }

1436: /*@C
1437:   PetscDSUseJacobianPreconditioner - Whether to construct a Jacobian preconditioner

1439:   Not collective

1441:   Input Parameters:
1442: + prob - The PetscDS
1443: - useJacPre - flag that enables construction of a Jacobian preconditioner

1445:   Level: intermediate

1447: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1448: @*/
1449: PetscErrorCode PetscDSUseJacobianPreconditioner(PetscDS prob, PetscBool useJacPre)
1450: {
1453:   prob->useJacPre = useJacPre;
1454:   return(0);
1455: }

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

1460:   Not collective

1462:   Input Parameter:
1463: . prob - The PetscDS

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

1468:   Level: intermediate

1470: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1471: @*/
1472: PetscErrorCode PetscDSHasJacobianPreconditioner(PetscDS ds, PetscBool *hasJacPre)
1473: {

1478:   *hasJacPre = PETSC_FALSE;
1479:   if (!ds->useJacPre) return(0);
1480:   PetscWeakFormHasJacobianPreconditioner(ds->wf, hasJacPre);
1481:   return(0);
1482: }

1484: /*@C
1485:   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.

1487:   Not collective

1489:   Input Parameters:
1490: + ds - 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: . numConstants - number of constant parameters
1527: . constants - constant parameters
1528: - g0 - output values at the current point

1530:   Level: intermediate

1532: .seealso: PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1533: @*/
1534: PetscErrorCode PetscDSGetJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g,
1535:                                   void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1536:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1537:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1538:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1539:                                   void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1540:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1541:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1542:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1543:                                   void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1544:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1545:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1546:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1547:                                   void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1548:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1549:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1550:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1551: {
1552:   PetscPointJac *tmp0, *tmp1, *tmp2, *tmp3;
1553:   PetscInt       n0, n1, n2, n3;

1558:   if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
1559:   if ((g < 0) || (g >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, ds->Nf);
1560:   PetscWeakFormGetJacobianPreconditioner(ds->wf, NULL, 0, f, g, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3);
1561:   *g0  = tmp0 ? tmp0[0] : NULL;
1562:   *g1  = tmp1 ? tmp1[0] : NULL;
1563:   *g2  = tmp2 ? tmp2[0] : NULL;
1564:   *g3  = tmp3 ? tmp3[0] : NULL;
1565:   return(0);
1566: }

1568: /*@C
1569:   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.

1571:   Not collective

1573:   Input Parameters:
1574: + ds - The PetscDS
1575: . f  - The test field number
1576: . g  - The field number
1577: . g0 - integrand for the test and basis function term
1578: . g1 - integrand for the test function and basis function gradient term
1579: . g2 - integrand for the test function gradient and basis function term
1580: - g3 - integrand for the test function gradient and basis function gradient term

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

1584:   \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

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

1588: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1589: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1590: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1591: $    PetscReal t, const PetscReal x[], PetscScalar g0[])

1593: + dim - the spatial dimension
1594: . Nf - the number of fields
1595: . uOff - the offset into u[] and u_t[] for each field
1596: . uOff_x - the offset into u_x[] for each field
1597: . u - each field evaluated at the current point
1598: . u_t - the time derivative of each field evaluated at the current point
1599: . u_x - the gradient of each field evaluated at the current point
1600: . aOff - the offset into a[] and a_t[] for each auxiliary field
1601: . aOff_x - the offset into a_x[] for each auxiliary field
1602: . a - each auxiliary field evaluated at the current point
1603: . a_t - the time derivative of each auxiliary field evaluated at the current point
1604: . a_x - the gradient of auxiliary each field evaluated at the current point
1605: . t - current time
1606: . u_tShift - the multiplier a for dF/dU_t
1607: . x - coordinates of the current point
1608: . numConstants - number of constant parameters
1609: . constants - constant parameters
1610: - g0 - output values at the current point

1612:   Level: intermediate

1614: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobian()
1615: @*/
1616: PetscErrorCode PetscDSSetJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g,
1617:                                   void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1618:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1619:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1620:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1621:                                   void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1622:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1623:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1624:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1625:                                   void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1626:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1627:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1628:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1629:                                   void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1630:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1631:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1632:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1633: {

1642:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1643:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
1644:   PetscWeakFormSetIndexJacobianPreconditioner(ds->wf, NULL, 0, f, g, 0, g0, 0, g1, 0, g2, 0, g3);
1645:   return(0);
1646: }

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

1651:   Not collective

1653:   Input Parameter:
1654: . ds - The PetscDS

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

1659:   Level: intermediate

1661: .seealso: PetscDSGetDynamicJacobian(), PetscDSSetDynamicJacobian(), PetscDSGetJacobian()
1662: @*/
1663: PetscErrorCode PetscDSHasDynamicJacobian(PetscDS ds, PetscBool *hasDynJac)
1664: {

1669:   PetscWeakFormHasDynamicJacobian(ds->wf, hasDynJac);
1670:   return(0);
1671: }

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

1676:   Not collective

1678:   Input Parameters:
1679: + ds - The PetscDS
1680: . f  - The test field number
1681: - g  - The field number

1683:   Output Parameters:
1684: + g0 - integrand for the test and basis function term
1685: . g1 - integrand for the test function and basis function gradient term
1686: . g2 - integrand for the test function gradient and basis function term
1687: - g3 - integrand for the test function gradient and basis function gradient term

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

1691:   \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

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

1695: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1696: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1697: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1698: $    PetscReal t, const PetscReal u_tShift, const PetscReal x[], PetscScalar g0[])

1700: + dim - the spatial dimension
1701: . Nf - the number of fields
1702: . uOff - the offset into u[] and u_t[] for each field
1703: . uOff_x - the offset into u_x[] for each field
1704: . u - each field evaluated at the current point
1705: . u_t - the time derivative of each field evaluated at the current point
1706: . u_x - the gradient of each field evaluated at the current point
1707: . aOff - the offset into a[] and a_t[] for each auxiliary field
1708: . aOff_x - the offset into a_x[] for each auxiliary field
1709: . a - each auxiliary field evaluated at the current point
1710: . a_t - the time derivative of each auxiliary field evaluated at the current point
1711: . a_x - the gradient of auxiliary each field evaluated at the current point
1712: . t - current time
1713: . u_tShift - the multiplier a for dF/dU_t
1714: . x - coordinates of the current point
1715: . numConstants - number of constant parameters
1716: . constants - constant parameters
1717: - g0 - output values at the current point

1719:   Level: intermediate

1721: .seealso: PetscDSSetJacobian()
1722: @*/
1723: PetscErrorCode PetscDSGetDynamicJacobian(PetscDS ds, PetscInt f, PetscInt g,
1724:                                          void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1725:                                                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1726:                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1727:                                                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1728:                                          void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1729:                                                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1730:                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1731:                                                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1732:                                          void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1733:                                                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1734:                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1735:                                                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1736:                                          void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1737:                                                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1738:                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1739:                                                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1740: {
1741:   PetscPointJac *tmp0, *tmp1, *tmp2, *tmp3;
1742:   PetscInt       n0, n1, n2, n3;

1747:   if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
1748:   if ((g < 0) || (g >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, ds->Nf);
1749:   PetscWeakFormGetDynamicJacobian(ds->wf, NULL, 0, f, g, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3);
1750:   *g0  = tmp0 ? tmp0[0] : NULL;
1751:   *g1  = tmp1 ? tmp1[0] : NULL;
1752:   *g2  = tmp2 ? tmp2[0] : NULL;
1753:   *g3  = tmp3 ? tmp3[0] : NULL;
1754:   return(0);
1755: }

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

1760:   Not collective

1762:   Input Parameters:
1763: + ds - The PetscDS
1764: . f  - The test field number
1765: . g  - The field number
1766: . g0 - integrand for the test and basis function term
1767: . g1 - integrand for the test function and basis function gradient term
1768: . g2 - integrand for the test function gradient and basis function term
1769: - g3 - integrand for the test function gradient and basis function gradient term

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

1773:   \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

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

1777: $ g0(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[], PetscScalar g0[])

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: . u_tShift - the multiplier a for dF/dU_t
1796: . x - coordinates of the current point
1797: . numConstants - number of constant parameters
1798: . constants - constant parameters
1799: - g0 - output values at the current point

1801:   Level: intermediate

1803: .seealso: PetscDSGetJacobian()
1804: @*/
1805: PetscErrorCode PetscDSSetDynamicJacobian(PetscDS ds, PetscInt f, PetscInt g,
1806:                                          void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1807:                                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1808:                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1809:                                                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1810:                                          void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1811:                                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1812:                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1813:                                                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1814:                                          void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1815:                                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1816:                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1817:                                                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1818:                                          void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1819:                                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1820:                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1821:                                                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1822: {

1831:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1832:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
1833:   PetscWeakFormSetIndexDynamicJacobian(ds->wf, NULL, 0, f, g, 0, g0, 0, g1, 0, g2, 0, g3);
1834:   return(0);
1835: }

1837: /*@C
1838:   PetscDSGetRiemannSolver - Returns the Riemann solver for the given field

1840:   Not collective

1842:   Input Arguments:
1843: + ds - The PetscDS object
1844: - f  - The field number

1846:   Output Argument:
1847: . r    - Riemann solver

1849:   Calling sequence for r:

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

1853: + dim  - The spatial dimension
1854: . Nf   - The number of fields
1855: . x    - The coordinates at a point on the interface
1856: . n    - The normal vector to the interface
1857: . uL   - The state vector to the left of the interface
1858: . uR   - The state vector to the right of the interface
1859: . flux - output array of flux through the interface
1860: . numConstants - number of constant parameters
1861: . constants - constant parameters
1862: - ctx  - optional user context

1864:   Level: intermediate

1866: .seealso: PetscDSSetRiemannSolver()
1867: @*/
1868: PetscErrorCode PetscDSGetRiemannSolver(PetscDS ds, PetscInt f,
1869:                                        void (**r)(PetscInt dim, PetscInt Nf, const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscInt numConstants, const PetscScalar constants[], PetscScalar flux[], void *ctx))
1870: {
1871:   PetscRiemannFunc *tmp;
1872:   PetscInt          n;
1873:   PetscErrorCode    ierr;

1878:   if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
1879:   PetscWeakFormGetRiemannSolver(ds->wf, NULL, 0, f, &n, &tmp);
1880:   *r   = tmp ? tmp[0] : NULL;
1881:   return(0);
1882: }

1884: /*@C
1885:   PetscDSSetRiemannSolver - Sets the Riemann solver for the given field

1887:   Not collective

1889:   Input Arguments:
1890: + ds - The PetscDS object
1891: . f  - The field number
1892: - r  - Riemann solver

1894:   Calling sequence for r:

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

1898: + dim  - The spatial dimension
1899: . Nf   - The number of fields
1900: . x    - The coordinates at a point on the interface
1901: . n    - The normal vector to the interface
1902: . uL   - The state vector to the left of the interface
1903: . uR   - The state vector to the right of the interface
1904: . flux - output array of flux through the interface
1905: . numConstants - number of constant parameters
1906: . constants - constant parameters
1907: - ctx  - optional user context

1909:   Level: intermediate

1911: .seealso: PetscDSGetRiemannSolver()
1912: @*/
1913: PetscErrorCode PetscDSSetRiemannSolver(PetscDS ds, PetscInt f,
1914:                                        void (*r)(PetscInt dim, PetscInt Nf, const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscInt numConstants, const PetscScalar constants[], PetscScalar flux[], void *ctx))
1915: {

1921:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1922:   PetscWeakFormSetIndexRiemannSolver(ds->wf, NULL, 0, f, 0, r);
1923:   return(0);
1924: }

1926: /*@C
1927:   PetscDSGetUpdate - Get the pointwise update function for a given field

1929:   Not collective

1931:   Input Parameters:
1932: + ds - The PetscDS
1933: - f  - The field number

1935:   Output Parameters:
1936: . update - update function

1938:   Note: The calling sequence for the callback update is given by:

1940: $ update(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1941: $        const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1942: $        const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1943: $        PetscReal t, const PetscReal x[], PetscScalar uNew[])

1945: + dim - the spatial dimension
1946: . Nf - the number of fields
1947: . uOff - the offset into u[] and u_t[] for each field
1948: . uOff_x - the offset into u_x[] for each field
1949: . u - each field evaluated at the current point
1950: . u_t - the time derivative of each field evaluated at the current point
1951: . u_x - the gradient of each field evaluated at the current point
1952: . aOff - the offset into a[] and a_t[] for each auxiliary field
1953: . aOff_x - the offset into a_x[] for each auxiliary field
1954: . a - each auxiliary field evaluated at the current point
1955: . a_t - the time derivative of each auxiliary field evaluated at the current point
1956: . a_x - the gradient of auxiliary each field evaluated at the current point
1957: . t - current time
1958: . x - coordinates of the current point
1959: - uNew - new value for field at the current point

1961:   Level: intermediate

1963: .seealso: PetscDSSetUpdate(), PetscDSSetResidual()
1964: @*/
1965: PetscErrorCode PetscDSGetUpdate(PetscDS ds, PetscInt f,
1966:                                   void (**update)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1967:                                                   const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1968:                                                   const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1969:                                                   PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar uNew[]))
1970: {
1973:   if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
1975:   return(0);
1976: }

1978: /*@C
1979:   PetscDSSetUpdate - Set the pointwise update function for a given field

1981:   Not collective

1983:   Input Parameters:
1984: + ds     - The PetscDS
1985: . f      - The field number
1986: - update - update function

1988:   Note: The calling sequence for the callback update is given by:

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

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

2011:   Level: intermediate

2013: .seealso: PetscDSGetResidual()
2014: @*/
2015: PetscErrorCode PetscDSSetUpdate(PetscDS ds, PetscInt f,
2016:                                 void (*update)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2017:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2018:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2019:                                                PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar uNew[]))
2020: {

2026:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2027:   PetscDSEnlarge_Static(ds, f+1);
2028:   ds->update[f] = update;
2029:   return(0);
2030: }

2032: PetscErrorCode PetscDSGetContext(PetscDS ds, PetscInt f, void **ctx)
2033: {
2036:   if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
2038:   *ctx = ds->ctx[f];
2039:   return(0);
2040: }

2042: PetscErrorCode PetscDSSetContext(PetscDS ds, PetscInt f, void *ctx)
2043: {

2048:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2049:   PetscDSEnlarge_Static(ds, f+1);
2050:   ds->ctx[f] = ctx;
2051:   return(0);
2052: }

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

2057:   Not collective

2059:   Input Parameters:
2060: + ds - The PetscDS
2061: - f  - The test field number

2063:   Output Parameters:
2064: + f0 - boundary integrand for the test function term
2065: - f1 - boundary integrand for the test function gradient term

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

2069:   \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

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

2073: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2074: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2075: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2076: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[])

2078: + dim - the spatial dimension
2079: . Nf - the number of fields
2080: . uOff - the offset into u[] and u_t[] for each field
2081: . uOff_x - the offset into u_x[] for each field
2082: . u - each field evaluated at the current point
2083: . u_t - the time derivative of each field evaluated at the current point
2084: . u_x - the gradient of each field evaluated at the current point
2085: . aOff - the offset into a[] and a_t[] for each auxiliary field
2086: . aOff_x - the offset into a_x[] for each auxiliary field
2087: . a - each auxiliary field evaluated at the current point
2088: . a_t - the time derivative of each auxiliary field evaluated at the current point
2089: . a_x - the gradient of auxiliary each field evaluated at the current point
2090: . t - current time
2091: . x - coordinates of the current point
2092: . n - unit normal at the current point
2093: . numConstants - number of constant parameters
2094: . constants - constant parameters
2095: - f0 - output values at the current point

2097:   Level: intermediate

2099: .seealso: PetscDSSetBdResidual()
2100: @*/
2101: PetscErrorCode PetscDSGetBdResidual(PetscDS ds, PetscInt f,
2102:                                     void (**f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2103:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2104:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2105:                                                 PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]),
2106:                                     void (**f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2107:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2108:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2109:                                                 PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
2110: {
2111:   PetscBdPointFunc *tmp0, *tmp1;
2112:   PetscInt          n0, n1;
2113:   PetscErrorCode    ierr;

2117:   if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
2118:   PetscWeakFormGetBdResidual(ds->wf, NULL, 0, f, &n0, &tmp0, &n1, &tmp1);
2119:   *f0  = tmp0 ? tmp0[0] : NULL;
2120:   *f1  = tmp1 ? tmp1[0] : NULL;
2121:   return(0);
2122: }

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

2127:   Not collective

2129:   Input Parameters:
2130: + ds - The PetscDS
2131: . f  - The test field number
2132: . f0 - boundary integrand for the test function term
2133: - f1 - boundary integrand for the test function gradient term

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

2137:   \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

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

2141: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2142: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2143: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2144: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[])

2146: + dim - the spatial dimension
2147: . Nf - the number of fields
2148: . uOff - the offset into u[] and u_t[] for each field
2149: . uOff_x - the offset into u_x[] for each field
2150: . u - each field evaluated at the current point
2151: . u_t - the time derivative of each field evaluated at the current point
2152: . u_x - the gradient of each field evaluated at the current point
2153: . aOff - the offset into a[] and a_t[] for each auxiliary field
2154: . aOff_x - the offset into a_x[] for each auxiliary field
2155: . a - each auxiliary field evaluated at the current point
2156: . a_t - the time derivative of each auxiliary field evaluated at the current point
2157: . a_x - the gradient of auxiliary each field evaluated at the current point
2158: . t - current time
2159: . x - coordinates of the current point
2160: . n - unit normal at the current point
2161: . numConstants - number of constant parameters
2162: . constants - constant parameters
2163: - f0 - output values at the current point

2165:   Level: intermediate

2167: .seealso: PetscDSGetBdResidual()
2168: @*/
2169: PetscErrorCode PetscDSSetBdResidual(PetscDS ds, PetscInt f,
2170:                                     void (*f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2171:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2172:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2173:                                                PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]),
2174:                                     void (*f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2175:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2176:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2177:                                                PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
2178: {

2183:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2184:   PetscWeakFormSetIndexBdResidual(ds->wf, NULL, 0, f, 0, f0, 0, f1);
2185:   return(0);
2186: }

2188: /*@
2189:   PetscDSHasBdJacobian - Signals that boundary Jacobian functions have been set

2191:   Not collective

2193:   Input Parameter:
2194: . ds - The PetscDS

2196:   Output Parameter:
2197: . hasBdJac - flag that pointwise function for the boundary Jacobian has been set

2199:   Level: intermediate

2201: .seealso: PetscDSHasJacobian(), PetscDSSetBdJacobian(), PetscDSGetBdJacobian()
2202: @*/
2203: PetscErrorCode PetscDSHasBdJacobian(PetscDS ds, PetscBool *hasBdJac)
2204: {

2210:   PetscWeakFormHasBdJacobian(ds->wf, hasBdJac);
2211:   return(0);
2212: }

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

2217:   Not collective

2219:   Input Parameters:
2220: + ds - The PetscDS
2221: . f  - The test field number
2222: - g  - The field number

2224:   Output Parameters:
2225: + g0 - integrand for the test and basis function term
2226: . g1 - integrand for the test function and basis function gradient term
2227: . g2 - integrand for the test function gradient and basis function term
2228: - g3 - integrand for the test function gradient and basis function gradient term

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

2232:   \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

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

2236: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2237: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2238: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2239: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar g0[])

2241: + dim - the spatial dimension
2242: . Nf - the number of fields
2243: . uOff - the offset into u[] and u_t[] for each field
2244: . uOff_x - the offset into u_x[] for each field
2245: . u - each field evaluated at the current point
2246: . u_t - the time derivative of each field evaluated at the current point
2247: . u_x - the gradient of each field evaluated at the current point
2248: . aOff - the offset into a[] and a_t[] for each auxiliary field
2249: . aOff_x - the offset into a_x[] for each auxiliary field
2250: . a - each auxiliary field evaluated at the current point
2251: . a_t - the time derivative of each auxiliary field evaluated at the current point
2252: . a_x - the gradient of auxiliary each field evaluated at the current point
2253: . t - current time
2254: . u_tShift - the multiplier a for dF/dU_t
2255: . x - coordinates of the current point
2256: . n - normal at the current point
2257: . numConstants - number of constant parameters
2258: . constants - constant parameters
2259: - g0 - output values at the current point

2261:   Level: intermediate

2263: .seealso: PetscDSSetBdJacobian()
2264: @*/
2265: PetscErrorCode PetscDSGetBdJacobian(PetscDS ds, PetscInt f, PetscInt g,
2266:                                     void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2267:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2268:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2269:                                                 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
2270:                                     void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2271:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2272:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2273:                                                 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
2274:                                     void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2275:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2276:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2277:                                                 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
2278:                                     void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2279:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2280:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2281:                                                 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
2282: {
2283:   PetscBdPointJac *tmp0, *tmp1, *tmp2, *tmp3;
2284:   PetscInt         n0, n1, n2, n3;
2285:   PetscErrorCode   ierr;

2289:   if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
2290:   if ((g < 0) || (g >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, ds->Nf);
2291:   PetscWeakFormGetBdJacobian(ds->wf, NULL, 0, f, g, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3);
2292:   *g0  = tmp0 ? tmp0[0] : NULL;
2293:   *g1  = tmp1 ? tmp1[0] : NULL;
2294:   *g2  = tmp2 ? tmp2[0] : NULL;
2295:   *g3  = tmp3 ? tmp3[0] : NULL;
2296:   return(0);
2297: }

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

2302:   Not collective

2304:   Input Parameters:
2305: + ds - The PetscDS
2306: . f  - The test field number
2307: . g  - The field number
2308: . g0 - integrand for the test and basis function term
2309: . g1 - integrand for the test function and basis function gradient term
2310: . g2 - integrand for the test function gradient and basis function term
2311: - g3 - integrand for the test function gradient and basis function gradient term

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

2315:   \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

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

2319: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2320: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2321: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2322: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar g0[])

2324: + dim - the spatial dimension
2325: . Nf - the number of fields
2326: . uOff - the offset into u[] and u_t[] for each field
2327: . uOff_x - the offset into u_x[] for each field
2328: . u - each field evaluated at the current point
2329: . u_t - the time derivative of each field evaluated at the current point
2330: . u_x - the gradient of each field evaluated at the current point
2331: . aOff - the offset into a[] and a_t[] for each auxiliary field
2332: . aOff_x - the offset into a_x[] for each auxiliary field
2333: . a - each auxiliary field evaluated at the current point
2334: . a_t - the time derivative of each auxiliary field evaluated at the current point
2335: . a_x - the gradient of auxiliary each field evaluated at the current point
2336: . t - current time
2337: . u_tShift - the multiplier a for dF/dU_t
2338: . x - coordinates of the current point
2339: . n - normal at the current point
2340: . numConstants - number of constant parameters
2341: . constants - constant parameters
2342: - g0 - output values at the current point

2344:   Level: intermediate

2346: .seealso: PetscDSGetBdJacobian()
2347: @*/
2348: PetscErrorCode PetscDSSetBdJacobian(PetscDS ds, PetscInt f, PetscInt g,
2349:                                     void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2350:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2351:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2352:                                                PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
2353:                                     void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2354:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2355:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2356:                                                PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
2357:                                     void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2358:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2359:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2360:                                                PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
2361:                                     void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2362:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2363:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2364:                                                PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
2365: {

2374:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2375:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
2376:   PetscWeakFormSetIndexBdJacobian(ds->wf, NULL, 0, f, g, 0, g0, 0, g1, 0, g2, 0, g3);
2377:   return(0);
2378: }

2380: /*@
2381:   PetscDSHasBdJacobianPreconditioner - Signals that boundary Jacobian preconditioner functions have been set

2383:   Not collective

2385:   Input Parameter:
2386: . ds - The PetscDS

2388:   Output Parameter:
2389: . hasBdJac - flag that pointwise function for the boundary Jacobian preconditioner has been set

2391:   Level: intermediate

2393: .seealso: PetscDSHasJacobian(), PetscDSSetBdJacobian(), PetscDSGetBdJacobian()
2394: @*/
2395: PetscErrorCode PetscDSHasBdJacobianPreconditioner(PetscDS ds, PetscBool *hasBdJacPre)
2396: {

2402:   PetscWeakFormHasBdJacobianPreconditioner(ds->wf, hasBdJacPre);
2403:   return(0);
2404: }

2406: /*@C
2407:   PetscDSGetBdJacobianPreconditioner - Get the pointwise boundary Jacobian preconditioner function for given test and basis field

2409:   Not collective

2411:   Input Parameters:
2412: + ds - The PetscDS
2413: . f  - The test field number
2414: - g  - The field number

2416:   Output Parameters:
2417: + g0 - integrand for the test and basis function term
2418: . g1 - integrand for the test function and basis function gradient term
2419: . g2 - integrand for the test function gradient and basis function term
2420: - g3 - integrand for the test function gradient and basis function gradient term

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

2424:   \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

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

2428: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2429: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2430: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2431: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])

2433: + dim - the spatial dimension
2434: . Nf - the number of fields
2435: . NfAux - the number of auxiliary fields
2436: . uOff - the offset into u[] and u_t[] for each field
2437: . uOff_x - the offset into u_x[] for each field
2438: . u - each field evaluated at the current point
2439: . u_t - the time derivative of each field evaluated at the current point
2440: . u_x - the gradient of each field evaluated at the current point
2441: . aOff - the offset into a[] and a_t[] for each auxiliary field
2442: . aOff_x - the offset into a_x[] for each auxiliary field
2443: . a - each auxiliary field evaluated at the current point
2444: . a_t - the time derivative of each auxiliary field evaluated at the current point
2445: . a_x - the gradient of auxiliary each field evaluated at the current point
2446: . t - current time
2447: . u_tShift - the multiplier a for dF/dU_t
2448: . x - coordinates of the current point
2449: . n - normal at the current point
2450: . numConstants - number of constant parameters
2451: . constants - constant parameters
2452: - g0 - output values at the current point

2454:   This is not yet available in Fortran.

2456:   Level: intermediate

2458: .seealso: PetscDSSetBdJacobianPreconditioner()
2459: @*/
2460: PetscErrorCode PetscDSGetBdJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g,
2461:                                                   void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2462:                                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2463:                                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2464:                                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
2465:                                                   void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2466:                                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2467:                                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2468:                                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
2469:                                                   void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2470:                                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2471:                                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2472:                                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
2473:                                                   void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2474:                                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2475:                                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2476:                                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
2477: {
2478:   PetscBdPointJac *tmp0, *tmp1, *tmp2, *tmp3;
2479:   PetscInt         n0, n1, n2, n3;
2480:   PetscErrorCode   ierr;

2484:   if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
2485:   if ((g < 0) || (g >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, ds->Nf);
2486:   PetscWeakFormGetBdJacobianPreconditioner(ds->wf, NULL, 0, f, g, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3);
2487:   *g0  = tmp0 ? tmp0[0] : NULL;
2488:   *g1  = tmp1 ? tmp1[0] : NULL;
2489:   *g2  = tmp2 ? tmp2[0] : NULL;
2490:   *g3  = tmp3 ? tmp3[0] : NULL;
2491:   return(0);
2492: }

2494: /*@C
2495:   PetscDSSetBdJacobianPreconditioner - Set the pointwise boundary Jacobian preconditioner function for given test and basis field

2497:   Not collective

2499:   Input Parameters:
2500: + ds - The PetscDS
2501: . f  - The test field number
2502: . g  - The field number
2503: . g0 - integrand for the test and basis function term
2504: . g1 - integrand for the test function and basis function gradient term
2505: . g2 - integrand for the test function gradient and basis function term
2506: - g3 - integrand for the test function gradient and basis function gradient term

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

2510:   \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

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

2514: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2515: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2516: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2517: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])

2519: + dim - the spatial dimension
2520: . Nf - the number of fields
2521: . NfAux - the number of auxiliary fields
2522: . uOff - the offset into u[] and u_t[] for each field
2523: . uOff_x - the offset into u_x[] for each field
2524: . u - each field evaluated at the current point
2525: . u_t - the time derivative of each field evaluated at the current point
2526: . u_x - the gradient of each field evaluated at the current point
2527: . aOff - the offset into a[] and a_t[] for each auxiliary field
2528: . aOff_x - the offset into a_x[] for each auxiliary field
2529: . a - each auxiliary field evaluated at the current point
2530: . a_t - the time derivative of each auxiliary field evaluated at the current point
2531: . a_x - the gradient of auxiliary each field evaluated at the current point
2532: . t - current time
2533: . u_tShift - the multiplier a for dF/dU_t
2534: . x - coordinates of the current point
2535: . n - normal at the current point
2536: . numConstants - number of constant parameters
2537: . constants - constant parameters
2538: - g0 - output values at the current point

2540:   This is not yet available in Fortran.

2542:   Level: intermediate

2544: .seealso: PetscDSGetBdJacobianPreconditioner()
2545: @*/
2546: PetscErrorCode PetscDSSetBdJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g,
2547:                                                   void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2548:                                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2549:                                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2550:                                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
2551:                                                   void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2552:                                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2553:                                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2554:                                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
2555:                                                   void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2556:                                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2557:                                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2558:                                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
2559:                                                   void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2560:                                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2561:                                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2562:                                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
2563: {

2572:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2573:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
2574:   PetscWeakFormSetIndexBdJacobianPreconditioner(ds->wf, NULL, 0, f, g, 0, g0, 0, g1, 0, g2, 0, g3);
2575:   return(0);
2576: }

2578: /*@C
2579:   PetscDSGetExactSolution - Get the pointwise exact solution function for a given test field

2581:   Not collective

2583:   Input Parameters:
2584: + prob - The PetscDS
2585: - f    - The test field number

2587:   Output Parameter:
2588: + exactSol - exact solution for the test field
2589: - exactCtx - exact solution context

2591:   Note: The calling sequence for the solution functions is given by:

2593: $ sol(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx)

2595: + dim - the spatial dimension
2596: . t - current time
2597: . x - coordinates of the current point
2598: . Nc - the number of field components
2599: . u - the solution field evaluated at the current point
2600: - ctx - a user context

2602:   Level: intermediate

2604: .seealso: PetscDSSetExactSolution(), PetscDSGetExactSolutionTimeDerivative()
2605: @*/
2606: PetscErrorCode PetscDSGetExactSolution(PetscDS prob, PetscInt f, PetscErrorCode (**sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void **ctx)
2607: {
2610:   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);
2613:   return(0);
2614: }

2616: /*@C
2617:   PetscDSSetExactSolution - Set the pointwise exact solution function for a given test field

2619:   Not collective

2621:   Input Parameters:
2622: + prob - The PetscDS
2623: . f    - The test field number
2624: . sol  - solution function for the test fields
2625: - ctx  - solution context or NULL

2627:   Note: The calling sequence for solution functions is given by:

2629: $ sol(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx)

2631: + dim - the spatial dimension
2632: . t - current time
2633: . x - coordinates of the current point
2634: . Nc - the number of field components
2635: . u - the solution field evaluated at the current point
2636: - ctx - a user context

2638:   Level: intermediate

2640: .seealso: PetscDSGetExactSolution()
2641: @*/
2642: PetscErrorCode PetscDSSetExactSolution(PetscDS prob, PetscInt f, PetscErrorCode (*sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void *ctx)
2643: {

2648:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2649:   PetscDSEnlarge_Static(prob, f+1);
2652:   return(0);
2653: }

2655: /*@C
2656:   PetscDSGetExactSolutionTimeDerivative - Get the pointwise time derivative of the exact solution function for a given test field

2658:   Not collective

2660:   Input Parameters:
2661: + prob - The PetscDS
2662: - f    - The test field number

2664:   Output Parameter:
2665: + exactSol - time derivative of the exact solution for the test field
2666: - exactCtx - time derivative of the exact solution context

2668:   Note: The calling sequence for the solution functions is given by:

2670: $ sol(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx)

2672: + dim - the spatial dimension
2673: . t - current time
2674: . x - coordinates of the current point
2675: . Nc - the number of field components
2676: . u - the solution field evaluated at the current point
2677: - ctx - a user context

2679:   Level: intermediate

2681: .seealso: PetscDSSetExactSolutionTimeDerivative(), PetscDSGetExactSolution()
2682: @*/
2683: PetscErrorCode PetscDSGetExactSolutionTimeDerivative(PetscDS prob, PetscInt f, PetscErrorCode (**sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void **ctx)
2684: {
2687:   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);
2690:   return(0);
2691: }

2693: /*@C
2694:   PetscDSSetExactSolutionTimeDerivative - Set the pointwise time derivative of the exact solution function for a given test field

2696:   Not collective

2698:   Input Parameters:
2699: + prob - The PetscDS
2700: . f    - The test field number
2701: . sol  - time derivative of the solution function for the test fields
2702: - ctx  - time derivative of the solution context or NULL

2704:   Note: The calling sequence for solution functions is given by:

2706: $ sol(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx)

2708: + dim - the spatial dimension
2709: . t - current time
2710: . x - coordinates of the current point
2711: . Nc - the number of field components
2712: . u - the solution field evaluated at the current point
2713: - ctx - a user context

2715:   Level: intermediate

2717: .seealso: PetscDSGetExactSolutionTimeDerivative(), PetscDSSetExactSolution()
2718: @*/
2719: PetscErrorCode PetscDSSetExactSolutionTimeDerivative(PetscDS prob, PetscInt f, PetscErrorCode (*sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void *ctx)
2720: {

2725:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2726:   PetscDSEnlarge_Static(prob, f+1);
2729:   return(0);
2730: }

2732: /*@C
2733:   PetscDSGetConstants - Returns the array of constants passed to point functions

2735:   Not collective

2737:   Input Parameter:
2738: . prob - The PetscDS object

2740:   Output Parameters:
2741: + numConstants - The number of constants
2742: - constants    - The array of constants, NULL if there are none

2744:   Level: intermediate

2746: .seealso: PetscDSSetConstants(), PetscDSCreate()
2747: @*/
2748: PetscErrorCode PetscDSGetConstants(PetscDS prob, PetscInt *numConstants, const PetscScalar *constants[])
2749: {
2754:   return(0);
2755: }

2757: /*@C
2758:   PetscDSSetConstants - Set the array of constants passed to point functions

2760:   Not collective

2762:   Input Parameters:
2763: + prob         - The PetscDS object
2764: . numConstants - The number of constants
2765: - constants    - The array of constants, NULL if there are none

2767:   Level: intermediate

2769: .seealso: PetscDSGetConstants(), PetscDSCreate()
2770: @*/
2771: PetscErrorCode PetscDSSetConstants(PetscDS prob, PetscInt numConstants, PetscScalar constants[])
2772: {

2777:   if (numConstants != prob->numConstants) {
2778:     PetscFree(prob->constants);
2779:     prob->numConstants = numConstants;
2780:     if (prob->numConstants) {
2781:       PetscMalloc1(prob->numConstants, &prob->constants);
2782:     } else {
2783:       prob->constants = NULL;
2784:     }
2785:   }
2786:   if (prob->numConstants) {
2788:     PetscArraycpy(prob->constants, constants, prob->numConstants);
2789:   }
2790:   return(0);
2791: }

2793: /*@
2794:   PetscDSGetFieldIndex - Returns the index of the given field

2796:   Not collective

2798:   Input Parameters:
2799: + prob - The PetscDS object
2800: - disc - The discretization object

2802:   Output Parameter:
2803: . f - The field number

2805:   Level: beginner

2807: .seealso: PetscGetDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
2808: @*/
2809: PetscErrorCode PetscDSGetFieldIndex(PetscDS prob, PetscObject disc, PetscInt *f)
2810: {
2811:   PetscInt g;

2816:   *f = -1;
2817:   for (g = 0; g < prob->Nf; ++g) {if (disc == prob->disc[g]) break;}
2818:   if (g == prob->Nf) SETERRQ(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Field not found in PetscDS.");
2819:   *f = g;
2820:   return(0);
2821: }

2823: /*@
2824:   PetscDSGetFieldSize - Returns the size of the given field in the full space basis

2826:   Not collective

2828:   Input Parameters:
2829: + prob - The PetscDS object
2830: - f - The field number

2832:   Output Parameter:
2833: . size - The size

2835:   Level: beginner

2837: .seealso: PetscDSGetFieldOffset(), PetscDSGetNumFields(), PetscDSCreate()
2838: @*/
2839: PetscErrorCode PetscDSGetFieldSize(PetscDS prob, PetscInt f, PetscInt *size)
2840: {

2846:   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);
2847:   PetscDSSetUp(prob);
2848:   *size = prob->Nb[f];
2849:   return(0);
2850: }

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

2855:   Not collective

2857:   Input Parameters:
2858: + prob - The PetscDS object
2859: - f - The field number

2861:   Output Parameter:
2862: . off - The offset

2864:   Level: beginner

2866: .seealso: PetscDSGetFieldSize(), PetscDSGetNumFields(), PetscDSCreate()
2867: @*/
2868: PetscErrorCode PetscDSGetFieldOffset(PetscDS prob, PetscInt f, PetscInt *off)
2869: {
2870:   PetscInt       size, g;

2876:   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);
2877:   *off = 0;
2878:   for (g = 0; g < f; ++g) {
2879:     PetscDSGetFieldSize(prob, g, &size);
2880:     *off += size;
2881:   }
2882:   return(0);
2883: }

2885: /*@
2886:   PetscDSGetDimensions - Returns the size of the approximation space for each field on an evaluation point

2888:   Not collective

2890:   Input Parameter:
2891: . prob - The PetscDS object

2893:   Output Parameter:
2894: . dimensions - The number of dimensions

2896:   Level: beginner

2898: .seealso: PetscDSGetComponentOffsets(), PetscDSGetNumFields(), PetscDSCreate()
2899: @*/
2900: PetscErrorCode PetscDSGetDimensions(PetscDS prob, PetscInt *dimensions[])
2901: {

2906:   PetscDSSetUp(prob);
2908:   *dimensions = prob->Nb;
2909:   return(0);
2910: }

2912: /*@
2913:   PetscDSGetComponents - Returns the number of components for each field on an evaluation point

2915:   Not collective

2917:   Input Parameter:
2918: . prob - The PetscDS object

2920:   Output Parameter:
2921: . components - The number of components

2923:   Level: beginner

2925: .seealso: PetscDSGetComponentOffsets(), PetscDSGetNumFields(), PetscDSCreate()
2926: @*/
2927: PetscErrorCode PetscDSGetComponents(PetscDS prob, PetscInt *components[])
2928: {

2933:   PetscDSSetUp(prob);
2935:   *components = prob->Nc;
2936:   return(0);
2937: }

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

2942:   Not collective

2944:   Input Parameters:
2945: + prob - The PetscDS object
2946: - f - The field number

2948:   Output Parameter:
2949: . off - The offset

2951:   Level: beginner

2953: .seealso: PetscDSGetNumFields(), PetscDSCreate()
2954: @*/
2955: PetscErrorCode PetscDSGetComponentOffset(PetscDS prob, PetscInt f, PetscInt *off)
2956: {
2960:   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);
2961:   *off = prob->off[f];
2962:   return(0);
2963: }

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

2968:   Not collective

2970:   Input Parameter:
2971: . prob - The PetscDS object

2973:   Output Parameter:
2974: . offsets - The offsets

2976:   Level: beginner

2978: .seealso: PetscDSGetNumFields(), PetscDSCreate()
2979: @*/
2980: PetscErrorCode PetscDSGetComponentOffsets(PetscDS prob, PetscInt *offsets[])
2981: {
2985:   *offsets = prob->off;
2986:   return(0);
2987: }

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

2992:   Not collective

2994:   Input Parameter:
2995: . prob - The PetscDS object

2997:   Output Parameter:
2998: . offsets - The offsets

3000:   Level: beginner

3002: .seealso: PetscDSGetNumFields(), PetscDSCreate()
3003: @*/
3004: PetscErrorCode PetscDSGetComponentDerivativeOffsets(PetscDS prob, PetscInt *offsets[])
3005: {
3009:   *offsets = prob->offDer;
3010:   return(0);
3011: }

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

3016:   Not collective

3018:   Input Parameter:
3019: . prob - The PetscDS object

3021:   Output Parameter:
3022: . T - The basis function and derivatives tabulation at quadrature points for each field

3024:   Level: intermediate

3026: .seealso: PetscDSCreate()
3027: @*/
3028: PetscErrorCode PetscDSGetTabulation(PetscDS prob, PetscTabulation *T[])
3029: {

3035:   PetscDSSetUp(prob);
3036:   *T = prob->T;
3037:   return(0);
3038: }

3040: /*@C
3041:   PetscDSGetFaceTabulation - Return the basis tabulation at quadrature points on the faces

3043:   Not collective

3045:   Input Parameter:
3046: . prob - The PetscDS object

3048:   Output Parameter:
3049: . Tf - The basis function and derviative tabulation on each local face at quadrature points for each and field

3051:   Level: intermediate

3053: .seealso: PetscDSGetTabulation(), PetscDSCreate()
3054: @*/
3055: PetscErrorCode PetscDSGetFaceTabulation(PetscDS prob, PetscTabulation *Tf[])
3056: {

3062:   PetscDSSetUp(prob);
3063:   *Tf = prob->Tf;
3064:   return(0);
3065: }

3067: PetscErrorCode PetscDSGetEvaluationArrays(PetscDS prob, PetscScalar **u, PetscScalar **u_t, PetscScalar **u_x)
3068: {

3073:   PetscDSSetUp(prob);
3077:   return(0);
3078: }

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

3086:   PetscDSSetUp(prob);
3093:   return(0);
3094: }

3096: PetscErrorCode PetscDSGetWorkspace(PetscDS prob, PetscReal **x, PetscScalar **basisReal, PetscScalar **basisDerReal, PetscScalar **testReal, PetscScalar **testDerReal)
3097: {

3102:   PetscDSSetUp(prob);
3108:   return(0);
3109: }

3111: /*@C
3112:   PetscDSAddBoundary - Add a boundary condition to the model. The pointwise functions are used to provide boundary values for essential boundary conditions. In FEM, they are acting upon by dual basis functionals to generate FEM coefficients which are fixed. Natural boundary conditions signal to PETSc that boundary integrals should be performaed, using the kernels from PetscDSSetBdResidual().

3114:   Collective on ds

3116:   Input Parameters:
3117: + ds          - The PetscDS object
3118: . type        - The type of condition, e.g. DM_BC_ESSENTIAL/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
3119: . name        - The BC name
3120: . labelname   - The label defining constrained points
3121: . field       - The field to constrain
3122: . numcomps    - The number of constrained field components (0 will constrain all fields)
3123: . comps       - An array of constrained component numbers
3124: . bcFunc      - A pointwise function giving boundary values
3125: . bcFunc_t    - A pointwise function giving the time derviative of the boundary values, or NULL
3126: . numids      - The number of DMLabel ids for constrained points
3127: . ids         - An array of ids for constrained points
3128: - ctx         - An optional user context for bcFunc

3130:   Options Database Keys:
3131: + -bc_<boundary name> <num> - Overrides the boundary ids
3132: - -bc_<boundary name>_comp <num> - Overrides the boundary components

3134:   Note:
3135:   Both bcFunc abd bcFunc_t will depend on the boundary condition type. If the type if DM_BC_ESSENTIAL, Then the calling sequence is:

3137: $ bcFunc(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar bcval[])

3139:   If the type is DM_BC_ESSENTIAL_FIELD or other _FIELD value, then the calling sequence is:

3141: $ bcFunc(PetscInt dim, PetscInt Nf, PetscInt NfAux,
3142: $        const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
3143: $        const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
3144: $        PetscReal time, const PetscReal x[], PetscScalar bcval[])

3146: + dim - the spatial dimension
3147: . Nf - the number of fields
3148: . uOff - the offset into u[] and u_t[] for each field
3149: . uOff_x - the offset into u_x[] for each field
3150: . u - each field evaluated at the current point
3151: . u_t - the time derivative of each field evaluated at the current point
3152: . u_x - the gradient of each field evaluated at the current point
3153: . aOff - the offset into a[] and a_t[] for each auxiliary field
3154: . aOff_x - the offset into a_x[] for each auxiliary field
3155: . a - each auxiliary field evaluated at the current point
3156: . a_t - the time derivative of each auxiliary field evaluated at the current point
3157: . a_x - the gradient of auxiliary each field evaluated at the current point
3158: . t - current time
3159: . x - coordinates of the current point
3160: . numConstants - number of constant parameters
3161: . constants - constant parameters
3162: - bcval - output values at the current point

3164:   Level: developer

3166: .seealso: PetscDSGetBoundary(), PetscDSSetResidual(), PetscDSSetBdResidual()
3167: @*/
3168: PetscErrorCode PetscDSAddBoundary(PetscDS ds, DMBoundaryConditionType type, const char name[], const char labelname[], PetscInt field, PetscInt numcomps, const PetscInt *comps, void (*bcFunc)(void), void (*bcFunc_t)(void), PetscInt numids, const PetscInt *ids, void *ctx)
3169: {
3170:   DSBoundary     b;

3179:   PetscNew(&b);
3180:   PetscStrallocpy(name, (char **) &b->name);
3181:   PetscStrallocpy(labelname, (char **) &b->labelname);
3182:   PetscMalloc1(numcomps, &b->comps);
3183:   if (numcomps) {PetscArraycpy(b->comps, comps, numcomps);}
3184:   PetscMalloc1(numids, &b->ids);
3185:   if (numids) {PetscArraycpy(b->ids, ids, numids);}
3186:   b->type            = type;
3187:   b->field           = field;
3188:   b->numcomps        = numcomps;
3189:   b->func            = bcFunc;
3190:   b->func_t          = bcFunc_t;
3191:   b->numids          = numids;
3192:   b->ctx             = ctx;
3193:   b->next            = ds->boundary;
3194:   ds->boundary       = b;
3195:   return(0);
3196: }

3198: /*@C
3199:   PetscDSUpdateBoundary - Change a boundary condition for the model. The pointwise functions are used to provide boundary values for essential boundary conditions. In FEM, they are acting upon by dual basis functionals to generate FEM coefficients which are fixed. Natural boundary conditions signal to PETSc that boundary integrals should be performaed, using the kernels from PetscDSSetBdResidual().

3201:   Input Parameters:
3202: + ds          - The PetscDS object
3203: . bd          - The boundary condition number
3204: . type        - The type of condition, e.g. DM_BC_ESSENTIAL/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
3205: . name        - The BC name
3206: . labelname   - The label defining constrained points
3207: . field       - The field to constrain
3208: . numcomps    - The number of constrained field components
3209: . comps       - An array of constrained component numbers
3210: . bcFunc      - A pointwise function giving boundary values
3211: . bcFunc_t    - A pointwise function giving the time derviative of the boundary values, or NULL
3212: . numids      - The number of DMLabel ids for constrained points
3213: . ids         - An array of ids for constrained points
3214: - ctx         - An optional user context for bcFunc

3216:   Note:
3217:   The boundary condition number is the order in which it was registered. The user can get the number of boundary conditions from PetscDSGetNumBoundary(). See PetscDSAddBoundary() for a description of the calling sequences for the callbacks.

3219:   Level: developer

3221: .seealso: PetscDSAddBoundary(), PetscDSGetBoundary(), PetscDSGetNumBoundary()
3222: @*/
3223: PetscErrorCode PetscDSUpdateBoundary(PetscDS ds, PetscInt bd, DMBoundaryConditionType type, const char name[], const char labelname[], PetscInt field, PetscInt numcomps, const PetscInt *comps, void (*bcFunc)(void), void (*bcFunc_t)(void), PetscInt numids, const PetscInt *ids, void *ctx)
3224: {
3225:   DSBoundary     b = ds->boundary;
3226:   PetscInt       n = 0;

3231:   while (b) {
3232:     if (n == bd) break;
3233:     b = b->next;
3234:     ++n;
3235:   }
3236:   if (!b) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %d is not in [0, %d)", bd, n);
3237:   if (name) {
3238:     PetscFree(b->name);
3239:     PetscStrallocpy(name, (char **) &b->name);
3240:   }
3241:   if (labelname) {
3242:     PetscFree(b->labelname);
3243:     PetscStrallocpy(labelname, (char **) &b->labelname);
3244:   }
3245:   if (numcomps >= 0 && numcomps != b->numcomps) {
3246:     b->numcomps = numcomps;
3247:     PetscFree(b->comps);
3248:     PetscMalloc1(numcomps, &b->comps);
3249:     if (numcomps) {PetscArraycpy(b->comps, comps, numcomps);}
3250:   }
3251:   if (numids >= 0 && numids != b->numids) {
3252:     b->numids = numids;
3253:     PetscFree(b->ids);
3254:     PetscMalloc1(numids, &b->ids);
3255:     if (numids) {PetscArraycpy(b->ids, ids, numids);}
3256:   }
3257:   b->type = type;
3258:   if (field >= 0) {b->field  = field;}
3259:   if (bcFunc)     {b->func   = bcFunc;}
3260:   if (bcFunc_t)   {b->func_t = bcFunc_t;}
3261:   if (ctx)        {b->ctx    = ctx;}
3262:   return(0);
3263: }

3265: /*@
3266:   PetscDSGetNumBoundary - Get the number of registered BC

3268:   Input Parameters:
3269: . ds - The PetscDS object

3271:   Output Parameters:
3272: . numBd - The number of BC

3274:   Level: intermediate

3276: .seealso: PetscDSAddBoundary(), PetscDSGetBoundary()
3277: @*/
3278: PetscErrorCode PetscDSGetNumBoundary(PetscDS ds, PetscInt *numBd)
3279: {
3280:   DSBoundary b = ds->boundary;

3285:   *numBd = 0;
3286:   while (b) {++(*numBd); b = b->next;}
3287:   return(0);
3288: }

3290: /*@C
3291:   PetscDSGetBoundary - Gets a boundary condition to the model

3293:   Input Parameters:
3294: + ds          - The PetscDS object
3295: - bd          - The BC number

3297:   Output Parameters:
3298: + type        - The type of condition, e.g. DM_BC_ESSENTIAL/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
3299: . name        - The BC name
3300: . labelname   - The label defining constrained points
3301: . field       - The field to constrain
3302: . numcomps    - The number of constrained field components
3303: . comps       - An array of constrained component numbers
3304: . bcFunc      - A pointwise function giving boundary values
3305: . bcFunc_t    - A pointwise function giving the time derviative of the boundary values
3306: . numids      - The number of DMLabel ids for constrained points
3307: . ids         - An array of ids for constrained points
3308: - ctx         - An optional user context for bcFunc

3310:   Options Database Keys:
3311: + -bc_<boundary name> <num> - Overrides the boundary ids
3312: - -bc_<boundary name>_comp <num> - Overrides the boundary components

3314:   Level: developer

3316: .seealso: PetscDSAddBoundary()
3317: @*/
3318: PetscErrorCode PetscDSGetBoundary(PetscDS ds, PetscInt bd, DMBoundaryConditionType *type, const char **name, const char **labelname, PetscInt *field, PetscInt *numcomps, const PetscInt **comps, void (**func)(void), void (**func_t)(void), PetscInt *numids, const PetscInt **ids, void **ctx)
3319: {
3320:   DSBoundary b    = ds->boundary;
3321:   PetscInt   n    = 0;

3325:   while (b) {
3326:     if (n == bd) break;
3327:     b = b->next;
3328:     ++n;
3329:   }
3330:   if (!b) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %d is not in [0, %d)", bd, n);
3331:   if (type) {
3333:     *type = b->type;
3334:   }
3335:   if (name) {
3337:     *name = b->name;
3338:   }
3339:   if (labelname) {
3341:     *labelname = b->labelname;
3342:   }
3343:   if (field) {
3345:     *field = b->field;
3346:   }
3347:   if (numcomps) {
3349:     *numcomps = b->numcomps;
3350:   }
3351:   if (comps) {
3353:     *comps = b->comps;
3354:   }
3355:   if (func) {
3357:     *func = b->func;
3358:   }
3359:   if (func_t) {
3361:     *func_t = b->func_t;
3362:   }
3363:   if (numids) {
3365:     *numids = b->numids;
3366:   }
3367:   if (ids) {
3369:     *ids = b->ids;
3370:   }
3371:   if (ctx) {
3373:     *ctx = b->ctx;
3374:   }
3375:   return(0);
3376: }

3378: /*@
3379:   PetscDSCopyBoundary - Copy all boundary condition objects to the new problem

3381:   Not collective

3383:   Input Parameters:
3384: + ds        - The source PetscDS object
3385: . numFields - The number of selected fields, or PETSC_DEFAULT for all fields
3386: - fields    - The selected fields, or NULL for all fields

3388:   Output Parameter:
3389: . newds     - The target PetscDS, now with a copy of the boundary conditions

3391:   Level: intermediate

3393: .seealso: PetscDSCopyEquations(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3394: @*/
3395: PetscErrorCode PetscDSCopyBoundary(PetscDS ds, PetscInt numFields, const PetscInt fields[], PetscDS newds)
3396: {
3397:   DSBoundary     b, next, *lastnext;

3403:   if (ds == newds) return(0);
3404:   next = newds->boundary;
3405:   while (next) {
3406:     DSBoundary b = next;

3408:     next = b->next;
3409:     PetscFree(b->comps);
3410:     PetscFree(b->ids);
3411:     PetscFree(b->name);
3412:     PetscFree(b->labelname);
3413:     PetscFree(b);
3414:   }
3415:   lastnext = &(newds->boundary);
3416:   for (b = ds->boundary; b; b = b->next) {
3417:     DSBoundary bNew;
3418:     PetscInt   fieldNew = -1;

3420:     if (numFields > 0 && fields) {
3421:       PetscInt f;

3423:       for (f = 0; f < numFields; ++f) if (b->field == fields[f]) break;
3424:       if (f == numFields) continue;
3425:       fieldNew = f;
3426:     }
3427:     PetscNew(&bNew);
3428:     bNew->numcomps = b->numcomps;
3429:     PetscMalloc1(bNew->numcomps, &bNew->comps);
3430:     PetscArraycpy(bNew->comps, b->comps, bNew->numcomps);
3431:     bNew->numids = b->numids;
3432:     PetscMalloc1(bNew->numids, &bNew->ids);
3433:     PetscArraycpy(bNew->ids, b->ids, bNew->numids);
3434:     PetscStrallocpy(b->labelname,(char **) &(bNew->labelname));
3435:     PetscStrallocpy(b->name,(char **) &(bNew->name));
3436:     bNew->ctx   = b->ctx;
3437:     bNew->type  = b->type;
3438:     bNew->field = fieldNew < 0 ? b->field : fieldNew;
3439:     bNew->func  = b->func;

3441:     *lastnext = bNew;
3442:     lastnext = &(bNew->next);
3443:   }
3444:   return(0);
3445: }

3447: /*@
3448:   PetscDSSelectDiscretizations - Copy discretizations to the new problem with different field layout

3450:   Not collective

3452:   Input Parameter:
3453: + prob - The PetscDS object
3454: . numFields - Number of new fields
3455: - fields - Old field number for each new field

3457:   Output Parameter:
3458: . newprob - The PetscDS copy

3460:   Level: intermediate

3462: .seealso: PetscDSSelectEquations(), PetscDSCopyBoundary(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3463: @*/
3464: PetscErrorCode PetscDSSelectDiscretizations(PetscDS prob, PetscInt numFields, const PetscInt fields[], PetscDS newprob)
3465: {
3466:   PetscInt       Nf, Nfn, fn;

3473:   PetscDSGetNumFields(prob, &Nf);
3474:   PetscDSGetNumFields(newprob, &Nfn);
3475:   if (numFields > Nfn) SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_SIZ, "Number of fields %D to transfer must not be greater then the total number of fields %D", numFields, Nfn);
3476:   for (fn = 0; fn < numFields; ++fn) {
3477:     const PetscInt f = fields ? fields[fn] : fn;
3478:     PetscObject    disc;

3480:     if (f >= Nf) continue;
3481:     PetscDSGetDiscretization(prob, f, &disc);
3482:     PetscDSSetDiscretization(newprob, fn, disc);
3483:   }
3484:   return(0);
3485: }

3487: /*@
3488:   PetscDSSelectEquations - Copy pointwise function pointers to the new problem with different field layout

3490:   Not collective

3492:   Input Parameter:
3493: + prob - The PetscDS object
3494: . numFields - Number of new fields
3495: - fields - Old field number for each new field

3497:   Output Parameter:
3498: . newprob - The PetscDS copy

3500:   Level: intermediate

3502: .seealso: PetscDSSelectDiscretizations(), PetscDSCopyBoundary(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3503: @*/
3504: PetscErrorCode PetscDSSelectEquations(PetscDS prob, PetscInt numFields, const PetscInt fields[], PetscDS newprob)
3505: {
3506:   PetscInt       Nf, Nfn, fn, gn;

3513:   PetscDSGetNumFields(prob, &Nf);
3514:   PetscDSGetNumFields(newprob, &Nfn);
3515:   if (numFields > Nfn) SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_SIZ, "Number of fields %D to transfer must not be greater then the total number of fields %D", numFields, Nfn);
3516:   for (fn = 0; fn < numFields; ++fn) {
3517:     const PetscInt   f = fields ? fields[fn] : fn;
3518:     PetscPointFunc   obj;
3519:     PetscPointFunc   f0, f1;
3520:     PetscBdPointFunc f0Bd, f1Bd;
3521:     PetscRiemannFunc r;

3523:     if (f >= Nf) continue;
3524:     PetscDSGetObjective(prob, f, &obj);
3525:     PetscDSGetResidual(prob, f, &f0, &f1);
3526:     PetscDSGetBdResidual(prob, f, &f0Bd, &f1Bd);
3527:     PetscDSGetRiemannSolver(prob, f, &r);
3528:     PetscDSSetObjective(newprob, fn, obj);
3529:     PetscDSSetResidual(newprob, fn, f0, f1);
3530:     PetscDSSetBdResidual(newprob, fn, f0Bd, f1Bd);
3531:     PetscDSSetRiemannSolver(newprob, fn, r);
3532:     for (gn = 0; gn < numFields; ++gn) {
3533:       const PetscInt  g = fields ? fields[gn] : gn;
3534:       PetscPointJac   g0, g1, g2, g3;
3535:       PetscPointJac   g0p, g1p, g2p, g3p;
3536:       PetscBdPointJac g0Bd, g1Bd, g2Bd, g3Bd;

3538:       if (g >= Nf) continue;
3539:       PetscDSGetJacobian(prob, f, g, &g0, &g1, &g2, &g3);
3540:       PetscDSGetJacobianPreconditioner(prob, f, g, &g0p, &g1p, &g2p, &g3p);
3541:       PetscDSGetBdJacobian(prob, f, g, &g0Bd, &g1Bd, &g2Bd, &g3Bd);
3542:       PetscDSSetJacobian(newprob, fn, gn, g0, g1, g2, g3);
3543:       PetscDSSetJacobianPreconditioner(newprob, fn, gn, g0p, g1p, g2p, g3p);
3544:       PetscDSSetBdJacobian(newprob, fn, gn, g0Bd, g1Bd, g2Bd, g3Bd);
3545:     }
3546:   }
3547:   return(0);
3548: }

3550: /*@
3551:   PetscDSCopyEquations - Copy all pointwise function pointers to the new problem

3553:   Not collective

3555:   Input Parameter:
3556: . prob - The PetscDS object

3558:   Output Parameter:
3559: . newprob - The PetscDS copy

3561:   Level: intermediate

3563: .seealso: PetscDSCopyBoundary(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3564: @*/
3565: PetscErrorCode PetscDSCopyEquations(PetscDS prob, PetscDS newprob)
3566: {
3567:   PetscInt       Nf, Ng;

3573:   PetscDSGetNumFields(prob, &Nf);
3574:   PetscDSGetNumFields(newprob, &Ng);
3575:   if (Nf != Ng) SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_SIZ, "Number of fields must match %D != %D", Nf, Ng);
3576:   PetscDSSelectEquations(prob, Nf, NULL, newprob);
3577:   return(0);
3578: }
3579: /*@
3580:   PetscDSCopyConstants - Copy all constants to the new problem

3582:   Not collective

3584:   Input Parameter:
3585: . prob - The PetscDS object

3587:   Output Parameter:
3588: . newprob - The PetscDS copy

3590:   Level: intermediate

3592: .seealso: PetscDSCopyBoundary(), PetscDSCopyEquations(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3593: @*/
3594: PetscErrorCode PetscDSCopyConstants(PetscDS prob, PetscDS newprob)
3595: {
3596:   PetscInt           Nc;
3597:   const PetscScalar *constants;
3598:   PetscErrorCode     ierr;

3603:   PetscDSGetConstants(prob, &Nc, &constants);
3604:   PetscDSSetConstants(newprob, Nc, (PetscScalar *) constants);
3605:   return(0);
3606: }

3608: PetscErrorCode PetscDSGetHeightSubspace(PetscDS prob, PetscInt height, PetscDS *subprob)
3609: {
3610:   PetscInt       dim, Nf, f;

3616:   if (height == 0) {*subprob = prob; return(0);}
3617:   PetscDSGetNumFields(prob, &Nf);
3618:   PetscDSGetSpatialDimension(prob, &dim);
3619:   if (height > dim) SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_OUTOFRANGE, "DS can only handle height in [0, %D], not %D", dim, height);
3620:   if (!prob->subprobs) {PetscCalloc1(dim, &prob->subprobs);}
3621:   if (!prob->subprobs[height-1]) {
3622:     PetscInt cdim;

3624:     PetscDSCreate(PetscObjectComm((PetscObject) prob), &prob->subprobs[height-1]);
3625:     PetscDSGetCoordinateDimension(prob, &cdim);
3626:     PetscDSSetCoordinateDimension(prob->subprobs[height-1], cdim);
3627:     for (f = 0; f < Nf; ++f) {
3628:       PetscFE      subfe;
3629:       PetscObject  obj;
3630:       PetscClassId id;

3632:       PetscDSGetDiscretization(prob, f, &obj);
3633:       PetscObjectGetClassId(obj, &id);
3634:       if (id == PETSCFE_CLASSID) {PetscFEGetHeightSubspace((PetscFE) obj, height, &subfe);}
3635:       else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unsupported discretization type for field %d", f);
3636:       PetscDSSetDiscretization(prob->subprobs[height-1], f, (PetscObject) subfe);
3637:     }
3638:   }
3639:   *subprob = prob->subprobs[height-1];
3640:   return(0);
3641: }

3643: PetscErrorCode PetscDSGetDiscType_Internal(PetscDS ds, PetscInt f, PetscDiscType *disctype)
3644: {
3645:   PetscObject    obj;
3646:   PetscClassId   id;
3647:   PetscInt       Nf;

3653:   *disctype = PETSC_DISC_NONE;
3654:   PetscDSGetNumFields(ds, &Nf);
3655:   if (f >= Nf) SETERRQ2(PetscObjectComm((PetscObject) ds), PETSC_ERR_ARG_SIZ, "Field %D must be in [0, %D)", f, Nf);
3656:   PetscDSGetDiscretization(ds, f, &obj);
3657:   if (obj) {
3658:     PetscObjectGetClassId(obj, &id);
3659:     if (id == PETSCFE_CLASSID) *disctype = PETSC_DISC_FE;
3660:     else                       *disctype = PETSC_DISC_FV;
3661:   }
3662:   return(0);
3663: }

3665: static PetscErrorCode PetscDSDestroy_Basic(PetscDS ds)
3666: {

3670:   PetscFree(ds->data);
3671:   return(0);
3672: }

3674: static PetscErrorCode PetscDSInitialize_Basic(PetscDS ds)
3675: {
3677:   ds->ops->setfromoptions = NULL;
3678:   ds->ops->setup          = NULL;
3679:   ds->ops->view           = NULL;
3680:   ds->ops->destroy        = PetscDSDestroy_Basic;
3681:   return(0);
3682: }

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

3687:   Level: intermediate

3689: .seealso: PetscDSType, PetscDSCreate(), PetscDSSetType()
3690: M*/

3692: PETSC_EXTERN PetscErrorCode PetscDSCreate_Basic(PetscDS ds)
3693: {
3694:   PetscDS_Basic *b;

3699:   PetscNewLog(ds, &b);
3700:   ds->data = b;

3702:   PetscDSInitialize_Basic(ds);
3703:   return(0);
3704: }