Actual source code: dtds.c

  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; No Fortran Support

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

 30:   Example Usage:
 31: .vb
 32:     PetscDSRegister("my_ds", MyPetscDSCreate);
 33: .ve

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

 45:   Level: advanced

 47:   Note:
 48:   `PetscDSRegister()` may be called multiple times to add several user-defined `PetscDSs`

 50: .seealso: `PetscDSType`, `PetscDS`, `PetscDSRegisterAll()`, `PetscDSRegisterDestroy()`
 51: @*/
 52: PetscErrorCode PetscDSRegister(const char sname[], PetscErrorCode (*function)(PetscDS))
 53: {
 54:   PetscFunctionBegin;
 55:   PetscCall(PetscFunctionListAdd(&PetscDSList, sname, function));
 56:   PetscFunctionReturn(PETSC_SUCCESS);
 57: }

 59: /*@
 60:   PetscDSSetType - Builds a particular `PetscDS`

 62:   Collective; No Fortran Support

 64:   Input Parameters:
 65: + prob - The `PetscDS` object
 66: - name - The `PetscDSType`

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

 71:   Level: intermediate

 73: .seealso: `PetscDSType`, `PetscDS`, `PetscDSGetType()`, `PetscDSCreate()`
 74: @*/
 75: PetscErrorCode PetscDSSetType(PetscDS prob, PetscDSType name)
 76: {
 77:   PetscErrorCode (*r)(PetscDS);
 78:   PetscBool match;

 80:   PetscFunctionBegin;
 82:   PetscCall(PetscObjectTypeCompare((PetscObject)prob, name, &match));
 83:   if (match) PetscFunctionReturn(PETSC_SUCCESS);

 85:   PetscCall(PetscDSRegisterAll());
 86:   PetscCall(PetscFunctionListFind(PetscDSList, name, &r));
 87:   PetscCheck(r, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscDS type: %s", name);

 89:   PetscTryTypeMethod(prob, destroy);
 90:   prob->ops->destroy = NULL;

 92:   PetscCall((*r)(prob));
 93:   PetscCall(PetscObjectChangeTypeName((PetscObject)prob, name));
 94:   PetscFunctionReturn(PETSC_SUCCESS);
 95: }

 97: /*@
 98:   PetscDSGetType - Gets the `PetscDSType` name (as a string) from the `PetscDS`

100:   Not Collective; No Fortran Support

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

105:   Output Parameter:
106: . name - The `PetscDSType` name

108:   Level: intermediate

110: .seealso: `PetscDSType`, `PetscDS`, `PetscDSSetType()`, `PetscDSCreate()`
111: @*/
112: PetscErrorCode PetscDSGetType(PetscDS prob, PetscDSType *name)
113: {
114:   PetscFunctionBegin;
116:   PetscAssertPointer(name, 2);
117:   PetscCall(PetscDSRegisterAll());
118:   *name = ((PetscObject)prob)->type_name;
119:   PetscFunctionReturn(PETSC_SUCCESS);
120: }

122: static PetscErrorCode PetscDSView_Ascii(PetscDS ds, PetscViewer viewer)
123: {
124:   PetscViewerFormat  format;
125:   const PetscScalar *constants;
126:   PetscInt           Nf, numConstants, f;

128:   PetscFunctionBegin;
129:   PetscCall(PetscDSGetNumFields(ds, &Nf));
130:   PetscCall(PetscViewerGetFormat(viewer, &format));
131:   PetscCall(PetscViewerASCIIPrintf(viewer, "Discrete System with %" PetscInt_FMT " fields\n", Nf));
132:   PetscCall(PetscViewerASCIIPushTab(viewer));
133:   PetscCall(PetscViewerASCIIPrintf(viewer, "  cell total dim %" PetscInt_FMT " total comp %" PetscInt_FMT "\n", ds->totDim, ds->totComp));
134:   if (ds->isCohesive) PetscCall(PetscViewerASCIIPrintf(viewer, "  cohesive cell\n"));
135:   for (f = 0; f < Nf; ++f) {
136:     DSBoundary      b;
137:     PetscObject     obj;
138:     PetscClassId    id;
139:     PetscQuadrature q;
140:     const char     *name;
141:     PetscInt        Nc, Nq, Nqc;

143:     PetscCall(PetscDSGetDiscretization(ds, f, &obj));
144:     PetscCall(PetscObjectGetClassId(obj, &id));
145:     PetscCall(PetscObjectGetName(obj, &name));
146:     PetscCall(PetscViewerASCIIPrintf(viewer, "Field %s", name ? name : "<unknown>"));
147:     PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
148:     if (id == PETSCFE_CLASSID) {
149:       PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
150:       PetscCall(PetscFEGetQuadrature((PetscFE)obj, &q));
151:       PetscCall(PetscViewerASCIIPrintf(viewer, " FEM"));
152:     } else if (id == PETSCFV_CLASSID) {
153:       PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc));
154:       PetscCall(PetscFVGetQuadrature((PetscFV)obj, &q));
155:       PetscCall(PetscViewerASCIIPrintf(viewer, " FVM"));
156:     } else SETERRQ(PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %" PetscInt_FMT, f);
157:     if (Nc > 1) PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT " components", Nc));
158:     else PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT " component ", Nc));
159:     if (ds->implicit[f]) PetscCall(PetscViewerASCIIPrintf(viewer, " (implicit)"));
160:     else PetscCall(PetscViewerASCIIPrintf(viewer, " (explicit)"));
161:     if (q) {
162:       PetscCall(PetscQuadratureGetData(q, NULL, &Nqc, &Nq, NULL, NULL));
163:       PetscCall(PetscViewerASCIIPrintf(viewer, " (Nq %" PetscInt_FMT " Nqc %" PetscInt_FMT ")", Nq, Nqc));
164:     }
165:     PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT "-jet", ds->jetDegree[f]));
166:     PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
167:     PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
168:     PetscCall(PetscViewerASCIIPushTab(viewer));
169:     if (id == PETSCFE_CLASSID) PetscCall(PetscFEView((PetscFE)obj, viewer));
170:     else if (id == PETSCFV_CLASSID) PetscCall(PetscFVView((PetscFV)obj, viewer));
171:     PetscCall(PetscViewerASCIIPopTab(viewer));

173:     for (b = ds->boundary; b; b = b->next) {
174:       char    *name;
175:       PetscInt c, i;

177:       if (b->field != f) continue;
178:       PetscCall(PetscViewerASCIIPushTab(viewer));
179:       PetscCall(PetscViewerASCIIPrintf(viewer, "Boundary %s (%s) %s\n", b->name, b->lname, DMBoundaryConditionTypes[b->type]));
180:       if (!b->Nc) {
181:         PetscCall(PetscViewerASCIIPrintf(viewer, "  all components\n"));
182:       } else {
183:         PetscCall(PetscViewerASCIIPrintf(viewer, "  components: "));
184:         PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
185:         for (c = 0; c < b->Nc; ++c) {
186:           if (c > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ", "));
187:           PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT, b->comps[c]));
188:         }
189:         PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
190:         PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
191:       }
192:       PetscCall(PetscViewerASCIIPrintf(viewer, "  values: "));
193:       PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
194:       for (i = 0; i < b->Nv; ++i) {
195:         if (i > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ", "));
196:         PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT, b->values[i]));
197:       }
198:       PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
199:       PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
200: #if defined(__clang__)
201:       PETSC_PRAGMA_DIAGNOSTIC_IGNORED_BEGIN("-Wformat-pedantic")
202: #elif defined(__GNUC__) || defined(__GNUG__)
203:       PETSC_PRAGMA_DIAGNOSTIC_IGNORED_BEGIN("-Wformat")
204: #endif
205:       if (b->func) {
206:         PetscCall(PetscDLAddr(b->func, &name));
207:         if (name) PetscCall(PetscViewerASCIIPrintf(viewer, "  func: %s\n", name));
208:         else PetscCall(PetscViewerASCIIPrintf(viewer, "  func: %p\n", b->func));
209:         PetscCall(PetscFree(name));
210:       }
211:       if (b->func_t) {
212:         PetscCall(PetscDLAddr(b->func_t, &name));
213:         if (name) PetscCall(PetscViewerASCIIPrintf(viewer, "  func_t: %s\n", name));
214:         else PetscCall(PetscViewerASCIIPrintf(viewer, "  func_t: %p\n", b->func_t));
215:         PetscCall(PetscFree(name));
216:       }
217:       PETSC_PRAGMA_DIAGNOSTIC_IGNORED_END()
218:       PetscCall(PetscWeakFormView(b->wf, viewer));
219:       PetscCall(PetscViewerASCIIPopTab(viewer));
220:     }
221:   }
222:   PetscCall(PetscDSGetConstants(ds, &numConstants, &constants));
223:   if (numConstants) {
224:     PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " constants\n", numConstants));
225:     PetscCall(PetscViewerASCIIPushTab(viewer));
226:     for (f = 0; f < numConstants; ++f) PetscCall(PetscViewerASCIIPrintf(viewer, "%g\n", (double)PetscRealPart(constants[f])));
227:     PetscCall(PetscViewerASCIIPopTab(viewer));
228:   }
229:   PetscCall(PetscWeakFormView(ds->wf, viewer));
230:   PetscCall(PetscViewerASCIIPopTab(viewer));
231:   PetscFunctionReturn(PETSC_SUCCESS);
232: }

234: /*@
235:   PetscDSViewFromOptions - View a `PetscDS` based on values in the options database

237:   Collective

239:   Input Parameters:
240: + A    - the `PetscDS` object
241: . obj  - Optional object that provides the options prefix used in the search
242: - name - command line option

244:   Level: intermediate

246: .seealso: `PetscDSType`, `PetscDS`, `PetscDSView()`, `PetscObjectViewFromOptions()`, `PetscDSCreate()`
247: @*/
248: PetscErrorCode PetscDSViewFromOptions(PetscDS A, PetscObject obj, const char name[])
249: {
250:   PetscFunctionBegin;
252:   PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
253:   PetscFunctionReturn(PETSC_SUCCESS);
254: }

256: /*@
257:   PetscDSView - Views a `PetscDS`

259:   Collective

261:   Input Parameters:
262: + prob - the `PetscDS` object to view
263: - v    - the viewer

265:   Level: developer

267: .seealso: `PetscDSType`, `PetscDS`, `PetscViewer`, `PetscDSDestroy()`, `PetscDSViewFromOptions()`
268: @*/
269: PetscErrorCode PetscDSView(PetscDS prob, PetscViewer v)
270: {
271:   PetscBool iascii;

273:   PetscFunctionBegin;
275:   if (!v) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)prob), &v));
277:   PetscCall(PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &iascii));
278:   if (iascii) PetscCall(PetscDSView_Ascii(prob, v));
279:   PetscTryTypeMethod(prob, view, v);
280:   PetscFunctionReturn(PETSC_SUCCESS);
281: }

283: /*@
284:   PetscDSSetFromOptions - sets parameters in a `PetscDS` from the options database

286:   Collective

288:   Input Parameter:
289: . prob - the `PetscDS` object to set options for

291:   Options Database Keys:
292: + -petscds_type <type>     - Set the `PetscDS` type
293: . -petscds_view <view opt> - View the `PetscDS`
294: . -petscds_jac_pre         - Turn formation of a separate Jacobian preconditioner on or off
295: . -bc_<name> <ids>         - Specify a list of label ids for a boundary condition
296: - -bc_<name>_comp <comps>  - Specify a list of field components to constrain for a boundary condition

298:   Level: intermediate

300: .seealso: `PetscDS`, `PetscDSView()`
301: @*/
302: PetscErrorCode PetscDSSetFromOptions(PetscDS prob)
303: {
304:   DSBoundary  b;
305:   const char *defaultType;
306:   char        name[256];
307:   PetscBool   flg;

309:   PetscFunctionBegin;
311:   if (!((PetscObject)prob)->type_name) {
312:     defaultType = PETSCDSBASIC;
313:   } else {
314:     defaultType = ((PetscObject)prob)->type_name;
315:   }
316:   PetscCall(PetscDSRegisterAll());

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

324:     PetscCall(PetscSNPrintf(optname, sizeof(optname), "-bc_%s", b->name));
325:     PetscCall(PetscMemzero(ids, sizeof(ids)));
326:     PetscCall(PetscOptionsIntArray(optname, "List of boundary IDs", "", ids, &len, &flg));
327:     if (flg) {
328:       b->Nv = len;
329:       PetscCall(PetscFree(b->values));
330:       PetscCall(PetscMalloc1(len, &b->values));
331:       PetscCall(PetscArraycpy(b->values, ids, len));
332:       PetscCall(PetscWeakFormRewriteKeys(b->wf, b->label, len, b->values));
333:     }
334:     len = 1024;
335:     PetscCall(PetscSNPrintf(optname, sizeof(optname), "-bc_%s_comp", b->name));
336:     PetscCall(PetscMemzero(ids, sizeof(ids)));
337:     PetscCall(PetscOptionsIntArray(optname, "List of boundary field components", "", ids, &len, &flg));
338:     if (flg) {
339:       b->Nc = len;
340:       PetscCall(PetscFree(b->comps));
341:       PetscCall(PetscMalloc1(len, &b->comps));
342:       PetscCall(PetscArraycpy(b->comps, ids, len));
343:     }
344:   }
345:   PetscCall(PetscOptionsFList("-petscds_type", "Discrete System", "PetscDSSetType", PetscDSList, defaultType, name, 256, &flg));
346:   if (flg) {
347:     PetscCall(PetscDSSetType(prob, name));
348:   } else if (!((PetscObject)prob)->type_name) {
349:     PetscCall(PetscDSSetType(prob, defaultType));
350:   }
351:   PetscCall(PetscOptionsBool("-petscds_jac_pre", "Discrete System", "PetscDSUseJacobianPreconditioner", prob->useJacPre, &prob->useJacPre, &flg));
352:   PetscCall(PetscOptionsBool("-petscds_force_quad", "Discrete System", "PetscDSSetForceQuad", prob->forceQuad, &prob->forceQuad, &flg));
353:   PetscCall(PetscOptionsInt("-petscds_print_integrate", "Discrete System", "", prob->printIntegrate, &prob->printIntegrate, NULL));
354:   PetscTryTypeMethod(prob, setfromoptions);
355:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
356:   PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)prob, PetscOptionsObject));
357:   PetscOptionsEnd();
358:   if (prob->Nf) PetscCall(PetscDSViewFromOptions(prob, NULL, "-petscds_view"));
359:   PetscFunctionReturn(PETSC_SUCCESS);
360: }

362: /*@
363:   PetscDSSetUp - Construct data structures for the `PetscDS`

365:   Collective

367:   Input Parameter:
368: . prob - the `PetscDS` object to setup

370:   Level: developer

372: .seealso: `PetscDS`, `PetscDSView()`, `PetscDSDestroy()`
373: @*/
374: PetscErrorCode PetscDSSetUp(PetscDS prob)
375: {
376:   const PetscInt Nf          = prob->Nf;
377:   PetscBool      hasH        = PETSC_FALSE;
378:   PetscInt       maxOrder[4] = {-2, -2, -2, -2};
379:   PetscInt       dim, dimEmbed, NbMax = 0, NcMax = 0, NqMax = 0, NsMax = 1, f;

381:   PetscFunctionBegin;
383:   if (prob->setup) PetscFunctionReturn(PETSC_SUCCESS);
384:   /* Calculate sizes */
385:   PetscCall(PetscDSGetSpatialDimension(prob, &dim));
386:   PetscCall(PetscDSGetCoordinateDimension(prob, &dimEmbed));
387:   prob->totDim = prob->totComp = 0;
388:   PetscCall(PetscMalloc2(Nf, &prob->Nc, Nf, &prob->Nb));
389:   PetscCall(PetscCalloc2(Nf + 1, &prob->off, Nf + 1, &prob->offDer));
390:   PetscCall(PetscCalloc6(Nf + 1, &prob->offCohesive[0], Nf + 1, &prob->offCohesive[1], Nf + 1, &prob->offCohesive[2], Nf + 1, &prob->offDerCohesive[0], Nf + 1, &prob->offDerCohesive[1], Nf + 1, &prob->offDerCohesive[2]));
391:   PetscCall(PetscMalloc2(Nf, &prob->T, Nf, &prob->Tf));
392:   if (prob->forceQuad) {
393:     // Note: This assumes we have one kind of cell at each dimension.
394:     //       We can fix this by having quadrature hold the celltype
395:     PetscQuadrature maxQuad[4] = {NULL, NULL, NULL, NULL};

397:     for (f = 0; f < Nf; ++f) {
398:       PetscObject     obj;
399:       PetscClassId    id;
400:       PetscQuadrature q = NULL, fq = NULL;
401:       PetscInt        dim = -1, order = -1, forder = -1;

403:       PetscCall(PetscDSGetDiscretization(prob, f, &obj));
404:       if (!obj) continue;
405:       PetscCall(PetscObjectGetClassId(obj, &id));
406:       if (id == PETSCFE_CLASSID) {
407:         PetscFE fe = (PetscFE)obj;

409:         PetscCall(PetscFEGetQuadrature(fe, &q));
410:         PetscCall(PetscFEGetFaceQuadrature(fe, &fq));
411:       } else if (id == PETSCFV_CLASSID) {
412:         PetscFV fv = (PetscFV)obj;

414:         PetscCall(PetscFVGetQuadrature(fv, &q));
415:       }
416:       if (q) {
417:         PetscCall(PetscQuadratureGetData(q, &dim, NULL, NULL, NULL, NULL));
418:         PetscCall(PetscQuadratureGetOrder(q, &order));
419:         if (order > maxOrder[dim]) {
420:           maxOrder[dim] = order;
421:           maxQuad[dim]  = q;
422:         }
423:       }
424:       if (fq) {
425:         PetscCall(PetscQuadratureGetData(fq, &dim, NULL, NULL, NULL, NULL));
426:         PetscCall(PetscQuadratureGetOrder(fq, &forder));
427:         if (forder > maxOrder[dim]) {
428:           maxOrder[dim] = forder;
429:           maxQuad[dim]  = fq;
430:         }
431:       }
432:     }
433:     for (f = 0; f < Nf; ++f) {
434:       PetscObject     obj;
435:       PetscClassId    id;
436:       PetscQuadrature q;
437:       PetscInt        dim;

439:       PetscCall(PetscDSGetDiscretization(prob, f, &obj));
440:       if (!obj) continue;
441:       PetscCall(PetscObjectGetClassId(obj, &id));
442:       if (id == PETSCFE_CLASSID) {
443:         PetscFE fe = (PetscFE)obj;

445:         PetscCall(PetscFEGetQuadrature(fe, &q));
446:         PetscCall(PetscQuadratureGetData(q, &dim, NULL, NULL, NULL, NULL));
447:         PetscCall(PetscFESetQuadrature(fe, maxQuad[dim]));
448:         PetscCall(PetscFESetFaceQuadrature(fe, dim ? maxQuad[dim - 1] : NULL));
449:       } else if (id == PETSCFV_CLASSID) {
450:         PetscFV fv = (PetscFV)obj;

452:         PetscCall(PetscFVGetQuadrature(fv, &q));
453:         PetscCall(PetscQuadratureGetData(q, &dim, NULL, NULL, NULL, NULL));
454:         PetscCall(PetscFVSetQuadrature(fv, maxQuad[dim]));
455:       }
456:     }
457:   }
458:   for (f = 0; f < Nf; ++f) {
459:     PetscObject     obj;
460:     PetscClassId    id;
461:     PetscQuadrature q  = NULL;
462:     PetscInt        Nq = 0, Nb, Nc;

464:     PetscCall(PetscDSGetDiscretization(prob, f, &obj));
465:     if (prob->jetDegree[f] > 1) hasH = PETSC_TRUE;
466:     if (!obj) {
467:       /* Empty mesh */
468:       Nb = Nc    = 0;
469:       prob->T[f] = prob->Tf[f] = NULL;
470:     } else {
471:       PetscCall(PetscObjectGetClassId(obj, &id));
472:       if (id == PETSCFE_CLASSID) {
473:         PetscFE fe = (PetscFE)obj;

475:         PetscCall(PetscFEGetQuadrature(fe, &q));
476:         {
477:           PetscQuadrature fq;
478:           PetscInt        dim, order;

480:           PetscCall(PetscQuadratureGetData(q, &dim, NULL, NULL, NULL, NULL));
481:           PetscCall(PetscQuadratureGetOrder(q, &order));
482:           if (maxOrder[dim] < 0) maxOrder[dim] = order;
483:           PetscCheck(order == maxOrder[dim], PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Field %" PetscInt_FMT " cell quadrature order %" PetscInt_FMT " != %" PetscInt_FMT " DS cell quadrature order", f, order, maxOrder[dim]);
484:           PetscCall(PetscFEGetFaceQuadrature(fe, &fq));
485:           if (fq) {
486:             PetscCall(PetscQuadratureGetData(fq, &dim, NULL, NULL, NULL, NULL));
487:             PetscCall(PetscQuadratureGetOrder(fq, &order));
488:             if (maxOrder[dim] < 0) maxOrder[dim] = order;
489:             PetscCheck(order == maxOrder[dim], PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Field %" PetscInt_FMT " face quadrature order %" PetscInt_FMT " != %" PetscInt_FMT " DS face quadrature order", f, order, maxOrder[dim]);
490:           }
491:         }
492:         PetscCall(PetscFEGetDimension(fe, &Nb));
493:         PetscCall(PetscFEGetNumComponents(fe, &Nc));
494:         PetscCall(PetscFEGetCellTabulation(fe, prob->jetDegree[f], &prob->T[f]));
495:         PetscCall(PetscFEGetFaceTabulation(fe, prob->jetDegree[f], &prob->Tf[f]));
496:       } else if (id == PETSCFV_CLASSID) {
497:         PetscFV fv = (PetscFV)obj;

499:         PetscCall(PetscFVGetQuadrature(fv, &q));
500:         PetscCall(PetscFVGetNumComponents(fv, &Nc));
501:         Nb = Nc;
502:         PetscCall(PetscFVGetCellTabulation(fv, &prob->T[f]));
503:         /* TODO: should PetscFV also have face tabulation? Otherwise there will be a null pointer in prob->basisFace */
504:       } else SETERRQ(PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %" PetscInt_FMT, f);
505:     }
506:     prob->Nc[f]                    = Nc;
507:     prob->Nb[f]                    = Nb;
508:     prob->off[f + 1]               = Nc + prob->off[f];
509:     prob->offDer[f + 1]            = Nc * dim + prob->offDer[f];
510:     prob->offCohesive[0][f + 1]    = (prob->cohesive[f] ? Nc : Nc * 2) + prob->offCohesive[0][f];
511:     prob->offDerCohesive[0][f + 1] = (prob->cohesive[f] ? Nc : Nc * 2) * dimEmbed + prob->offDerCohesive[0][f];
512:     prob->offCohesive[1][f]        = (prob->cohesive[f] ? 0 : Nc) + prob->offCohesive[0][f];
513:     prob->offDerCohesive[1][f]     = (prob->cohesive[f] ? 0 : Nc) * dimEmbed + prob->offDerCohesive[0][f];
514:     prob->offCohesive[2][f + 1]    = (prob->cohesive[f] ? Nc : Nc * 2) + prob->offCohesive[2][f];
515:     prob->offDerCohesive[2][f + 1] = (prob->cohesive[f] ? Nc : Nc * 2) * dimEmbed + prob->offDerCohesive[2][f];
516:     if (q) PetscCall(PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL));
517:     NqMax = PetscMax(NqMax, Nq);
518:     NbMax = PetscMax(NbMax, Nb);
519:     NcMax = PetscMax(NcMax, Nc);
520:     prob->totDim += Nb;
521:     prob->totComp += Nc;
522:     /* There are two faces for all fields on a cohesive cell, except for cohesive fields */
523:     if (prob->isCohesive && !prob->cohesive[f]) prob->totDim += Nb;
524:   }
525:   prob->offCohesive[1][Nf]    = prob->offCohesive[0][Nf];
526:   prob->offDerCohesive[1][Nf] = prob->offDerCohesive[0][Nf];
527:   /* Allocate works space */
528:   NsMax = 2; /* A non-cohesive discretizations can be used on a cohesive cell, so we need this extra workspace for all DS */
529:   PetscCall(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));
530:   PetscCall(PetscMalloc5(dimEmbed, &prob->x, NbMax * NcMax, &prob->basisReal, NbMax * NcMax * dimEmbed, &prob->basisDerReal, NbMax * NcMax, &prob->testReal, NbMax * NcMax * dimEmbed, &prob->testDerReal));
531:   PetscCall(PetscMalloc6(NsMax * NqMax * NcMax, &prob->f0, NsMax * NqMax * NcMax * dimEmbed, &prob->f1, NsMax * NsMax * NqMax * NcMax * NcMax, &prob->g0, NsMax * NsMax * NqMax * NcMax * NcMax * dimEmbed, &prob->g1, NsMax * NsMax * NqMax * NcMax * NcMax * dimEmbed,
532:                          &prob->g2, NsMax * NsMax * NqMax * NcMax * NcMax * dimEmbed * dimEmbed, &prob->g3));
533:   PetscTryTypeMethod(prob, setup);
534:   prob->setup = PETSC_TRUE;
535:   PetscFunctionReturn(PETSC_SUCCESS);
536: }

538: static PetscErrorCode PetscDSDestroyStructs_Static(PetscDS prob)
539: {
540:   PetscFunctionBegin;
541:   PetscCall(PetscFree2(prob->Nc, prob->Nb));
542:   PetscCall(PetscFree2(prob->off, prob->offDer));
543:   PetscCall(PetscFree6(prob->offCohesive[0], prob->offCohesive[1], prob->offCohesive[2], prob->offDerCohesive[0], prob->offDerCohesive[1], prob->offDerCohesive[2]));
544:   PetscCall(PetscFree2(prob->T, prob->Tf));
545:   PetscCall(PetscFree3(prob->u, prob->u_t, prob->u_x));
546:   PetscCall(PetscFree5(prob->x, prob->basisReal, prob->basisDerReal, prob->testReal, prob->testDerReal));
547:   PetscCall(PetscFree6(prob->f0, prob->f1, prob->g0, prob->g1, prob->g2, prob->g3));
548:   PetscFunctionReturn(PETSC_SUCCESS);
549: }

551: static PetscErrorCode PetscDSEnlarge_Static(PetscDS prob, PetscInt NfNew)
552: {
553:   PetscObject         *tmpd;
554:   PetscBool           *tmpi;
555:   PetscInt            *tmpk;
556:   PetscBool           *tmpc;
557:   PetscPointFunc      *tmpup;
558:   PetscSimplePointFn **tmpexactSol, **tmpexactSol_t, **tmplowerBound, **tmpupperBound;
559:   void               **tmpexactCtx, **tmpexactCtx_t, **tmplowerCtx, **tmpupperCtx;
560:   void               **tmpctx;
561:   PetscInt             Nf = prob->Nf, f;

563:   PetscFunctionBegin;
564:   if (Nf >= NfNew) PetscFunctionReturn(PETSC_SUCCESS);
565:   prob->setup = PETSC_FALSE;
566:   PetscCall(PetscDSDestroyStructs_Static(prob));
567:   PetscCall(PetscMalloc4(NfNew, &tmpd, NfNew, &tmpi, NfNew, &tmpc, NfNew, &tmpk));
568:   for (f = 0; f < Nf; ++f) {
569:     tmpd[f] = prob->disc[f];
570:     tmpi[f] = prob->implicit[f];
571:     tmpc[f] = prob->cohesive[f];
572:     tmpk[f] = prob->jetDegree[f];
573:   }
574:   for (f = Nf; f < NfNew; ++f) {
575:     tmpd[f] = NULL;
576:     tmpi[f] = PETSC_TRUE, tmpc[f] = PETSC_FALSE;
577:     tmpk[f] = 1;
578:   }
579:   PetscCall(PetscFree4(prob->disc, prob->implicit, prob->cohesive, prob->jetDegree));
580:   PetscCall(PetscWeakFormSetNumFields(prob->wf, NfNew));
581:   prob->Nf        = NfNew;
582:   prob->disc      = tmpd;
583:   prob->implicit  = tmpi;
584:   prob->cohesive  = tmpc;
585:   prob->jetDegree = tmpk;
586:   PetscCall(PetscCalloc2(NfNew, &tmpup, NfNew, &tmpctx));
587:   for (f = 0; f < Nf; ++f) tmpup[f] = prob->update[f];
588:   for (f = 0; f < Nf; ++f) tmpctx[f] = prob->ctx[f];
589:   for (f = Nf; f < NfNew; ++f) tmpup[f] = NULL;
590:   for (f = Nf; f < NfNew; ++f) tmpctx[f] = NULL;
591:   PetscCall(PetscFree2(prob->update, prob->ctx));
592:   prob->update = tmpup;
593:   prob->ctx    = tmpctx;
594:   PetscCall(PetscCalloc4(NfNew, &tmpexactSol, NfNew, &tmpexactCtx, NfNew, &tmpexactSol_t, NfNew, &tmpexactCtx_t));
595:   PetscCall(PetscCalloc4(NfNew, &tmplowerBound, NfNew, &tmplowerCtx, NfNew, &tmpupperBound, NfNew, &tmpupperCtx));
596:   for (f = 0; f < Nf; ++f) tmpexactSol[f] = prob->exactSol[f];
597:   for (f = 0; f < Nf; ++f) tmpexactCtx[f] = prob->exactCtx[f];
598:   for (f = 0; f < Nf; ++f) tmpexactSol_t[f] = prob->exactSol_t[f];
599:   for (f = 0; f < Nf; ++f) tmpexactCtx_t[f] = prob->exactCtx_t[f];
600:   for (f = 0; f < Nf; ++f) tmplowerBound[f] = prob->lowerBound[f];
601:   for (f = 0; f < Nf; ++f) tmplowerCtx[f] = prob->lowerCtx[f];
602:   for (f = 0; f < Nf; ++f) tmpupperBound[f] = prob->upperBound[f];
603:   for (f = 0; f < Nf; ++f) tmpupperCtx[f] = prob->upperCtx[f];
604:   for (f = Nf; f < NfNew; ++f) tmpexactSol[f] = NULL;
605:   for (f = Nf; f < NfNew; ++f) tmpexactCtx[f] = NULL;
606:   for (f = Nf; f < NfNew; ++f) tmpexactSol_t[f] = NULL;
607:   for (f = Nf; f < NfNew; ++f) tmpexactCtx_t[f] = NULL;
608:   for (f = Nf; f < NfNew; ++f) tmplowerBound[f] = NULL;
609:   for (f = Nf; f < NfNew; ++f) tmplowerCtx[f] = NULL;
610:   for (f = Nf; f < NfNew; ++f) tmpupperBound[f] = NULL;
611:   for (f = Nf; f < NfNew; ++f) tmpupperCtx[f] = NULL;
612:   PetscCall(PetscFree4(prob->exactSol, prob->exactCtx, prob->exactSol_t, prob->exactCtx_t));
613:   PetscCall(PetscFree4(prob->lowerBound, prob->lowerCtx, prob->upperBound, prob->upperCtx));
614:   prob->exactSol   = tmpexactSol;
615:   prob->exactCtx   = tmpexactCtx;
616:   prob->exactSol_t = tmpexactSol_t;
617:   prob->exactCtx_t = tmpexactCtx_t;
618:   prob->lowerBound = tmplowerBound;
619:   prob->lowerCtx   = tmplowerCtx;
620:   prob->upperBound = tmpupperBound;
621:   prob->upperCtx   = tmpupperCtx;
622:   PetscFunctionReturn(PETSC_SUCCESS);
623: }

625: /*@
626:   PetscDSDestroy - Destroys a `PetscDS` object

628:   Collective

630:   Input Parameter:
631: . ds - the `PetscDS` object to destroy

633:   Level: developer

635: .seealso: `PetscDSView()`
636: @*/
637: PetscErrorCode PetscDSDestroy(PetscDS *ds)
638: {
639:   PetscInt f;

641:   PetscFunctionBegin;
642:   if (!*ds) PetscFunctionReturn(PETSC_SUCCESS);

645:   if (--((PetscObject)*ds)->refct > 0) {
646:     *ds = NULL;
647:     PetscFunctionReturn(PETSC_SUCCESS);
648:   }
649:   ((PetscObject)*ds)->refct = 0;
650:   if ((*ds)->subprobs) {
651:     PetscInt dim, d;

653:     PetscCall(PetscDSGetSpatialDimension(*ds, &dim));
654:     for (d = 0; d < dim; ++d) PetscCall(PetscDSDestroy(&(*ds)->subprobs[d]));
655:   }
656:   PetscCall(PetscFree((*ds)->subprobs));
657:   PetscCall(PetscDSDestroyStructs_Static(*ds));
658:   for (f = 0; f < (*ds)->Nf; ++f) PetscCall(PetscObjectDereference((*ds)->disc[f]));
659:   PetscCall(PetscFree4((*ds)->disc, (*ds)->implicit, (*ds)->cohesive, (*ds)->jetDegree));
660:   PetscCall(PetscWeakFormDestroy(&(*ds)->wf));
661:   PetscCall(PetscFree2((*ds)->update, (*ds)->ctx));
662:   PetscCall(PetscFree4((*ds)->exactSol, (*ds)->exactCtx, (*ds)->exactSol_t, (*ds)->exactCtx_t));
663:   PetscCall(PetscFree4((*ds)->lowerBound, (*ds)->lowerCtx, (*ds)->upperBound, (*ds)->upperCtx));
664:   PetscTryTypeMethod(*ds, destroy);
665:   PetscCall(PetscDSDestroyBoundary(*ds));
666:   PetscCall(PetscFree((*ds)->constants));
667:   for (PetscInt c = 0; c < DM_NUM_POLYTOPES; ++c) {
668:     const PetscInt Na = DMPolytopeTypeGetNumArrangements((DMPolytopeType)c);
669:     if ((*ds)->quadPerm[c])
670:       for (PetscInt o = 0; o < Na; ++o) PetscCall(ISDestroy(&(*ds)->quadPerm[c][o]));
671:     PetscCall(PetscFree((*ds)->quadPerm[c]));
672:     (*ds)->quadPerm[c] = NULL;
673:   }
674:   PetscCall(PetscHeaderDestroy(ds));
675:   PetscFunctionReturn(PETSC_SUCCESS);
676: }

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

681:   Collective

683:   Input Parameter:
684: . comm - The communicator for the `PetscDS` object

686:   Output Parameter:
687: . ds - The `PetscDS` object

689:   Level: beginner

691: .seealso: `PetscDS`, `PetscDSSetType()`, `PETSCDSBASIC`, `PetscDSType`
692: @*/
693: PetscErrorCode PetscDSCreate(MPI_Comm comm, PetscDS *ds)
694: {
695:   PetscDS p;

697:   PetscFunctionBegin;
698:   PetscAssertPointer(ds, 2);
699:   PetscCall(PetscDSInitializePackage());

701:   PetscCall(PetscHeaderCreate(p, PETSCDS_CLASSID, "PetscDS", "Discrete System", "PetscDS", comm, PetscDSDestroy, PetscDSView));
702:   p->Nf               = 0;
703:   p->setup            = PETSC_FALSE;
704:   p->numConstants     = 0;
705:   p->numFuncConstants = 3; // Row and col fields, cell size
706:   p->dimEmbed         = -1;
707:   p->useJacPre        = PETSC_TRUE;
708:   p->forceQuad        = PETSC_TRUE;
709:   PetscCall(PetscMalloc1(p->numConstants + p->numFuncConstants, &p->constants));
710:   PetscCall(PetscWeakFormCreate(comm, &p->wf));
711:   PetscCall(PetscArrayzero(p->quadPerm, DM_NUM_POLYTOPES));
712:   *ds = p;
713:   PetscFunctionReturn(PETSC_SUCCESS);
714: }

716: /*@
717:   PetscDSGetNumFields - Returns the number of fields in the `PetscDS`

719:   Not Collective

721:   Input Parameter:
722: . prob - The `PetscDS` object

724:   Output Parameter:
725: . Nf - The number of fields

727:   Level: beginner

729: .seealso: `PetscDS`, `PetscDSGetSpatialDimension()`, `PetscDSCreate()`
730: @*/
731: PetscErrorCode PetscDSGetNumFields(PetscDS prob, PetscInt *Nf)
732: {
733:   PetscFunctionBegin;
735:   PetscAssertPointer(Nf, 2);
736:   *Nf = prob->Nf;
737:   PetscFunctionReturn(PETSC_SUCCESS);
738: }

740: /*@
741:   PetscDSGetSpatialDimension - Returns the spatial dimension of the `PetscDS`, meaning the topological dimension of the discretizations

743:   Not Collective

745:   Input Parameter:
746: . prob - The `PetscDS` object

748:   Output Parameter:
749: . dim - The spatial dimension

751:   Level: beginner

753: .seealso: `PetscDS`, `PetscDSGetCoordinateDimension()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
754: @*/
755: PetscErrorCode PetscDSGetSpatialDimension(PetscDS prob, PetscInt *dim)
756: {
757:   PetscFunctionBegin;
759:   PetscAssertPointer(dim, 2);
760:   *dim = 0;
761:   if (prob->Nf) {
762:     PetscObject  obj;
763:     PetscClassId id;

765:     PetscCall(PetscDSGetDiscretization(prob, 0, &obj));
766:     if (obj) {
767:       PetscCall(PetscObjectGetClassId(obj, &id));
768:       if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetSpatialDimension((PetscFE)obj, dim));
769:       else if (id == PETSCFV_CLASSID) PetscCall(PetscFVGetSpatialDimension((PetscFV)obj, dim));
770:       else SETERRQ(PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", 0);
771:     }
772:   }
773:   PetscFunctionReturn(PETSC_SUCCESS);
774: }

776: /*@
777:   PetscDSGetCoordinateDimension - Returns the coordinate dimension of the `PetscDS`, meaning the dimension of the space into which the discretiaztions are embedded

779:   Not Collective

781:   Input Parameter:
782: . prob - The `PetscDS` object

784:   Output Parameter:
785: . dimEmbed - The coordinate dimension

787:   Level: beginner

789: .seealso: `PetscDS`, `PetscDSSetCoordinateDimension()`, `PetscDSGetSpatialDimension()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
790: @*/
791: PetscErrorCode PetscDSGetCoordinateDimension(PetscDS prob, PetscInt *dimEmbed)
792: {
793:   PetscFunctionBegin;
795:   PetscAssertPointer(dimEmbed, 2);
796:   PetscCheck(prob->dimEmbed >= 0, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONGSTATE, "No coordinate dimension set for this DS");
797:   *dimEmbed = prob->dimEmbed;
798:   PetscFunctionReturn(PETSC_SUCCESS);
799: }

801: /*@
802:   PetscDSSetCoordinateDimension - Set the coordinate dimension of the `PetscDS`, meaning the dimension of the space into which the discretiaztions are embedded

804:   Logically Collective

806:   Input Parameters:
807: + prob     - The `PetscDS` object
808: - dimEmbed - The coordinate dimension

810:   Level: beginner

812: .seealso: `PetscDS`, `PetscDSGetCoordinateDimension()`, `PetscDSGetSpatialDimension()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
813: @*/
814: PetscErrorCode PetscDSSetCoordinateDimension(PetscDS prob, PetscInt dimEmbed)
815: {
816:   PetscFunctionBegin;
818:   PetscCheck(dimEmbed >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Coordinate dimension must be non-negative, not %" PetscInt_FMT, dimEmbed);
819:   prob->dimEmbed = dimEmbed;
820:   PetscFunctionReturn(PETSC_SUCCESS);
821: }

823: /*@
824:   PetscDSGetForceQuad - Returns the flag to force matching quadratures among the field discretizations

826:   Not collective

828:   Input Parameter:
829: . ds - The `PetscDS` object

831:   Output Parameter:
832: . forceQuad - The flag

834:   Level: intermediate

836: .seealso: `PetscDS`, `PetscDSSetForceQuad()`, `PetscDSGetDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
837: @*/
838: PetscErrorCode PetscDSGetForceQuad(PetscDS ds, PetscBool *forceQuad)
839: {
840:   PetscFunctionBegin;
842:   PetscAssertPointer(forceQuad, 2);
843:   *forceQuad = ds->forceQuad;
844:   PetscFunctionReturn(PETSC_SUCCESS);
845: }

847: /*@
848:   PetscDSSetForceQuad - Set the flag to force matching quadratures among the field discretizations

850:   Logically collective on ds

852:   Input Parameters:
853: + ds        - The `PetscDS` object
854: - forceQuad - The flag

856:   Level: intermediate

858: .seealso: `PetscDS`, `PetscDSGetForceQuad()`, `PetscDSGetDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
859: @*/
860: PetscErrorCode PetscDSSetForceQuad(PetscDS ds, PetscBool forceQuad)
861: {
862:   PetscFunctionBegin;
864:   ds->forceQuad = forceQuad;
865:   PetscFunctionReturn(PETSC_SUCCESS);
866: }

868: /*@
869:   PetscDSIsCohesive - Returns the flag indicating that this `PetscDS` is for a cohesive cell

871:   Not Collective

873:   Input Parameter:
874: . ds - The `PetscDS` object

876:   Output Parameter:
877: . isCohesive - The flag

879:   Level: developer

881: .seealso: `PetscDS`, `PetscDSGetNumCohesive()`, `PetscDSGetCohesive()`, `PetscDSSetCohesive()`, `PetscDSCreate()`
882: @*/
883: PetscErrorCode PetscDSIsCohesive(PetscDS ds, PetscBool *isCohesive)
884: {
885:   PetscFunctionBegin;
887:   PetscAssertPointer(isCohesive, 2);
888:   *isCohesive = ds->isCohesive;
889:   PetscFunctionReturn(PETSC_SUCCESS);
890: }

892: /*@
893:   PetscDSGetNumCohesive - Returns the number of cohesive fields, meaning those defined on the interior of a cohesive cell

895:   Not Collective

897:   Input Parameter:
898: . ds - The `PetscDS` object

900:   Output Parameter:
901: . numCohesive - The number of cohesive fields

903:   Level: developer

905: .seealso: `PetscDS`, `PetscDSSetCohesive()`, `PetscDSCreate()`
906: @*/
907: PetscErrorCode PetscDSGetNumCohesive(PetscDS ds, PetscInt *numCohesive)
908: {
909:   PetscInt f;

911:   PetscFunctionBegin;
913:   PetscAssertPointer(numCohesive, 2);
914:   *numCohesive = 0;
915:   for (f = 0; f < ds->Nf; ++f) *numCohesive += ds->cohesive[f] ? 1 : 0;
916:   PetscFunctionReturn(PETSC_SUCCESS);
917: }

919: /*@
920:   PetscDSGetCohesive - Returns the flag indicating that a field is cohesive, meaning it is defined on the interior of a cohesive cell

922:   Not Collective

924:   Input Parameters:
925: + ds - The `PetscDS` object
926: - f  - The field index

928:   Output Parameter:
929: . isCohesive - The flag

931:   Level: developer

933: .seealso: `PetscDS`, `PetscDSSetCohesive()`, `PetscDSIsCohesive()`, `PetscDSCreate()`
934: @*/
935: PetscErrorCode PetscDSGetCohesive(PetscDS ds, PetscInt f, PetscBool *isCohesive)
936: {
937:   PetscFunctionBegin;
939:   PetscAssertPointer(isCohesive, 3);
940:   PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
941:   *isCohesive = ds->cohesive[f];
942:   PetscFunctionReturn(PETSC_SUCCESS);
943: }

945: /*@
946:   PetscDSSetCohesive - Set the flag indicating that a field is cohesive, meaning it is defined on the interior of a cohesive cell

948:   Not Collective

950:   Input Parameters:
951: + ds         - The `PetscDS` object
952: . f          - The field index
953: - isCohesive - The flag for a cohesive field

955:   Level: developer

957: .seealso: `PetscDS`, `PetscDSGetCohesive()`, `PetscDSIsCohesive()`, `PetscDSCreate()`
958: @*/
959: PetscErrorCode PetscDSSetCohesive(PetscDS ds, PetscInt f, PetscBool isCohesive)
960: {
961:   PetscInt i;

963:   PetscFunctionBegin;
965:   PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
966:   ds->cohesive[f] = isCohesive;
967:   ds->isCohesive  = PETSC_FALSE;
968:   for (i = 0; i < ds->Nf; ++i) ds->isCohesive = ds->isCohesive || ds->cohesive[f] ? PETSC_TRUE : PETSC_FALSE;
969:   PetscFunctionReturn(PETSC_SUCCESS);
970: }

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

975:   Not Collective

977:   Input Parameter:
978: . prob - The `PetscDS` object

980:   Output Parameter:
981: . dim - The total problem dimension

983:   Level: beginner

985: .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
986: @*/
987: PetscErrorCode PetscDSGetTotalDimension(PetscDS prob, PetscInt *dim)
988: {
989:   PetscFunctionBegin;
991:   PetscCall(PetscDSSetUp(prob));
992:   PetscAssertPointer(dim, 2);
993:   *dim = prob->totDim;
994:   PetscFunctionReturn(PETSC_SUCCESS);
995: }

997: /*@
998:   PetscDSGetTotalComponents - Returns the total number of components in this system

1000:   Not Collective

1002:   Input Parameter:
1003: . prob - The `PetscDS` object

1005:   Output Parameter:
1006: . Nc - The total number of components

1008:   Level: beginner

1010: .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1011: @*/
1012: PetscErrorCode PetscDSGetTotalComponents(PetscDS prob, PetscInt *Nc)
1013: {
1014:   PetscFunctionBegin;
1016:   PetscCall(PetscDSSetUp(prob));
1017:   PetscAssertPointer(Nc, 2);
1018:   *Nc = prob->totComp;
1019:   PetscFunctionReturn(PETSC_SUCCESS);
1020: }

1022: /*@
1023:   PetscDSGetDiscretization - Returns the discretization object for the given field

1025:   Not Collective

1027:   Input Parameters:
1028: + prob - The `PetscDS` object
1029: - f    - The field number

1031:   Output Parameter:
1032: . disc - The discretization object

1034:   Level: beginner

1036: .seealso: `PetscDS`, `PetscFE`, `PetscFV`, `PetscDSSetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1037: @*/
1038: PetscErrorCode PetscDSGetDiscretization(PetscDS prob, PetscInt f, PetscObject *disc)
1039: {
1040:   PetscFunctionBeginHot;
1042:   PetscAssertPointer(disc, 3);
1043:   PetscCheck(!(f < 0) && !(f >= prob->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, prob->Nf);
1044:   *disc = prob->disc[f];
1045:   PetscFunctionReturn(PETSC_SUCCESS);
1046: }

1048: /*@
1049:   PetscDSSetDiscretization - Sets the discretization object for the given field

1051:   Not Collective

1053:   Input Parameters:
1054: + prob - The `PetscDS` object
1055: . f    - The field number
1056: - disc - The discretization object

1058:   Level: beginner

1060: .seealso: `PetscDS`, `PetscFE`, `PetscFV`, `PetscDSGetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1061: @*/
1062: PetscErrorCode PetscDSSetDiscretization(PetscDS prob, PetscInt f, PetscObject disc)
1063: {
1064:   PetscFunctionBegin;
1066:   if (disc) PetscAssertPointer(disc, 3);
1067:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
1068:   PetscCall(PetscDSEnlarge_Static(prob, f + 1));
1069:   PetscCall(PetscObjectDereference(prob->disc[f]));
1070:   prob->disc[f] = disc;
1071:   PetscCall(PetscObjectReference(disc));
1072:   if (disc) {
1073:     PetscClassId id;

1075:     PetscCall(PetscObjectGetClassId(disc, &id));
1076:     if (id == PETSCFE_CLASSID) {
1077:       PetscCall(PetscDSSetImplicit(prob, f, PETSC_TRUE));
1078:     } else if (id == PETSCFV_CLASSID) {
1079:       PetscCall(PetscDSSetImplicit(prob, f, PETSC_FALSE));
1080:     }
1081:     PetscCall(PetscDSSetJetDegree(prob, f, 1));
1082:   }
1083:   PetscFunctionReturn(PETSC_SUCCESS);
1084: }

1086: /*@
1087:   PetscDSGetWeakForm - Returns the weak form object

1089:   Not Collective

1091:   Input Parameter:
1092: . ds - The `PetscDS` object

1094:   Output Parameter:
1095: . wf - The weak form object

1097:   Level: beginner

1099: .seealso: `PetscWeakForm`, `PetscDSSetWeakForm()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1100: @*/
1101: PetscErrorCode PetscDSGetWeakForm(PetscDS ds, PetscWeakForm *wf)
1102: {
1103:   PetscFunctionBegin;
1105:   PetscAssertPointer(wf, 2);
1106:   *wf = ds->wf;
1107:   PetscFunctionReturn(PETSC_SUCCESS);
1108: }

1110: /*@
1111:   PetscDSSetWeakForm - Sets the weak form object

1113:   Not Collective

1115:   Input Parameters:
1116: + ds - The `PetscDS` object
1117: - wf - The weak form object

1119:   Level: beginner

1121: .seealso: `PetscWeakForm`, `PetscDSGetWeakForm()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1122: @*/
1123: PetscErrorCode PetscDSSetWeakForm(PetscDS ds, PetscWeakForm wf)
1124: {
1125:   PetscFunctionBegin;
1128:   PetscCall(PetscObjectDereference((PetscObject)ds->wf));
1129:   ds->wf = wf;
1130:   PetscCall(PetscObjectReference((PetscObject)wf));
1131:   PetscCall(PetscWeakFormSetNumFields(wf, ds->Nf));
1132:   PetscFunctionReturn(PETSC_SUCCESS);
1133: }

1135: /*@
1136:   PetscDSAddDiscretization - Adds a discretization object

1138:   Not Collective

1140:   Input Parameters:
1141: + prob - The `PetscDS` object
1142: - disc - The boundary discretization object

1144:   Level: beginner

1146: .seealso: `PetscWeakForm`, `PetscDSGetDiscretization()`, `PetscDSSetDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1147: @*/
1148: PetscErrorCode PetscDSAddDiscretization(PetscDS prob, PetscObject disc)
1149: {
1150:   PetscFunctionBegin;
1151:   PetscCall(PetscDSSetDiscretization(prob, prob->Nf, disc));
1152:   PetscFunctionReturn(PETSC_SUCCESS);
1153: }

1155: /*@
1156:   PetscDSGetQuadrature - Returns the quadrature, which must agree for all fields in the `PetscDS`

1158:   Not Collective

1160:   Input Parameter:
1161: . prob - The `PetscDS` object

1163:   Output Parameter:
1164: . q - The quadrature object

1166:   Level: intermediate

1168: .seealso: `PetscDS`, `PetscQuadrature`, `PetscDSSetImplicit()`, `PetscDSSetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1169: @*/
1170: PetscErrorCode PetscDSGetQuadrature(PetscDS prob, PetscQuadrature *q)
1171: {
1172:   PetscObject  obj;
1173:   PetscClassId id;

1175:   PetscFunctionBegin;
1176:   *q = NULL;
1177:   if (!prob->Nf) PetscFunctionReturn(PETSC_SUCCESS);
1178:   PetscCall(PetscDSGetDiscretization(prob, 0, &obj));
1179:   PetscCall(PetscObjectGetClassId(obj, &id));
1180:   if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetQuadrature((PetscFE)obj, q));
1181:   else if (id == PETSCFV_CLASSID) PetscCall(PetscFVGetQuadrature((PetscFV)obj, q));
1182:   else SETERRQ(PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", 0);
1183:   PetscFunctionReturn(PETSC_SUCCESS);
1184: }

1186: /*@
1187:   PetscDSGetImplicit - Returns the flag for implicit solve for this field. This is just a guide for `TSARKIMEX`

1189:   Not Collective

1191:   Input Parameters:
1192: + prob - The `PetscDS` object
1193: - f    - The field number

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

1198:   Level: developer

1200: .seealso: `TSARKIMEX`, `PetscDS`, `PetscDSSetImplicit()`, `PetscDSSetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1201: @*/
1202: PetscErrorCode PetscDSGetImplicit(PetscDS prob, PetscInt f, PetscBool *implicit)
1203: {
1204:   PetscFunctionBegin;
1206:   PetscAssertPointer(implicit, 3);
1207:   PetscCheck(!(f < 0) && !(f >= prob->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, prob->Nf);
1208:   *implicit = prob->implicit[f];
1209:   PetscFunctionReturn(PETSC_SUCCESS);
1210: }

1212: /*@
1213:   PetscDSSetImplicit - Set the flag for implicit solve for this field. This is just a guide for `TSARKIMEX`

1215:   Not Collective

1217:   Input Parameters:
1218: + prob     - The `PetscDS` object
1219: . f        - The field number
1220: - implicit - The flag indicating what kind of solve to use for this field

1222:   Level: developer

1224: .seealso: `TSARKIMEX`, `PetscDSGetImplicit()`, `PetscDSSetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1225: @*/
1226: PetscErrorCode PetscDSSetImplicit(PetscDS prob, PetscInt f, PetscBool implicit)
1227: {
1228:   PetscFunctionBegin;
1230:   PetscCheck(!(f < 0) && !(f >= prob->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, prob->Nf);
1231:   prob->implicit[f] = implicit;
1232:   PetscFunctionReturn(PETSC_SUCCESS);
1233: }

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

1238:   Not Collective

1240:   Input Parameters:
1241: + ds - The `PetscDS` object
1242: - f  - The field number

1244:   Output Parameter:
1245: . k - The highest derivative we need to tabulate

1247:   Level: developer

1249: .seealso: `PetscDS`, `PetscDSSetJetDegree()`, `PetscDSSetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1250: @*/
1251: PetscErrorCode PetscDSGetJetDegree(PetscDS ds, PetscInt f, PetscInt *k)
1252: {
1253:   PetscFunctionBegin;
1255:   PetscAssertPointer(k, 3);
1256:   PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
1257:   *k = ds->jetDegree[f];
1258:   PetscFunctionReturn(PETSC_SUCCESS);
1259: }

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

1264:   Not Collective

1266:   Input Parameters:
1267: + ds - The `PetscDS` object
1268: . f  - The field number
1269: - k  - The highest derivative we need to tabulate

1271:   Level: developer

1273: .seealso: `PetscDS`, `PetscDSGetJetDegree()`, `PetscDSSetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1274: @*/
1275: PetscErrorCode PetscDSSetJetDegree(PetscDS ds, PetscInt f, PetscInt k)
1276: {
1277:   PetscFunctionBegin;
1279:   PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
1280:   ds->jetDegree[f] = k;
1281:   PetscFunctionReturn(PETSC_SUCCESS);
1282: }

1284: /*@C
1285:   PetscDSGetObjective - Get the pointwise objective function for a given test field

1287:   Not Collective

1289:   Input Parameters:
1290: + ds - The `PetscDS`
1291: - f  - The test field number

1293:   Output Parameter:
1294: . obj - integrand for the test function term

1296:   Calling sequence of `obj`:
1297: + dim          - the coordinate dimension
1298: . Nf           - the number of fields
1299: . NfAux        - the number of auxiliary fields
1300: . uOff         - the offset into u[] and u_t[] for each field
1301: . uOff_x       - the offset into u_x[] for each field
1302: . u            - each field evaluated at the current point
1303: . u_t          - the time derivative of each field evaluated at the current point
1304: . u_x          - the gradient of each field evaluated at the current point
1305: . aOff         - the offset into a[] and a_t[] for each auxiliary field
1306: . aOff_x       - the offset into a_x[] for each auxiliary field
1307: . a            - each auxiliary field evaluated at the current point
1308: . a_t          - the time derivative of each auxiliary field evaluated at the current point
1309: . a_x          - the gradient of auxiliary each field evaluated at the current point
1310: . t            - current time
1311: . x            - coordinates of the current point
1312: . numConstants - number of constant parameters
1313: . constants    - constant parameters
1314: - obj          - output values at the current point

1316:   Level: intermediate

1318:   Note:
1319:   We are using a first order FEM model for the weak form\: $  \int_\Omega \phi\,\mathrm{obj}(u, u_t, \nabla u, x, t)$

1321: .seealso: `PetscDS`, `PetscDSSetObjective()`, `PetscDSGetResidual()`
1322: @*/
1323: PetscErrorCode PetscDSGetObjective(PetscDS ds, PetscInt f, void (**obj)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[]))
1324: {
1325:   PetscPointFunc *tmp;
1326:   PetscInt        n;

1328:   PetscFunctionBegin;
1330:   PetscAssertPointer(obj, 3);
1331:   PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
1332:   PetscCall(PetscWeakFormGetObjective(ds->wf, NULL, 0, f, 0, &n, &tmp));
1333:   *obj = tmp ? tmp[0] : NULL;
1334:   PetscFunctionReturn(PETSC_SUCCESS);
1335: }

1337: /*@C
1338:   PetscDSSetObjective - Set the pointwise objective function for a given test field

1340:   Not Collective

1342:   Input Parameters:
1343: + ds  - The `PetscDS`
1344: . f   - The test field number
1345: - obj - integrand for the test function term

1347:   Calling sequence of `obj`:
1348: + dim          - the coordinate dimension
1349: . Nf           - the number of fields
1350: . NfAux        - the number of auxiliary fields
1351: . uOff         - the offset into u[] and u_t[] for each field
1352: . uOff_x       - the offset into u_x[] for each field
1353: . u            - each field evaluated at the current point
1354: . u_t          - the time derivative of each field evaluated at the current point
1355: . u_x          - the gradient of each field evaluated at the current point
1356: . aOff         - the offset into a[] and a_t[] for each auxiliary field
1357: . aOff_x       - the offset into a_x[] for each auxiliary field
1358: . a            - each auxiliary field evaluated at the current point
1359: . a_t          - the time derivative of each auxiliary field evaluated at the current point
1360: . a_x          - the gradient of auxiliary each field evaluated at the current point
1361: . t            - current time
1362: . x            - coordinates of the current point
1363: . numConstants - number of constant parameters
1364: . constants    - constant parameters
1365: - obj          - output values at the current point

1367:   Level: intermediate

1369:   Note:
1370:   We are using a first order FEM model for the weak form\: $  \int_\Omega \phi\,\mathrm{obj}(u, u_t, \nabla u, x, t)$

1372: .seealso: `PetscDS`, `PetscDSGetObjective()`, `PetscDSSetResidual()`
1373: @*/
1374: PetscErrorCode PetscDSSetObjective(PetscDS ds, PetscInt f, void (*obj)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[]))
1375: {
1376:   PetscFunctionBegin;
1379:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
1380:   PetscCall(PetscWeakFormSetIndexObjective(ds->wf, NULL, 0, f, 0, 0, obj));
1381:   PetscFunctionReturn(PETSC_SUCCESS);
1382: }

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

1387:   Not Collective

1389:   Input Parameters:
1390: + ds - The `PetscDS`
1391: - f  - The test field number

1393:   Output Parameters:
1394: + f0 - integrand for the test function term
1395: - f1 - integrand for the test function gradient term

1397:   Calling sequence of `f0`:
1398: + dim          - the coordinate dimension
1399: . Nf           - the number of fields
1400: . NfAux        - the number of auxiliary fields
1401: . uOff         - the offset into u[] and u_t[] for each field
1402: . uOff_x       - the offset into u_x[] for each field
1403: . u            - each field evaluated at the current point
1404: . u_t          - the time derivative of each field evaluated at the current point
1405: . u_x          - the gradient of each field evaluated at the current point
1406: . aOff         - the offset into a[] and a_t[] for each auxiliary field
1407: . aOff_x       - the offset into a_x[] for each auxiliary field
1408: . a            - each auxiliary field evaluated at the current point
1409: . a_t          - the time derivative of each auxiliary field evaluated at the current point
1410: . a_x          - the gradient of auxiliary each field evaluated at the current point
1411: . t            - current time
1412: . x            - coordinates of the current point
1413: . numConstants - number of constant parameters
1414: . constants    - constant parameters
1415: - f0           - output values at the current point

1417:   Level: intermediate

1419:   Note:
1420:   `f1` has an identical form and is omitted for brevity.

1422:   We are using a first order FEM model for the weak form\: $  \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)$

1424: .seealso: `PetscDS`, `PetscDSSetResidual()`
1425: @*/
1426: PetscErrorCode PetscDSGetResidual(PetscDS ds, PetscInt f, void (**f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]), void (**f1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))
1427: {
1428:   PetscPointFunc *tmp0, *tmp1;
1429:   PetscInt        n0, n1;

1431:   PetscFunctionBegin;
1433:   PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
1434:   PetscCall(PetscWeakFormGetResidual(ds->wf, NULL, 0, f, 0, &n0, &tmp0, &n1, &tmp1));
1435:   *f0 = tmp0 ? tmp0[0] : NULL;
1436:   *f1 = tmp1 ? tmp1[0] : NULL;
1437:   PetscFunctionReturn(PETSC_SUCCESS);
1438: }

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

1443:   Not Collective

1445:   Input Parameters:
1446: + ds - The `PetscDS`
1447: . f  - The test field number
1448: . f0 - integrand for the test function term
1449: - f1 - integrand for the test function gradient term

1451:   Calling sequence of `f0`:
1452: + dim          - the coordinate dimension
1453: . Nf           - the number of fields
1454: . NfAux        - the number of auxiliary fields
1455: . uOff         - the offset into u[] and u_t[] for each field
1456: . uOff_x       - the offset into u_x[] for each field
1457: . u            - each field evaluated at the current point
1458: . u_t          - the time derivative of each field evaluated at the current point
1459: . u_x          - the gradient of each field evaluated at the current point
1460: . aOff         - the offset into a[] and a_t[] for each auxiliary field
1461: . aOff_x       - the offset into a_x[] for each auxiliary field
1462: . a            - each auxiliary field evaluated at the current point
1463: . a_t          - the time derivative of each auxiliary field evaluated at the current point
1464: . a_x          - the gradient of auxiliary each field evaluated at the current point
1465: . t            - current time
1466: . x            - coordinates of the current point
1467: . numConstants - number of constant parameters
1468: . constants    - constant parameters
1469: - f0           - output values at the current point

1471:   Level: intermediate

1473:   Note:
1474:   `f1` has an identical form and is omitted for brevity.

1476:   We are using a first order FEM model for the weak form\: $  \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)$

1478: .seealso: `PetscDS`, `PetscDSGetResidual()`
1479: @*/
1480: PetscErrorCode PetscDSSetResidual(PetscDS ds, PetscInt f, void (*f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]), void (*f1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))
1481: {
1482:   PetscFunctionBegin;
1486:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
1487:   PetscCall(PetscWeakFormSetIndexResidual(ds->wf, NULL, 0, f, 0, 0, f0, 0, f1));
1488:   PetscFunctionReturn(PETSC_SUCCESS);
1489: }

1491: /*@C
1492:   PetscDSGetRHSResidual - Get the pointwise RHS residual function for explicit timestepping for a given test field

1494:   Not Collective

1496:   Input Parameters:
1497: + ds - The `PetscDS`
1498: - f  - The test field number

1500:   Output Parameters:
1501: + f0 - integrand for the test function term
1502: - f1 - integrand for the test function gradient term

1504:   Calling sequence of `f0`:
1505: + dim          - the coordinate dimension
1506: . Nf           - the number of fields
1507: . NfAux        - the number of auxiliary fields
1508: . uOff         - the offset into u[] and u_t[] for each field
1509: . uOff_x       - the offset into u_x[] for each field
1510: . u            - each field evaluated at the current point
1511: . u_t          - the time derivative of each field evaluated at the current point
1512: . u_x          - the gradient of each field evaluated at the current point
1513: . aOff         - the offset into a[] and a_t[] for each auxiliary field
1514: . aOff_x       - the offset into a_x[] for each auxiliary field
1515: . a            - each auxiliary field evaluated at the current point
1516: . a_t          - the time derivative of each auxiliary field evaluated at the current point
1517: . a_x          - the gradient of auxiliary each field evaluated at the current point
1518: . t            - current time
1519: . x            - coordinates of the current point
1520: . numConstants - number of constant parameters
1521: . constants    - constant parameters
1522: - f0           - output values at the current point

1524:   Level: intermediate

1526:   Note:
1527:   `f1` has an identical form and is omitted for brevity.

1529:   We are using a first order FEM model for the weak form\: $ \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)$

1531: .seealso: `PetscDS`, `PetscDSSetRHSResidual()`
1532: @*/
1533: PetscErrorCode PetscDSGetRHSResidual(PetscDS ds, PetscInt f, void (**f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]), void (**f1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))
1534: {
1535:   PetscPointFunc *tmp0, *tmp1;
1536:   PetscInt        n0, n1;

1538:   PetscFunctionBegin;
1540:   PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
1541:   PetscCall(PetscWeakFormGetResidual(ds->wf, NULL, 0, f, 100, &n0, &tmp0, &n1, &tmp1));
1542:   *f0 = tmp0 ? tmp0[0] : NULL;
1543:   *f1 = tmp1 ? tmp1[0] : NULL;
1544:   PetscFunctionReturn(PETSC_SUCCESS);
1545: }

1547: /*@C
1548:   PetscDSSetRHSResidual - Set the pointwise residual function for explicit timestepping for a given test field

1550:   Not Collective

1552:   Input Parameters:
1553: + ds - The `PetscDS`
1554: . f  - The test field number
1555: . f0 - integrand for the test function term
1556: - f1 - integrand for the test function gradient term

1558:   Calling sequence for the callbacks `f0`:
1559: + dim          - the coordinate dimension
1560: . Nf           - the number of fields
1561: . NfAux        - the number of auxiliary fields
1562: . uOff         - the offset into u[] and u_t[] for each field
1563: . uOff_x       - the offset into u_x[] for each field
1564: . u            - each field evaluated at the current point
1565: . u_t          - the time derivative of each field evaluated at the current point
1566: . u_x          - the gradient of each field evaluated at the current point
1567: . aOff         - the offset into a[] and a_t[] for each auxiliary field
1568: . aOff_x       - the offset into a_x[] for each auxiliary field
1569: . a            - each auxiliary field evaluated at the current point
1570: . a_t          - the time derivative of each auxiliary field evaluated at the current point
1571: . a_x          - the gradient of auxiliary each field evaluated at the current point
1572: . t            - current time
1573: . x            - coordinates of the current point
1574: . numConstants - number of constant parameters
1575: . constants    - constant parameters
1576: - f0           - output values at the current point

1578:   Level: intermediate

1580:   Note:
1581:   `f1` has an identical form and is omitted for brevity.

1583:   We are using a first order FEM model for the weak form\: $ \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)$

1585: .seealso: `PetscDS`, `PetscDSGetResidual()`
1586: @*/
1587: PetscErrorCode PetscDSSetRHSResidual(PetscDS ds, PetscInt f, void (*f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]), void (*f1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))
1588: {
1589:   PetscFunctionBegin;
1593:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
1594:   PetscCall(PetscWeakFormSetIndexResidual(ds->wf, NULL, 0, f, 100, 0, f0, 0, f1));
1595:   PetscFunctionReturn(PETSC_SUCCESS);
1596: }

1598: /*@
1599:   PetscDSHasJacobian - Checks that the Jacobian functions have been set

1601:   Not Collective

1603:   Input Parameter:
1604: . ds - The `PetscDS`

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

1609:   Level: intermediate

1611: .seealso: `PetscDS`, `PetscDSGetJacobianPreconditioner()`, `PetscDSSetJacobianPreconditioner()`, `PetscDSGetJacobian()`
1612: @*/
1613: PetscErrorCode PetscDSHasJacobian(PetscDS ds, PetscBool *hasJac)
1614: {
1615:   PetscFunctionBegin;
1617:   PetscCall(PetscWeakFormHasJacobian(ds->wf, hasJac));
1618:   PetscFunctionReturn(PETSC_SUCCESS);
1619: }

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

1624:   Not Collective

1626:   Input Parameters:
1627: + ds - The `PetscDS`
1628: . f  - The test field number
1629: - g  - The field number

1631:   Output Parameters:
1632: + g0 - integrand for the test and basis function term
1633: . g1 - integrand for the test function and basis function gradient term
1634: . g2 - integrand for the test function gradient and basis function term
1635: - g3 - integrand for the test function gradient and basis function gradient term

1637:   Calling sequence of `g0`:
1638: + dim          - the coordinate dimension
1639: . Nf           - the number of fields
1640: . NfAux        - the number of auxiliary fields
1641: . uOff         - the offset into u[] and u_t[] for each field
1642: . uOff_x       - the offset into u_x[] for each field
1643: . u            - each field evaluated at the current point
1644: . u_t          - the time derivative of each field evaluated at the current point
1645: . u_x          - the gradient of each field evaluated at the current point
1646: . aOff         - the offset into a[] and a_t[] for each auxiliary field
1647: . aOff_x       - the offset into a_x[] for each auxiliary field
1648: . a            - each auxiliary field evaluated at the current point
1649: . a_t          - the time derivative of each auxiliary field evaluated at the current point
1650: . a_x          - the gradient of auxiliary each field evaluated at the current point
1651: . t            - current time
1652: . u_tShift     - the multiplier a for dF/dU_t
1653: . x            - coordinates of the current point
1654: . numConstants - number of constant parameters
1655: . constants    - constant parameters
1656: - g0           - output values at the current point

1658:   Level: intermediate

1660:   Note:
1661:   `g1`, `g2`, and `g3` have identical calling sequences to `g0` and are omitted for brevity.

1663:   We are using a first order FEM model for the weak form\:

1665:   $$
1666:   \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
1667:   + \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
1668:   $$

1670: .seealso: `PetscDS`, `PetscDSSetJacobian()`
1671: @*/
1672: PetscErrorCode PetscDSGetJacobian(PetscDS ds, PetscInt f, PetscInt g, void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]), void (**g1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (**g2)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (**g3)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))
1673: {
1674:   PetscPointJac *tmp0, *tmp1, *tmp2, *tmp3;
1675:   PetscInt       n0, n1, n2, n3;

1677:   PetscFunctionBegin;
1679:   PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
1680:   PetscCheck(!(g < 0) && !(g >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", g, ds->Nf);
1681:   PetscCall(PetscWeakFormGetJacobian(ds->wf, NULL, 0, f, g, 0, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3));
1682:   *g0 = tmp0 ? tmp0[0] : NULL;
1683:   *g1 = tmp1 ? tmp1[0] : NULL;
1684:   *g2 = tmp2 ? tmp2[0] : NULL;
1685:   *g3 = tmp3 ? tmp3[0] : NULL;
1686:   PetscFunctionReturn(PETSC_SUCCESS);
1687: }

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

1692:   Not Collective

1694:   Input Parameters:
1695: + ds - The `PetscDS`
1696: . f  - The test field number
1697: . g  - The field number
1698: . g0 - integrand for the test and basis function term
1699: . g1 - integrand for the test function and basis function gradient term
1700: . g2 - integrand for the test function gradient and basis function term
1701: - g3 - integrand for the test function gradient and basis function gradient term

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

1724:   Level: intermediate

1726:   Note:
1727:   `g1`, `g2`, and `g3` have identical calling sequences to `g0` and are omitted for brevity.

1729:   We are using a first order FEM model for the weak form\:

1731:   $$
1732:   \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
1733:   + \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
1734:   $$

1736: .seealso: `PetscDS`, `PetscDSGetJacobian()`
1737: @*/
1738: PetscErrorCode PetscDSSetJacobian(PetscDS ds, PetscInt f, PetscInt g, void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]), void (*g1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (*g2)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (*g3)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))
1739: {
1740:   PetscFunctionBegin;
1746:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
1747:   PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", g);
1748:   PetscCall(PetscWeakFormSetIndexJacobian(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3));
1749:   PetscFunctionReturn(PETSC_SUCCESS);
1750: }

1752: /*@
1753:   PetscDSUseJacobianPreconditioner - Set whether to construct a Jacobian preconditioner

1755:   Not Collective

1757:   Input Parameters:
1758: + prob      - The `PetscDS`
1759: - useJacPre - flag that enables construction of a Jacobian preconditioner

1761:   Level: intermediate

1763:   Developer Note:
1764:   Should be called `PetscDSSetUseJacobianPreconditioner()`

1766: .seealso: `PetscDS`, `PetscDSGetJacobianPreconditioner()`, `PetscDSSetJacobianPreconditioner()`, `PetscDSGetJacobian()`
1767: @*/
1768: PetscErrorCode PetscDSUseJacobianPreconditioner(PetscDS prob, PetscBool useJacPre)
1769: {
1770:   PetscFunctionBegin;
1772:   prob->useJacPre = useJacPre;
1773:   PetscFunctionReturn(PETSC_SUCCESS);
1774: }

1776: /*@
1777:   PetscDSHasJacobianPreconditioner - Checks if a Jacobian preconditioner matrix has been set

1779:   Not Collective

1781:   Input Parameter:
1782: . ds - The `PetscDS`

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

1787:   Level: intermediate

1789: .seealso: `PetscDS`, `PetscDSGetJacobianPreconditioner()`, `PetscDSSetJacobianPreconditioner()`, `PetscDSGetJacobian()`
1790: @*/
1791: PetscErrorCode PetscDSHasJacobianPreconditioner(PetscDS ds, PetscBool *hasJacPre)
1792: {
1793:   PetscFunctionBegin;
1795:   *hasJacPre = PETSC_FALSE;
1796:   if (!ds->useJacPre) PetscFunctionReturn(PETSC_SUCCESS);
1797:   PetscCall(PetscWeakFormHasJacobianPreconditioner(ds->wf, hasJacPre));
1798:   PetscFunctionReturn(PETSC_SUCCESS);
1799: }

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

1805:   Not Collective

1807:   Input Parameters:
1808: + ds - The `PetscDS`
1809: . f  - The test field number
1810: - g  - The field number

1812:   Output Parameters:
1813: + g0 - integrand for the test and basis function term
1814: . g1 - integrand for the test function and basis function gradient term
1815: . g2 - integrand for the test function gradient and basis function term
1816: - g3 - integrand for the test function gradient and basis function gradient term

1818:   Calling sequence of `g0`:
1819: + dim          - the coordinate dimension
1820: . Nf           - the number of fields
1821: . NfAux        - the number of auxiliary fields
1822: . uOff         - the offset into u[] and u_t[] for each field
1823: . uOff_x       - the offset into u_x[] for each field
1824: . u            - each field evaluated at the current point
1825: . u_t          - the time derivative of each field evaluated at the current point
1826: . u_x          - the gradient of each field evaluated at the current point
1827: . aOff         - the offset into a[] and a_t[] for each auxiliary field
1828: . aOff_x       - the offset into a_x[] for each auxiliary field
1829: . a            - each auxiliary field evaluated at the current point
1830: . a_t          - the time derivative of each auxiliary field evaluated at the current point
1831: . a_x          - the gradient of auxiliary each field evaluated at the current point
1832: . t            - current time
1833: . u_tShift     - the multiplier a for dF/dU_t
1834: . x            - coordinates of the current point
1835: . numConstants - number of constant parameters
1836: . constants    - constant parameters
1837: - g0           - output values at the current point

1839:   Level: intermediate

1841:   Note:
1842:   `g1`, `g2`, and `g3` have identical calling sequences to `g0` and are omitted for brevity.
1843:   We are using a first order FEM model for the weak form\:

1845:   $$
1846:   \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
1847:   + \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
1848:   $$

1850: .seealso: `PetscDS`, `PetscDSSetJacobianPreconditioner()`, `PetscDSGetJacobian()`
1851: @*/
1852: PetscErrorCode PetscDSGetJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g, void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]), void (**g1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (**g2)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (**g3)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))
1853: {
1854:   PetscPointJac *tmp0, *tmp1, *tmp2, *tmp3;
1855:   PetscInt       n0, n1, n2, n3;

1857:   PetscFunctionBegin;
1859:   PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
1860:   PetscCheck(!(g < 0) && !(g >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", g, ds->Nf);
1861:   PetscCall(PetscWeakFormGetJacobianPreconditioner(ds->wf, NULL, 0, f, g, 0, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3));
1862:   *g0 = tmp0 ? tmp0[0] : NULL;
1863:   *g1 = tmp1 ? tmp1[0] : NULL;
1864:   *g2 = tmp2 ? tmp2[0] : NULL;
1865:   *g3 = tmp3 ? tmp3[0] : NULL;
1866:   PetscFunctionReturn(PETSC_SUCCESS);
1867: }

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

1873:   Not Collective

1875:   Input Parameters:
1876: + ds - The `PetscDS`
1877: . f  - The test field number
1878: . g  - The field number
1879: . g0 - integrand for the test and basis function term
1880: . g1 - integrand for the test function and basis function gradient term
1881: . g2 - integrand for the test function gradient and basis function term
1882: - g3 - integrand for the test function gradient and basis function gradient term

1884:   Calling sequence of `g0`:
1885: + dim          - the coordinate dimension
1886: . Nf           - the number of fields
1887: . NfAux        - the number of auxiliary fields
1888: . uOff         - the offset into u[] and u_t[] for each field
1889: . uOff_x       - the offset into u_x[] for each field
1890: . u            - each field evaluated at the current point
1891: . u_t          - the time derivative of each field evaluated at the current point
1892: . u_x          - the gradient of each field evaluated at the current point
1893: . aOff         - the offset into a[] and a_t[] for each auxiliary field
1894: . aOff_x       - the offset into a_x[] for each auxiliary field
1895: . a            - each auxiliary field evaluated at the current point
1896: . a_t          - the time derivative of each auxiliary field evaluated at the current point
1897: . a_x          - the gradient of auxiliary each field evaluated at the current point
1898: . t            - current time
1899: . u_tShift     - the multiplier a for dF/dU_t
1900: . x            - coordinates of the current point
1901: . numConstants - number of constant parameters
1902: . constants    - constant parameters
1903: - g0           - output values at the current point

1905:   Level: intermediate

1907:   Note:
1908:   `g1`, `g2`, and `g3` have identical calling sequences to `g0` and are omitted for brevity.

1910:   We are using a first order FEM model for the weak form\:

1912:   $$
1913:   \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
1914:   + \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
1915:   $$

1917: .seealso: `PetscDS`, `PetscDSGetJacobianPreconditioner()`, `PetscDSSetJacobian()`
1918: @*/
1919: PetscErrorCode PetscDSSetJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g, void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]), void (*g1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (*g2)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (*g3)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))
1920: {
1921:   PetscFunctionBegin;
1927:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
1928:   PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", g);
1929:   PetscCall(PetscWeakFormSetIndexJacobianPreconditioner(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3));
1930:   PetscFunctionReturn(PETSC_SUCCESS);
1931: }

1933: /*@
1934:   PetscDSHasDynamicJacobian - Signals that a dynamic Jacobian, dF/du_t, has been set

1936:   Not Collective

1938:   Input Parameter:
1939: . ds - The `PetscDS`

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

1944:   Level: intermediate

1946: .seealso: `PetscDS`, `PetscDSGetDynamicJacobian()`, `PetscDSSetDynamicJacobian()`, `PetscDSGetJacobian()`
1947: @*/
1948: PetscErrorCode PetscDSHasDynamicJacobian(PetscDS ds, PetscBool *hasDynJac)
1949: {
1950:   PetscFunctionBegin;
1952:   PetscCall(PetscWeakFormHasDynamicJacobian(ds->wf, hasDynJac));
1953:   PetscFunctionReturn(PETSC_SUCCESS);
1954: }

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

1959:   Not Collective

1961:   Input Parameters:
1962: + ds - The `PetscDS`
1963: . f  - The test field number
1964: - g  - The field number

1966:   Output Parameters:
1967: + g0 - integrand for the test and basis function term
1968: . g1 - integrand for the test function and basis function gradient term
1969: . g2 - integrand for the test function gradient and basis function term
1970: - g3 - integrand for the test function gradient and basis function gradient term

1972:   Calling sequence of `g0`:
1973: + dim          - the coordinate dimension
1974: . Nf           - the number of fields
1975: . NfAux        - the number of auxiliary fields
1976: . uOff         - the offset into u[] and u_t[] for each field
1977: . uOff_x       - the offset into u_x[] for each field
1978: . u            - each field evaluated at the current point
1979: . u_t          - the time derivative of each field evaluated at the current point
1980: . u_x          - the gradient of each field evaluated at the current point
1981: . aOff         - the offset into a[] and a_t[] for each auxiliary field
1982: . aOff_x       - the offset into a_x[] for each auxiliary field
1983: . a            - each auxiliary field evaluated at the current point
1984: . a_t          - the time derivative of each auxiliary field evaluated at the current point
1985: . a_x          - the gradient of auxiliary each field evaluated at the current point
1986: . t            - current time
1987: . u_tShift     - the multiplier a for dF/dU_t
1988: . x            - coordinates of the current point
1989: . numConstants - number of constant parameters
1990: . constants    - constant parameters
1991: - g0           - output values at the current point

1993:   Level: intermediate

1995:   Note:
1996:   `g1`, `g2`, and `g3` have identical calling sequences to `g0` and are omitted for brevity.

1998:   We are using a first order FEM model for the weak form\:

2000:   $$
2001:   \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
2002:   + \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
2003:   $$

2005: .seealso: `PetscDS`, `PetscDSSetJacobian()`
2006: @*/
2007: PetscErrorCode PetscDSGetDynamicJacobian(PetscDS ds, PetscInt f, PetscInt g, void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]), void (**g1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (**g2)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (**g3)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))
2008: {
2009:   PetscPointJac *tmp0, *tmp1, *tmp2, *tmp3;
2010:   PetscInt       n0, n1, n2, n3;

2012:   PetscFunctionBegin;
2014:   PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
2015:   PetscCheck(!(g < 0) && !(g >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", g, ds->Nf);
2016:   PetscCall(PetscWeakFormGetDynamicJacobian(ds->wf, NULL, 0, f, g, 0, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3));
2017:   *g0 = tmp0 ? tmp0[0] : NULL;
2018:   *g1 = tmp1 ? tmp1[0] : NULL;
2019:   *g2 = tmp2 ? tmp2[0] : NULL;
2020:   *g3 = tmp3 ? tmp3[0] : NULL;
2021:   PetscFunctionReturn(PETSC_SUCCESS);
2022: }

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

2027:   Not Collective

2029:   Input Parameters:
2030: + ds - The `PetscDS`
2031: . f  - The test field number
2032: . g  - The field number
2033: . g0 - integrand for the test and basis function term
2034: . g1 - integrand for the test function and basis function gradient term
2035: . g2 - integrand for the test function gradient and basis function term
2036: - g3 - integrand for the test function gradient and basis function gradient term

2038:   Calling sequence of `g0`:
2039: + dim          - the coordinate dimension
2040: . Nf           - the number of fields
2041: . NfAux        - the number of auxiliary fields
2042: . uOff         - the offset into u[] and u_t[] for each field
2043: . uOff_x       - the offset into u_x[] for each field
2044: . u            - each field evaluated at the current point
2045: . u_t          - the time derivative of each field evaluated at the current point
2046: . u_x          - the gradient of each field evaluated at the current point
2047: . aOff         - the offset into a[] and a_t[] for each auxiliary field
2048: . aOff_x       - the offset into a_x[] for each auxiliary field
2049: . a            - each auxiliary field evaluated at the current point
2050: . a_t          - the time derivative of each auxiliary field evaluated at the current point
2051: . a_x          - the gradient of auxiliary each field evaluated at the current point
2052: . t            - current time
2053: . u_tShift     - the multiplier a for dF/dU_t
2054: . x            - coordinates of the current point
2055: . numConstants - number of constant parameters
2056: . constants    - constant parameters
2057: - g0           - output values at the current point

2059:   Level: intermediate

2061:   Note:
2062:   `g1`, `g2`, and `g3` have identical calling sequences to `g0` and are omitted for brevity.

2064:   We are using a first order FEM model for the weak form\:

2066:   $$
2067:   \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
2068:   + \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
2069:   $$

2071: .seealso: `PetscDS`, `PetscDSGetJacobian()`
2072: @*/
2073: PetscErrorCode PetscDSSetDynamicJacobian(PetscDS ds, PetscInt f, PetscInt g, void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]), void (*g1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (*g2)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (*g3)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))
2074: {
2075:   PetscFunctionBegin;
2081:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2082:   PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", g);
2083:   PetscCall(PetscWeakFormSetIndexDynamicJacobian(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3));
2084:   PetscFunctionReturn(PETSC_SUCCESS);
2085: }

2087: /*@C
2088:   PetscDSGetRiemannSolver - Returns the Riemann solver for the given field

2090:   Not Collective

2092:   Input Parameters:
2093: + ds - The `PetscDS` object
2094: - f  - The field number

2096:   Output Parameter:
2097: . r - Riemann solver

2099:   Calling sequence of `r`:
2100: + dim          - the coordinate dimension
2101: . Nf           - The number of fields
2102: . x            - The coordinates at a point on the interface
2103: . n            - The normal vector to the interface
2104: . uL           - The state vector to the left of the interface
2105: . uR           - The state vector to the right of the interface
2106: . flux         - output array of flux through the interface
2107: . numConstants - number of constant parameters
2108: . constants    - constant parameters
2109: - ctx          - optional user context

2111:   Level: intermediate

2113: .seealso: `PetscDS`, `PetscDSSetRiemannSolver()`
2114: @*/
2115: PetscErrorCode PetscDSGetRiemannSolver(PetscDS ds, PetscInt f, 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))
2116: {
2117:   PetscRiemannFunc *tmp;
2118:   PetscInt          n;

2120:   PetscFunctionBegin;
2122:   PetscAssertPointer(r, 3);
2123:   PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
2124:   PetscCall(PetscWeakFormGetRiemannSolver(ds->wf, NULL, 0, f, 0, &n, &tmp));
2125:   *r = tmp ? tmp[0] : NULL;
2126:   PetscFunctionReturn(PETSC_SUCCESS);
2127: }

2129: /*@C
2130:   PetscDSSetRiemannSolver - Sets the Riemann solver for the given field

2132:   Not Collective

2134:   Input Parameters:
2135: + ds - The `PetscDS` object
2136: . f  - The field number
2137: - r  - Riemann solver

2139:   Calling sequence of `r`:
2140: + dim          - the coordinate dimension
2141: . Nf           - The number of fields
2142: . x            - The coordinates at a point on the interface
2143: . n            - The normal vector to the interface
2144: . uL           - The state vector to the left of the interface
2145: . uR           - The state vector to the right of the interface
2146: . flux         - output array of flux through the interface
2147: . numConstants - number of constant parameters
2148: . constants    - constant parameters
2149: - ctx          - optional user context

2151:   Level: intermediate

2153: .seealso: `PetscDS`, `PetscDSGetRiemannSolver()`
2154: @*/
2155: PetscErrorCode PetscDSSetRiemannSolver(PetscDS ds, PetscInt f, 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))
2156: {
2157:   PetscFunctionBegin;
2160:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2161:   PetscCall(PetscWeakFormSetIndexRiemannSolver(ds->wf, NULL, 0, f, 0, 0, r));
2162:   PetscFunctionReturn(PETSC_SUCCESS);
2163: }

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

2168:   Not Collective

2170:   Input Parameters:
2171: + ds - The `PetscDS`
2172: - f  - The field number

2174:   Output Parameter:
2175: . update - update function

2177:   Calling sequence of `update`:
2178: + dim          - the coordinate dimension
2179: . Nf           - the number of fields
2180: . NfAux        - the number of auxiliary fields
2181: . uOff         - the offset into u[] and u_t[] for each field
2182: . uOff_x       - the offset into u_x[] for each field
2183: . u            - each field evaluated at the current point
2184: . u_t          - the time derivative of each field evaluated at the current point
2185: . u_x          - the gradient of each field evaluated at the current point
2186: . aOff         - the offset into a[] and a_t[] for each auxiliary field
2187: . aOff_x       - the offset into a_x[] for each auxiliary field
2188: . a            - each auxiliary field evaluated at the current point
2189: . a_t          - the time derivative of each auxiliary field evaluated at the current point
2190: . a_x          - the gradient of auxiliary each field evaluated at the current point
2191: . t            - current time
2192: . x            - coordinates of the current point
2193: . numConstants - number of constant parameters
2194: . constants    - constant parameters
2195: - uNew         - new value for field at the current point

2197:   Level: intermediate

2199: .seealso: `PetscDS`, `PetscDSSetUpdate()`, `PetscDSSetResidual()`
2200: @*/
2201: PetscErrorCode PetscDSGetUpdate(PetscDS ds, PetscInt f, void (**update)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar uNew[]))
2202: {
2203:   PetscFunctionBegin;
2205:   PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
2206:   if (update) {
2207:     PetscAssertPointer(update, 3);
2208:     *update = ds->update[f];
2209:   }
2210:   PetscFunctionReturn(PETSC_SUCCESS);
2211: }

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

2216:   Not Collective

2218:   Input Parameters:
2219: + ds     - The `PetscDS`
2220: . f      - The field number
2221: - update - update function

2223:   Calling sequence of `update`:
2224: + dim          - the coordinate dimension
2225: . Nf           - the number of fields
2226: . NfAux        - the number of auxiliary fields
2227: . uOff         - the offset into u[] and u_t[] for each field
2228: . uOff_x       - the offset into u_x[] for each field
2229: . u            - each field evaluated at the current point
2230: . u_t          - the time derivative of each field evaluated at the current point
2231: . u_x          - the gradient of each field evaluated at the current point
2232: . aOff         - the offset into a[] and a_t[] for each auxiliary field
2233: . aOff_x       - the offset into a_x[] for each auxiliary field
2234: . a            - each auxiliary field evaluated at the current point
2235: . a_t          - the time derivative of each auxiliary field evaluated at the current point
2236: . a_x          - the gradient of auxiliary each field evaluated at the current point
2237: . t            - current time
2238: . x            - coordinates of the current point
2239: . numConstants - number of constant parameters
2240: . constants    - constant parameters
2241: - uNew         - new field values at the current point

2243:   Level: intermediate

2245: .seealso: `PetscDS`, `PetscDSGetResidual()`
2246: @*/
2247: PetscErrorCode PetscDSSetUpdate(PetscDS ds, PetscInt f, void (*update)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar uNew[]))
2248: {
2249:   PetscFunctionBegin;
2252:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2253:   PetscCall(PetscDSEnlarge_Static(ds, f + 1));
2254:   ds->update[f] = update;
2255:   PetscFunctionReturn(PETSC_SUCCESS);
2256: }

2258: PetscErrorCode PetscDSGetContext(PetscDS ds, PetscInt f, void *ctx)
2259: {
2260:   PetscFunctionBegin;
2262:   PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
2263:   PetscAssertPointer(ctx, 3);
2264:   *(void **)ctx = ds->ctx[f];
2265:   PetscFunctionReturn(PETSC_SUCCESS);
2266: }

2268: PetscErrorCode PetscDSSetContext(PetscDS ds, PetscInt f, void *ctx)
2269: {
2270:   PetscFunctionBegin;
2272:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2273:   PetscCall(PetscDSEnlarge_Static(ds, f + 1));
2274:   ds->ctx[f] = ctx;
2275:   PetscFunctionReturn(PETSC_SUCCESS);
2276: }

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

2281:   Not Collective

2283:   Input Parameters:
2284: + ds - The PetscDS
2285: - f  - The test field number

2287:   Output Parameters:
2288: + f0 - boundary integrand for the test function term
2289: - f1 - boundary integrand for the test function gradient term

2291:   Calling sequence of `f0`:
2292: + dim          - the coordinate dimension
2293: . Nf           - the number of fields
2294: . NfAux        - the number of auxiliary fields
2295: . uOff         - the offset into u[] and u_t[] for each field
2296: . uOff_x       - the offset into u_x[] for each field
2297: . u            - each field evaluated at the current point
2298: . u_t          - the time derivative of each field evaluated at the current point
2299: . u_x          - the gradient of each field evaluated at the current point
2300: . aOff         - the offset into a[] and a_t[] for each auxiliary field
2301: . aOff_x       - the offset into a_x[] for each auxiliary field
2302: . a            - each auxiliary field evaluated at the current point
2303: . a_t          - the time derivative of each auxiliary field evaluated at the current point
2304: . a_x          - the gradient of auxiliary each field evaluated at the current point
2305: . t            - current time
2306: . x            - coordinates of the current point
2307: . n            - unit normal at the current point
2308: . numConstants - number of constant parameters
2309: . constants    - constant parameters
2310: - f0           - output values at the current point

2312:   Level: intermediate

2314:   Note:
2315:   The calling sequence of `f1` is identical, and therefore omitted for brevity.

2317:   We are using a first order FEM model for the weak form\:

2319:   $$
2320:   \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
2321:   $$

2323: .seealso: `PetscDS`, `PetscDSSetBdResidual()`
2324: @*/
2325: PetscErrorCode PetscDSGetBdResidual(PetscDS ds, PetscInt f, void (**f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]), void (**f1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))
2326: {
2327:   PetscBdPointFunc *tmp0, *tmp1;
2328:   PetscInt          n0, n1;

2330:   PetscFunctionBegin;
2332:   PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
2333:   PetscCall(PetscWeakFormGetBdResidual(ds->wf, NULL, 0, f, 0, &n0, &tmp0, &n1, &tmp1));
2334:   *f0 = tmp0 ? tmp0[0] : NULL;
2335:   *f1 = tmp1 ? tmp1[0] : NULL;
2336:   PetscFunctionReturn(PETSC_SUCCESS);
2337: }

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

2342:   Not Collective

2344:   Input Parameters:
2345: + ds - The `PetscDS`
2346: . f  - The test field number
2347: . f0 - boundary integrand for the test function term
2348: - f1 - boundary integrand for the test function gradient term

2350:   Calling sequence of `f0`:
2351: + dim          - the coordinate dimension
2352: . Nf           - the number of fields
2353: . NfAux        - the number of auxiliary fields
2354: . uOff         - the offset into u[] and u_t[] for each field
2355: . uOff_x       - the offset into u_x[] for each field
2356: . u            - each field evaluated at the current point
2357: . u_t          - the time derivative of each field evaluated at the current point
2358: . u_x          - the gradient of each field evaluated at the current point
2359: . aOff         - the offset into a[] and a_t[] for each auxiliary field
2360: . aOff_x       - the offset into a_x[] for each auxiliary field
2361: . a            - each auxiliary field evaluated at the current point
2362: . a_t          - the time derivative of each auxiliary field evaluated at the current point
2363: . a_x          - the gradient of auxiliary each field evaluated at the current point
2364: . t            - current time
2365: . x            - coordinates of the current point
2366: . n            - unit normal at the current point
2367: . numConstants - number of constant parameters
2368: . constants    - constant parameters
2369: - f0           - output values at the current point

2371:   Level: intermediate

2373:   Note:
2374:   The calling sequence of `f1` is identical, and therefore omitted for brevity.

2376:   We are using a first order FEM model for the weak form\:

2378:   $$
2379:   \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
2380:   $$

2382: .seealso: `PetscDS`, `PetscDSGetBdResidual()`
2383: @*/
2384: PetscErrorCode PetscDSSetBdResidual(PetscDS ds, PetscInt f, void (*f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]), void (*f1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))
2385: {
2386:   PetscFunctionBegin;
2388:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2389:   PetscCall(PetscWeakFormSetIndexBdResidual(ds->wf, NULL, 0, f, 0, 0, f0, 0, f1));
2390:   PetscFunctionReturn(PETSC_SUCCESS);
2391: }

2393: /*@
2394:   PetscDSHasBdJacobian - Indicates that boundary Jacobian functions have been set

2396:   Not Collective

2398:   Input Parameter:
2399: . ds - The `PetscDS`

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

2404:   Level: intermediate

2406: .seealso: `PetscDS`, `PetscDSHasJacobian()`, `PetscDSSetBdJacobian()`, `PetscDSGetBdJacobian()`
2407: @*/
2408: PetscErrorCode PetscDSHasBdJacobian(PetscDS ds, PetscBool *hasBdJac)
2409: {
2410:   PetscFunctionBegin;
2412:   PetscAssertPointer(hasBdJac, 2);
2413:   PetscCall(PetscWeakFormHasBdJacobian(ds->wf, hasBdJac));
2414:   PetscFunctionReturn(PETSC_SUCCESS);
2415: }

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

2420:   Not Collective

2422:   Input Parameters:
2423: + ds - The `PetscDS`
2424: . f  - The test field number
2425: - g  - The field number

2427:   Output Parameters:
2428: + g0 - integrand for the test and basis function term
2429: . g1 - integrand for the test function and basis function gradient term
2430: . g2 - integrand for the test function gradient and basis function term
2431: - g3 - integrand for the test function gradient and basis function gradient term

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

2455:   Level: intermediate

2457:   Note:
2458:   `g1`, `g2`, and `g3` have identical calling sequences to `g0` and are omitted for brevity.

2460:   We are using a first order FEM model for the weak form\:

2462:   $$
2463:   \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
2464:   + \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
2465:   $$

2467: .seealso: `PetscDS`, `PetscDSSetBdJacobian()`
2468: @*/
2469: PetscErrorCode PetscDSGetBdJacobian(PetscDS ds, PetscInt f, PetscInt g, void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]), void (**g1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (**g2)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (**g3)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))
2470: {
2471:   PetscBdPointJac *tmp0, *tmp1, *tmp2, *tmp3;
2472:   PetscInt         n0, n1, n2, n3;

2474:   PetscFunctionBegin;
2476:   PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
2477:   PetscCheck(!(g < 0) && !(g >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", g, ds->Nf);
2478:   PetscCall(PetscWeakFormGetBdJacobian(ds->wf, NULL, 0, f, g, 0, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3));
2479:   *g0 = tmp0 ? tmp0[0] : NULL;
2480:   *g1 = tmp1 ? tmp1[0] : NULL;
2481:   *g2 = tmp2 ? tmp2[0] : NULL;
2482:   *g3 = tmp3 ? tmp3[0] : NULL;
2483:   PetscFunctionReturn(PETSC_SUCCESS);
2484: }

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

2489:   Not Collective

2491:   Input Parameters:
2492: + ds - The PetscDS
2493: . f  - The test field number
2494: . g  - The field number
2495: . g0 - integrand for the test and basis function term
2496: . g1 - integrand for the test function and basis function gradient term
2497: . g2 - integrand for the test function gradient and basis function term
2498: - g3 - integrand for the test function gradient and basis function gradient term

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

2522:   Level: intermediate

2524:   Note:
2525:   `g1`, `g2`, and `g3` have identical calling sequences to `g0` and are omitted for brevity.

2527:   We are using a first order FEM model for the weak form\:

2529:   $$
2530:   \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
2531:   + \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
2532:   $$

2534: .seealso: `PetscDS`, `PetscDSGetBdJacobian()`
2535: @*/
2536: PetscErrorCode PetscDSSetBdJacobian(PetscDS ds, PetscInt f, PetscInt g, void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]), void (*g1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (*g2)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (*g3)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))
2537: {
2538:   PetscFunctionBegin;
2544:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2545:   PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", g);
2546:   PetscCall(PetscWeakFormSetIndexBdJacobian(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3));
2547:   PetscFunctionReturn(PETSC_SUCCESS);
2548: }

2550: /*@
2551:   PetscDSHasBdJacobianPreconditioner - Signals that boundary Jacobian preconditioner functions have been set

2553:   Not Collective

2555:   Input Parameter:
2556: . ds - The `PetscDS`

2558:   Output Parameter:
2559: . hasBdJacPre - flag that pointwise function for the boundary Jacobian preconditioner has been set

2561:   Level: intermediate

2563: .seealso: `PetscDS`, `PetscDSHasJacobian()`, `PetscDSSetBdJacobian()`, `PetscDSGetBdJacobian()`
2564: @*/
2565: PetscErrorCode PetscDSHasBdJacobianPreconditioner(PetscDS ds, PetscBool *hasBdJacPre)
2566: {
2567:   PetscFunctionBegin;
2569:   PetscAssertPointer(hasBdJacPre, 2);
2570:   PetscCall(PetscWeakFormHasBdJacobianPreconditioner(ds->wf, hasBdJacPre));
2571:   PetscFunctionReturn(PETSC_SUCCESS);
2572: }

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

2577:   Not Collective; No Fortran Support

2579:   Input Parameters:
2580: + ds - The `PetscDS`
2581: . f  - The test field number
2582: - g  - The field number

2584:   Output Parameters:
2585: + g0 - integrand for the test and basis function term
2586: . g1 - integrand for the test function and basis function gradient term
2587: . g2 - integrand for the test function gradient and basis function term
2588: - g3 - integrand for the test function gradient and basis function gradient term

2590:   Calling sequence of `g0`:
2591: + dim          - the coordinate dimension
2592: . Nf           - the number of fields
2593: . NfAux        - the number of auxiliary fields
2594: . uOff         - the offset into u[] and u_t[] for each field
2595: . uOff_x       - the offset into u_x[] for each field
2596: . u            - each field evaluated at the current point
2597: . u_t          - the time derivative of each field evaluated at the current point
2598: . u_x          - the gradient of each field evaluated at the current point
2599: . aOff         - the offset into a[] and a_t[] for each auxiliary field
2600: . aOff_x       - the offset into a_x[] for each auxiliary field
2601: . a            - each auxiliary field evaluated at the current point
2602: . a_t          - the time derivative of each auxiliary field evaluated at the current point
2603: . a_x          - the gradient of auxiliary each field evaluated at the current point
2604: . t            - current time
2605: . u_tShift     - the multiplier a for dF/dU_t
2606: . x            - coordinates of the current point
2607: . n            - normal at the current point
2608: . numConstants - number of constant parameters
2609: . constants    - constant parameters
2610: - g0           - output values at the current point

2612:   Level: intermediate

2614:   Note:
2615:   `g1`, `g2`, and `g3` have identical calling sequences to `g0` and are omitted for brevity.

2617:   We are using a first order FEM model for the weak form\:

2619:   $$
2620:   \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
2621:   + \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
2622:   $$

2624: .seealso: `PetscDS`, `PetscDSSetBdJacobianPreconditioner()`
2625: @*/
2626: PetscErrorCode PetscDSGetBdJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g, void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]), void (**g1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (**g2)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (**g3)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))
2627: {
2628:   PetscBdPointJac *tmp0, *tmp1, *tmp2, *tmp3;
2629:   PetscInt         n0, n1, n2, n3;

2631:   PetscFunctionBegin;
2633:   PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
2634:   PetscCheck(!(g < 0) && !(g >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", g, ds->Nf);
2635:   PetscCall(PetscWeakFormGetBdJacobianPreconditioner(ds->wf, NULL, 0, f, g, 0, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3));
2636:   *g0 = tmp0 ? tmp0[0] : NULL;
2637:   *g1 = tmp1 ? tmp1[0] : NULL;
2638:   *g2 = tmp2 ? tmp2[0] : NULL;
2639:   *g3 = tmp3 ? tmp3[0] : NULL;
2640:   PetscFunctionReturn(PETSC_SUCCESS);
2641: }

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

2646:   Not Collective; No Fortran Support

2648:   Input Parameters:
2649: + ds - The `PetscDS`
2650: . f  - The test field number
2651: . g  - The field number
2652: . g0 - integrand for the test and basis function term
2653: . g1 - integrand for the test function and basis function gradient term
2654: . g2 - integrand for the test function gradient and basis function term
2655: - g3 - integrand for the test function gradient and basis function gradient term

2657:   Calling sequence of `g0':
2658: + dim          - the coordinate dimension
2659: . Nf           - the number of fields
2660: . NfAux        - the number of auxiliary fields
2661: . uOff         - the offset into u[] and u_t[] for each field
2662: . uOff_x       - the offset into u_x[] for each field
2663: . u            - each field evaluated at the current point
2664: . u_t          - the time derivative of each field evaluated at the current point
2665: . u_x          - the gradient of each field evaluated at the current point
2666: . aOff         - the offset into a[] and a_t[] for each auxiliary field
2667: . aOff_x       - the offset into a_x[] for each auxiliary field
2668: . a            - each auxiliary field evaluated at the current point
2669: . a_t          - the time derivative of each auxiliary field evaluated at the current point
2670: . a_x          - the gradient of auxiliary each field evaluated at the current point
2671: . t            - current time
2672: . u_tShift     - the multiplier a for dF/dU_t
2673: . x            - coordinates of the current point
2674: . n            - normal at the current point
2675: . numConstants - number of constant parameters
2676: . constants    - constant parameters
2677: - g0           - output values at the current point

2679:   Level: intermediate

2681:   Note:
2682:   `g1`, `g2`, and `g3` have identical calling sequences to `g0` and are omitted for brevity.

2684:   We are using a first order FEM model for the weak form\:

2686:   $$
2687:   \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
2688:   + \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
2689:   $$

2691: .seealso: `PetscDS`, `PetscDSGetBdJacobianPreconditioner()`
2692: @*/
2693: PetscErrorCode PetscDSSetBdJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g, void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]), void (*g1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (*g2)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (*g3)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))
2694: {
2695:   PetscFunctionBegin;
2701:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2702:   PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", g);
2703:   PetscCall(PetscWeakFormSetIndexBdJacobianPreconditioner(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3));
2704:   PetscFunctionReturn(PETSC_SUCCESS);
2705: }

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

2710:   Not Collective

2712:   Input Parameters:
2713: + prob - The PetscDS
2714: - f    - The test field number

2716:   Output Parameters:
2717: + sol - exact solution for the test field
2718: - ctx - exact solution context

2720:   Calling sequence of `exactSol`:
2721: + dim - the coordinate dimension
2722: . t   - current time
2723: . x   - coordinates of the current point
2724: . Nc  - the number of field components
2725: . u   - the solution field evaluated at the current point
2726: - ctx - a user context

2728:   Level: intermediate

2730: .seealso: `PetscDS`, `PetscDSSetExactSolution()`, `PetscDSGetExactSolutionTimeDerivative()`
2731: @*/
2732: PetscErrorCode PetscDSGetExactSolution(PetscDS prob, PetscInt f, PetscErrorCode (**sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void **ctx)
2733: {
2734:   PetscFunctionBegin;
2736:   PetscCheck(!(f < 0) && !(f >= prob->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, prob->Nf);
2737:   if (sol) {
2738:     PetscAssertPointer(sol, 3);
2739:     *sol = prob->exactSol[f];
2740:   }
2741:   if (ctx) {
2742:     PetscAssertPointer(ctx, 4);
2743:     *ctx = prob->exactCtx[f];
2744:   }
2745:   PetscFunctionReturn(PETSC_SUCCESS);
2746: }

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

2751:   Not Collective

2753:   Input Parameters:
2754: + prob - The `PetscDS`
2755: . f    - The test field number
2756: . sol  - solution function for the test fields
2757: - ctx  - solution context or `NULL`

2759:   Calling sequence of `sol`:
2760: + dim - the coordinate dimension
2761: . t   - current time
2762: . x   - coordinates of the current point
2763: . Nc  - the number of field components
2764: . u   - the solution field evaluated at the current point
2765: - ctx - a user context

2767:   Level: intermediate

2769: .seealso: `PetscDS`, `PetscDSGetExactSolution()`
2770: @*/
2771: PetscErrorCode PetscDSSetExactSolution(PetscDS prob, PetscInt f, PetscErrorCode (*sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void *ctx)
2772: {
2773:   PetscFunctionBegin;
2775:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2776:   PetscCall(PetscDSEnlarge_Static(prob, f + 1));
2777:   if (sol) {
2779:     prob->exactSol[f] = sol;
2780:   }
2781:   if (ctx) {
2783:     prob->exactCtx[f] = ctx;
2784:   }
2785:   PetscFunctionReturn(PETSC_SUCCESS);
2786: }

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

2791:   Not Collective

2793:   Input Parameters:
2794: + prob - The `PetscDS`
2795: - f    - The test field number

2797:   Output Parameters:
2798: + sol - time derivative of the exact solution for the test field
2799: - ctx - time derivative of the exact solution context

2801:   Calling sequence of `exactSol`:
2802: + dim - the coordinate dimension
2803: . t   - current time
2804: . x   - coordinates of the current point
2805: . Nc  - the number of field components
2806: . u   - the solution field evaluated at the current point
2807: - ctx - a user context

2809:   Level: intermediate

2811: .seealso: `PetscDS`, `PetscDSSetExactSolutionTimeDerivative()`, `PetscDSGetExactSolution()`
2812: @*/
2813: PetscErrorCode PetscDSGetExactSolutionTimeDerivative(PetscDS prob, PetscInt f, PetscErrorCode (**sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void **ctx)
2814: {
2815:   PetscFunctionBegin;
2817:   PetscCheck(!(f < 0) && !(f >= prob->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, prob->Nf);
2818:   if (sol) {
2819:     PetscAssertPointer(sol, 3);
2820:     *sol = prob->exactSol_t[f];
2821:   }
2822:   if (ctx) {
2823:     PetscAssertPointer(ctx, 4);
2824:     *ctx = prob->exactCtx_t[f];
2825:   }
2826:   PetscFunctionReturn(PETSC_SUCCESS);
2827: }

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

2832:   Not Collective

2834:   Input Parameters:
2835: + prob - The `PetscDS`
2836: . f    - The test field number
2837: . sol  - time derivative of the solution function for the test fields
2838: - ctx  - time derivative of the solution context or `NULL`

2840:   Calling sequence of `sol`:
2841: + dim - the coordinate dimension
2842: . t   - current time
2843: . x   - coordinates of the current point
2844: . Nc  - the number of field components
2845: . u   - the solution field evaluated at the current point
2846: - ctx - a user context

2848:   Level: intermediate

2850: .seealso: `PetscDS`, `PetscDSGetExactSolutionTimeDerivative()`, `PetscDSSetExactSolution()`
2851: @*/
2852: PetscErrorCode PetscDSSetExactSolutionTimeDerivative(PetscDS prob, PetscInt f, PetscErrorCode (*sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void *ctx)
2853: {
2854:   PetscFunctionBegin;
2856:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2857:   PetscCall(PetscDSEnlarge_Static(prob, f + 1));
2858:   if (sol) {
2860:     prob->exactSol_t[f] = sol;
2861:   }
2862:   if (ctx) {
2864:     prob->exactCtx_t[f] = ctx;
2865:   }
2866:   PetscFunctionReturn(PETSC_SUCCESS);
2867: }

2869: /*@C
2870:   PetscDSGetLowerBound - Get the pointwise lower bound function for a given field

2872:   Not Collective

2874:   Input Parameters:
2875: + ds - The PetscDS
2876: - f  - The field number

2878:   Output Parameters:
2879: + lb  - lower bound for the field
2880: - ctx - lower bound context

2882:   Calling sequence of `lb`:
2883: + dim - the coordinate dimension
2884: . t   - current time
2885: . x   - coordinates of the current point
2886: . Nc  - the number of field components
2887: . u   - the lower bound evaluated at the current point
2888: - ctx - a user context

2890:   Level: intermediate

2892: .seealso: `PetscDS`, `PetscDSSetLowerBound()`, `PetscDSGetUpperBound()`, `PetscDSGetExactSolution()`
2893: @*/
2894: PetscErrorCode PetscDSGetLowerBound(PetscDS ds, PetscInt f, PetscErrorCode (**lb)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void **ctx)
2895: {
2896:   PetscFunctionBegin;
2898:   PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
2899:   if (lb) {
2900:     PetscAssertPointer(lb, 3);
2901:     *lb = ds->lowerBound[f];
2902:   }
2903:   if (ctx) {
2904:     PetscAssertPointer(ctx, 4);
2905:     *ctx = ds->lowerCtx[f];
2906:   }
2907:   PetscFunctionReturn(PETSC_SUCCESS);
2908: }

2910: /*@C
2911:   PetscDSSetLowerBound - Set the pointwise lower bound function for a given field

2913:   Not Collective

2915:   Input Parameters:
2916: + ds  - The `PetscDS`
2917: . f   - The field number
2918: . lb  - solution function for the test fields
2919: - ctx - solution context or `NULL`

2921:   Calling sequence of `lb`:
2922: + dim - the coordinate dimension
2923: . t   - current time
2924: . x   - coordinates of the current point
2925: . Nc  - the number of field components
2926: . u   - the lower bound evaluated at the current point
2927: - ctx - a user context

2929:   Level: intermediate

2931: .seealso: `PetscDS`, `PetscDSGetLowerBound()`, `PetscDSGetUpperBound()`, `PetscDSGetExactSolution()`
2932: @*/
2933: PetscErrorCode PetscDSSetLowerBound(PetscDS ds, PetscInt f, PetscErrorCode (*lb)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void *ctx)
2934: {
2935:   PetscFunctionBegin;
2937:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2938:   PetscCall(PetscDSEnlarge_Static(ds, f + 1));
2939:   if (lb) {
2941:     ds->lowerBound[f] = lb;
2942:   }
2943:   if (ctx) {
2945:     ds->lowerCtx[f] = ctx;
2946:   }
2947:   PetscFunctionReturn(PETSC_SUCCESS);
2948: }

2950: /*@C
2951:   PetscDSGetUpperBound - Get the pointwise upper bound function for a given field

2953:   Not Collective

2955:   Input Parameters:
2956: + ds - The PetscDS
2957: - f  - The field number

2959:   Output Parameters:
2960: + ub  - upper bound for the field
2961: - ctx - upper bound context

2963:   Calling sequence of `ub`:
2964: + dim - the coordinate dimension
2965: . t   - current time
2966: . x   - coordinates of the current point
2967: . Nc  - the number of field components
2968: . u   - the upper bound evaluated at the current point
2969: - ctx - a user context

2971:   Level: intermediate

2973: .seealso: `PetscDS`, `PetscDSSetUpperBound()`, `PetscDSGetLowerBound()`, `PetscDSGetExactSolution()`
2974: @*/
2975: PetscErrorCode PetscDSGetUpperBound(PetscDS ds, PetscInt f, PetscErrorCode (**ub)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void **ctx)
2976: {
2977:   PetscFunctionBegin;
2979:   PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
2980:   if (ub) {
2981:     PetscAssertPointer(ub, 3);
2982:     *ub = ds->upperBound[f];
2983:   }
2984:   if (ctx) {
2985:     PetscAssertPointer(ctx, 4);
2986:     *ctx = ds->upperCtx[f];
2987:   }
2988:   PetscFunctionReturn(PETSC_SUCCESS);
2989: }

2991: /*@C
2992:   PetscDSSetUpperBound - Set the pointwise upper bound function for a given field

2994:   Not Collective

2996:   Input Parameters:
2997: + ds  - The `PetscDS`
2998: . f   - The field number
2999: . ub  - solution function for the test fields
3000: - ctx - solution context or `NULL`

3002:   Calling sequence of `ub`:
3003: + dim - the coordinate dimension
3004: . t   - current time
3005: . x   - coordinates of the current point
3006: . Nc  - the number of field components
3007: . u   - the upper bound evaluated at the current point
3008: - ctx - a user context

3010:   Level: intermediate

3012: .seealso: `PetscDS`, `PetscDSGetUpperBound()`, `PetscDSGetLowerBound()`, `PetscDSGetExactSolution()`
3013: @*/
3014: PetscErrorCode PetscDSSetUpperBound(PetscDS ds, PetscInt f, PetscErrorCode (*ub)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void *ctx)
3015: {
3016:   PetscFunctionBegin;
3018:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
3019:   PetscCall(PetscDSEnlarge_Static(ds, f + 1));
3020:   if (ub) {
3022:     ds->upperBound[f] = ub;
3023:   }
3024:   if (ctx) {
3026:     ds->upperCtx[f] = ctx;
3027:   }
3028:   PetscFunctionReturn(PETSC_SUCCESS);
3029: }

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

3034:   Not Collective

3036:   Input Parameter:
3037: . ds - The `PetscDS` object

3039:   Output Parameters:
3040: + numConstants - The number of constants, or pass in `NULL` if not required
3041: - constants    - The array of constants, `NULL` if there are none

3043:   Level: intermediate

3045: .seealso: `PetscDS`, `PetscDSSetConstants()`, `PetscDSCreate()`
3046: @*/
3047: PetscErrorCode PetscDSGetConstants(PetscDS ds, PeOp PetscInt *numConstants, PeOp const PetscScalar *constants[])
3048: {
3049:   PetscFunctionBegin;
3051:   if (numConstants) {
3052:     PetscAssertPointer(numConstants, 2);
3053:     *numConstants = ds->numConstants;
3054:   }
3055:   if (constants) {
3056:     PetscAssertPointer(constants, 3);
3057:     *constants = ds->constants;
3058:   }
3059:   PetscFunctionReturn(PETSC_SUCCESS);
3060: }

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

3065:   Not Collective

3067:   Input Parameters:
3068: + ds           - The `PetscDS` object
3069: . numConstants - The number of constants
3070: - constants    - The array of constants, `NULL` if there are none

3072:   Level: intermediate

3074: .seealso: `PetscDS`, `PetscDSGetConstants()`, `PetscDSCreate()`
3075: @*/
3076: PetscErrorCode PetscDSSetConstants(PetscDS ds, PetscInt numConstants, PetscScalar constants[])
3077: {
3078:   PetscFunctionBegin;
3080:   if (numConstants != ds->numConstants) {
3081:     PetscCall(PetscFree(ds->constants));
3082:     ds->numConstants = numConstants;
3083:     PetscCall(PetscMalloc1(ds->numConstants + ds->numFuncConstants, &ds->constants));
3084:   }
3085:   if (ds->numConstants) {
3086:     PetscAssertPointer(constants, 3);
3087:     PetscCall(PetscArraycpy(ds->constants, constants, ds->numConstants));
3088:   }
3089:   PetscFunctionReturn(PETSC_SUCCESS);
3090: }

3092: /*@C
3093:   PetscDSSetIntegrationParameters - Set the parameters for a particular integration

3095:   Not Collective

3097:   Input Parameters:
3098: + ds     - The `PetscDS` object
3099: . fieldI - The test field for a given point function, or PETSC_DETERMINE
3100: - fieldJ - The basis field for a given point function, or PETSC_DETERMINE

3102:   Level: intermediate

3104: .seealso: `PetscDS`, `PetscDSSetConstants()`, `PetscDSGetConstants()`, `PetscDSCreate()`
3105: @*/
3106: PetscErrorCode PetscDSSetIntegrationParameters(PetscDS ds, PetscInt fieldI, PetscInt fieldJ)
3107: {
3108:   PetscFunctionBegin;
3110:   ds->constants[ds->numConstants]     = fieldI;
3111:   ds->constants[ds->numConstants + 1] = fieldJ;
3112:   PetscFunctionReturn(PETSC_SUCCESS);
3113: }

3115: /*@C
3116:   PetscDSSetCellParameters - Set the parameters for a particular cell

3118:   Not Collective

3120:   Input Parameters:
3121: + ds     - The `PetscDS` object
3122: - volume - The cell volume

3124:   Level: intermediate

3126: .seealso: `PetscDS`, `PetscDSSetConstants()`, `PetscDSGetConstants()`, `PetscDSCreate()`
3127: @*/
3128: PetscErrorCode PetscDSSetCellParameters(PetscDS ds, PetscReal volume)
3129: {
3130:   PetscFunctionBegin;
3132:   ds->constants[ds->numConstants + 2] = volume;
3133:   PetscFunctionReturn(PETSC_SUCCESS);
3134: }

3136: /*@
3137:   PetscDSGetFieldIndex - Returns the index of the given field

3139:   Not Collective

3141:   Input Parameters:
3142: + prob - The `PetscDS` object
3143: - disc - The discretization object

3145:   Output Parameter:
3146: . f - The field number

3148:   Level: beginner

3150: .seealso: `PetscDS`, `PetscGetDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3151: @*/
3152: PetscErrorCode PetscDSGetFieldIndex(PetscDS prob, PetscObject disc, PetscInt *f)
3153: {
3154:   PetscInt g;

3156:   PetscFunctionBegin;
3158:   PetscAssertPointer(f, 3);
3159:   *f = -1;
3160:   for (g = 0; g < prob->Nf; ++g) {
3161:     if (disc == prob->disc[g]) break;
3162:   }
3163:   PetscCheck(g != prob->Nf, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONG, "Field not found in PetscDS.");
3164:   *f = g;
3165:   PetscFunctionReturn(PETSC_SUCCESS);
3166: }

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

3171:   Not Collective

3173:   Input Parameters:
3174: + prob - The `PetscDS` object
3175: - f    - The field number

3177:   Output Parameter:
3178: . size - The size

3180:   Level: beginner

3182: .seealso: `PetscDS`, `PetscDSGetFieldOffset()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3183: @*/
3184: PetscErrorCode PetscDSGetFieldSize(PetscDS prob, PetscInt f, PetscInt *size)
3185: {
3186:   PetscFunctionBegin;
3188:   PetscAssertPointer(size, 3);
3189:   PetscCheck(!(f < 0) && !(f >= prob->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, prob->Nf);
3190:   PetscCall(PetscDSSetUp(prob));
3191:   *size = prob->Nb[f];
3192:   PetscFunctionReturn(PETSC_SUCCESS);
3193: }

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

3198:   Not Collective

3200:   Input Parameters:
3201: + prob - The `PetscDS` object
3202: - f    - The field number

3204:   Output Parameter:
3205: . off - The offset

3207:   Level: beginner

3209: .seealso: `PetscDS`, `PetscDSGetFieldSize()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3210: @*/
3211: PetscErrorCode PetscDSGetFieldOffset(PetscDS prob, PetscInt f, PetscInt *off)
3212: {
3213:   PetscInt size, g;

3215:   PetscFunctionBegin;
3217:   PetscAssertPointer(off, 3);
3218:   PetscCheck(!(f < 0) && !(f >= prob->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, prob->Nf);
3219:   *off = 0;
3220:   for (g = 0; g < f; ++g) {
3221:     PetscCall(PetscDSGetFieldSize(prob, g, &size));
3222:     *off += size;
3223:   }
3224:   PetscFunctionReturn(PETSC_SUCCESS);
3225: }

3227: /*@
3228:   PetscDSGetFieldOffsetCohesive - Returns the offset of the given field in the full space basis on a cohesive cell

3230:   Not Collective

3232:   Input Parameters:
3233: + ds - The `PetscDS` object
3234: - f  - The field number

3236:   Output Parameter:
3237: . off - The offset

3239:   Level: beginner

3241: .seealso: `PetscDS`, `PetscDSGetFieldSize()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3242: @*/
3243: PetscErrorCode PetscDSGetFieldOffsetCohesive(PetscDS ds, PetscInt f, PetscInt *off)
3244: {
3245:   PetscInt size, g;

3247:   PetscFunctionBegin;
3249:   PetscAssertPointer(off, 3);
3250:   PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
3251:   *off = 0;
3252:   for (g = 0; g < f; ++g) {
3253:     PetscBool cohesive;

3255:     PetscCall(PetscDSGetCohesive(ds, g, &cohesive));
3256:     PetscCall(PetscDSGetFieldSize(ds, g, &size));
3257:     *off += cohesive ? size : size * 2;
3258:   }
3259:   PetscFunctionReturn(PETSC_SUCCESS);
3260: }

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

3265:   Not Collective

3267:   Input Parameter:
3268: . prob - The `PetscDS` object

3270:   Output Parameter:
3271: . dimensions - The number of dimensions

3273:   Level: beginner

3275: .seealso: `PetscDS`, `PetscDSGetComponentOffsets()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3276: @*/
3277: PetscErrorCode PetscDSGetDimensions(PetscDS prob, PetscInt *dimensions[])
3278: {
3279:   PetscFunctionBegin;
3281:   PetscCall(PetscDSSetUp(prob));
3282:   PetscAssertPointer(dimensions, 2);
3283:   *dimensions = prob->Nb;
3284:   PetscFunctionReturn(PETSC_SUCCESS);
3285: }

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

3290:   Not Collective

3292:   Input Parameter:
3293: . prob - The `PetscDS` object

3295:   Output Parameter:
3296: . components - The number of components

3298:   Level: beginner

3300: .seealso: `PetscDS`, `PetscDSGetComponentOffsets()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3301: @*/
3302: PetscErrorCode PetscDSGetComponents(PetscDS prob, PetscInt *components[])
3303: {
3304:   PetscFunctionBegin;
3306:   PetscCall(PetscDSSetUp(prob));
3307:   PetscAssertPointer(components, 2);
3308:   *components = prob->Nc;
3309:   PetscFunctionReturn(PETSC_SUCCESS);
3310: }

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

3315:   Not Collective

3317:   Input Parameters:
3318: + prob - The `PetscDS` object
3319: - f    - The field number

3321:   Output Parameter:
3322: . off - The offset

3324:   Level: beginner

3326: .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3327: @*/
3328: PetscErrorCode PetscDSGetComponentOffset(PetscDS prob, PetscInt f, PetscInt *off)
3329: {
3330:   PetscFunctionBegin;
3332:   PetscAssertPointer(off, 3);
3333:   PetscCheck(!(f < 0) && !(f >= prob->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, prob->Nf);
3334:   PetscCall(PetscDSSetUp(prob));
3335:   *off = prob->off[f];
3336:   PetscFunctionReturn(PETSC_SUCCESS);
3337: }

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

3342:   Not Collective

3344:   Input Parameter:
3345: . prob - The `PetscDS` object

3347:   Output Parameter:
3348: . offsets - The offsets

3350:   Level: beginner

3352: .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3353: @*/
3354: PetscErrorCode PetscDSGetComponentOffsets(PetscDS prob, PetscInt *offsets[])
3355: {
3356:   PetscFunctionBegin;
3358:   PetscAssertPointer(offsets, 2);
3359:   PetscCall(PetscDSSetUp(prob));
3360:   *offsets = prob->off;
3361:   PetscFunctionReturn(PETSC_SUCCESS);
3362: }

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

3367:   Not Collective

3369:   Input Parameter:
3370: . prob - The `PetscDS` object

3372:   Output Parameter:
3373: . offsets - The offsets

3375:   Level: beginner

3377: .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3378: @*/
3379: PetscErrorCode PetscDSGetComponentDerivativeOffsets(PetscDS prob, PetscInt *offsets[])
3380: {
3381:   PetscFunctionBegin;
3383:   PetscAssertPointer(offsets, 2);
3384:   PetscCall(PetscDSSetUp(prob));
3385:   *offsets = prob->offDer;
3386:   PetscFunctionReturn(PETSC_SUCCESS);
3387: }

3389: /*@
3390:   PetscDSGetComponentOffsetsCohesive - Returns the offset of each field on an evaluation point

3392:   Not Collective

3394:   Input Parameters:
3395: + ds - The `PetscDS` object
3396: - s  - The cohesive side, 0 for negative, 1 for positive, 2 for cohesive

3398:   Output Parameter:
3399: . offsets - The offsets

3401:   Level: beginner

3403: .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3404: @*/
3405: PetscErrorCode PetscDSGetComponentOffsetsCohesive(PetscDS ds, PetscInt s, PetscInt *offsets[])
3406: {
3407:   PetscFunctionBegin;
3409:   PetscAssertPointer(offsets, 3);
3410:   PetscCheck(ds->isCohesive, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cohesive offsets are only valid for a cohesive DS");
3411:   PetscCheck(!(s < 0) && !(s > 2), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cohesive side %" PetscInt_FMT " is not in [0, 2]", s);
3412:   PetscCall(PetscDSSetUp(ds));
3413:   *offsets = ds->offCohesive[s];
3414:   PetscFunctionReturn(PETSC_SUCCESS);
3415: }

3417: /*@
3418:   PetscDSGetComponentDerivativeOffsetsCohesive - Returns the offset of each field derivative on an evaluation point

3420:   Not Collective

3422:   Input Parameters:
3423: + ds - The `PetscDS` object
3424: - s  - The cohesive side, 0 for negative, 1 for positive, 2 for cohesive

3426:   Output Parameter:
3427: . offsets - The offsets

3429:   Level: beginner

3431: .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3432: @*/
3433: PetscErrorCode PetscDSGetComponentDerivativeOffsetsCohesive(PetscDS ds, PetscInt s, PetscInt *offsets[])
3434: {
3435:   PetscFunctionBegin;
3437:   PetscAssertPointer(offsets, 3);
3438:   PetscCheck(ds->isCohesive, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cohesive offsets are only valid for a cohesive DS");
3439:   PetscCheck(!(s < 0) && !(s > 2), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cohesive side %" PetscInt_FMT " is not in [0, 2]", s);
3440:   PetscCall(PetscDSSetUp(ds));
3441:   *offsets = ds->offDerCohesive[s];
3442:   PetscFunctionReturn(PETSC_SUCCESS);
3443: }

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

3448:   Not Collective

3450:   Input Parameter:
3451: . prob - The `PetscDS` object

3453:   Output Parameter:
3454: . T - The basis function and derivatives tabulation at quadrature points for each field, see `PetscTabulation` for its details

3456:   Level: intermediate

3458:   Note:
3459:   The tabulation is only valid so long as the `PetscDS` has not be destroyed. There is no `PetscDSRestoreTabulation()` in C.

3461:   Fortran Note:
3462:   Use the declaration
3463: .vb
3464:   PetscTabulation, pointer :: tab(:)
3465: .ve
3466:   and access the values using, for example,
3467: .vb
3468:   tab(i)%ptr%K
3469:   tab(i)%ptr%T(j)%ptr
3470: .ve
3471:   where $ i = 1, 2, ..., Nf $ and $ j = 1, 2, ..., tab(i)%ptr%K+1 $.

3473:   Use `PetscDSRestoreTabulation()` to restore the array

3475:   Developer Note:
3476:   The Fortran language syntax does not directly support arrays of pointers, the '%ptr' notation allows mimicking their use in Fortran.

3478: .seealso: `PetscDS`, `PetscTabulation`, `PetscDSCreate()`
3479: @*/
3480: PetscErrorCode PetscDSGetTabulation(PetscDS prob, PetscTabulation *T[]) PeNS
3481: {
3482:   PetscFunctionBegin;
3484:   PetscAssertPointer(T, 2);
3485:   PetscCall(PetscDSSetUp(prob));
3486:   *T = prob->T;
3487:   PetscFunctionReturn(PETSC_SUCCESS);
3488: }

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

3493:   Not Collective

3495:   Input Parameter:
3496: . prob - The `PetscDS` object

3498:   Output Parameter:
3499: . Tf - The basis function and derivative tabulation on each local face at quadrature points for each field

3501:   Level: intermediate

3503:   Note:
3504:   The tabulation is only valid so long as the `PetscDS` has not be destroyed. There is no `PetscDSRestoreFaceTabulation()` in C.

3506: .seealso: `PetscTabulation`, `PetscDS`, `PetscDSGetTabulation()`, `PetscDSCreate()`
3507: @*/
3508: PetscErrorCode PetscDSGetFaceTabulation(PetscDS prob, PetscTabulation *Tf[])
3509: {
3510:   PetscFunctionBegin;
3512:   PetscAssertPointer(Tf, 2);
3513:   PetscCall(PetscDSSetUp(prob));
3514:   *Tf = prob->Tf;
3515:   PetscFunctionReturn(PETSC_SUCCESS);
3516: }

3518: PetscErrorCode PetscDSGetEvaluationArrays(PetscDS prob, PetscScalar *u[], PetscScalar *u_t[], PetscScalar *u_x[])
3519: {
3520:   PetscFunctionBegin;
3522:   PetscCall(PetscDSSetUp(prob));
3523:   if (u) {
3524:     PetscAssertPointer(u, 2);
3525:     *u = prob->u;
3526:   }
3527:   if (u_t) {
3528:     PetscAssertPointer(u_t, 3);
3529:     *u_t = prob->u_t;
3530:   }
3531:   if (u_x) {
3532:     PetscAssertPointer(u_x, 4);
3533:     *u_x = prob->u_x;
3534:   }
3535:   PetscFunctionReturn(PETSC_SUCCESS);
3536: }

3538: PetscErrorCode PetscDSGetWeakFormArrays(PetscDS prob, PetscScalar *f0[], PetscScalar *f1[], PetscScalar *g0[], PetscScalar *g1[], PetscScalar *g2[], PetscScalar *g3[])
3539: {
3540:   PetscFunctionBegin;
3542:   PetscCall(PetscDSSetUp(prob));
3543:   if (f0) {
3544:     PetscAssertPointer(f0, 2);
3545:     *f0 = prob->f0;
3546:   }
3547:   if (f1) {
3548:     PetscAssertPointer(f1, 3);
3549:     *f1 = prob->f1;
3550:   }
3551:   if (g0) {
3552:     PetscAssertPointer(g0, 4);
3553:     *g0 = prob->g0;
3554:   }
3555:   if (g1) {
3556:     PetscAssertPointer(g1, 5);
3557:     *g1 = prob->g1;
3558:   }
3559:   if (g2) {
3560:     PetscAssertPointer(g2, 6);
3561:     *g2 = prob->g2;
3562:   }
3563:   if (g3) {
3564:     PetscAssertPointer(g3, 7);
3565:     *g3 = prob->g3;
3566:   }
3567:   PetscFunctionReturn(PETSC_SUCCESS);
3568: }

3570: PetscErrorCode PetscDSGetWorkspace(PetscDS prob, PetscReal **x, PetscScalar **basisReal, PetscScalar **basisDerReal, PetscScalar **testReal, PetscScalar **testDerReal)
3571: {
3572:   PetscFunctionBegin;
3574:   PetscCall(PetscDSSetUp(prob));
3575:   if (x) {
3576:     PetscAssertPointer(x, 2);
3577:     *x = prob->x;
3578:   }
3579:   if (basisReal) {
3580:     PetscAssertPointer(basisReal, 3);
3581:     *basisReal = prob->basisReal;
3582:   }
3583:   if (basisDerReal) {
3584:     PetscAssertPointer(basisDerReal, 4);
3585:     *basisDerReal = prob->basisDerReal;
3586:   }
3587:   if (testReal) {
3588:     PetscAssertPointer(testReal, 5);
3589:     *testReal = prob->testReal;
3590:   }
3591:   if (testDerReal) {
3592:     PetscAssertPointer(testDerReal, 6);
3593:     *testDerReal = prob->testDerReal;
3594:   }
3595:   PetscFunctionReturn(PETSC_SUCCESS);
3596: }

3598: /*@C
3599:   PetscDSAddBoundary - Add a boundary condition to the model.

3601:   Collective

3603:   Input Parameters:
3604: + ds       - The PetscDS object
3605: . type     - The type of condition, e.g. `DM_BC_ESSENTIAL`/`DM_BC_ESSENTIAL_FIELD` (Dirichlet), or `DM_BC_NATURAL` (Neumann)
3606: . name     - The BC name
3607: . label    - The label defining constrained points
3608: . Nv       - The number of `DMLabel` values for constrained points
3609: . values   - An array of label values for constrained points
3610: . field    - The field to constrain
3611: . Nc       - The number of constrained field components (0 will constrain all fields)
3612: . comps    - An array of constrained component numbers
3613: . bcFunc   - A pointwise function giving boundary values
3614: . bcFunc_t - A pointwise function giving the time derivative of the boundary values, or NULL
3615: - ctx      - An optional user context for bcFunc

3617:   Output Parameter:
3618: . bd - The boundary number

3620:   Options Database Keys:
3621: + -bc_<boundary name> <num>      - Overrides the boundary ids
3622: - -bc_<boundary name>_comp <num> - Overrides the boundary components

3624:   Level: developer

3626:   Note:
3627:   Both `bcFunc` and `bcFunc_t` will depend on the boundary condition type. If the type if `DM_BC_ESSENTIAL`, then the calling sequence is\:

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

3631:   If the type is `DM_BC_ESSENTIAL_FIELD` or other _FIELD value, then the calling sequence is\:
3632: .vb
3633:   void bcFunc(PetscInt dim, PetscInt Nf, PetscInt NfAux,
3634:               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
3635:               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
3636:               PetscReal time, const PetscReal x[], PetscScalar bcval[])
3637: .ve
3638: + dim - the coordinate dimension
3639: . Nf - the number of fields
3640: . uOff - the offset into u[] and u_t[] for each field
3641: . uOff_x - the offset into u_x[] for each field
3642: . u - each field evaluated at the current point
3643: . u_t - the time derivative of each field evaluated at the current point
3644: . u_x - the gradient of each field evaluated at the current point
3645: . aOff - the offset into a[] and a_t[] for each auxiliary field
3646: . aOff_x - the offset into a_x[] for each auxiliary field
3647: . a - each auxiliary field evaluated at the current point
3648: . a_t - the time derivative of each auxiliary field evaluated at the current point
3649: . a_x - the gradient of auxiliary each field evaluated at the current point
3650: . t - current time
3651: . x - coordinates of the current point
3652: . numConstants - number of constant parameters
3653: . constants - constant parameters
3654: - bcval - output values at the current point

3656:   Notes:
3657:   The pointwise functions are used to provide boundary values for essential boundary
3658:   conditions. In FEM, they are acting upon by dual basis functionals to generate FEM
3659:   coefficients which are fixed. Natural boundary conditions signal to PETSc that boundary
3660:   integrals should be performed, using the kernels from `PetscDSSetBdResidual()`.

3662: .seealso: `PetscDS`, `PetscWeakForm`, `DMLabel`, `DMBoundaryConditionType`, `PetscDSAddBoundaryByName()`, `PetscDSGetBoundary()`, `PetscDSSetResidual()`, `PetscDSSetBdResidual()`
3663: @*/
3664: PetscErrorCode PetscDSAddBoundary(PetscDS ds, DMBoundaryConditionType type, const char name[], DMLabel label, PetscInt Nv, const PetscInt values[], PetscInt field, PetscInt Nc, const PetscInt comps[], void (*bcFunc)(void), void (*bcFunc_t)(void), void *ctx, PetscInt *bd)
3665: {
3666:   DSBoundary  head = ds->boundary, b;
3667:   PetscInt    n    = 0;
3668:   const char *lname;

3670:   PetscFunctionBegin;
3673:   PetscAssertPointer(name, 3);
3678:   PetscCheck(field >= 0 && field < ds->Nf, PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_OUTOFRANGE, "Field %" PetscInt_FMT " is not in [0, %" PetscInt_FMT ")", field, ds->Nf);
3679:   if (Nc > 0) {
3680:     PetscInt *fcomps;
3681:     PetscInt  c;

3683:     PetscCall(PetscDSGetComponents(ds, &fcomps));
3684:     PetscCheck(Nc <= fcomps[field], PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_OUTOFRANGE, "Number of constrained components %" PetscInt_FMT " > %" PetscInt_FMT " components for field %" PetscInt_FMT, Nc, fcomps[field], field);
3685:     for (c = 0; c < Nc; ++c) {
3686:       PetscCheck(comps[c] >= 0 && comps[c] < fcomps[field], PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_OUTOFRANGE, "Constrained component[%" PetscInt_FMT "] %" PetscInt_FMT " not in [0, %" PetscInt_FMT ") components for field %" PetscInt_FMT, c, comps[c], fcomps[field], field);
3687:     }
3688:   }
3689:   PetscCall(PetscNew(&b));
3690:   PetscCall(PetscStrallocpy(name, (char **)&b->name));
3691:   PetscCall(PetscWeakFormCreate(PETSC_COMM_SELF, &b->wf));
3692:   PetscCall(PetscWeakFormSetNumFields(b->wf, ds->Nf));
3693:   PetscCall(PetscMalloc1(Nv, &b->values));
3694:   if (Nv) PetscCall(PetscArraycpy(b->values, values, Nv));
3695:   PetscCall(PetscMalloc1(Nc, &b->comps));
3696:   if (Nc) PetscCall(PetscArraycpy(b->comps, comps, Nc));
3697:   PetscCall(PetscObjectGetName((PetscObject)label, &lname));
3698:   PetscCall(PetscStrallocpy(lname, (char **)&b->lname));
3699:   b->type   = type;
3700:   b->label  = label;
3701:   b->Nv     = Nv;
3702:   b->field  = field;
3703:   b->Nc     = Nc;
3704:   b->func   = bcFunc;
3705:   b->func_t = bcFunc_t;
3706:   b->ctx    = ctx;
3707:   b->next   = NULL;
3708:   /* Append to linked list so that we can preserve the order */
3709:   if (!head) ds->boundary = b;
3710:   while (head) {
3711:     if (!head->next) {
3712:       head->next = b;
3713:       head       = b;
3714:     }
3715:     head = head->next;
3716:     ++n;
3717:   }
3718:   if (bd) {
3719:     PetscAssertPointer(bd, 13);
3720:     *bd = n;
3721:   }
3722:   PetscFunctionReturn(PETSC_SUCCESS);
3723: }

3725: // PetscClangLinter pragma ignore: -fdoc-section-header-unknown
3726: /*@C
3727:   PetscDSAddBoundaryByName - Add a boundary condition to the model.

3729:   Collective

3731:   Input Parameters:
3732: + ds       - The `PetscDS` object
3733: . type     - The type of condition, e.g. `DM_BC_ESSENTIAL`/`DM_BC_ESSENTIAL_FIELD` (Dirichlet), or `DM_BC_NATURAL` (Neumann)
3734: . name     - The BC name
3735: . lname    - The naem of the label defining constrained points
3736: . Nv       - The number of `DMLabel` values for constrained points
3737: . values   - An array of label values for constrained points
3738: . field    - The field to constrain
3739: . Nc       - The number of constrained field components (0 will constrain all fields)
3740: . comps    - An array of constrained component numbers
3741: . bcFunc   - A pointwise function giving boundary values
3742: . bcFunc_t - A pointwise function giving the time derivative of the boundary values, or NULL
3743: - ctx      - An optional user context for bcFunc

3745:   Output Parameter:
3746: . bd - The boundary number

3748:   Options Database Keys:
3749: + -bc_<boundary name> <num>      - Overrides the boundary ids
3750: - -bc_<boundary name>_comp <num> - Overrides the boundary components

3752:   Calling Sequence of `bcFunc` and `bcFunc_t`:
3753:   If the type is `DM_BC_ESSENTIAL`
3754: .vb
3755:   void bcFunc(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar bcval[])
3756: .ve
3757:   If the type is `DM_BC_ESSENTIAL_FIELD` or other _FIELD value,
3758: .vb
3759:   void bcFunc(PetscInt dim, PetscInt Nf, PetscInt NfAux,
3760:               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
3761:               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
3762:               PetscReal time, const PetscReal x[], PetscScalar bcval[])
3763: .ve
3764: + dim - the coordinate dimension
3765: . Nf - the number of fields
3766: . uOff - the offset into u[] and u_t[] for each field
3767: . uOff_x - the offset into u_x[] for each field
3768: . u - each field evaluated at the current point
3769: . u_t - the time derivative of each field evaluated at the current point
3770: . u_x - the gradient of each field evaluated at the current point
3771: . aOff - the offset into a[] and a_t[] for each auxiliary field
3772: . aOff_x - the offset into a_x[] for each auxiliary field
3773: . a - each auxiliary field evaluated at the current point
3774: . a_t - the time derivative of each auxiliary field evaluated at the current point
3775: . a_x - the gradient of auxiliary each field evaluated at the current point
3776: . t - current time
3777: . x - coordinates of the current point
3778: . numConstants - number of constant parameters
3779: . constants - constant parameters
3780: - bcval - output values at the current point

3782:   Level: developer

3784:   Notes:
3785:   The pointwise functions are used to provide boundary values for essential boundary
3786:   conditions. In FEM, they are acting upon by dual basis functionals to generate FEM
3787:   coefficients which are fixed. Natural boundary conditions signal to PETSc that boundary
3788:   integrals should be performed, using the kernels from `PetscDSSetBdResidual()`.

3790:   This function should only be used with `DMFOREST` currently, since labels cannot be defined before the underlying `DMPLEX` is built.

3792: .seealso: `PetscDS`, `PetscWeakForm`, `DMLabel`, `DMBoundaryConditionType`, `PetscDSAddBoundary()`, `PetscDSGetBoundary()`, `PetscDSSetResidual()`, `PetscDSSetBdResidual()`
3793: @*/
3794: PetscErrorCode PetscDSAddBoundaryByName(PetscDS ds, DMBoundaryConditionType type, const char name[], const char lname[], PetscInt Nv, const PetscInt values[], PetscInt field, PetscInt Nc, const PetscInt comps[], void (*bcFunc)(void), void (*bcFunc_t)(void), void *ctx, PetscInt *bd)
3795: {
3796:   DSBoundary head = ds->boundary, b;
3797:   PetscInt   n    = 0;

3799:   PetscFunctionBegin;
3802:   PetscAssertPointer(name, 3);
3803:   PetscAssertPointer(lname, 4);
3807:   PetscCall(PetscNew(&b));
3808:   PetscCall(PetscStrallocpy(name, (char **)&b->name));
3809:   PetscCall(PetscWeakFormCreate(PETSC_COMM_SELF, &b->wf));
3810:   PetscCall(PetscWeakFormSetNumFields(b->wf, ds->Nf));
3811:   PetscCall(PetscMalloc1(Nv, &b->values));
3812:   if (Nv) PetscCall(PetscArraycpy(b->values, values, Nv));
3813:   PetscCall(PetscMalloc1(Nc, &b->comps));
3814:   if (Nc) PetscCall(PetscArraycpy(b->comps, comps, Nc));
3815:   PetscCall(PetscStrallocpy(lname, (char **)&b->lname));
3816:   b->type   = type;
3817:   b->label  = NULL;
3818:   b->Nv     = Nv;
3819:   b->field  = field;
3820:   b->Nc     = Nc;
3821:   b->func   = bcFunc;
3822:   b->func_t = bcFunc_t;
3823:   b->ctx    = ctx;
3824:   b->next   = NULL;
3825:   /* Append to linked list so that we can preserve the order */
3826:   if (!head) ds->boundary = b;
3827:   while (head) {
3828:     if (!head->next) {
3829:       head->next = b;
3830:       head       = b;
3831:     }
3832:     head = head->next;
3833:     ++n;
3834:   }
3835:   if (bd) {
3836:     PetscAssertPointer(bd, 13);
3837:     *bd = n;
3838:   }
3839:   PetscFunctionReturn(PETSC_SUCCESS);
3840: }

3842: /*@C
3843:   PetscDSUpdateBoundary - Change a boundary condition for the model.

3845:   Input Parameters:
3846: + ds       - The `PetscDS` object
3847: . bd       - The boundary condition number
3848: . type     - The type of condition, e.g. `DM_BC_ESSENTIAL`/`DM_BC_ESSENTIAL_FIELD` (Dirichlet), or `DM_BC_NATURAL` (Neumann)
3849: . name     - The BC name
3850: . label    - The label defining constrained points
3851: . Nv       - The number of `DMLabel` ids for constrained points
3852: . values   - An array of ids for constrained points
3853: . field    - The field to constrain
3854: . Nc       - The number of constrained field components
3855: . comps    - An array of constrained component numbers
3856: . bcFunc   - A pointwise function giving boundary values
3857: . bcFunc_t - A pointwise function giving the time derivative of the boundary values, or NULL
3858: - ctx      - An optional user context for bcFunc

3860:   Level: developer

3862:   Notes:
3863:   The pointwise functions are used to provide boundary values for essential boundary
3864:   conditions. In FEM, they are acting upon by dual basis functionals to generate FEM
3865:   coefficients which are fixed. Natural boundary conditions signal to PETSc that boundary
3866:   integrals should be performed, using the kernels from `PetscDSSetBdResidual()`.

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

3871: .seealso: `PetscDS`, `PetscWeakForm`, `DMBoundaryConditionType`, `PetscDSAddBoundary()`, `PetscDSGetBoundary()`, `PetscDSGetNumBoundary()`, `DMLabel`
3872: @*/
3873: PetscErrorCode PetscDSUpdateBoundary(PetscDS ds, PetscInt bd, DMBoundaryConditionType type, const char name[], DMLabel label, PetscInt Nv, const PetscInt values[], PetscInt field, PetscInt Nc, const PetscInt comps[], void (*bcFunc)(void), void (*bcFunc_t)(void), void *ctx)
3874: {
3875:   DSBoundary b = ds->boundary;
3876:   PetscInt   n = 0;

3878:   PetscFunctionBegin;
3880:   while (b) {
3881:     if (n == bd) break;
3882:     b = b->next;
3883:     ++n;
3884:   }
3885:   PetscCheck(b, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %" PetscInt_FMT " is not in [0, %" PetscInt_FMT ")", bd, n);
3886:   if (name) {
3887:     PetscCall(PetscFree(b->name));
3888:     PetscCall(PetscStrallocpy(name, (char **)&b->name));
3889:   }
3890:   b->type = type;
3891:   if (label) {
3892:     const char *name;

3894:     b->label = label;
3895:     PetscCall(PetscFree(b->lname));
3896:     PetscCall(PetscObjectGetName((PetscObject)label, &name));
3897:     PetscCall(PetscStrallocpy(name, (char **)&b->lname));
3898:   }
3899:   if (Nv >= 0) {
3900:     b->Nv = Nv;
3901:     PetscCall(PetscFree(b->values));
3902:     PetscCall(PetscMalloc1(Nv, &b->values));
3903:     if (Nv) PetscCall(PetscArraycpy(b->values, values, Nv));
3904:   }
3905:   if (field >= 0) b->field = field;
3906:   if (Nc >= 0) {
3907:     b->Nc = Nc;
3908:     PetscCall(PetscFree(b->comps));
3909:     PetscCall(PetscMalloc1(Nc, &b->comps));
3910:     if (Nc) PetscCall(PetscArraycpy(b->comps, comps, Nc));
3911:   }
3912:   if (bcFunc) b->func = bcFunc;
3913:   if (bcFunc_t) b->func_t = bcFunc_t;
3914:   if (ctx) b->ctx = ctx;
3915:   PetscFunctionReturn(PETSC_SUCCESS);
3916: }

3918: /*@
3919:   PetscDSGetNumBoundary - Get the number of registered BC

3921:   Input Parameter:
3922: . ds - The `PetscDS` object

3924:   Output Parameter:
3925: . numBd - The number of BC

3927:   Level: intermediate

3929: .seealso: `PetscDS`, `PetscDSAddBoundary()`, `PetscDSGetBoundary()`
3930: @*/
3931: PetscErrorCode PetscDSGetNumBoundary(PetscDS ds, PetscInt *numBd)
3932: {
3933:   DSBoundary b = ds->boundary;

3935:   PetscFunctionBegin;
3937:   PetscAssertPointer(numBd, 2);
3938:   *numBd = 0;
3939:   while (b) {
3940:     ++(*numBd);
3941:     b = b->next;
3942:   }
3943:   PetscFunctionReturn(PETSC_SUCCESS);
3944: }

3946: /*@C
3947:   PetscDSGetBoundary - Gets a boundary condition to the model

3949:   Input Parameters:
3950: + ds - The `PetscDS` object
3951: - bd - The BC number

3953:   Output Parameters:
3954: + wf     - The `PetscWeakForm` holding the pointwise functions
3955: . type   - The type of condition, e.g. `DM_BC_ESSENTIAL`/`DM_BC_ESSENTIAL_FIELD` (Dirichlet), or `DM_BC_NATURAL` (Neumann)
3956: . name   - The BC name
3957: . label  - The label defining constrained points
3958: . Nv     - The number of `DMLabel` ids for constrained points
3959: . values - An array of ids for constrained points
3960: . field  - The field to constrain
3961: . Nc     - The number of constrained field components
3962: . comps  - An array of constrained component numbers
3963: . func   - A pointwise function giving boundary values
3964: . func_t - A pointwise function giving the time derivative of the boundary values
3965: - ctx    - An optional user context for bcFunc

3967:   Options Database Keys:
3968: + -bc_<boundary name> <num>      - Overrides the boundary ids
3969: - -bc_<boundary name>_comp <num> - Overrides the boundary components

3971:   Level: developer

3973: .seealso: `PetscDS`, `PetscWeakForm`, `DMBoundaryConditionType`, `PetscDSAddBoundary()`, `DMLabel`
3974: @*/
3975: PetscErrorCode PetscDSGetBoundary(PetscDS ds, PetscInt bd, PetscWeakForm *wf, DMBoundaryConditionType *type, const char *name[], DMLabel *label, PetscInt *Nv, const PetscInt *values[], PetscInt *field, PetscInt *Nc, const PetscInt *comps[], void (**func)(void), void (**func_t)(void), void **ctx)
3976: {
3977:   DSBoundary b = ds->boundary;
3978:   PetscInt   n = 0;

3980:   PetscFunctionBegin;
3982:   while (b) {
3983:     if (n == bd) break;
3984:     b = b->next;
3985:     ++n;
3986:   }
3987:   PetscCheck(b, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %" PetscInt_FMT " is not in [0, %" PetscInt_FMT ")", bd, n);
3988:   if (wf) {
3989:     PetscAssertPointer(wf, 3);
3990:     *wf = b->wf;
3991:   }
3992:   if (type) {
3993:     PetscAssertPointer(type, 4);
3994:     *type = b->type;
3995:   }
3996:   if (name) {
3997:     PetscAssertPointer(name, 5);
3998:     *name = b->name;
3999:   }
4000:   if (label) {
4001:     PetscAssertPointer(label, 6);
4002:     *label = b->label;
4003:   }
4004:   if (Nv) {
4005:     PetscAssertPointer(Nv, 7);
4006:     *Nv = b->Nv;
4007:   }
4008:   if (values) {
4009:     PetscAssertPointer(values, 8);
4010:     *values = b->values;
4011:   }
4012:   if (field) {
4013:     PetscAssertPointer(field, 9);
4014:     *field = b->field;
4015:   }
4016:   if (Nc) {
4017:     PetscAssertPointer(Nc, 10);
4018:     *Nc = b->Nc;
4019:   }
4020:   if (comps) {
4021:     PetscAssertPointer(comps, 11);
4022:     *comps = b->comps;
4023:   }
4024:   if (func) {
4025:     PetscAssertPointer(func, 12);
4026:     *func = b->func;
4027:   }
4028:   if (func_t) {
4029:     PetscAssertPointer(func_t, 13);
4030:     *func_t = b->func_t;
4031:   }
4032:   if (ctx) {
4033:     PetscAssertPointer(ctx, 14);
4034:     *ctx = b->ctx;
4035:   }
4036:   PetscFunctionReturn(PETSC_SUCCESS);
4037: }

4039: /*@
4040:   PetscDSUpdateBoundaryLabels - Update `DMLabel` in each boundary condition using the label name and the input `DM`

4042:   Not Collective

4044:   Input Parameters:
4045: + ds - The source `PetscDS` object
4046: - dm - The `DM` holding labels

4048:   Level: intermediate

4050: .seealso: `PetscDS`, `DMBoundary`, `DM`, `PetscDSCopyBoundary()`, `PetscDSCreate()`, `DMGetLabel()`
4051: @*/
4052: PetscErrorCode PetscDSUpdateBoundaryLabels(PetscDS ds, DM dm)
4053: {
4054:   DSBoundary b;

4056:   PetscFunctionBegin;
4059:   for (b = ds->boundary; b; b = b->next) {
4060:     if (b->lname) PetscCall(DMGetLabel(dm, b->lname, &b->label));
4061:   }
4062:   PetscFunctionReturn(PETSC_SUCCESS);
4063: }

4065: static PetscErrorCode DSBoundaryDuplicate_Internal(DSBoundary b, DSBoundary *bNew)
4066: {
4067:   PetscFunctionBegin;
4068:   PetscCall(PetscNew(bNew));
4069:   PetscCall(PetscWeakFormCreate(PETSC_COMM_SELF, &(*bNew)->wf));
4070:   PetscCall(PetscWeakFormCopy(b->wf, (*bNew)->wf));
4071:   PetscCall(PetscStrallocpy(b->name, (char **)&((*bNew)->name)));
4072:   PetscCall(PetscStrallocpy(b->lname, (char **)&((*bNew)->lname)));
4073:   (*bNew)->type  = b->type;
4074:   (*bNew)->label = b->label;
4075:   (*bNew)->Nv    = b->Nv;
4076:   PetscCall(PetscMalloc1(b->Nv, &(*bNew)->values));
4077:   PetscCall(PetscArraycpy((*bNew)->values, b->values, b->Nv));
4078:   (*bNew)->field = b->field;
4079:   (*bNew)->Nc    = b->Nc;
4080:   PetscCall(PetscMalloc1(b->Nc, &(*bNew)->comps));
4081:   PetscCall(PetscArraycpy((*bNew)->comps, b->comps, b->Nc));
4082:   (*bNew)->func   = b->func;
4083:   (*bNew)->func_t = b->func_t;
4084:   (*bNew)->ctx    = b->ctx;
4085:   PetscFunctionReturn(PETSC_SUCCESS);
4086: }

4088: /*@
4089:   PetscDSCopyBoundary - Copy all boundary condition objects to the new problem

4091:   Not Collective

4093:   Input Parameters:
4094: + ds        - The source `PetscDS` object
4095: . numFields - The number of selected fields, or `PETSC_DEFAULT` for all fields
4096: - fields    - The selected fields, or NULL for all fields

4098:   Output Parameter:
4099: . newds - The target `PetscDS`, now with a copy of the boundary conditions

4101:   Level: intermediate

4103: .seealso: `PetscDS`, `DMBoundary`, `PetscDSCopyEquations()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
4104: @*/
4105: PetscErrorCode PetscDSCopyBoundary(PetscDS ds, PetscInt numFields, const PetscInt fields[], PetscDS newds)
4106: {
4107:   DSBoundary b, *lastnext;

4109:   PetscFunctionBegin;
4112:   if (ds == newds) PetscFunctionReturn(PETSC_SUCCESS);
4113:   PetscCall(PetscDSDestroyBoundary(newds));
4114:   lastnext = &newds->boundary;
4115:   for (b = ds->boundary; b; b = b->next) {
4116:     DSBoundary bNew;
4117:     PetscInt   fieldNew = -1;

4119:     if (numFields > 0 && fields) {
4120:       PetscInt f;

4122:       for (f = 0; f < numFields; ++f)
4123:         if (b->field == fields[f]) break;
4124:       if (f == numFields) continue;
4125:       fieldNew = f;
4126:     }
4127:     PetscCall(DSBoundaryDuplicate_Internal(b, &bNew));
4128:     bNew->field = fieldNew < 0 ? b->field : fieldNew;
4129:     *lastnext   = bNew;
4130:     lastnext    = &bNew->next;
4131:   }
4132:   PetscFunctionReturn(PETSC_SUCCESS);
4133: }

4135: /*@
4136:   PetscDSDestroyBoundary - Remove all `DMBoundary` objects from the `PetscDS`

4138:   Not Collective

4140:   Input Parameter:
4141: . ds - The `PetscDS` object

4143:   Level: intermediate

4145: .seealso: `PetscDS`, `DMBoundary`, `PetscDSCopyBoundary()`, `PetscDSCopyEquations()`
4146: @*/
4147: PetscErrorCode PetscDSDestroyBoundary(PetscDS ds)
4148: {
4149:   DSBoundary next = ds->boundary;

4151:   PetscFunctionBegin;
4152:   while (next) {
4153:     DSBoundary b = next;

4155:     next = b->next;
4156:     PetscCall(PetscWeakFormDestroy(&b->wf));
4157:     PetscCall(PetscFree(b->name));
4158:     PetscCall(PetscFree(b->lname));
4159:     PetscCall(PetscFree(b->values));
4160:     PetscCall(PetscFree(b->comps));
4161:     PetscCall(PetscFree(b));
4162:   }
4163:   PetscFunctionReturn(PETSC_SUCCESS);
4164: }

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

4169:   Not Collective

4171:   Input Parameters:
4172: + prob      - The `PetscDS` object
4173: . numFields - Number of new fields
4174: . fields    - Old field number for each new field
4175: . minDegree - Minimum degree for a discretization, or `PETSC_DETERMINE` for no limit
4176: - maxDegree - Maximum degree for a discretization, or `PETSC_DETERMINE` for no limit

4178:   Output Parameter:
4179: . newprob - The `PetscDS` copy

4181:   Level: intermediate

4183: .seealso: `PetscDS`, `PetscDSSelectEquations()`, `PetscDSCopyBoundary()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
4184: @*/
4185: PetscErrorCode PetscDSSelectDiscretizations(PetscDS prob, PetscInt numFields, const PetscInt fields[], PetscInt minDegree, PetscInt maxDegree, PetscDS newprob)
4186: {
4187:   PetscInt Nf, Nfn, fn;

4189:   PetscFunctionBegin;
4191:   if (fields) PetscAssertPointer(fields, 3);
4193:   PetscCall(PetscDSGetNumFields(prob, &Nf));
4194:   PetscCall(PetscDSGetNumFields(newprob, &Nfn));
4195:   numFields = numFields < 0 ? Nf : numFields;
4196:   for (fn = 0; fn < numFields; ++fn) {
4197:     const PetscInt f = fields ? fields[fn] : fn;
4198:     PetscObject    disc;
4199:     PetscClassId   id;

4201:     if (f >= Nf) continue;
4202:     PetscCall(PetscDSGetDiscretization(prob, f, &disc));
4203:     PetscCallContinue(PetscObjectGetClassId(disc, &id));
4204:     if (id == PETSCFE_CLASSID) {
4205:       PetscFE fe;

4207:       PetscCall(PetscFELimitDegree((PetscFE)disc, minDegree, maxDegree, &fe));
4208:       PetscCall(PetscDSSetDiscretization(newprob, fn, (PetscObject)fe));
4209:       PetscCall(PetscFEDestroy(&fe));
4210:     } else {
4211:       PetscCall(PetscDSSetDiscretization(newprob, fn, disc));
4212:     }
4213:   }
4214:   PetscFunctionReturn(PETSC_SUCCESS);
4215: }

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

4220:   Not Collective

4222:   Input Parameters:
4223: + prob      - The `PetscDS` object
4224: . numFields - Number of new fields
4225: - fields    - Old field number for each new field

4227:   Output Parameter:
4228: . newprob - The `PetscDS` copy

4230:   Level: intermediate

4232: .seealso: `PetscDS`, `PetscDSSelectDiscretizations()`, `PetscDSCopyBoundary()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
4233: @*/
4234: PetscErrorCode PetscDSSelectEquations(PetscDS prob, PetscInt numFields, const PetscInt fields[], PetscDS newprob)
4235: {
4236:   PetscInt Nf, Nfn, fn, gn;

4238:   PetscFunctionBegin;
4240:   if (fields) PetscAssertPointer(fields, 3);
4242:   PetscCall(PetscDSGetNumFields(prob, &Nf));
4243:   PetscCall(PetscDSGetNumFields(newprob, &Nfn));
4244:   PetscCheck(numFields <= Nfn, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_SIZ, "Number of fields %" PetscInt_FMT " to transfer must not be greater then the total number of fields %" PetscInt_FMT, numFields, Nfn);
4245:   for (fn = 0; fn < numFields; ++fn) {
4246:     const PetscInt   f = fields ? fields[fn] : fn;
4247:     PetscPointFunc   obj;
4248:     PetscPointFunc   f0, f1;
4249:     PetscBdPointFunc f0Bd, f1Bd;
4250:     PetscRiemannFunc r;

4252:     if (f >= Nf) continue;
4253:     PetscCall(PetscDSGetObjective(prob, f, &obj));
4254:     PetscCall(PetscDSGetResidual(prob, f, &f0, &f1));
4255:     PetscCall(PetscDSGetBdResidual(prob, f, &f0Bd, &f1Bd));
4256:     PetscCall(PetscDSGetRiemannSolver(prob, f, &r));
4257:     PetscCall(PetscDSSetObjective(newprob, fn, obj));
4258:     PetscCall(PetscDSSetResidual(newprob, fn, f0, f1));
4259:     PetscCall(PetscDSSetBdResidual(newprob, fn, f0Bd, f1Bd));
4260:     PetscCall(PetscDSSetRiemannSolver(newprob, fn, r));
4261:     for (gn = 0; gn < numFields; ++gn) {
4262:       const PetscInt  g = fields ? fields[gn] : gn;
4263:       PetscPointJac   g0, g1, g2, g3;
4264:       PetscPointJac   g0p, g1p, g2p, g3p;
4265:       PetscBdPointJac g0Bd, g1Bd, g2Bd, g3Bd;

4267:       if (g >= Nf) continue;
4268:       PetscCall(PetscDSGetJacobian(prob, f, g, &g0, &g1, &g2, &g3));
4269:       PetscCall(PetscDSGetJacobianPreconditioner(prob, f, g, &g0p, &g1p, &g2p, &g3p));
4270:       PetscCall(PetscDSGetBdJacobian(prob, f, g, &g0Bd, &g1Bd, &g2Bd, &g3Bd));
4271:       PetscCall(PetscDSSetJacobian(newprob, fn, gn, g0, g1, g2, g3));
4272:       PetscCall(PetscDSSetJacobianPreconditioner(newprob, fn, gn, g0p, g1p, g2p, g3p));
4273:       PetscCall(PetscDSSetBdJacobian(newprob, fn, gn, g0Bd, g1Bd, g2Bd, g3Bd));
4274:     }
4275:   }
4276:   PetscFunctionReturn(PETSC_SUCCESS);
4277: }

4279: /*@
4280:   PetscDSCopyEquations - Copy all pointwise function pointers to another `PetscDS`

4282:   Not Collective

4284:   Input Parameter:
4285: . prob - The `PetscDS` object

4287:   Output Parameter:
4288: . newprob - The `PetscDS` copy

4290:   Level: intermediate

4292: .seealso: `PetscDS`, `PetscDSCopyBoundary()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
4293: @*/
4294: PetscErrorCode PetscDSCopyEquations(PetscDS prob, PetscDS newprob)
4295: {
4296:   PetscWeakForm wf, newwf;
4297:   PetscInt      Nf, Ng;

4299:   PetscFunctionBegin;
4302:   PetscCall(PetscDSGetNumFields(prob, &Nf));
4303:   PetscCall(PetscDSGetNumFields(newprob, &Ng));
4304:   PetscCheck(Nf == Ng, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_SIZ, "Number of fields must match %" PetscInt_FMT " != %" PetscInt_FMT, Nf, Ng);
4305:   PetscCall(PetscDSGetWeakForm(prob, &wf));
4306:   PetscCall(PetscDSGetWeakForm(newprob, &newwf));
4307:   PetscCall(PetscWeakFormCopy(wf, newwf));
4308:   PetscFunctionReturn(PETSC_SUCCESS);
4309: }

4311: /*@
4312:   PetscDSCopyConstants - Copy all constants to another `PetscDS`

4314:   Not Collective

4316:   Input Parameter:
4317: . prob - The `PetscDS` object

4319:   Output Parameter:
4320: . newprob - The `PetscDS` copy

4322:   Level: intermediate

4324: .seealso: `PetscDS`, `PetscDSCopyBoundary()`, `PetscDSCopyEquations()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
4325: @*/
4326: PetscErrorCode PetscDSCopyConstants(PetscDS prob, PetscDS newprob)
4327: {
4328:   PetscInt           Nc;
4329:   const PetscScalar *constants;

4331:   PetscFunctionBegin;
4334:   PetscCall(PetscDSGetConstants(prob, &Nc, &constants));
4335:   PetscCall(PetscDSSetConstants(newprob, Nc, (PetscScalar *)constants));
4336:   PetscFunctionReturn(PETSC_SUCCESS);
4337: }

4339: /*@
4340:   PetscDSCopyExactSolutions - Copy all exact solutions to another `PetscDS`

4342:   Not Collective

4344:   Input Parameter:
4345: . ds - The `PetscDS` object

4347:   Output Parameter:
4348: . newds - The `PetscDS` copy

4350:   Level: intermediate

4352: .seealso: `PetscDS`, `PetscDSCopyBoundary()`, `PetscDSCopyEquations()`, `PetscDSCopyBounds()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
4353: @*/
4354: PetscErrorCode PetscDSCopyExactSolutions(PetscDS ds, PetscDS newds)
4355: {
4356:   PetscSimplePointFn *sol;
4357:   void               *ctx;
4358:   PetscInt            Nf, f;

4360:   PetscFunctionBegin;
4363:   PetscCall(PetscDSGetNumFields(ds, &Nf));
4364:   for (f = 0; f < Nf; ++f) {
4365:     PetscCall(PetscDSGetExactSolution(ds, f, &sol, &ctx));
4366:     PetscCall(PetscDSSetExactSolution(newds, f, sol, ctx));
4367:     PetscCall(PetscDSGetExactSolutionTimeDerivative(ds, f, &sol, &ctx));
4368:     PetscCall(PetscDSSetExactSolutionTimeDerivative(newds, f, sol, ctx));
4369:   }
4370:   PetscFunctionReturn(PETSC_SUCCESS);
4371: }

4373: /*@
4374:   PetscDSCopyBounds - Copy lower and upper solution bounds to another `PetscDS`

4376:   Not Collective

4378:   Input Parameter:
4379: . ds - The `PetscDS` object

4381:   Output Parameter:
4382: . newds - The `PetscDS` copy

4384:   Level: intermediate

4386: .seealso: `PetscDS`, `PetscDSCopyBoundary()`, `PetscDSCopyEquations()`, `PetscDSCopyExactSolutions()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
4387: @*/
4388: PetscErrorCode PetscDSCopyBounds(PetscDS ds, PetscDS newds)
4389: {
4390:   PetscSimplePointFn *bound;
4391:   void               *ctx;
4392:   PetscInt            Nf, f;

4394:   PetscFunctionBegin;
4397:   PetscCall(PetscDSGetNumFields(ds, &Nf));
4398:   for (f = 0; f < Nf; ++f) {
4399:     PetscCall(PetscDSGetLowerBound(ds, f, &bound, &ctx));
4400:     PetscCall(PetscDSSetLowerBound(newds, f, bound, ctx));
4401:     PetscCall(PetscDSGetUpperBound(ds, f, &bound, &ctx));
4402:     PetscCall(PetscDSSetUpperBound(newds, f, bound, ctx));
4403:   }
4404:   PetscFunctionReturn(PETSC_SUCCESS);
4405: }

4407: PetscErrorCode PetscDSCopy(PetscDS ds, PetscInt minDegree, PetscInt maxDegree, DM dmNew, PetscDS dsNew)
4408: {
4409:   DSBoundary b;
4410:   PetscInt   cdim, Nf, f, d;
4411:   PetscBool  isCohesive;
4412:   void      *ctx;

4414:   PetscFunctionBegin;
4415:   PetscCall(PetscDSCopyConstants(ds, dsNew));
4416:   PetscCall(PetscDSCopyExactSolutions(ds, dsNew));
4417:   PetscCall(PetscDSCopyBounds(ds, dsNew));
4418:   PetscCall(PetscDSSelectDiscretizations(ds, PETSC_DETERMINE, NULL, minDegree, maxDegree, dsNew));
4419:   PetscCall(PetscDSCopyEquations(ds, dsNew));
4420:   PetscCall(PetscDSGetNumFields(ds, &Nf));
4421:   for (f = 0; f < Nf; ++f) {
4422:     PetscCall(PetscDSGetContext(ds, f, &ctx));
4423:     PetscCall(PetscDSSetContext(dsNew, f, ctx));
4424:     PetscCall(PetscDSGetCohesive(ds, f, &isCohesive));
4425:     PetscCall(PetscDSSetCohesive(dsNew, f, isCohesive));
4426:     PetscCall(PetscDSGetJetDegree(ds, f, &d));
4427:     PetscCall(PetscDSSetJetDegree(dsNew, f, d));
4428:   }
4429:   if (Nf) {
4430:     PetscCall(PetscDSGetCoordinateDimension(ds, &cdim));
4431:     PetscCall(PetscDSSetCoordinateDimension(dsNew, cdim));
4432:   }
4433:   PetscCall(PetscDSCopyBoundary(ds, PETSC_DETERMINE, NULL, dsNew));
4434:   for (b = dsNew->boundary; b; b = b->next) {
4435:     PetscCall(DMGetLabel(dmNew, b->lname, &b->label));
4436:     /* Do not check if label exists here, since p4est calls this for the reference tree which does not have the labels */
4437:     //PetscCheck(b->label,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Label %s missing in new DM", name);
4438:   }
4439:   PetscFunctionReturn(PETSC_SUCCESS);
4440: }

4442: PetscErrorCode PetscDSGetHeightSubspace(PetscDS prob, PetscInt height, PetscDS *subprob)
4443: {
4444:   PetscInt dim, Nf, f;

4446:   PetscFunctionBegin;
4448:   PetscAssertPointer(subprob, 3);
4449:   if (height == 0) {
4450:     *subprob = prob;
4451:     PetscFunctionReturn(PETSC_SUCCESS);
4452:   }
4453:   PetscCall(PetscDSGetNumFields(prob, &Nf));
4454:   PetscCall(PetscDSGetSpatialDimension(prob, &dim));
4455:   PetscCheck(height <= dim, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_OUTOFRANGE, "DS can only handle height in [0, %" PetscInt_FMT "], not %" PetscInt_FMT, dim, height);
4456:   if (!prob->subprobs) PetscCall(PetscCalloc1(dim, &prob->subprobs));
4457:   if (!prob->subprobs[height - 1]) {
4458:     PetscInt cdim;

4460:     PetscCall(PetscDSCreate(PetscObjectComm((PetscObject)prob), &prob->subprobs[height - 1]));
4461:     PetscCall(PetscDSGetCoordinateDimension(prob, &cdim));
4462:     PetscCall(PetscDSSetCoordinateDimension(prob->subprobs[height - 1], cdim));
4463:     for (f = 0; f < Nf; ++f) {
4464:       PetscFE      subfe;
4465:       PetscObject  obj;
4466:       PetscClassId id;

4468:       PetscCall(PetscDSGetDiscretization(prob, f, &obj));
4469:       PetscCall(PetscObjectGetClassId(obj, &id));
4470:       if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetHeightSubspace((PetscFE)obj, height, &subfe));
4471:       else SETERRQ(PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONG, "Unsupported discretization type for field %" PetscInt_FMT, f);
4472:       PetscCall(PetscDSSetDiscretization(prob->subprobs[height - 1], f, (PetscObject)subfe));
4473:     }
4474:   }
4475:   *subprob = prob->subprobs[height - 1];
4476:   PetscFunctionReturn(PETSC_SUCCESS);
4477: }

4479: PetscErrorCode PetscDSPermuteQuadPoint(PetscDS ds, PetscInt ornt, PetscInt field, PetscInt q, PetscInt *qperm)
4480: {
4481:   IS              permIS;
4482:   PetscQuadrature quad;
4483:   DMPolytopeType  ct;
4484:   const PetscInt *perm;
4485:   PetscInt        Na, Nq;

4487:   PetscFunctionBeginHot;
4488:   PetscCall(PetscFEGetQuadrature((PetscFE)ds->disc[field], &quad));
4489:   PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL));
4490:   PetscCall(PetscQuadratureGetCellType(quad, &ct));
4491:   PetscCheck(q >= 0 && q < Nq, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Quadrature point %" PetscInt_FMT " is not in [0, %" PetscInt_FMT ")", q, Nq);
4492:   Na = DMPolytopeTypeGetNumArrangements(ct) / 2;
4493:   PetscCheck(ornt >= -Na && ornt < Na, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Orientation %" PetscInt_FMT " of %s is not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", ornt, DMPolytopeTypes[ct], -Na, Na);
4494:   if (!ds->quadPerm[(PetscInt)ct]) PetscCall(PetscQuadratureComputePermutations(quad, NULL, &ds->quadPerm[(PetscInt)ct]));
4495:   permIS = ds->quadPerm[(PetscInt)ct][ornt + Na];
4496:   PetscCall(ISGetIndices(permIS, &perm));
4497:   *qperm = perm[q];
4498:   PetscCall(ISRestoreIndices(permIS, &perm));
4499:   PetscFunctionReturn(PETSC_SUCCESS);
4500: }

4502: PetscErrorCode PetscDSGetDiscType_Internal(PetscDS ds, PetscInt f, PetscDiscType *disctype)
4503: {
4504:   PetscObject  obj;
4505:   PetscClassId id;
4506:   PetscInt     Nf;

4508:   PetscFunctionBegin;
4510:   PetscAssertPointer(disctype, 3);
4511:   *disctype = PETSC_DISC_NONE;
4512:   PetscCall(PetscDSGetNumFields(ds, &Nf));
4513:   PetscCheck(f < Nf, PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_SIZ, "Field %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, Nf);
4514:   PetscCall(PetscDSGetDiscretization(ds, f, &obj));
4515:   if (obj) {
4516:     PetscCall(PetscObjectGetClassId(obj, &id));
4517:     if (id == PETSCFE_CLASSID) *disctype = PETSC_DISC_FE;
4518:     else *disctype = PETSC_DISC_FV;
4519:   }
4520:   PetscFunctionReturn(PETSC_SUCCESS);
4521: }

4523: static PetscErrorCode PetscDSDestroy_Basic(PetscDS ds)
4524: {
4525:   PetscFunctionBegin;
4526:   PetscCall(PetscFree(ds->data));
4527:   PetscFunctionReturn(PETSC_SUCCESS);
4528: }

4530: static PetscErrorCode PetscDSInitialize_Basic(PetscDS ds)
4531: {
4532:   PetscFunctionBegin;
4533:   ds->ops->setfromoptions = NULL;
4534:   ds->ops->setup          = NULL;
4535:   ds->ops->view           = NULL;
4536:   ds->ops->destroy        = PetscDSDestroy_Basic;
4537:   PetscFunctionReturn(PETSC_SUCCESS);
4538: }

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

4543:   Level: intermediate

4545: .seealso: `PetscDSType`, `PetscDSCreate()`, `PetscDSSetType()`
4546: M*/

4548: PETSC_EXTERN PetscErrorCode PetscDSCreate_Basic(PetscDS ds)
4549: {
4550:   PetscDS_Basic *b;

4552:   PetscFunctionBegin;
4554:   PetscCall(PetscNew(&b));
4555:   ds->data = b;

4557:   PetscCall(PetscDSInitialize_Basic(ds));
4558:   PetscFunctionReturn(PETSC_SUCCESS);
4559: }