Actual source code: dtds.c

petsc-main 2021-03-01
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:   PetscViewerASCIIPopTab(viewer);
224:   return(0);
225: }

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

230:    Collective on PetscDS

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

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

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

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

253:   Collective on prob

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

259:   Level: developer

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

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

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

281:   Collective on prob

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

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

293:   Level: developer

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

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

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

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

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

358:   Collective on prob

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

363:   Level: developer

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

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

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

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

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

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

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

455: static PetscErrorCode PetscDSEnlarge_Static(PetscDS prob, PetscInt NfNew)
456: {
457:   PetscObject      *tmpd;
458:   PetscBool        *tmpi;
459:   PetscInt         *tmpk;
460:   PetscPointFunc   *tmpobj, *tmpf, *tmpup;
461:   PetscPointJac    *tmpg, *tmpgp, *tmpgt;
462:   PetscBdPointFunc *tmpfbd;
463:   PetscBdPointJac  *tmpgbd, *tmpgpbd;
464:   PetscRiemannFunc *tmpr;
465:   PetscSimplePointFunc *tmpexactSol,  *tmpexactSol_t;
466:   void                **tmpexactCtx, **tmpexactCtx_t;
467:   void            **tmpctx;
468:   PetscInt          Nf = prob->Nf, f;
469:   PetscErrorCode    ierr;

472:   if (Nf >= NfNew) return(0);
473:   prob->setup = PETSC_FALSE;
474:   PetscDSDestroyStructs_Static(prob);
475:   PetscMalloc3(NfNew, &tmpd, NfNew, &tmpi, NfNew, &tmpk);
476:   for (f = 0; f < Nf; ++f) {tmpd[f] = prob->disc[f]; tmpi[f] = prob->implicit[f]; tmpk[f] = prob->jetDegree[f];}
477:   for (f = Nf; f < NfNew; ++f) {tmpd[f] = NULL; tmpi[f] = PETSC_TRUE; tmpk[f] = 1;}
478:   PetscFree3(prob->disc, prob->implicit, prob->jetDegree);
479:   prob->Nf        = NfNew;
480:   prob->disc      = tmpd;
481:   prob->implicit  = tmpi;
482:   prob->jetDegree = tmpk;
483:   PetscCalloc7(NfNew, &tmpobj, NfNew*2, &tmpf, NfNew*NfNew*4, &tmpg, NfNew*NfNew*4, &tmpgp, NfNew*NfNew*4, &tmpgt, NfNew, &tmpr, NfNew, &tmpctx);
484:   PetscCalloc1(NfNew, &tmpup);
485:   for (f = 0; f < Nf; ++f) tmpobj[f] = prob->obj[f];
486:   for (f = 0; f < Nf*2; ++f) tmpf[f] = prob->f[f];
487:   for (f = 0; f < Nf*Nf*4; ++f) tmpg[f] = prob->g[f];
488:   for (f = 0; f < Nf*Nf*4; ++f) tmpgp[f] = prob->gp[f];
489:   for (f = 0; f < Nf; ++f) tmpr[f] = prob->r[f];
490:   for (f = 0; f < Nf; ++f) tmpup[f] = prob->update[f];
491:   for (f = 0; f < Nf; ++f) tmpctx[f] = prob->ctx[f];
492:   for (f = Nf; f < NfNew; ++f) tmpobj[f] = NULL;
493:   for (f = Nf*2; f < NfNew*2; ++f) tmpf[f] = NULL;
494:   for (f = Nf*Nf*4; f < NfNew*NfNew*4; ++f) tmpg[f] = NULL;
495:   for (f = Nf*Nf*4; f < NfNew*NfNew*4; ++f) tmpgp[f] = NULL;
496:   for (f = Nf*Nf*4; f < NfNew*NfNew*4; ++f) tmpgt[f] = NULL;
497:   for (f = Nf; f < NfNew; ++f) tmpr[f] = NULL;
498:   for (f = Nf; f < NfNew; ++f) tmpup[f] = NULL;
499:   for (f = Nf; f < NfNew; ++f) tmpctx[f] = NULL;
500:   PetscFree7(prob->obj, prob->f, prob->g, prob->gp, prob->gt, prob->r, prob->ctx);
501:   PetscFree(prob->update);
502:   prob->obj = tmpobj;
503:   prob->f   = tmpf;
504:   prob->g   = tmpg;
505:   prob->gp  = tmpgp;
506:   prob->gt  = tmpgt;
507:   prob->r   = tmpr;
508:   prob->update = tmpup;
509:   prob->ctx = tmpctx;
510:   PetscCalloc7(NfNew*2, &tmpfbd, NfNew*NfNew*4, &tmpgbd, NfNew*NfNew*4, &tmpgpbd, NfNew, &tmpexactSol, NfNew, &tmpexactCtx, NfNew, &tmpexactSol_t, NfNew, &tmpexactCtx_t);
511:   for (f = 0; f < Nf*2; ++f) tmpfbd[f] = prob->fBd[f];
512:   for (f = 0; f < Nf*Nf*4; ++f) tmpgbd[f] = prob->gBd[f];
513:   for (f = 0; f < Nf*Nf*4; ++f) tmpgpbd[f] = prob->gpBd[f];
514:   for (f = 0; f < Nf; ++f) tmpexactSol[f] = prob->exactSol[f];
515:   for (f = 0; f < Nf; ++f) tmpexactCtx[f] = prob->exactCtx[f];
516:   for (f = 0; f < Nf; ++f) tmpexactSol_t[f] = prob->exactSol_t[f];
517:   for (f = 0; f < Nf; ++f) tmpexactCtx_t[f] = prob->exactCtx_t[f];
518:   for (f = Nf*2; f < NfNew*2; ++f) tmpfbd[f] = NULL;
519:   for (f = Nf*Nf*4; f < NfNew*NfNew*4; ++f) tmpgbd[f] = NULL;
520:   for (f = Nf*Nf*4; f < NfNew*NfNew*4; ++f) tmpgpbd[f] = NULL;
521:   for (f = Nf; f < NfNew; ++f) tmpexactSol[f] = NULL;
522:   for (f = Nf; f < NfNew; ++f) tmpexactCtx[f] = NULL;
523:   for (f = Nf; f < NfNew; ++f) tmpexactSol_t[f] = NULL;
524:   for (f = Nf; f < NfNew; ++f) tmpexactCtx_t[f] = NULL;
525:   PetscFree7(prob->fBd, prob->gBd, prob->gpBd, prob->exactSol, prob->exactCtx, prob->exactSol_t, prob->exactCtx_t);
526:   prob->fBd = tmpfbd;
527:   prob->gBd = tmpgbd;
528:   prob->gpBd = tmpgpbd;
529:   prob->exactSol = tmpexactSol;
530:   prob->exactCtx = tmpexactCtx;
531:   prob->exactSol_t = tmpexactSol_t;
532:   prob->exactCtx_t = tmpexactCtx_t;
533:   return(0);
534: }

536: /*@
537:   PetscDSDestroy - Destroys a PetscDS object

539:   Collective on prob

541:   Input Parameter:
542: . prob - the PetscDS object to destroy

544:   Level: developer

546: .seealso PetscDSView()
547: @*/
548: PetscErrorCode PetscDSDestroy(PetscDS *prob)
549: {
550:   PetscInt       f;
551:   DSBoundary     next;

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

558:   if (--((PetscObject)(*prob))->refct > 0) {*prob = NULL; return(0);}
559:   ((PetscObject) (*prob))->refct = 0;
560:   if ((*prob)->subprobs) {
561:     PetscInt dim, d;

563:     PetscDSGetSpatialDimension(*prob, &dim);
564:     for (d = 0; d < dim; ++d) {PetscDSDestroy(&(*prob)->subprobs[d]);}
565:   }
566:   PetscFree((*prob)->subprobs);
567:   PetscDSDestroyStructs_Static(*prob);
568:   for (f = 0; f < (*prob)->Nf; ++f) {
569:     PetscObjectDereference((*prob)->disc[f]);
570:   }
571:   PetscFree3((*prob)->disc, (*prob)->implicit, (*prob)->jetDegree);
572:   PetscFree7((*prob)->obj,(*prob)->f,(*prob)->g,(*prob)->gp,(*prob)->gt,(*prob)->r,(*prob)->ctx);
573:   PetscFree((*prob)->update);
574:   PetscFree7((*prob)->fBd,(*prob)->gBd,(*prob)->gpBd,(*prob)->exactSol,(*prob)->exactCtx,(*prob)->exactSol_t,(*prob)->exactCtx_t);
575:   if ((*prob)->ops->destroy) {(*(*prob)->ops->destroy)(*prob);}
576:   next = (*prob)->boundary;
577:   while (next) {
578:     DSBoundary b = next;

580:     next = b->next;
581:     PetscFree(b->comps);
582:     PetscFree(b->ids);
583:     PetscFree(b->name);
584:     PetscFree(b->labelname);
585:     PetscFree(b);
586:   }
587:   PetscFree((*prob)->constants);
588:   PetscHeaderDestroy(prob);
589:   return(0);
590: }

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

595:   Collective

597:   Input Parameter:
598: . comm - The communicator for the PetscDS object

600:   Output Parameter:
601: . prob - The PetscDS object

603:   Level: beginner

605: .seealso: PetscDSSetType(), PETSCDSBASIC
606: @*/
607: PetscErrorCode PetscDSCreate(MPI_Comm comm, PetscDS *prob)
608: {
609:   PetscDS   p;

614:   *prob  = NULL;
615:   PetscDSInitializePackage();

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

619:   p->Nf           = 0;
620:   p->setup        = PETSC_FALSE;
621:   p->numConstants = 0;
622:   p->constants    = NULL;
623:   p->dimEmbed     = -1;
624:   p->useJacPre    = PETSC_TRUE;

626:   *prob = p;
627:   return(0);
628: }

630: /*@
631:   PetscDSGetNumFields - Returns the number of fields in the DS

633:   Not collective

635:   Input Parameter:
636: . prob - The PetscDS object

638:   Output Parameter:
639: . Nf - The number of fields

641:   Level: beginner

643: .seealso: PetscDSGetSpatialDimension(), PetscDSCreate()
644: @*/
645: PetscErrorCode PetscDSGetNumFields(PetscDS prob, PetscInt *Nf)
646: {
650:   *Nf = prob->Nf;
651:   return(0);
652: }

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

657:   Not collective

659:   Input Parameter:
660: . prob - The PetscDS object

662:   Output Parameter:
663: . dim - The spatial dimension

665:   Level: beginner

667: .seealso: PetscDSGetCoordinateDimension(), PetscDSGetNumFields(), PetscDSCreate()
668: @*/
669: PetscErrorCode PetscDSGetSpatialDimension(PetscDS prob, PetscInt *dim)
670: {

676:   *dim = 0;
677:   if (prob->Nf) {
678:     PetscObject  obj;
679:     PetscClassId id;

681:     PetscDSGetDiscretization(prob, 0, &obj);
682:     if (obj) {
683:       PetscObjectGetClassId(obj, &id);
684:       if (id == PETSCFE_CLASSID)      {PetscFEGetSpatialDimension((PetscFE) obj, dim);}
685:       else if (id == PETSCFV_CLASSID) {PetscFVGetSpatialDimension((PetscFV) obj, dim);}
686:       else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", 0);
687:     }
688:   }
689:   return(0);
690: }

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

695:   Not collective

697:   Input Parameter:
698: . prob - The PetscDS object

700:   Output Parameter:
701: . dimEmbed - The coordinate dimension

703:   Level: beginner

705: .seealso: PetscDSSetCoordinateDimension(), PetscDSGetSpatialDimension(), PetscDSGetNumFields(), PetscDSCreate()
706: @*/
707: PetscErrorCode PetscDSGetCoordinateDimension(PetscDS prob, PetscInt *dimEmbed)
708: {
712:   if (prob->dimEmbed < 0) SETERRQ(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONGSTATE, "No coordinate dimension set for this DS");
713:   *dimEmbed = prob->dimEmbed;
714:   return(0);
715: }

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

720:   Logically collective on prob

722:   Input Parameters:
723: + prob - The PetscDS object
724: - dimEmbed - The coordinate dimension

726:   Level: beginner

728: .seealso: PetscDSGetCoordinateDimension(), PetscDSGetSpatialDimension(), PetscDSGetNumFields(), PetscDSCreate()
729: @*/
730: PetscErrorCode PetscDSSetCoordinateDimension(PetscDS prob, PetscInt dimEmbed)
731: {
734:   if (dimEmbed < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Coordinate dimension must be non-negative, not %D", dimEmbed);
735:   prob->dimEmbed = dimEmbed;
736:   return(0);
737: }

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

742:   Not collective

744:   Input Parameter:
745: . prob - The PetscDS object

747:   Output Parameter:
748: . isHybrid - The flag

750:   Level: developer

752: .seealso: PetscDSSetHybrid(), PetscDSCreate()
753: @*/
754: PetscErrorCode PetscDSGetHybrid(PetscDS prob, PetscBool *isHybrid)
755: {
759:   *isHybrid = prob->isHybrid;
760:   return(0);
761: }

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

766:   Not collective

768:   Input Parameters:
769: + prob - The PetscDS object
770: - isHybrid - The flag

772:   Level: developer

774: .seealso: PetscDSGetHybrid(), PetscDSCreate()
775: @*/
776: PetscErrorCode PetscDSSetHybrid(PetscDS prob, PetscBool isHybrid)
777: {
780:   prob->isHybrid = isHybrid;
781:   return(0);
782: }

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

787:   Not collective

789:   Input Parameter:
790: . prob - The PetscDS object

792:   Output Parameter:
793: . dim - The total problem dimension

795:   Level: beginner

797: .seealso: PetscDSGetNumFields(), PetscDSCreate()
798: @*/
799: PetscErrorCode PetscDSGetTotalDimension(PetscDS prob, PetscInt *dim)
800: {

805:   PetscDSSetUp(prob);
807:   *dim = prob->totDim;
808:   return(0);
809: }

811: /*@
812:   PetscDSGetTotalComponents - Returns the total number of components in this system

814:   Not collective

816:   Input Parameter:
817: . prob - The PetscDS object

819:   Output Parameter:
820: . dim - The total number of components

822:   Level: beginner

824: .seealso: PetscDSGetNumFields(), PetscDSCreate()
825: @*/
826: PetscErrorCode PetscDSGetTotalComponents(PetscDS prob, PetscInt *Nc)
827: {

832:   PetscDSSetUp(prob);
834:   *Nc = prob->totComp;
835:   return(0);
836: }

838: /*@
839:   PetscDSGetDiscretization - Returns the discretization object for the given field

841:   Not collective

843:   Input Parameters:
844: + prob - The PetscDS object
845: - f - The field number

847:   Output Parameter:
848: . disc - The discretization object

850:   Level: beginner

852: .seealso: PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
853: @*/
854: PetscErrorCode PetscDSGetDiscretization(PetscDS prob, PetscInt f, PetscObject *disc)
855: {
859:   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);
860:   *disc = prob->disc[f];
861:   return(0);
862: }

864: /*@
865:   PetscDSSetDiscretization - Sets the discretization object for the given field

867:   Not collective

869:   Input Parameters:
870: + prob - The PetscDS object
871: . f - The field number
872: - disc - The discretization object

874:   Level: beginner

876: .seealso: PetscDSGetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
877: @*/
878: PetscErrorCode PetscDSSetDiscretization(PetscDS prob, PetscInt f, PetscObject disc)
879: {

885:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
886:   PetscDSEnlarge_Static(prob, f+1);
887:   PetscObjectDereference(prob->disc[f]);
888:   prob->disc[f] = disc;
889:   PetscObjectReference(disc);
890:   if (disc) {
891:     PetscClassId id;

893:     PetscObjectGetClassId(disc, &id);
894:     if (id == PETSCFE_CLASSID) {
895:       PetscDSSetImplicit(prob, f, PETSC_TRUE);
896:     } else if (id == PETSCFV_CLASSID) {
897:       PetscDSSetImplicit(prob, f, PETSC_FALSE);
898:     }
899:     PetscDSSetJetDegree(prob, f, 1);
900:   }
901:   return(0);
902: }

904: /*@
905:   PetscDSAddDiscretization - Adds a discretization object

907:   Not collective

909:   Input Parameters:
910: + prob - The PetscDS object
911: - disc - The boundary discretization object

913:   Level: beginner

915: .seealso: PetscDSGetDiscretization(), PetscDSSetDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
916: @*/
917: PetscErrorCode PetscDSAddDiscretization(PetscDS prob, PetscObject disc)
918: {

922:   PetscDSSetDiscretization(prob, prob->Nf, disc);
923:   return(0);
924: }

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

929:   Not collective

931:   Input Parameter:
932: . prob - The PetscDS object

934:   Output Parameter:
935: . q - The quadrature object

937: Level: intermediate

939: .seealso: PetscDSSetImplicit(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
940: @*/
941: PetscErrorCode PetscDSGetQuadrature(PetscDS prob, PetscQuadrature *q)
942: {
943:   PetscObject    obj;
944:   PetscClassId   id;

948:   *q = NULL;
949:   if (!prob->Nf) return(0);
950:   PetscDSGetDiscretization(prob, 0, &obj);
951:   PetscObjectGetClassId(obj, &id);
952:   if      (id == PETSCFE_CLASSID) {PetscFEGetQuadrature((PetscFE) obj, q);}
953:   else if (id == PETSCFV_CLASSID) {PetscFVGetQuadrature((PetscFV) obj, q);}
954:   else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", 0);
955:   return(0);
956: }

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

961:   Not collective

963:   Input Parameters:
964: + prob - The PetscDS object
965: - f - The field number

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

970:   Level: developer

972: .seealso: PetscDSSetImplicit(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
973: @*/
974: PetscErrorCode PetscDSGetImplicit(PetscDS prob, PetscInt f, PetscBool *implicit)
975: {
979:   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);
980:   *implicit = prob->implicit[f];
981:   return(0);
982: }

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

987:   Not collective

989:   Input Parameters:
990: + prob - The PetscDS object
991: . f - The field number
992: - implicit - The flag indicating what kind of solve to use for this field

994:   Level: developer

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

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

1010:   Not collective

1012:   Input Parameters:
1013: + ds - The PetscDS object
1014: - f  - The field number

1016:   Output Parameter:
1017: . k  - The highest derivative we need to tabulate

1019:   Level: developer

1021: .seealso: PetscDSSetJetDegree(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
1022: @*/
1023: PetscErrorCode PetscDSGetJetDegree(PetscDS ds, PetscInt f, PetscInt *k)
1024: {
1028:   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);
1029:   *k = ds->jetDegree[f];
1030:   return(0);
1031: }

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

1036:   Not collective

1038:   Input Parameters:
1039: + ds - The PetscDS object
1040: . f  - The field number
1041: - k  - The highest derivative we need to tabulate

1043:   Level: developer

1045: .seealso: PetscDSGetJetDegree(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
1046: @*/
1047: PetscErrorCode PetscDSSetJetDegree(PetscDS ds, PetscInt f, PetscInt k)
1048: {
1051:   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);
1052:   ds->jetDegree[f] = k;
1053:   return(0);
1054: }

1056: PetscErrorCode PetscDSGetObjective(PetscDS prob, PetscInt f,
1057:                                    void (**obj)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1058:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1059:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1060:                                                 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[]))
1061: {
1065:   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);
1066:   *obj = prob->obj[f];
1067:   return(0);
1068: }

1070: PetscErrorCode PetscDSSetObjective(PetscDS prob, PetscInt f,
1071:                                    void (*obj)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1072:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1073:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1074:                                                PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[]))
1075: {

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

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

1090:   Not collective

1092:   Input Parameters:
1093: + prob - The PetscDS
1094: - f    - The test field number

1096:   Output Parameters:
1097: + f0 - integrand for the test function term
1098: - f1 - integrand for the test function gradient term

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

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

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

1106: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1107: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1108: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1109: $    PetscReal t, const PetscReal x[], PetscScalar f0[])

1111: + dim - the spatial dimension
1112: . Nf - the number of fields
1113: . uOff - the offset into u[] and u_t[] for each field
1114: . uOff_x - the offset into u_x[] for each field
1115: . u - each field evaluated at the current point
1116: . u_t - the time derivative of each field evaluated at the current point
1117: . u_x - the gradient of each field evaluated at the current point
1118: . aOff - the offset into a[] and a_t[] for each auxiliary field
1119: . aOff_x - the offset into a_x[] for each auxiliary field
1120: . a - each auxiliary field evaluated at the current point
1121: . a_t - the time derivative of each auxiliary field evaluated at the current point
1122: . a_x - the gradient of auxiliary each field evaluated at the current point
1123: . t - current time
1124: . x - coordinates of the current point
1125: . numConstants - number of constant parameters
1126: . constants - constant parameters
1127: - f0 - output values at the current point

1129:   Level: intermediate

1131: .seealso: PetscDSSetResidual()
1132: @*/
1133: PetscErrorCode PetscDSGetResidual(PetscDS prob, PetscInt f,
1134:                                   void (**f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1135:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1136:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1137:                                               PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]),
1138:                                   void (**f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1139:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1140:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1141:                                               PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
1142: {
1145:   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);
1148:   return(0);
1149: }

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

1154:   Not collective

1156:   Input Parameters:
1157: + prob - The PetscDS
1158: . f    - The test field number
1159: . f0 - integrand for the test function term
1160: - f1 - integrand for the test function gradient term

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

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

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

1168: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1169: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1170: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1171: $    PetscReal t, const PetscReal x[], PetscScalar f0[])

1173: + dim - the spatial dimension
1174: . Nf - the number of fields
1175: . uOff - the offset into u[] and u_t[] for each field
1176: . uOff_x - the offset into u_x[] for each field
1177: . u - each field evaluated at the current point
1178: . u_t - the time derivative of each field evaluated at the current point
1179: . u_x - the gradient of each field evaluated at the current point
1180: . aOff - the offset into a[] and a_t[] for each auxiliary field
1181: . aOff_x - the offset into a_x[] for each auxiliary field
1182: . a - each auxiliary field evaluated at the current point
1183: . a_t - the time derivative of each auxiliary field evaluated at the current point
1184: . a_x - the gradient of auxiliary each field evaluated at the current point
1185: . t - current time
1186: . x - coordinates of the current point
1187: . numConstants - number of constant parameters
1188: . constants - constant parameters
1189: - f0 - output values at the current point

1191:   Level: intermediate

1193: .seealso: PetscDSGetResidual()
1194: @*/
1195: PetscErrorCode PetscDSSetResidual(PetscDS prob, PetscInt f,
1196:                                   void (*f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1197:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1198:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1199:                                              PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]),
1200:                                   void (*f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1201:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1202:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1203:                                              PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
1204: {

1211:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1212:   PetscDSEnlarge_Static(prob, f+1);
1213:   prob->f[f*2+0] = f0;
1214:   prob->f[f*2+1] = f1;
1215:   return(0);
1216: }

1218: /*@C
1219:   PetscDSHasJacobian - Signals that Jacobian functions have been set

1221:   Not collective

1223:   Input Parameter:
1224: . prob - The PetscDS

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

1229:   Level: intermediate

1231: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1232: @*/
1233: PetscErrorCode PetscDSHasJacobian(PetscDS prob, PetscBool *hasJac)
1234: {
1235:   PetscInt f, g, h;

1239:   *hasJac = PETSC_FALSE;
1240:   for (f = 0; f < prob->Nf; ++f) {
1241:     for (g = 0; g < prob->Nf; ++g) {
1242:       for (h = 0; h < 4; ++h) {
1243:         if (prob->g[(f*prob->Nf + g)*4+h]) *hasJac = PETSC_TRUE;
1244:       }
1245:     }
1246:   }
1247:   return(0);
1248: }

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

1253:   Not collective

1255:   Input Parameters:
1256: + prob - The PetscDS
1257: . f    - The test field number
1258: - g    - The field number

1260:   Output Parameters:
1261: + g0 - integrand for the test and basis function term
1262: . g1 - integrand for the test function and basis function gradient term
1263: . g2 - integrand for the test function gradient and basis function term
1264: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

1272: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1273: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1274: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1275: $    PetscReal t, const PetscReal u_tShift, const PetscReal x[], PetscScalar g0[])

1277: + dim - the spatial dimension
1278: . Nf - the number of fields
1279: . uOff - the offset into u[] and u_t[] for each field
1280: . uOff_x - the offset into u_x[] for each field
1281: . u - each field evaluated at the current point
1282: . u_t - the time derivative of each field evaluated at the current point
1283: . u_x - the gradient of each field evaluated at the current point
1284: . aOff - the offset into a[] and a_t[] for each auxiliary field
1285: . aOff_x - the offset into a_x[] for each auxiliary field
1286: . a - each auxiliary field evaluated at the current point
1287: . a_t - the time derivative of each auxiliary field evaluated at the current point
1288: . a_x - the gradient of auxiliary each field evaluated at the current point
1289: . t - current time
1290: . u_tShift - the multiplier a for dF/dU_t
1291: . x - coordinates of the current point
1292: . numConstants - number of constant parameters
1293: . constants - constant parameters
1294: - g0 - output values at the current point

1296:   Level: intermediate

1298: .seealso: PetscDSSetJacobian()
1299: @*/
1300: PetscErrorCode PetscDSGetJacobian(PetscDS prob, PetscInt f, PetscInt g,
1301:                                   void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1302:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1303:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1304:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1305:                                   void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1306:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1307:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1308:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1309:                                   void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1310:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1311:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1312:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1313:                                   void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1314:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1315:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1316:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1317: {
1320:   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);
1321:   if ((g < 0) || (g >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, prob->Nf);
1326:   return(0);
1327: }

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

1332:   Not collective

1334:   Input Parameters:
1335: + prob - The PetscDS
1336: . f    - The test field number
1337: . g    - The field number
1338: . g0 - integrand for the test and basis function term
1339: . g1 - integrand for the test function and basis function gradient term
1340: . g2 - integrand for the test function gradient and basis function term
1341: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

1349: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1350: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1351: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1352: $    PetscReal t, const PetscReal x[], PetscScalar g0[])

1354: + dim - the spatial dimension
1355: . Nf - the number of fields
1356: . uOff - the offset into u[] and u_t[] for each field
1357: . uOff_x - the offset into u_x[] for each field
1358: . u - each field evaluated at the current point
1359: . u_t - the time derivative of each field evaluated at the current point
1360: . u_x - the gradient of each field evaluated at the current point
1361: . aOff - the offset into a[] and a_t[] for each auxiliary field
1362: . aOff_x - the offset into a_x[] for each auxiliary field
1363: . a - each auxiliary field evaluated at the current point
1364: . a_t - the time derivative of each auxiliary field evaluated at the current point
1365: . a_x - the gradient of auxiliary each field evaluated at the current point
1366: . t - current time
1367: . u_tShift - the multiplier a for dF/dU_t
1368: . x - coordinates of the current point
1369: . numConstants - number of constant parameters
1370: . constants - constant parameters
1371: - g0 - output values at the current point

1373:   Level: intermediate

1375: .seealso: PetscDSGetJacobian()
1376: @*/
1377: PetscErrorCode PetscDSSetJacobian(PetscDS prob, PetscInt f, PetscInt g,
1378:                                   void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1379:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1380:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1381:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1382:                                   void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1383:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1384:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1385:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1386:                                   void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1387:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1388:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1389:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1390:                                   void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1391:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1392:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1393:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1394: {

1403:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1404:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
1405:   PetscDSEnlarge_Static(prob, PetscMax(f, g)+1);
1406:   prob->g[(f*prob->Nf + g)*4+0] = g0;
1407:   prob->g[(f*prob->Nf + g)*4+1] = g1;
1408:   prob->g[(f*prob->Nf + g)*4+2] = g2;
1409:   prob->g[(f*prob->Nf + g)*4+3] = g3;
1410:   return(0);
1411: }

1413: /*@C
1414:   PetscDSUseJacobianPreconditioner - Whether to construct a Jacobian preconditioner

1416:   Not collective

1418:   Input Parameters:
1419: + prob - The PetscDS
1420: - useJacPre - flag that enables construction of a Jacobian preconditioner

1422:   Level: intermediate

1424: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1425: @*/
1426: PetscErrorCode PetscDSUseJacobianPreconditioner(PetscDS prob, PetscBool useJacPre)
1427: {
1430:   prob->useJacPre = useJacPre;
1431:   return(0);
1432: }

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

1437:   Not collective

1439:   Input Parameter:
1440: . prob - The PetscDS

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

1445:   Level: intermediate

1447: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1448: @*/
1449: PetscErrorCode PetscDSHasJacobianPreconditioner(PetscDS prob, PetscBool *hasJacPre)
1450: {
1451:   PetscInt f, g, h;

1455:   *hasJacPre = PETSC_FALSE;
1456:   if (!prob->useJacPre) return(0);
1457:   for (f = 0; f < prob->Nf; ++f) {
1458:     for (g = 0; g < prob->Nf; ++g) {
1459:       for (h = 0; h < 4; ++h) {
1460:         if (prob->gp[(f*prob->Nf + g)*4+h]) *hasJacPre = PETSC_TRUE;
1461:       }
1462:     }
1463:   }
1464:   return(0);
1465: }

1467: /*@C
1468:   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.

1470:   Not collective

1472:   Input Parameters:
1473: + prob - The PetscDS
1474: . f    - The test field number
1475: - g    - The field number

1477:   Output Parameters:
1478: + g0 - integrand for the test and basis function term
1479: . g1 - integrand for the test function and basis function gradient term
1480: . g2 - integrand for the test function gradient and basis function term
1481: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

1489: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1490: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1491: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1492: $    PetscReal t, const PetscReal u_tShift, const PetscReal x[], PetscScalar g0[])

1494: + dim - the spatial dimension
1495: . Nf - the number of fields
1496: . uOff - the offset into u[] and u_t[] for each field
1497: . uOff_x - the offset into u_x[] for each field
1498: . u - each field evaluated at the current point
1499: . u_t - the time derivative of each field evaluated at the current point
1500: . u_x - the gradient of each field evaluated at the current point
1501: . aOff - the offset into a[] and a_t[] for each auxiliary field
1502: . aOff_x - the offset into a_x[] for each auxiliary field
1503: . a - each auxiliary field evaluated at the current point
1504: . a_t - the time derivative of each auxiliary field evaluated at the current point
1505: . a_x - the gradient of auxiliary each field evaluated at the current point
1506: . t - current time
1507: . u_tShift - the multiplier a for dF/dU_t
1508: . x - coordinates of the current point
1509: . numConstants - number of constant parameters
1510: . constants - constant parameters
1511: - g0 - output values at the current point

1513:   Level: intermediate

1515: .seealso: PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1516: @*/
1517: PetscErrorCode PetscDSGetJacobianPreconditioner(PetscDS prob, PetscInt f, PetscInt g,
1518:                                   void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1519:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1520:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1521:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1522:                                   void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1523:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1524:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1525:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1526:                                   void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1527:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1528:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1529:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1530:                                   void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1531:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1532:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1533:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1534: {
1537:   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);
1538:   if ((g < 0) || (g >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, prob->Nf);
1543:   return(0);
1544: }

1546: /*@C
1547:   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.

1549:   Not collective

1551:   Input Parameters:
1552: + prob - The PetscDS
1553: . f    - The test field number
1554: . g    - The field number
1555: . g0 - integrand for the test and basis function term
1556: . g1 - integrand for the test function and basis function gradient term
1557: . g2 - integrand for the test function gradient and basis function term
1558: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

1566: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1567: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1568: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1569: $    PetscReal t, const PetscReal x[], PetscScalar g0[])

1571: + dim - the spatial dimension
1572: . Nf - the number of fields
1573: . uOff - the offset into u[] and u_t[] for each field
1574: . uOff_x - the offset into u_x[] for each field
1575: . u - each field evaluated at the current point
1576: . u_t - the time derivative of each field evaluated at the current point
1577: . u_x - the gradient of each field evaluated at the current point
1578: . aOff - the offset into a[] and a_t[] for each auxiliary field
1579: . aOff_x - the offset into a_x[] for each auxiliary field
1580: . a - each auxiliary field evaluated at the current point
1581: . a_t - the time derivative of each auxiliary field evaluated at the current point
1582: . a_x - the gradient of auxiliary each field evaluated at the current point
1583: . t - current time
1584: . u_tShift - the multiplier a for dF/dU_t
1585: . x - coordinates of the current point
1586: . numConstants - number of constant parameters
1587: . constants - constant parameters
1588: - g0 - output values at the current point

1590:   Level: intermediate

1592: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobian()
1593: @*/
1594: PetscErrorCode PetscDSSetJacobianPreconditioner(PetscDS prob, PetscInt f, PetscInt g,
1595:                                   void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1596:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1597:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1598:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1599:                                   void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1600:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1601:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1602:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1603:                                   void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1604:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1605:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1606:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1607:                                   void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1608:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1609:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1610:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1611: {

1620:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1621:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
1622:   PetscDSEnlarge_Static(prob, PetscMax(f, g)+1);
1623:   prob->gp[(f*prob->Nf + g)*4+0] = g0;
1624:   prob->gp[(f*prob->Nf + g)*4+1] = g1;
1625:   prob->gp[(f*prob->Nf + g)*4+2] = g2;
1626:   prob->gp[(f*prob->Nf + g)*4+3] = g3;
1627:   return(0);
1628: }

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

1633:   Not collective

1635:   Input Parameter:
1636: . prob - The PetscDS

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

1641:   Level: intermediate

1643: .seealso: PetscDSGetDynamicJacobian(), PetscDSSetDynamicJacobian(), PetscDSGetJacobian()
1644: @*/
1645: PetscErrorCode PetscDSHasDynamicJacobian(PetscDS prob, PetscBool *hasDynJac)
1646: {
1647:   PetscInt f, g, h;

1651:   *hasDynJac = PETSC_FALSE;
1652:   for (f = 0; f < prob->Nf; ++f) {
1653:     for (g = 0; g < prob->Nf; ++g) {
1654:       for (h = 0; h < 4; ++h) {
1655:         if (prob->gt[(f*prob->Nf + g)*4+h]) *hasDynJac = PETSC_TRUE;
1656:       }
1657:     }
1658:   }
1659:   return(0);
1660: }

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

1665:   Not collective

1667:   Input Parameters:
1668: + prob - The PetscDS
1669: . f    - The test field number
1670: - g    - The field number

1672:   Output Parameters:
1673: + g0 - integrand for the test and basis function term
1674: . g1 - integrand for the test function and basis function gradient term
1675: . g2 - integrand for the test function gradient and basis function term
1676: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

1684: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1685: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1686: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1687: $    PetscReal t, const PetscReal u_tShift, const PetscReal x[], PetscScalar g0[])

1689: + dim - the spatial dimension
1690: . Nf - the number of fields
1691: . uOff - the offset into u[] and u_t[] for each field
1692: . uOff_x - the offset into u_x[] for each field
1693: . u - each field evaluated at the current point
1694: . u_t - the time derivative of each field evaluated at the current point
1695: . u_x - the gradient of each field evaluated at the current point
1696: . aOff - the offset into a[] and a_t[] for each auxiliary field
1697: . aOff_x - the offset into a_x[] for each auxiliary field
1698: . a - each auxiliary field evaluated at the current point
1699: . a_t - the time derivative of each auxiliary field evaluated at the current point
1700: . a_x - the gradient of auxiliary each field evaluated at the current point
1701: . t - current time
1702: . u_tShift - the multiplier a for dF/dU_t
1703: . x - coordinates of the current point
1704: . numConstants - number of constant parameters
1705: . constants - constant parameters
1706: - g0 - output values at the current point

1708:   Level: intermediate

1710: .seealso: PetscDSSetJacobian()
1711: @*/
1712: PetscErrorCode PetscDSGetDynamicJacobian(PetscDS prob, PetscInt f, PetscInt g,
1713:                                          void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1714:                                                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1715:                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1716:                                                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1717:                                          void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1718:                                                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1719:                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1720:                                                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1721:                                          void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1722:                                                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1723:                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1724:                                                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1725:                                          void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1726:                                                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1727:                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1728:                                                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1729: {
1732:   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);
1733:   if ((g < 0) || (g >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, prob->Nf);
1738:   return(0);
1739: }

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

1744:   Not collective

1746:   Input Parameters:
1747: + prob - The PetscDS
1748: . f    - The test field number
1749: . g    - The field number
1750: . g0 - integrand for the test and basis function term
1751: . g1 - integrand for the test function and basis function gradient term
1752: . g2 - integrand for the test function gradient and basis function term
1753: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

1761: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1762: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1763: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1764: $    PetscReal t, const PetscReal x[], PetscScalar g0[])

1766: + dim - the spatial dimension
1767: . Nf - the number of fields
1768: . uOff - the offset into u[] and u_t[] for each field
1769: . uOff_x - the offset into u_x[] for each field
1770: . u - each field evaluated at the current point
1771: . u_t - the time derivative of each field evaluated at the current point
1772: . u_x - the gradient of each field evaluated at the current point
1773: . aOff - the offset into a[] and a_t[] for each auxiliary field
1774: . aOff_x - the offset into a_x[] for each auxiliary field
1775: . a - each auxiliary field evaluated at the current point
1776: . a_t - the time derivative of each auxiliary field evaluated at the current point
1777: . a_x - the gradient of auxiliary each field evaluated at the current point
1778: . t - current time
1779: . u_tShift - the multiplier a for dF/dU_t
1780: . x - coordinates of the current point
1781: . numConstants - number of constant parameters
1782: . constants - constant parameters
1783: - g0 - output values at the current point

1785:   Level: intermediate

1787: .seealso: PetscDSGetJacobian()
1788: @*/
1789: PetscErrorCode PetscDSSetDynamicJacobian(PetscDS prob, PetscInt f, PetscInt g,
1790:                                          void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1791:                                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1792:                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1793:                                                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1794:                                          void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1795:                                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1796:                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1797:                                                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1798:                                          void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1799:                                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1800:                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1801:                                                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1802:                                          void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1803:                                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1804:                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1805:                                                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1806: {

1815:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1816:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
1817:   PetscDSEnlarge_Static(prob, PetscMax(f, g)+1);
1818:   prob->gt[(f*prob->Nf + g)*4+0] = g0;
1819:   prob->gt[(f*prob->Nf + g)*4+1] = g1;
1820:   prob->gt[(f*prob->Nf + g)*4+2] = g2;
1821:   prob->gt[(f*prob->Nf + g)*4+3] = g3;
1822:   return(0);
1823: }

1825: /*@C
1826:   PetscDSGetRiemannSolver - Returns the Riemann solver for the given field

1828:   Not collective

1830:   Input Arguments:
1831: + prob - The PetscDS object
1832: - f    - The field number

1834:   Output Argument:
1835: . r    - Riemann solver

1837:   Calling sequence for r:

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

1841: + dim  - The spatial dimension
1842: . Nf   - The number of fields
1843: . x    - The coordinates at a point on the interface
1844: . n    - The normal vector to the interface
1845: . uL   - The state vector to the left of the interface
1846: . uR   - The state vector to the right of the interface
1847: . flux - output array of flux through the interface
1848: . numConstants - number of constant parameters
1849: . constants - constant parameters
1850: - ctx  - optional user context

1852:   Level: intermediate

1854: .seealso: PetscDSSetRiemannSolver()
1855: @*/
1856: PetscErrorCode PetscDSGetRiemannSolver(PetscDS prob, PetscInt f,
1857:                                        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))
1858: {
1861:   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);
1863:   *r = prob->r[f];
1864:   return(0);
1865: }

1867: /*@C
1868:   PetscDSSetRiemannSolver - Sets the Riemann solver for the given field

1870:   Not collective

1872:   Input Arguments:
1873: + prob - The PetscDS object
1874: . f    - The field number
1875: - r    - Riemann solver

1877:   Calling sequence for r:

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

1881: + dim  - The spatial dimension
1882: . Nf   - The number of fields
1883: . x    - The coordinates at a point on the interface
1884: . n    - The normal vector to the interface
1885: . uL   - The state vector to the left of the interface
1886: . uR   - The state vector to the right of the interface
1887: . flux - output array of flux through the interface
1888: . numConstants - number of constant parameters
1889: . constants - constant parameters
1890: - ctx  - optional user context

1892:   Level: intermediate

1894: .seealso: PetscDSGetRiemannSolver()
1895: @*/
1896: PetscErrorCode PetscDSSetRiemannSolver(PetscDS prob, PetscInt f,
1897:                                        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))
1898: {

1904:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1905:   PetscDSEnlarge_Static(prob, f+1);
1906:   prob->r[f] = r;
1907:   return(0);
1908: }

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

1913:   Not collective

1915:   Input Parameters:
1916: + prob - The PetscDS
1917: - f    - The field number

1919:   Output Parameters:
1920: . update - update function

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

1924: $ update(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1925: $        const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1926: $        const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1927: $        PetscReal t, const PetscReal x[], PetscScalar uNew[])

1929: + dim - the spatial dimension
1930: . Nf - the number of fields
1931: . uOff - the offset into u[] and u_t[] for each field
1932: . uOff_x - the offset into u_x[] for each field
1933: . u - each field evaluated at the current point
1934: . u_t - the time derivative of each field evaluated at the current point
1935: . u_x - the gradient of each field evaluated at the current point
1936: . aOff - the offset into a[] and a_t[] for each auxiliary field
1937: . aOff_x - the offset into a_x[] for each auxiliary field
1938: . a - each auxiliary field evaluated at the current point
1939: . a_t - the time derivative of each auxiliary field evaluated at the current point
1940: . a_x - the gradient of auxiliary each field evaluated at the current point
1941: . t - current time
1942: . x - coordinates of the current point
1943: - uNew - new value for field at the current point

1945:   Level: intermediate

1947: .seealso: PetscDSSetUpdate(), PetscDSSetResidual()
1948: @*/
1949: PetscErrorCode PetscDSGetUpdate(PetscDS prob, PetscInt f,
1950:                                   void (**update)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1951:                                                   const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1952:                                                   const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1953:                                                   PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar uNew[]))
1954: {
1957:   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
1959:   return(0);
1960: }

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

1965:   Not collective

1967:   Input Parameters:
1968: + prob   - The PetscDS
1969: . f      - The field number
1970: - update - update function

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

1974: $ update(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1975: $        const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1976: $        const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1977: $        PetscReal t, const PetscReal x[], PetscScalar uNew[])

1979: + dim - the spatial dimension
1980: . Nf - the number of fields
1981: . uOff - the offset into u[] and u_t[] for each field
1982: . uOff_x - the offset into u_x[] for each field
1983: . u - each field evaluated at the current point
1984: . u_t - the time derivative of each field evaluated at the current point
1985: . u_x - the gradient of each field evaluated at the current point
1986: . aOff - the offset into a[] and a_t[] for each auxiliary field
1987: . aOff_x - the offset into a_x[] for each auxiliary field
1988: . a - each auxiliary field evaluated at the current point
1989: . a_t - the time derivative of each auxiliary field evaluated at the current point
1990: . a_x - the gradient of auxiliary each field evaluated at the current point
1991: . t - current time
1992: . x - coordinates of the current point
1993: - uNew - new field values at the current point

1995:   Level: intermediate

1997: .seealso: PetscDSGetResidual()
1998: @*/
1999: PetscErrorCode PetscDSSetUpdate(PetscDS prob, PetscInt f,
2000:                                 void (*update)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2001:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2002:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2003:                                                PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar uNew[]))
2004: {

2010:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2011:   PetscDSEnlarge_Static(prob, f+1);
2012:   prob->update[f] = update;
2013:   return(0);
2014: }

2016: PetscErrorCode PetscDSGetContext(PetscDS prob, PetscInt f, void **ctx)
2017: {
2020:   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);
2022:   *ctx = prob->ctx[f];
2023:   return(0);
2024: }

2026: PetscErrorCode PetscDSSetContext(PetscDS prob, PetscInt f, void *ctx)
2027: {

2032:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2033:   PetscDSEnlarge_Static(prob, f+1);
2034:   prob->ctx[f] = ctx;
2035:   return(0);
2036: }

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

2041:   Not collective

2043:   Input Parameters:
2044: + prob - The PetscDS
2045: - f    - The test field number

2047:   Output Parameters:
2048: + f0 - boundary integrand for the test function term
2049: - f1 - boundary integrand for the test function gradient term

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

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

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

2057: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2058: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2059: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2060: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[])

2062: + dim - the spatial dimension
2063: . Nf - the number of fields
2064: . uOff - the offset into u[] and u_t[] for each field
2065: . uOff_x - the offset into u_x[] for each field
2066: . u - each field evaluated at the current point
2067: . u_t - the time derivative of each field evaluated at the current point
2068: . u_x - the gradient of each field evaluated at the current point
2069: . aOff - the offset into a[] and a_t[] for each auxiliary field
2070: . aOff_x - the offset into a_x[] for each auxiliary field
2071: . a - each auxiliary field evaluated at the current point
2072: . a_t - the time derivative of each auxiliary field evaluated at the current point
2073: . a_x - the gradient of auxiliary each field evaluated at the current point
2074: . t - current time
2075: . x - coordinates of the current point
2076: . n - unit normal at the current point
2077: . numConstants - number of constant parameters
2078: . constants - constant parameters
2079: - f0 - output values at the current point

2081:   Level: intermediate

2083: .seealso: PetscDSSetBdResidual()
2084: @*/
2085: PetscErrorCode PetscDSGetBdResidual(PetscDS prob, PetscInt f,
2086:                                     void (**f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2087:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2088:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2089:                                                 PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]),
2090:                                     void (**f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2091:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2092:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2093:                                                 PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
2094: {
2097:   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);
2100:   return(0);
2101: }

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

2106:   Not collective

2108:   Input Parameters:
2109: + prob - The PetscDS
2110: . f    - The test field number
2111: . f0 - boundary integrand for the test function term
2112: - f1 - boundary integrand for the test function gradient term

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

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

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

2120: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2121: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2122: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2123: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[])

2125: + dim - the spatial dimension
2126: . Nf - the number of fields
2127: . uOff - the offset into u[] and u_t[] for each field
2128: . uOff_x - the offset into u_x[] for each field
2129: . u - each field evaluated at the current point
2130: . u_t - the time derivative of each field evaluated at the current point
2131: . u_x - the gradient of each field evaluated at the current point
2132: . aOff - the offset into a[] and a_t[] for each auxiliary field
2133: . aOff_x - the offset into a_x[] for each auxiliary field
2134: . a - each auxiliary field evaluated at the current point
2135: . a_t - the time derivative of each auxiliary field evaluated at the current point
2136: . a_x - the gradient of auxiliary each field evaluated at the current point
2137: . t - current time
2138: . x - coordinates of the current point
2139: . n - unit normal at the current point
2140: . numConstants - number of constant parameters
2141: . constants - constant parameters
2142: - f0 - output values at the current point

2144:   Level: intermediate

2146: .seealso: PetscDSGetBdResidual()
2147: @*/
2148: PetscErrorCode PetscDSSetBdResidual(PetscDS prob, PetscInt f,
2149:                                     void (*f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2150:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2151:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2152:                                                PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]),
2153:                                     void (*f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2154:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2155:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2156:                                                PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
2157: {

2162:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2163:   PetscDSEnlarge_Static(prob, f+1);
2166:   return(0);
2167: }

2169: /*@
2170:   PetscDSHasBdJacobian - Signals that boundary Jacobian functions have been set

2172:   Not collective

2174:   Input Parameter:
2175: . prob - The PetscDS

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

2180:   Level: intermediate

2182: .seealso: PetscDSHasJacobian(), PetscDSSetBdJacobian(), PetscDSGetBdJacobian()
2183: @*/
2184: PetscErrorCode PetscDSHasBdJacobian(PetscDS prob, PetscBool *hasBdJac)
2185: {
2186:   PetscInt f, g, h;

2190:   *hasBdJac = PETSC_FALSE;
2191:   for (f = 0; f < prob->Nf; ++f) {
2192:     for (g = 0; g < prob->Nf; ++g) {
2193:       for (h = 0; h < 4; ++h) {
2194:         if (prob->gBd[(f*prob->Nf + g)*4+h]) *hasBdJac = PETSC_TRUE;
2195:       }
2196:     }
2197:   }
2198:   return(0);
2199: }

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

2204:   Not collective

2206:   Input Parameters:
2207: + prob - The PetscDS
2208: . f    - The test field number
2209: - g    - The field number

2211:   Output Parameters:
2212: + g0 - integrand for the test and basis function term
2213: . g1 - integrand for the test function and basis function gradient term
2214: . g2 - integrand for the test function gradient and basis function term
2215: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

2223: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2224: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2225: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2226: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar g0[])

2228: + dim - the spatial dimension
2229: . Nf - the number of fields
2230: . uOff - the offset into u[] and u_t[] for each field
2231: . uOff_x - the offset into u_x[] for each field
2232: . u - each field evaluated at the current point
2233: . u_t - the time derivative of each field evaluated at the current point
2234: . u_x - the gradient of each field evaluated at the current point
2235: . aOff - the offset into a[] and a_t[] for each auxiliary field
2236: . aOff_x - the offset into a_x[] for each auxiliary field
2237: . a - each auxiliary field evaluated at the current point
2238: . a_t - the time derivative of each auxiliary field evaluated at the current point
2239: . a_x - the gradient of auxiliary each field evaluated at the current point
2240: . t - current time
2241: . u_tShift - the multiplier a for dF/dU_t
2242: . x - coordinates of the current point
2243: . n - normal at the current point
2244: . numConstants - number of constant parameters
2245: . constants - constant parameters
2246: - g0 - output values at the current point

2248:   Level: intermediate

2250: .seealso: PetscDSSetBdJacobian()
2251: @*/
2252: PetscErrorCode PetscDSGetBdJacobian(PetscDS prob, PetscInt f, PetscInt g,
2253:                                     void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2254:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2255:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2256:                                                 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
2257:                                     void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2258:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2259:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2260:                                                 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
2261:                                     void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2262:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2263:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2264:                                                 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
2265:                                     void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2266:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2267:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2268:                                                 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
2269: {
2272:   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);
2273:   if ((g < 0) || (g >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, prob->Nf);
2278:   return(0);
2279: }

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

2284:   Not collective

2286:   Input Parameters:
2287: + prob - The PetscDS
2288: . f    - The test field number
2289: . g    - The field number
2290: . g0 - integrand for the test and basis function term
2291: . g1 - integrand for the test function and basis function gradient term
2292: . g2 - integrand for the test function gradient and basis function term
2293: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

2301: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2302: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2303: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2304: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar g0[])

2306: + dim - the spatial dimension
2307: . Nf - the number of fields
2308: . uOff - the offset into u[] and u_t[] for each field
2309: . uOff_x - the offset into u_x[] for each field
2310: . u - each field evaluated at the current point
2311: . u_t - the time derivative of each field evaluated at the current point
2312: . u_x - the gradient of each field evaluated at the current point
2313: . aOff - the offset into a[] and a_t[] for each auxiliary field
2314: . aOff_x - the offset into a_x[] for each auxiliary field
2315: . a - each auxiliary field evaluated at the current point
2316: . a_t - the time derivative of each auxiliary field evaluated at the current point
2317: . a_x - the gradient of auxiliary each field evaluated at the current point
2318: . t - current time
2319: . u_tShift - the multiplier a for dF/dU_t
2320: . x - coordinates of the current point
2321: . n - normal at the current point
2322: . numConstants - number of constant parameters
2323: . constants - constant parameters
2324: - g0 - output values at the current point

2326:   Level: intermediate

2328: .seealso: PetscDSGetBdJacobian()
2329: @*/
2330: PetscErrorCode PetscDSSetBdJacobian(PetscDS prob, PetscInt f, PetscInt g,
2331:                                     void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2332:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2333:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2334:                                                PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
2335:                                     void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2336:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2337:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2338:                                                PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
2339:                                     void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2340:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2341:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2342:                                                PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
2343:                                     void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2344:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2345:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2346:                                                PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
2347: {

2356:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2357:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
2358:   PetscDSEnlarge_Static(prob, PetscMax(f, g)+1);
2359:   prob->gBd[(f*prob->Nf + g)*4+0] = g0;
2360:   prob->gBd[(f*prob->Nf + g)*4+1] = g1;
2361:   prob->gBd[(f*prob->Nf + g)*4+2] = g2;
2362:   prob->gBd[(f*prob->Nf + g)*4+3] = g3;
2363:   return(0);
2364: }

2366: /*@
2367:   PetscDSHasBdJacobianPreconditioner - Signals that boundary Jacobian preconditioner functions have been set

2369:   Not collective

2371:   Input Parameter:
2372: . prob - The PetscDS

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

2377:   Level: intermediate

2379: .seealso: PetscDSHasJacobian(), PetscDSSetBdJacobian(), PetscDSGetBdJacobian()
2380: @*/
2381: PetscErrorCode PetscDSHasBdJacobianPreconditioner(PetscDS prob, PetscBool *hasBdJacPre)
2382: {
2383:   PetscInt f, g, h;

2387:   *hasBdJacPre = PETSC_FALSE;
2388:   for (f = 0; f < prob->Nf; ++f) {
2389:     for (g = 0; g < prob->Nf; ++g) {
2390:       for (h = 0; h < 4; ++h) {
2391:         if (prob->gpBd[(f*prob->Nf + g)*4+h]) *hasBdJacPre = PETSC_TRUE;
2392:       }
2393:     }
2394:   }
2395:   return(0);
2396: }

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

2401:   Not collective

2403:   Input Parameters:
2404: + prob - The PetscDS
2405: . f    - The test field number
2406: - g    - The field number

2408:   Output Parameters:
2409: + g0 - integrand for the test and basis function term
2410: . g1 - integrand for the test function and basis function gradient term
2411: . g2 - integrand for the test function gradient and basis function term
2412: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

2420: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2421: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2422: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2423: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])

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

2446:   This is not yet available in Fortran.

2448:   Level: intermediate

2450: .seealso: PetscDSSetBdJacobianPreconditioner()
2451: @*/
2452: PetscErrorCode PetscDSGetBdJacobianPreconditioner(PetscDS prob, PetscInt f, PetscInt g,
2453:                                                   void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2454:                                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2455:                                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2456:                                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
2457:                                                   void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2458:                                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2459:                                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2460:                                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
2461:                                                   void (**g2)(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 g2[]),
2465:                                                   void (**g3)(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 g3[]))
2469: {
2472:   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);
2473:   if ((g < 0) || (g >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, prob->Nf);
2478:   return(0);
2479: }

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

2484:   Not collective

2486:   Input Parameters:
2487: + prob - The PetscDS
2488: . f    - The test field number
2489: . g    - The field number
2490: . g0 - integrand for the test and basis function term
2491: . g1 - integrand for the test function and basis function gradient term
2492: . g2 - integrand for the test function gradient and basis function term
2493: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

2501: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2502: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2503: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2504: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])

2506: + dim - the spatial dimension
2507: . Nf - the number of fields
2508: . NfAux - the number of auxiliary fields
2509: . uOff - the offset into u[] and u_t[] for each field
2510: . uOff_x - the offset into u_x[] for each field
2511: . u - each field evaluated at the current point
2512: . u_t - the time derivative of each field evaluated at the current point
2513: . u_x - the gradient of each field evaluated at the current point
2514: . aOff - the offset into a[] and a_t[] for each auxiliary field
2515: . aOff_x - the offset into a_x[] for each auxiliary field
2516: . a - each auxiliary field evaluated at the current point
2517: . a_t - the time derivative of each auxiliary field evaluated at the current point
2518: . a_x - the gradient of auxiliary each field evaluated at the current point
2519: . t - current time
2520: . u_tShift - the multiplier a for dF/dU_t
2521: . x - coordinates of the current point
2522: . n - normal at the current point
2523: . numConstants - number of constant parameters
2524: . constants - constant parameters
2525: - g0 - output values at the current point

2527:   This is not yet available in Fortran.

2529:   Level: intermediate

2531: .seealso: PetscDSGetBdJacobianPreconditioner()
2532: @*/
2533: PetscErrorCode PetscDSSetBdJacobianPreconditioner(PetscDS prob, PetscInt f, PetscInt g,
2534:                                                   void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2535:                                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2536:                                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2537:                                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
2538:                                                   void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2539:                                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2540:                                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2541:                                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
2542:                                                   void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2543:                                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2544:                                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2545:                                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
2546:                                                   void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2547:                                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2548:                                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2549:                                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
2550: {

2559:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2560:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
2561:   PetscDSEnlarge_Static(prob, PetscMax(f, g)+1);
2562:   prob->gpBd[(f*prob->Nf + g)*4+0] = g0;
2563:   prob->gpBd[(f*prob->Nf + g)*4+1] = g1;
2564:   prob->gpBd[(f*prob->Nf + g)*4+2] = g2;
2565:   prob->gpBd[(f*prob->Nf + g)*4+3] = g3;
2566:   return(0);
2567: }

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

2572:   Not collective

2574:   Input Parameters:
2575: + prob - The PetscDS
2576: - f    - The test field number

2578:   Output Parameter:
2579: + exactSol - exact solution for the test field
2580: - exactCtx - exact solution context

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

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

2586: + dim - the spatial dimension
2587: . t - current time
2588: . x - coordinates of the current point
2589: . Nc - the number of field components
2590: . u - the solution field evaluated at the current point
2591: - ctx - a user context

2593:   Level: intermediate

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

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

2610:   Not collective

2612:   Input Parameters:
2613: + prob - The PetscDS
2614: . f    - The test field number
2615: . sol  - solution function for the test fields
2616: - ctx  - solution context or NULL

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

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

2622: + dim - the spatial dimension
2623: . t - current time
2624: . x - coordinates of the current point
2625: . Nc - the number of field components
2626: . u - the solution field evaluated at the current point
2627: - ctx - a user context

2629:   Level: intermediate

2631: .seealso: PetscDSGetExactSolution()
2632: @*/
2633: PetscErrorCode PetscDSSetExactSolution(PetscDS prob, PetscInt f, PetscErrorCode (*sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void *ctx)
2634: {

2639:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2640:   PetscDSEnlarge_Static(prob, f+1);
2643:   return(0);
2644: }

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

2649:   Not collective

2651:   Input Parameters:
2652: + prob - The PetscDS
2653: - f    - The test field number

2655:   Output Parameter:
2656: + exactSol - time derivative of the exact solution for the test field
2657: - exactCtx - time derivative of the exact solution context

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

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

2663: + dim - the spatial dimension
2664: . t - current time
2665: . x - coordinates of the current point
2666: . Nc - the number of field components
2667: . u - the solution field evaluated at the current point
2668: - ctx - a user context

2670:   Level: intermediate

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

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

2687:   Not collective

2689:   Input Parameters:
2690: + prob - The PetscDS
2691: . f    - The test field number
2692: . sol  - time derivative of the solution function for the test fields
2693: - ctx  - time derivative of the solution context or NULL

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

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

2699: + dim - the spatial dimension
2700: . t - current time
2701: . x - coordinates of the current point
2702: . Nc - the number of field components
2703: . u - the solution field evaluated at the current point
2704: - ctx - a user context

2706:   Level: intermediate

2708: .seealso: PetscDSGetExactSolutionTimeDerivative(), PetscDSSetExactSolution()
2709: @*/
2710: PetscErrorCode PetscDSSetExactSolutionTimeDerivative(PetscDS prob, PetscInt f, PetscErrorCode (*sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void *ctx)
2711: {

2716:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2717:   PetscDSEnlarge_Static(prob, f+1);
2720:   return(0);
2721: }

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

2726:   Not collective

2728:   Input Parameter:
2729: . prob - The PetscDS object

2731:   Output Parameters:
2732: + numConstants - The number of constants
2733: - constants    - The array of constants, NULL if there are none

2735:   Level: intermediate

2737: .seealso: PetscDSSetConstants(), PetscDSCreate()
2738: @*/
2739: PetscErrorCode PetscDSGetConstants(PetscDS prob, PetscInt *numConstants, const PetscScalar *constants[])
2740: {
2745:   return(0);
2746: }

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

2751:   Not collective

2753:   Input Parameters:
2754: + prob         - The PetscDS object
2755: . numConstants - The number of constants
2756: - constants    - The array of constants, NULL if there are none

2758:   Level: intermediate

2760: .seealso: PetscDSGetConstants(), PetscDSCreate()
2761: @*/
2762: PetscErrorCode PetscDSSetConstants(PetscDS prob, PetscInt numConstants, PetscScalar constants[])
2763: {

2768:   if (numConstants != prob->numConstants) {
2769:     PetscFree(prob->constants);
2770:     prob->numConstants = numConstants;
2771:     if (prob->numConstants) {
2772:       PetscMalloc1(prob->numConstants, &prob->constants);
2773:     } else {
2774:       prob->constants = NULL;
2775:     }
2776:   }
2777:   if (prob->numConstants) {
2779:     PetscArraycpy(prob->constants, constants, prob->numConstants);
2780:   }
2781:   return(0);
2782: }

2784: /*@
2785:   PetscDSGetFieldIndex - Returns the index of the given field

2787:   Not collective

2789:   Input Parameters:
2790: + prob - The PetscDS object
2791: - disc - The discretization object

2793:   Output Parameter:
2794: . f - The field number

2796:   Level: beginner

2798: .seealso: PetscGetDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
2799: @*/
2800: PetscErrorCode PetscDSGetFieldIndex(PetscDS prob, PetscObject disc, PetscInt *f)
2801: {
2802:   PetscInt g;

2807:   *f = -1;
2808:   for (g = 0; g < prob->Nf; ++g) {if (disc == prob->disc[g]) break;}
2809:   if (g == prob->Nf) SETERRQ(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Field not found in PetscDS.");
2810:   *f = g;
2811:   return(0);
2812: }

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

2817:   Not collective

2819:   Input Parameters:
2820: + prob - The PetscDS object
2821: - f - The field number

2823:   Output Parameter:
2824: . size - The size

2826:   Level: beginner

2828: .seealso: PetscDSGetFieldOffset(), PetscDSGetNumFields(), PetscDSCreate()
2829: @*/
2830: PetscErrorCode PetscDSGetFieldSize(PetscDS prob, PetscInt f, PetscInt *size)
2831: {

2837:   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);
2838:   PetscDSSetUp(prob);
2839:   *size = prob->Nb[f];
2840:   return(0);
2841: }

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

2846:   Not collective

2848:   Input Parameters:
2849: + prob - The PetscDS object
2850: - f - The field number

2852:   Output Parameter:
2853: . off - The offset

2855:   Level: beginner

2857: .seealso: PetscDSGetFieldSize(), PetscDSGetNumFields(), PetscDSCreate()
2858: @*/
2859: PetscErrorCode PetscDSGetFieldOffset(PetscDS prob, PetscInt f, PetscInt *off)
2860: {
2861:   PetscInt       size, g;

2867:   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);
2868:   *off = 0;
2869:   for (g = 0; g < f; ++g) {
2870:     PetscDSGetFieldSize(prob, g, &size);
2871:     *off += size;
2872:   }
2873:   return(0);
2874: }

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

2879:   Not collective

2881:   Input Parameter:
2882: . prob - The PetscDS object

2884:   Output Parameter:
2885: . dimensions - The number of dimensions

2887:   Level: beginner

2889: .seealso: PetscDSGetComponentOffsets(), PetscDSGetNumFields(), PetscDSCreate()
2890: @*/
2891: PetscErrorCode PetscDSGetDimensions(PetscDS prob, PetscInt *dimensions[])
2892: {

2897:   PetscDSSetUp(prob);
2899:   *dimensions = prob->Nb;
2900:   return(0);
2901: }

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

2906:   Not collective

2908:   Input Parameter:
2909: . prob - The PetscDS object

2911:   Output Parameter:
2912: . components - The number of components

2914:   Level: beginner

2916: .seealso: PetscDSGetComponentOffsets(), PetscDSGetNumFields(), PetscDSCreate()
2917: @*/
2918: PetscErrorCode PetscDSGetComponents(PetscDS prob, PetscInt *components[])
2919: {

2924:   PetscDSSetUp(prob);
2926:   *components = prob->Nc;
2927:   return(0);
2928: }

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

2933:   Not collective

2935:   Input Parameters:
2936: + prob - The PetscDS object
2937: - f - The field number

2939:   Output Parameter:
2940: . off - The offset

2942:   Level: beginner

2944: .seealso: PetscDSGetNumFields(), PetscDSCreate()
2945: @*/
2946: PetscErrorCode PetscDSGetComponentOffset(PetscDS prob, PetscInt f, PetscInt *off)
2947: {
2951:   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);
2952:   *off = prob->off[f];
2953:   return(0);
2954: }

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

2959:   Not collective

2961:   Input Parameter:
2962: . prob - The PetscDS object

2964:   Output Parameter:
2965: . offsets - The offsets

2967:   Level: beginner

2969: .seealso: PetscDSGetNumFields(), PetscDSCreate()
2970: @*/
2971: PetscErrorCode PetscDSGetComponentOffsets(PetscDS prob, PetscInt *offsets[])
2972: {
2976:   *offsets = prob->off;
2977:   return(0);
2978: }

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

2983:   Not collective

2985:   Input Parameter:
2986: . prob - The PetscDS object

2988:   Output Parameter:
2989: . offsets - The offsets

2991:   Level: beginner

2993: .seealso: PetscDSGetNumFields(), PetscDSCreate()
2994: @*/
2995: PetscErrorCode PetscDSGetComponentDerivativeOffsets(PetscDS prob, PetscInt *offsets[])
2996: {
3000:   *offsets = prob->offDer;
3001:   return(0);
3002: }

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

3007:   Not collective

3009:   Input Parameter:
3010: . prob - The PetscDS object

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

3015:   Level: intermediate

3017: .seealso: PetscDSCreate()
3018: @*/
3019: PetscErrorCode PetscDSGetTabulation(PetscDS prob, PetscTabulation *T[])
3020: {

3026:   PetscDSSetUp(prob);
3027:   *T = prob->T;
3028:   return(0);
3029: }

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

3034:   Not collective

3036:   Input Parameter:
3037: . prob - The PetscDS object

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

3042:   Level: intermediate

3044: .seealso: PetscDSGetTabulation(), PetscDSCreate()
3045: @*/
3046: PetscErrorCode PetscDSGetFaceTabulation(PetscDS prob, PetscTabulation *Tf[])
3047: {

3053:   PetscDSSetUp(prob);
3054:   *Tf = prob->Tf;
3055:   return(0);
3056: }

3058: PetscErrorCode PetscDSGetEvaluationArrays(PetscDS prob, PetscScalar **u, PetscScalar **u_t, PetscScalar **u_x)
3059: {

3064:   PetscDSSetUp(prob);
3068:   return(0);
3069: }

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

3077:   PetscDSSetUp(prob);
3084:   return(0);
3085: }

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

3093:   PetscDSSetUp(prob);
3099:   return(0);
3100: }

3102: /*@C
3103:   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().

3105:   Collective on ds

3107:   Input Parameters:
3108: + ds          - The PetscDS object
3109: . type        - The type of condition, e.g. DM_BC_ESSENTIAL/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
3110: . name        - The BC name
3111: . labelname   - The label defining constrained points
3112: . field       - The field to constrain
3113: . numcomps    - The number of constrained field components (0 will constrain all fields)
3114: . comps       - An array of constrained component numbers
3115: . bcFunc      - A pointwise function giving boundary values
3116: . bcFunc_t    - A pointwise function giving the time derviative of the boundary values, or NULL
3117: . numids      - The number of DMLabel ids for constrained points
3118: . ids         - An array of ids for constrained points
3119: - ctx         - An optional user context for bcFunc

3121:   Options Database Keys:
3122: + -bc_<boundary name> <num> - Overrides the boundary ids
3123: - -bc_<boundary name>_comp <num> - Overrides the boundary components

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

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

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

3132: $ bcFunc(PetscInt dim, PetscInt Nf, PetscInt NfAux,
3133: $        const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
3134: $        const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
3135: $        PetscReal time, const PetscReal x[], PetscScalar bcval[])

3137: + dim - the spatial dimension
3138: . Nf - the number of fields
3139: . uOff - the offset into u[] and u_t[] for each field
3140: . uOff_x - the offset into u_x[] for each field
3141: . u - each field evaluated at the current point
3142: . u_t - the time derivative of each field evaluated at the current point
3143: . u_x - the gradient of each field evaluated at the current point
3144: . aOff - the offset into a[] and a_t[] for each auxiliary field
3145: . aOff_x - the offset into a_x[] for each auxiliary field
3146: . a - each auxiliary field evaluated at the current point
3147: . a_t - the time derivative of each auxiliary field evaluated at the current point
3148: . a_x - the gradient of auxiliary each field evaluated at the current point
3149: . t - current time
3150: . x - coordinates of the current point
3151: . numConstants - number of constant parameters
3152: . constants - constant parameters
3153: - bcval - output values at the current point

3155:   Level: developer

3157: .seealso: PetscDSGetBoundary(), PetscDSSetResidual(), PetscDSSetBdResidual()
3158: @*/
3159: 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)
3160: {
3161:   DSBoundary     b;

3170:   PetscNew(&b);
3171:   PetscStrallocpy(name, (char **) &b->name);
3172:   PetscStrallocpy(labelname, (char **) &b->labelname);
3173:   PetscMalloc1(numcomps, &b->comps);
3174:   if (numcomps) {PetscArraycpy(b->comps, comps, numcomps);}
3175:   PetscMalloc1(numids, &b->ids);
3176:   if (numids) {PetscArraycpy(b->ids, ids, numids);}
3177:   b->type            = type;
3178:   b->field           = field;
3179:   b->numcomps        = numcomps;
3180:   b->func            = bcFunc;
3181:   b->func_t          = bcFunc_t;
3182:   b->numids          = numids;
3183:   b->ctx             = ctx;
3184:   b->next            = ds->boundary;
3185:   ds->boundary       = b;
3186:   return(0);
3187: }

3189: /*@C
3190:   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().

3192:   Input Parameters:
3193: + ds          - The PetscDS object
3194: . bd          - The boundary condition number
3195: . type        - The type of condition, e.g. DM_BC_ESSENTIAL/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
3196: . name        - The BC name
3197: . labelname   - The label defining constrained points
3198: . field       - The field to constrain
3199: . numcomps    - The number of constrained field components
3200: . comps       - An array of constrained component numbers
3201: . bcFunc      - A pointwise function giving boundary values
3202: . bcFunc_t    - A pointwise function giving the time derviative of the boundary values, or NULL
3203: . numids      - The number of DMLabel ids for constrained points
3204: . ids         - An array of ids for constrained points
3205: - ctx         - An optional user context for bcFunc

3207:   Note:
3208:   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.

3210:   Level: developer

3212: .seealso: PetscDSAddBoundary(), PetscDSGetBoundary(), PetscDSGetNumBoundary()
3213: @*/
3214: 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)
3215: {
3216:   DSBoundary     b = ds->boundary;
3217:   PetscInt       n = 0;

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

3256: /*@
3257:   PetscDSGetNumBoundary - Get the number of registered BC

3259:   Input Parameters:
3260: . ds - The PetscDS object

3262:   Output Parameters:
3263: . numBd - The number of BC

3265:   Level: intermediate

3267: .seealso: PetscDSAddBoundary(), PetscDSGetBoundary()
3268: @*/
3269: PetscErrorCode PetscDSGetNumBoundary(PetscDS ds, PetscInt *numBd)
3270: {
3271:   DSBoundary b = ds->boundary;

3276:   *numBd = 0;
3277:   while (b) {++(*numBd); b = b->next;}
3278:   return(0);
3279: }

3281: /*@C
3282:   PetscDSGetBoundary - Gets a boundary condition to the model

3284:   Input Parameters:
3285: + ds          - The PetscDS object
3286: - bd          - The BC number

3288:   Output Parameters:
3289: + type        - The type of condition, e.g. DM_BC_ESSENTIAL/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
3290: . name        - The BC name
3291: . labelname   - The label defining constrained points
3292: . field       - The field to constrain
3293: . numcomps    - The number of constrained field components
3294: . comps       - An array of constrained component numbers
3295: . bcFunc      - A pointwise function giving boundary values
3296: . bcFunc_t    - A pointwise function giving the time derviative of the boundary values
3297: . numids      - The number of DMLabel ids for constrained points
3298: . ids         - An array of ids for constrained points
3299: - ctx         - An optional user context for bcFunc

3301:   Options Database Keys:
3302: + -bc_<boundary name> <num> - Overrides the boundary ids
3303: - -bc_<boundary name>_comp <num> - Overrides the boundary components

3305:   Level: developer

3307: .seealso: PetscDSAddBoundary()
3308: @*/
3309: 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)
3310: {
3311:   DSBoundary b    = ds->boundary;
3312:   PetscInt   n    = 0;

3316:   while (b) {
3317:     if (n == bd) break;
3318:     b = b->next;
3319:     ++n;
3320:   }
3321:   if (!b) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %d is not in [0, %d)", bd, n);
3322:   if (type) {
3324:     *type = b->type;
3325:   }
3326:   if (name) {
3328:     *name = b->name;
3329:   }
3330:   if (labelname) {
3332:     *labelname = b->labelname;
3333:   }
3334:   if (field) {
3336:     *field = b->field;
3337:   }
3338:   if (numcomps) {
3340:     *numcomps = b->numcomps;
3341:   }
3342:   if (comps) {
3344:     *comps = b->comps;
3345:   }
3346:   if (func) {
3348:     *func = b->func;
3349:   }
3350:   if (func_t) {
3352:     *func_t = b->func_t;
3353:   }
3354:   if (numids) {
3356:     *numids = b->numids;
3357:   }
3358:   if (ids) {
3360:     *ids = b->ids;
3361:   }
3362:   if (ctx) {
3364:     *ctx = b->ctx;
3365:   }
3366:   return(0);
3367: }

3369: /*@
3370:   PetscDSCopyBoundary - Copy all boundary condition objects to the new problem

3372:   Not collective

3374:   Input Parameters:
3375: + ds        - The source PetscDS object
3376: . numFields - The number of selected fields, or PETSC_DEFAULT for all fields
3377: - fields    - The selected fields, or NULL for all fields

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

3382:   Level: intermediate

3384: .seealso: PetscDSCopyEquations(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3385: @*/
3386: PetscErrorCode PetscDSCopyBoundary(PetscDS ds, PetscInt numFields, const PetscInt fields[], PetscDS newds)
3387: {
3388:   DSBoundary     b, next, *lastnext;

3394:   if (ds == newds) return(0);
3395:   next = newds->boundary;
3396:   while (next) {
3397:     DSBoundary b = next;

3399:     next = b->next;
3400:     PetscFree(b->comps);
3401:     PetscFree(b->ids);
3402:     PetscFree(b->name);
3403:     PetscFree(b->labelname);
3404:     PetscFree(b);
3405:   }
3406:   lastnext = &(newds->boundary);
3407:   for (b = ds->boundary; b; b = b->next) {
3408:     DSBoundary bNew;
3409:     PetscInt   fieldNew = -1;

3411:     if (numFields > 0 && fields) {
3412:       PetscInt f;

3414:       for (f = 0; f < numFields; ++f) if (b->field == fields[f]) break;
3415:       if (f == numFields) continue;
3416:       fieldNew = f;
3417:     }
3418:     PetscNew(&bNew);
3419:     bNew->numcomps = b->numcomps;
3420:     PetscMalloc1(bNew->numcomps, &bNew->comps);
3421:     PetscArraycpy(bNew->comps, b->comps, bNew->numcomps);
3422:     bNew->numids = b->numids;
3423:     PetscMalloc1(bNew->numids, &bNew->ids);
3424:     PetscArraycpy(bNew->ids, b->ids, bNew->numids);
3425:     PetscStrallocpy(b->labelname,(char **) &(bNew->labelname));
3426:     PetscStrallocpy(b->name,(char **) &(bNew->name));
3427:     bNew->ctx   = b->ctx;
3428:     bNew->type  = b->type;
3429:     bNew->field = fieldNew < 0 ? b->field : fieldNew;
3430:     bNew->func  = b->func;

3432:     *lastnext = bNew;
3433:     lastnext = &(bNew->next);
3434:   }
3435:   return(0);
3436: }

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

3441:   Not collective

3443:   Input Parameter:
3444: + prob - The PetscDS object
3445: . numFields - Number of new fields
3446: - fields - Old field number for each new field

3448:   Output Parameter:
3449: . newprob - The PetscDS copy

3451:   Level: intermediate

3453: .seealso: PetscDSSelectEquations(), PetscDSCopyBoundary(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3454: @*/
3455: PetscErrorCode PetscDSSelectDiscretizations(PetscDS prob, PetscInt numFields, const PetscInt fields[], PetscDS newprob)
3456: {
3457:   PetscInt       Nf, Nfn, fn;

3464:   PetscDSGetNumFields(prob, &Nf);
3465:   PetscDSGetNumFields(newprob, &Nfn);
3466:   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);
3467:   for (fn = 0; fn < numFields; ++fn) {
3468:     const PetscInt f = fields ? fields[fn] : fn;
3469:     PetscObject    disc;

3471:     if (f >= Nf) continue;
3472:     PetscDSGetDiscretization(prob, f, &disc);
3473:     PetscDSSetDiscretization(newprob, fn, disc);
3474:   }
3475:   return(0);
3476: }

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

3481:   Not collective

3483:   Input Parameter:
3484: + prob - The PetscDS object
3485: . numFields - Number of new fields
3486: - fields - Old field number for each new field

3488:   Output Parameter:
3489: . newprob - The PetscDS copy

3491:   Level: intermediate

3493: .seealso: PetscDSSelectDiscretizations(), PetscDSCopyBoundary(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3494: @*/
3495: PetscErrorCode PetscDSSelectEquations(PetscDS prob, PetscInt numFields, const PetscInt fields[], PetscDS newprob)
3496: {
3497:   PetscInt       Nf, Nfn, fn, gn;

3504:   PetscDSGetNumFields(prob, &Nf);
3505:   PetscDSGetNumFields(newprob, &Nfn);
3506:   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);
3507:   for (fn = 0; fn < numFields; ++fn) {
3508:     const PetscInt   f = fields ? fields[fn] : fn;
3509:     PetscPointFunc   obj;
3510:     PetscPointFunc   f0, f1;
3511:     PetscBdPointFunc f0Bd, f1Bd;
3512:     PetscRiemannFunc r;

3514:     if (f >= Nf) continue;
3515:     PetscDSGetObjective(prob, f, &obj);
3516:     PetscDSGetResidual(prob, f, &f0, &f1);
3517:     PetscDSGetBdResidual(prob, f, &f0Bd, &f1Bd);
3518:     PetscDSGetRiemannSolver(prob, f, &r);
3519:     PetscDSSetObjective(newprob, fn, obj);
3520:     PetscDSSetResidual(newprob, fn, f0, f1);
3521:     PetscDSSetBdResidual(newprob, fn, f0Bd, f1Bd);
3522:     PetscDSSetRiemannSolver(newprob, fn, r);
3523:     for (gn = 0; gn < numFields; ++gn) {
3524:       const PetscInt  g = fields ? fields[gn] : gn;
3525:       PetscPointJac   g0, g1, g2, g3;
3526:       PetscPointJac   g0p, g1p, g2p, g3p;
3527:       PetscBdPointJac g0Bd, g1Bd, g2Bd, g3Bd;

3529:       if (g >= Nf) continue;
3530:       PetscDSGetJacobian(prob, f, g, &g0, &g1, &g2, &g3);
3531:       PetscDSGetJacobianPreconditioner(prob, f, g, &g0p, &g1p, &g2p, &g3p);
3532:       PetscDSGetBdJacobian(prob, f, g, &g0Bd, &g1Bd, &g2Bd, &g3Bd);
3533:       PetscDSSetJacobian(newprob, fn, gn, g0, g1, g2, g3);
3534:       PetscDSSetJacobianPreconditioner(newprob, fn, gn, g0p, g1p, g2p, g3p);
3535:       PetscDSSetBdJacobian(newprob, fn, gn, g0Bd, g1Bd, g2Bd, g3Bd);
3536:     }
3537:   }
3538:   return(0);
3539: }

3541: /*@
3542:   PetscDSCopyEquations - Copy all pointwise function pointers to the new problem

3544:   Not collective

3546:   Input Parameter:
3547: . prob - The PetscDS object

3549:   Output Parameter:
3550: . newprob - The PetscDS copy

3552:   Level: intermediate

3554: .seealso: PetscDSCopyBoundary(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3555: @*/
3556: PetscErrorCode PetscDSCopyEquations(PetscDS prob, PetscDS newprob)
3557: {
3558:   PetscInt       Nf, Ng;

3564:   PetscDSGetNumFields(prob, &Nf);
3565:   PetscDSGetNumFields(newprob, &Ng);
3566:   if (Nf != Ng) SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_SIZ, "Number of fields must match %D != %D", Nf, Ng);
3567:   PetscDSSelectEquations(prob, Nf, NULL, newprob);
3568:   return(0);
3569: }
3570: /*@
3571:   PetscDSCopyConstants - Copy all constants to the new problem

3573:   Not collective

3575:   Input Parameter:
3576: . prob - The PetscDS object

3578:   Output Parameter:
3579: . newprob - The PetscDS copy

3581:   Level: intermediate

3583: .seealso: PetscDSCopyBoundary(), PetscDSCopyEquations(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3584: @*/
3585: PetscErrorCode PetscDSCopyConstants(PetscDS prob, PetscDS newprob)
3586: {
3587:   PetscInt           Nc;
3588:   const PetscScalar *constants;
3589:   PetscErrorCode     ierr;

3594:   PetscDSGetConstants(prob, &Nc, &constants);
3595:   PetscDSSetConstants(newprob, Nc, (PetscScalar *) constants);
3596:   return(0);
3597: }

3599: PetscErrorCode PetscDSGetHeightSubspace(PetscDS prob, PetscInt height, PetscDS *subprob)
3600: {
3601:   PetscInt       dim, Nf, f;

3607:   if (height == 0) {*subprob = prob; return(0);}
3608:   PetscDSGetNumFields(prob, &Nf);
3609:   PetscDSGetSpatialDimension(prob, &dim);
3610:   if (height > dim) SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_OUTOFRANGE, "DS can only handle height in [0, %D], not %D", dim, height);
3611:   if (!prob->subprobs) {PetscCalloc1(dim, &prob->subprobs);}
3612:   if (!prob->subprobs[height-1]) {
3613:     PetscInt cdim;

3615:     PetscDSCreate(PetscObjectComm((PetscObject) prob), &prob->subprobs[height-1]);
3616:     PetscDSGetCoordinateDimension(prob, &cdim);
3617:     PetscDSSetCoordinateDimension(prob->subprobs[height-1], cdim);
3618:     for (f = 0; f < Nf; ++f) {
3619:       PetscFE      subfe;
3620:       PetscObject  obj;
3621:       PetscClassId id;

3623:       PetscDSGetDiscretization(prob, f, &obj);
3624:       PetscObjectGetClassId(obj, &id);
3625:       if (id == PETSCFE_CLASSID) {PetscFEGetHeightSubspace((PetscFE) obj, height, &subfe);}
3626:       else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unsupported discretization type for field %d", f);
3627:       PetscDSSetDiscretization(prob->subprobs[height-1], f, (PetscObject) subfe);
3628:     }
3629:   }
3630:   *subprob = prob->subprobs[height-1];
3631:   return(0);
3632: }

3634: PetscErrorCode PetscDSGetDiscType_Internal(PetscDS ds, PetscInt f, PetscDiscType *disctype)
3635: {
3636:   PetscObject    obj;
3637:   PetscClassId   id;
3638:   PetscInt       Nf;

3644:   *disctype = PETSC_DISC_NONE;
3645:   PetscDSGetNumFields(ds, &Nf);
3646:   if (f >= Nf) SETERRQ2(PetscObjectComm((PetscObject) ds), PETSC_ERR_ARG_SIZ, "Field %D must be in [0, %D)", f, Nf);
3647:   PetscDSGetDiscretization(ds, f, &obj);
3648:   if (obj) {
3649:     PetscObjectGetClassId(obj, &id);
3650:     if (id == PETSCFE_CLASSID) *disctype = PETSC_DISC_FE;
3651:     else                       *disctype = PETSC_DISC_FV;
3652:   }
3653:   return(0);
3654: }

3656: static PetscErrorCode PetscDSDestroy_Basic(PetscDS prob)
3657: {
3658:   PetscErrorCode      ierr;

3661:   PetscFree(prob->data);
3662:   return(0);
3663: }

3665: static PetscErrorCode PetscDSInitialize_Basic(PetscDS prob)
3666: {
3668:   prob->ops->setfromoptions = NULL;
3669:   prob->ops->setup          = NULL;
3670:   prob->ops->view           = NULL;
3671:   prob->ops->destroy        = PetscDSDestroy_Basic;
3672:   return(0);
3673: }

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

3678:   Level: intermediate

3680: .seealso: PetscDSType, PetscDSCreate(), PetscDSSetType()
3681: M*/

3683: PETSC_EXTERN PetscErrorCode PetscDSCreate_Basic(PetscDS prob)
3684: {
3685:   PetscDS_Basic *b;
3686:   PetscErrorCode      ierr;

3690:   PetscNewLog(prob, &b);
3691:   prob->data = b;

3693:   PetscDSInitialize_Basic(prob);
3694:   return(0);
3695: }