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;
559:   void               **tmpexactCtx, **tmpexactCtx_t;
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:   for (f = 0; f < Nf; ++f) tmpexactSol[f] = prob->exactSol[f];
596:   for (f = 0; f < Nf; ++f) tmpexactCtx[f] = prob->exactCtx[f];
597:   for (f = 0; f < Nf; ++f) tmpexactSol_t[f] = prob->exactSol_t[f];
598:   for (f = 0; f < Nf; ++f) tmpexactCtx_t[f] = prob->exactCtx_t[f];
599:   for (f = Nf; f < NfNew; ++f) tmpexactSol[f] = NULL;
600:   for (f = Nf; f < NfNew; ++f) tmpexactCtx[f] = NULL;
601:   for (f = Nf; f < NfNew; ++f) tmpexactSol_t[f] = NULL;
602:   for (f = Nf; f < NfNew; ++f) tmpexactCtx_t[f] = NULL;
603:   PetscCall(PetscFree4(prob->exactSol, prob->exactCtx, prob->exactSol_t, prob->exactCtx_t));
604:   prob->exactSol   = tmpexactSol;
605:   prob->exactCtx   = tmpexactCtx;
606:   prob->exactSol_t = tmpexactSol_t;
607:   prob->exactCtx_t = tmpexactCtx_t;
608:   PetscFunctionReturn(PETSC_SUCCESS);
609: }

611: /*@
612:   PetscDSDestroy - Destroys a `PetscDS` object

614:   Collective

616:   Input Parameter:
617: . ds - the `PetscDS` object to destroy

619:   Level: developer

621: .seealso: `PetscDSView()`
622: @*/
623: PetscErrorCode PetscDSDestroy(PetscDS *ds)
624: {
625:   PetscInt f;

627:   PetscFunctionBegin;
628:   if (!*ds) PetscFunctionReturn(PETSC_SUCCESS);

631:   if (--((PetscObject)*ds)->refct > 0) {
632:     *ds = NULL;
633:     PetscFunctionReturn(PETSC_SUCCESS);
634:   }
635:   ((PetscObject)*ds)->refct = 0;
636:   if ((*ds)->subprobs) {
637:     PetscInt dim, d;

639:     PetscCall(PetscDSGetSpatialDimension(*ds, &dim));
640:     for (d = 0; d < dim; ++d) PetscCall(PetscDSDestroy(&(*ds)->subprobs[d]));
641:   }
642:   PetscCall(PetscFree((*ds)->subprobs));
643:   PetscCall(PetscDSDestroyStructs_Static(*ds));
644:   for (f = 0; f < (*ds)->Nf; ++f) PetscCall(PetscObjectDereference((*ds)->disc[f]));
645:   PetscCall(PetscFree4((*ds)->disc, (*ds)->implicit, (*ds)->cohesive, (*ds)->jetDegree));
646:   PetscCall(PetscWeakFormDestroy(&(*ds)->wf));
647:   PetscCall(PetscFree2((*ds)->update, (*ds)->ctx));
648:   PetscCall(PetscFree4((*ds)->exactSol, (*ds)->exactCtx, (*ds)->exactSol_t, (*ds)->exactCtx_t));
649:   PetscTryTypeMethod(*ds, destroy);
650:   PetscCall(PetscDSDestroyBoundary(*ds));
651:   PetscCall(PetscFree((*ds)->constants));
652:   for (PetscInt c = 0; c < DM_NUM_POLYTOPES; ++c) {
653:     const PetscInt Na = DMPolytopeTypeGetNumArrangements((DMPolytopeType)c);
654:     if ((*ds)->quadPerm[c])
655:       for (PetscInt o = 0; o < Na; ++o) PetscCall(ISDestroy(&(*ds)->quadPerm[c][o]));
656:     PetscCall(PetscFree((*ds)->quadPerm[c]));
657:     (*ds)->quadPerm[c] = NULL;
658:   }
659:   PetscCall(PetscHeaderDestroy(ds));
660:   PetscFunctionReturn(PETSC_SUCCESS);
661: }

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

666:   Collective

668:   Input Parameter:
669: . comm - The communicator for the `PetscDS` object

671:   Output Parameter:
672: . ds - The `PetscDS` object

674:   Level: beginner

676: .seealso: `PetscDS`, `PetscDSSetType()`, `PETSCDSBASIC`, `PetscDSType`
677: @*/
678: PetscErrorCode PetscDSCreate(MPI_Comm comm, PetscDS *ds)
679: {
680:   PetscDS p;

682:   PetscFunctionBegin;
683:   PetscAssertPointer(ds, 2);
684:   PetscCall(PetscDSInitializePackage());

686:   PetscCall(PetscHeaderCreate(p, PETSCDS_CLASSID, "PetscDS", "Discrete System", "PetscDS", comm, PetscDSDestroy, PetscDSView));
687:   p->Nf               = 0;
688:   p->setup            = PETSC_FALSE;
689:   p->numConstants     = 0;
690:   p->numFuncConstants = 3; // Row and col fields, cell size
691:   p->dimEmbed         = -1;
692:   p->useJacPre        = PETSC_TRUE;
693:   p->forceQuad        = PETSC_TRUE;
694:   PetscCall(PetscMalloc1(p->numConstants + p->numFuncConstants, &p->constants));
695:   PetscCall(PetscWeakFormCreate(comm, &p->wf));
696:   PetscCall(PetscArrayzero(p->quadPerm, DM_NUM_POLYTOPES));
697:   *ds = p;
698:   PetscFunctionReturn(PETSC_SUCCESS);
699: }

701: /*@
702:   PetscDSGetNumFields - Returns the number of fields in the `PetscDS`

704:   Not Collective

706:   Input Parameter:
707: . prob - The `PetscDS` object

709:   Output Parameter:
710: . Nf - The number of fields

712:   Level: beginner

714: .seealso: `PetscDS`, `PetscDSGetSpatialDimension()`, `PetscDSCreate()`
715: @*/
716: PetscErrorCode PetscDSGetNumFields(PetscDS prob, PetscInt *Nf)
717: {
718:   PetscFunctionBegin;
720:   PetscAssertPointer(Nf, 2);
721:   *Nf = prob->Nf;
722:   PetscFunctionReturn(PETSC_SUCCESS);
723: }

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

728:   Not Collective

730:   Input Parameter:
731: . prob - The `PetscDS` object

733:   Output Parameter:
734: . dim - The spatial dimension

736:   Level: beginner

738: .seealso: `PetscDS`, `PetscDSGetCoordinateDimension()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
739: @*/
740: PetscErrorCode PetscDSGetSpatialDimension(PetscDS prob, PetscInt *dim)
741: {
742:   PetscFunctionBegin;
744:   PetscAssertPointer(dim, 2);
745:   *dim = 0;
746:   if (prob->Nf) {
747:     PetscObject  obj;
748:     PetscClassId id;

750:     PetscCall(PetscDSGetDiscretization(prob, 0, &obj));
751:     if (obj) {
752:       PetscCall(PetscObjectGetClassId(obj, &id));
753:       if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetSpatialDimension((PetscFE)obj, dim));
754:       else if (id == PETSCFV_CLASSID) PetscCall(PetscFVGetSpatialDimension((PetscFV)obj, dim));
755:       else SETERRQ(PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", 0);
756:     }
757:   }
758:   PetscFunctionReturn(PETSC_SUCCESS);
759: }

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

764:   Not Collective

766:   Input Parameter:
767: . prob - The `PetscDS` object

769:   Output Parameter:
770: . dimEmbed - The coordinate dimension

772:   Level: beginner

774: .seealso: `PetscDS`, `PetscDSSetCoordinateDimension()`, `PetscDSGetSpatialDimension()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
775: @*/
776: PetscErrorCode PetscDSGetCoordinateDimension(PetscDS prob, PetscInt *dimEmbed)
777: {
778:   PetscFunctionBegin;
780:   PetscAssertPointer(dimEmbed, 2);
781:   PetscCheck(prob->dimEmbed >= 0, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONGSTATE, "No coordinate dimension set for this DS");
782:   *dimEmbed = prob->dimEmbed;
783:   PetscFunctionReturn(PETSC_SUCCESS);
784: }

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

789:   Logically Collective

791:   Input Parameters:
792: + prob     - The `PetscDS` object
793: - dimEmbed - The coordinate dimension

795:   Level: beginner

797: .seealso: `PetscDS`, `PetscDSGetCoordinateDimension()`, `PetscDSGetSpatialDimension()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
798: @*/
799: PetscErrorCode PetscDSSetCoordinateDimension(PetscDS prob, PetscInt dimEmbed)
800: {
801:   PetscFunctionBegin;
803:   PetscCheck(dimEmbed >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Coordinate dimension must be non-negative, not %" PetscInt_FMT, dimEmbed);
804:   prob->dimEmbed = dimEmbed;
805:   PetscFunctionReturn(PETSC_SUCCESS);
806: }

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

811:   Not collective

813:   Input Parameter:
814: . ds - The `PetscDS` object

816:   Output Parameter:
817: . forceQuad - The flag

819:   Level: intermediate

821: .seealso: `PetscDS`, `PetscDSSetForceQuad()`, `PetscDSGetDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
822: @*/
823: PetscErrorCode PetscDSGetForceQuad(PetscDS ds, PetscBool *forceQuad)
824: {
825:   PetscFunctionBegin;
827:   PetscAssertPointer(forceQuad, 2);
828:   *forceQuad = ds->forceQuad;
829:   PetscFunctionReturn(PETSC_SUCCESS);
830: }

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

835:   Logically collective on ds

837:   Input Parameters:
838: + ds        - The `PetscDS` object
839: - forceQuad - The flag

841:   Level: intermediate

843: .seealso: `PetscDS`, `PetscDSGetForceQuad()`, `PetscDSGetDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
844: @*/
845: PetscErrorCode PetscDSSetForceQuad(PetscDS ds, PetscBool forceQuad)
846: {
847:   PetscFunctionBegin;
849:   ds->forceQuad = forceQuad;
850:   PetscFunctionReturn(PETSC_SUCCESS);
851: }

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

856:   Not Collective

858:   Input Parameter:
859: . ds - The `PetscDS` object

861:   Output Parameter:
862: . isCohesive - The flag

864:   Level: developer

866: .seealso: `PetscDS`, `PetscDSGetNumCohesive()`, `PetscDSGetCohesive()`, `PetscDSSetCohesive()`, `PetscDSCreate()`
867: @*/
868: PetscErrorCode PetscDSIsCohesive(PetscDS ds, PetscBool *isCohesive)
869: {
870:   PetscFunctionBegin;
872:   PetscAssertPointer(isCohesive, 2);
873:   *isCohesive = ds->isCohesive;
874:   PetscFunctionReturn(PETSC_SUCCESS);
875: }

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

880:   Not Collective

882:   Input Parameter:
883: . ds - The `PetscDS` object

885:   Output Parameter:
886: . numCohesive - The number of cohesive fields

888:   Level: developer

890: .seealso: `PetscDS`, `PetscDSSetCohesive()`, `PetscDSCreate()`
891: @*/
892: PetscErrorCode PetscDSGetNumCohesive(PetscDS ds, PetscInt *numCohesive)
893: {
894:   PetscInt f;

896:   PetscFunctionBegin;
898:   PetscAssertPointer(numCohesive, 2);
899:   *numCohesive = 0;
900:   for (f = 0; f < ds->Nf; ++f) *numCohesive += ds->cohesive[f] ? 1 : 0;
901:   PetscFunctionReturn(PETSC_SUCCESS);
902: }

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

907:   Not Collective

909:   Input Parameters:
910: + ds - The `PetscDS` object
911: - f  - The field index

913:   Output Parameter:
914: . isCohesive - The flag

916:   Level: developer

918: .seealso: `PetscDS`, `PetscDSSetCohesive()`, `PetscDSIsCohesive()`, `PetscDSCreate()`
919: @*/
920: PetscErrorCode PetscDSGetCohesive(PetscDS ds, PetscInt f, PetscBool *isCohesive)
921: {
922:   PetscFunctionBegin;
924:   PetscAssertPointer(isCohesive, 3);
925:   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);
926:   *isCohesive = ds->cohesive[f];
927:   PetscFunctionReturn(PETSC_SUCCESS);
928: }

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

933:   Not Collective

935:   Input Parameters:
936: + ds         - The `PetscDS` object
937: . f          - The field index
938: - isCohesive - The flag for a cohesive field

940:   Level: developer

942: .seealso: `PetscDS`, `PetscDSGetCohesive()`, `PetscDSIsCohesive()`, `PetscDSCreate()`
943: @*/
944: PetscErrorCode PetscDSSetCohesive(PetscDS ds, PetscInt f, PetscBool isCohesive)
945: {
946:   PetscInt i;

948:   PetscFunctionBegin;
950:   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);
951:   ds->cohesive[f] = isCohesive;
952:   ds->isCohesive  = PETSC_FALSE;
953:   for (i = 0; i < ds->Nf; ++i) ds->isCohesive = ds->isCohesive || ds->cohesive[f] ? PETSC_TRUE : PETSC_FALSE;
954:   PetscFunctionReturn(PETSC_SUCCESS);
955: }

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

960:   Not Collective

962:   Input Parameter:
963: . prob - The `PetscDS` object

965:   Output Parameter:
966: . dim - The total problem dimension

968:   Level: beginner

970: .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
971: @*/
972: PetscErrorCode PetscDSGetTotalDimension(PetscDS prob, PetscInt *dim)
973: {
974:   PetscFunctionBegin;
976:   PetscCall(PetscDSSetUp(prob));
977:   PetscAssertPointer(dim, 2);
978:   *dim = prob->totDim;
979:   PetscFunctionReturn(PETSC_SUCCESS);
980: }

982: /*@
983:   PetscDSGetTotalComponents - Returns the total number of components in this system

985:   Not Collective

987:   Input Parameter:
988: . prob - The `PetscDS` object

990:   Output Parameter:
991: . Nc - The total number of components

993:   Level: beginner

995: .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
996: @*/
997: PetscErrorCode PetscDSGetTotalComponents(PetscDS prob, PetscInt *Nc)
998: {
999:   PetscFunctionBegin;
1001:   PetscCall(PetscDSSetUp(prob));
1002:   PetscAssertPointer(Nc, 2);
1003:   *Nc = prob->totComp;
1004:   PetscFunctionReturn(PETSC_SUCCESS);
1005: }

1007: /*@
1008:   PetscDSGetDiscretization - Returns the discretization object for the given field

1010:   Not Collective

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

1016:   Output Parameter:
1017: . disc - The discretization object

1019:   Level: beginner

1021: .seealso: `PetscDS`, `PetscFE`, `PetscFV`, `PetscDSSetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1022: @*/
1023: PetscErrorCode PetscDSGetDiscretization(PetscDS prob, PetscInt f, PetscObject *disc)
1024: {
1025:   PetscFunctionBeginHot;
1027:   PetscAssertPointer(disc, 3);
1028:   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);
1029:   *disc = prob->disc[f];
1030:   PetscFunctionReturn(PETSC_SUCCESS);
1031: }

1033: /*@
1034:   PetscDSSetDiscretization - Sets the discretization object for the given field

1036:   Not Collective

1038:   Input Parameters:
1039: + prob - The `PetscDS` object
1040: . f    - The field number
1041: - disc - The discretization object

1043:   Level: beginner

1045: .seealso: `PetscDS`, `PetscFE`, `PetscFV`, `PetscDSGetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1046: @*/
1047: PetscErrorCode PetscDSSetDiscretization(PetscDS prob, PetscInt f, PetscObject disc)
1048: {
1049:   PetscFunctionBegin;
1051:   if (disc) PetscAssertPointer(disc, 3);
1052:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
1053:   PetscCall(PetscDSEnlarge_Static(prob, f + 1));
1054:   PetscCall(PetscObjectDereference(prob->disc[f]));
1055:   prob->disc[f] = disc;
1056:   PetscCall(PetscObjectReference(disc));
1057:   if (disc) {
1058:     PetscClassId id;

1060:     PetscCall(PetscObjectGetClassId(disc, &id));
1061:     if (id == PETSCFE_CLASSID) {
1062:       PetscCall(PetscDSSetImplicit(prob, f, PETSC_TRUE));
1063:     } else if (id == PETSCFV_CLASSID) {
1064:       PetscCall(PetscDSSetImplicit(prob, f, PETSC_FALSE));
1065:     }
1066:     PetscCall(PetscDSSetJetDegree(prob, f, 1));
1067:   }
1068:   PetscFunctionReturn(PETSC_SUCCESS);
1069: }

1071: /*@
1072:   PetscDSGetWeakForm - Returns the weak form object

1074:   Not Collective

1076:   Input Parameter:
1077: . ds - The `PetscDS` object

1079:   Output Parameter:
1080: . wf - The weak form object

1082:   Level: beginner

1084: .seealso: `PetscWeakForm`, `PetscDSSetWeakForm()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1085: @*/
1086: PetscErrorCode PetscDSGetWeakForm(PetscDS ds, PetscWeakForm *wf)
1087: {
1088:   PetscFunctionBegin;
1090:   PetscAssertPointer(wf, 2);
1091:   *wf = ds->wf;
1092:   PetscFunctionReturn(PETSC_SUCCESS);
1093: }

1095: /*@
1096:   PetscDSSetWeakForm - Sets the weak form object

1098:   Not Collective

1100:   Input Parameters:
1101: + ds - The `PetscDS` object
1102: - wf - The weak form object

1104:   Level: beginner

1106: .seealso: `PetscWeakForm`, `PetscDSGetWeakForm()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1107: @*/
1108: PetscErrorCode PetscDSSetWeakForm(PetscDS ds, PetscWeakForm wf)
1109: {
1110:   PetscFunctionBegin;
1113:   PetscCall(PetscObjectDereference((PetscObject)ds->wf));
1114:   ds->wf = wf;
1115:   PetscCall(PetscObjectReference((PetscObject)wf));
1116:   PetscCall(PetscWeakFormSetNumFields(wf, ds->Nf));
1117:   PetscFunctionReturn(PETSC_SUCCESS);
1118: }

1120: /*@
1121:   PetscDSAddDiscretization - Adds a discretization object

1123:   Not Collective

1125:   Input Parameters:
1126: + prob - The `PetscDS` object
1127: - disc - The boundary discretization object

1129:   Level: beginner

1131: .seealso: `PetscWeakForm`, `PetscDSGetDiscretization()`, `PetscDSSetDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1132: @*/
1133: PetscErrorCode PetscDSAddDiscretization(PetscDS prob, PetscObject disc)
1134: {
1135:   PetscFunctionBegin;
1136:   PetscCall(PetscDSSetDiscretization(prob, prob->Nf, disc));
1137:   PetscFunctionReturn(PETSC_SUCCESS);
1138: }

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

1143:   Not Collective

1145:   Input Parameter:
1146: . prob - The `PetscDS` object

1148:   Output Parameter:
1149: . q - The quadrature object

1151:   Level: intermediate

1153: .seealso: `PetscDS`, `PetscQuadrature`, `PetscDSSetImplicit()`, `PetscDSSetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1154: @*/
1155: PetscErrorCode PetscDSGetQuadrature(PetscDS prob, PetscQuadrature *q)
1156: {
1157:   PetscObject  obj;
1158:   PetscClassId id;

1160:   PetscFunctionBegin;
1161:   *q = NULL;
1162:   if (!prob->Nf) PetscFunctionReturn(PETSC_SUCCESS);
1163:   PetscCall(PetscDSGetDiscretization(prob, 0, &obj));
1164:   PetscCall(PetscObjectGetClassId(obj, &id));
1165:   if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetQuadrature((PetscFE)obj, q));
1166:   else if (id == PETSCFV_CLASSID) PetscCall(PetscFVGetQuadrature((PetscFV)obj, q));
1167:   else SETERRQ(PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", 0);
1168:   PetscFunctionReturn(PETSC_SUCCESS);
1169: }

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

1174:   Not Collective

1176:   Input Parameters:
1177: + prob - The `PetscDS` object
1178: - f    - The field number

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

1183:   Level: developer

1185: .seealso: `TSARKIMEX`, `PetscDS`, `PetscDSSetImplicit()`, `PetscDSSetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1186: @*/
1187: PetscErrorCode PetscDSGetImplicit(PetscDS prob, PetscInt f, PetscBool *implicit)
1188: {
1189:   PetscFunctionBegin;
1191:   PetscAssertPointer(implicit, 3);
1192:   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);
1193:   *implicit = prob->implicit[f];
1194:   PetscFunctionReturn(PETSC_SUCCESS);
1195: }

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

1200:   Not Collective

1202:   Input Parameters:
1203: + prob     - The `PetscDS` object
1204: . f        - The field number
1205: - implicit - The flag indicating what kind of solve to use for this field

1207:   Level: developer

1209: .seealso: `TSARKIMEX`, `PetscDSGetImplicit()`, `PetscDSSetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1210: @*/
1211: PetscErrorCode PetscDSSetImplicit(PetscDS prob, PetscInt f, PetscBool implicit)
1212: {
1213:   PetscFunctionBegin;
1215:   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);
1216:   prob->implicit[f] = implicit;
1217:   PetscFunctionReturn(PETSC_SUCCESS);
1218: }

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

1223:   Not Collective

1225:   Input Parameters:
1226: + ds - The `PetscDS` object
1227: - f  - The field number

1229:   Output Parameter:
1230: . k - The highest derivative we need to tabulate

1232:   Level: developer

1234: .seealso: `PetscDS`, `PetscDSSetJetDegree()`, `PetscDSSetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1235: @*/
1236: PetscErrorCode PetscDSGetJetDegree(PetscDS ds, PetscInt f, PetscInt *k)
1237: {
1238:   PetscFunctionBegin;
1240:   PetscAssertPointer(k, 3);
1241:   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);
1242:   *k = ds->jetDegree[f];
1243:   PetscFunctionReturn(PETSC_SUCCESS);
1244: }

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

1249:   Not Collective

1251:   Input Parameters:
1252: + ds - The `PetscDS` object
1253: . f  - The field number
1254: - k  - The highest derivative we need to tabulate

1256:   Level: developer

1258: .seealso: `PetscDS`, `PetscDSGetJetDegree()`, `PetscDSSetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1259: @*/
1260: PetscErrorCode PetscDSSetJetDegree(PetscDS ds, PetscInt f, PetscInt k)
1261: {
1262:   PetscFunctionBegin;
1264:   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);
1265:   ds->jetDegree[f] = k;
1266:   PetscFunctionReturn(PETSC_SUCCESS);
1267: }

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

1272:   Not Collective

1274:   Input Parameters:
1275: + ds - The `PetscDS`
1276: - f  - The test field number

1278:   Output Parameter:
1279: . obj - integrand for the test function term

1281:   Calling sequence of `obj`:
1282: + dim          - the spatial dimension
1283: . Nf           - the number of fields
1284: . NfAux        - the number of auxiliary fields
1285: . uOff         - the offset into u[] and u_t[] for each field
1286: . uOff_x       - the offset into u_x[] for each field
1287: . u            - each field evaluated at the current point
1288: . u_t          - the time derivative of each field evaluated at the current point
1289: . u_x          - the gradient of each field evaluated at the current point
1290: . aOff         - the offset into a[] and a_t[] for each auxiliary field
1291: . aOff_x       - the offset into a_x[] for each auxiliary field
1292: . a            - each auxiliary field evaluated at the current point
1293: . a_t          - the time derivative of each auxiliary field evaluated at the current point
1294: . a_x          - the gradient of auxiliary each field evaluated at the current point
1295: . t            - current time
1296: . x            - coordinates of the current point
1297: . numConstants - number of constant parameters
1298: . constants    - constant parameters
1299: - obj          - output values at the current point

1301:   Level: intermediate

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

1306: .seealso: `PetscDS`, `PetscDSSetObjective()`, `PetscDSGetResidual()`
1307: @*/
1308: 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[]))
1309: {
1310:   PetscPointFunc *tmp;
1311:   PetscInt        n;

1313:   PetscFunctionBegin;
1315:   PetscAssertPointer(obj, 3);
1316:   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);
1317:   PetscCall(PetscWeakFormGetObjective(ds->wf, NULL, 0, f, 0, &n, &tmp));
1318:   *obj = tmp ? tmp[0] : NULL;
1319:   PetscFunctionReturn(PETSC_SUCCESS);
1320: }

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

1325:   Not Collective

1327:   Input Parameters:
1328: + ds  - The `PetscDS`
1329: . f   - The test field number
1330: - obj - integrand for the test function term

1332:   Calling sequence of `obj`:
1333: + dim          - the spatial dimension
1334: . Nf           - the number of fields
1335: . NfAux        - the number of auxiliary fields
1336: . uOff         - the offset into u[] and u_t[] for each field
1337: . uOff_x       - the offset into u_x[] for each field
1338: . u            - each field evaluated at the current point
1339: . u_t          - the time derivative of each field evaluated at the current point
1340: . u_x          - the gradient of each field evaluated at the current point
1341: . aOff         - the offset into a[] and a_t[] for each auxiliary field
1342: . aOff_x       - the offset into a_x[] for each auxiliary field
1343: . a            - each auxiliary field evaluated at the current point
1344: . a_t          - the time derivative of each auxiliary field evaluated at the current point
1345: . a_x          - the gradient of auxiliary each field evaluated at the current point
1346: . t            - current time
1347: . x            - coordinates of the current point
1348: . numConstants - number of constant parameters
1349: . constants    - constant parameters
1350: - obj          - output values at the current point

1352:   Level: intermediate

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

1357: .seealso: `PetscDS`, `PetscDSGetObjective()`, `PetscDSSetResidual()`
1358: @*/
1359: 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[]))
1360: {
1361:   PetscFunctionBegin;
1364:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
1365:   PetscCall(PetscWeakFormSetIndexObjective(ds->wf, NULL, 0, f, 0, 0, obj));
1366:   PetscFunctionReturn(PETSC_SUCCESS);
1367: }

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

1372:   Not Collective

1374:   Input Parameters:
1375: + ds - The `PetscDS`
1376: - f  - The test field number

1378:   Output Parameters:
1379: + f0 - integrand for the test function term
1380: - f1 - integrand for the test function gradient term

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

1402:   Level: intermediate

1404:   Note:
1405:   `f1` has an identical form and is omitted for brevity.

1407:   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)$

1409: .seealso: `PetscDS`, `PetscDSSetResidual()`
1410: @*/
1411: 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[]))
1412: {
1413:   PetscPointFunc *tmp0, *tmp1;
1414:   PetscInt        n0, n1;

1416:   PetscFunctionBegin;
1418:   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);
1419:   PetscCall(PetscWeakFormGetResidual(ds->wf, NULL, 0, f, 0, &n0, &tmp0, &n1, &tmp1));
1420:   *f0 = tmp0 ? tmp0[0] : NULL;
1421:   *f1 = tmp1 ? tmp1[0] : NULL;
1422:   PetscFunctionReturn(PETSC_SUCCESS);
1423: }

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

1428:   Not Collective

1430:   Input Parameters:
1431: + ds - The `PetscDS`
1432: . f  - The test field number
1433: . f0 - integrand for the test function term
1434: - f1 - integrand for the test function gradient term

1436:   Calling sequence of `f0`:
1437: + dim          - the spatial dimension
1438: . Nf           - the number of fields
1439: . NfAux        - the number of auxiliary fields
1440: . uOff         - the offset into u[] and u_t[] for each field
1441: . uOff_x       - the offset into u_x[] for each field
1442: . u            - each field evaluated at the current point
1443: . u_t          - the time derivative of each field evaluated at the current point
1444: . u_x          - the gradient of each field evaluated at the current point
1445: . aOff         - the offset into a[] and a_t[] for each auxiliary field
1446: . aOff_x       - the offset into a_x[] for each auxiliary field
1447: . a            - each auxiliary field evaluated at the current point
1448: . a_t          - the time derivative of each auxiliary field evaluated at the current point
1449: . a_x          - the gradient of auxiliary each field evaluated at the current point
1450: . t            - current time
1451: . x            - coordinates of the current point
1452: . numConstants - number of constant parameters
1453: . constants    - constant parameters
1454: - f0           - output values at the current point

1456:   Level: intermediate

1458:   Note:
1459:   `f1` has an identical form and is omitted for brevity.

1461:   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)$

1463: .seealso: `PetscDS`, `PetscDSGetResidual()`
1464: @*/
1465: 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[]))
1466: {
1467:   PetscFunctionBegin;
1471:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
1472:   PetscCall(PetscWeakFormSetIndexResidual(ds->wf, NULL, 0, f, 0, 0, f0, 0, f1));
1473:   PetscFunctionReturn(PETSC_SUCCESS);
1474: }

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

1479:   Not Collective

1481:   Input Parameters:
1482: + ds - The `PetscDS`
1483: - f  - The test field number

1485:   Output Parameters:
1486: + f0 - integrand for the test function term
1487: - f1 - integrand for the test function gradient term

1489:   Calling sequence of `f0`:
1490: + dim          - the spatial dimension
1491: . Nf           - the number of fields
1492: . NfAux        - the number of auxiliary fields
1493: . uOff         - the offset into u[] and u_t[] for each field
1494: . uOff_x       - the offset into u_x[] for each field
1495: . u            - each field evaluated at the current point
1496: . u_t          - the time derivative of each field evaluated at the current point
1497: . u_x          - the gradient of each field evaluated at the current point
1498: . aOff         - the offset into a[] and a_t[] for each auxiliary field
1499: . aOff_x       - the offset into a_x[] for each auxiliary field
1500: . a            - each auxiliary field evaluated at the current point
1501: . a_t          - the time derivative of each auxiliary field evaluated at the current point
1502: . a_x          - the gradient of auxiliary each field evaluated at the current point
1503: . t            - current time
1504: . x            - coordinates of the current point
1505: . numConstants - number of constant parameters
1506: . constants    - constant parameters
1507: - f0           - output values at the current point

1509:   Level: intermediate

1511:   Note:
1512:   `f1` has an identical form and is omitted for brevity.

1514:   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)$

1516: .seealso: `PetscDS`, `PetscDSSetRHSResidual()`
1517: @*/
1518: 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[]))
1519: {
1520:   PetscPointFunc *tmp0, *tmp1;
1521:   PetscInt        n0, n1;

1523:   PetscFunctionBegin;
1525:   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);
1526:   PetscCall(PetscWeakFormGetResidual(ds->wf, NULL, 0, f, 100, &n0, &tmp0, &n1, &tmp1));
1527:   *f0 = tmp0 ? tmp0[0] : NULL;
1528:   *f1 = tmp1 ? tmp1[0] : NULL;
1529:   PetscFunctionReturn(PETSC_SUCCESS);
1530: }

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

1535:   Not Collective

1537:   Input Parameters:
1538: + ds - The `PetscDS`
1539: . f  - The test field number
1540: . f0 - integrand for the test function term
1541: - f1 - integrand for the test function gradient term

1543:   Calling sequence for the callbacks `f0`:
1544: + dim          - the spatial dimension
1545: . Nf           - the number of fields
1546: . NfAux        - the number of auxiliary fields
1547: . uOff         - the offset into u[] and u_t[] for each field
1548: . uOff_x       - the offset into u_x[] for each field
1549: . u            - each field evaluated at the current point
1550: . u_t          - the time derivative of each field evaluated at the current point
1551: . u_x          - the gradient of each field evaluated at the current point
1552: . aOff         - the offset into a[] and a_t[] for each auxiliary field
1553: . aOff_x       - the offset into a_x[] for each auxiliary field
1554: . a            - each auxiliary field evaluated at the current point
1555: . a_t          - the time derivative of each auxiliary field evaluated at the current point
1556: . a_x          - the gradient of auxiliary each field evaluated at the current point
1557: . t            - current time
1558: . x            - coordinates of the current point
1559: . numConstants - number of constant parameters
1560: . constants    - constant parameters
1561: - f0           - output values at the current point

1563:   Level: intermediate

1565:   Note:
1566:   `f1` has an identical form and is omitted for brevity.

1568:   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)$

1570: .seealso: `PetscDS`, `PetscDSGetResidual()`
1571: @*/
1572: 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[]))
1573: {
1574:   PetscFunctionBegin;
1578:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
1579:   PetscCall(PetscWeakFormSetIndexResidual(ds->wf, NULL, 0, f, 100, 0, f0, 0, f1));
1580:   PetscFunctionReturn(PETSC_SUCCESS);
1581: }

1583: /*@
1584:   PetscDSHasJacobian - Checks that the Jacobian functions have been set

1586:   Not Collective

1588:   Input Parameter:
1589: . ds - The `PetscDS`

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

1594:   Level: intermediate

1596: .seealso: `PetscDS`, `PetscDSGetJacobianPreconditioner()`, `PetscDSSetJacobianPreconditioner()`, `PetscDSGetJacobian()`
1597: @*/
1598: PetscErrorCode PetscDSHasJacobian(PetscDS ds, PetscBool *hasJac)
1599: {
1600:   PetscFunctionBegin;
1602:   PetscCall(PetscWeakFormHasJacobian(ds->wf, hasJac));
1603:   PetscFunctionReturn(PETSC_SUCCESS);
1604: }

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

1609:   Not Collective

1611:   Input Parameters:
1612: + ds - The `PetscDS`
1613: . f  - The test field number
1614: - g  - The field number

1616:   Output Parameters:
1617: + g0 - integrand for the test and basis function term
1618: . g1 - integrand for the test function and basis function gradient term
1619: . g2 - integrand for the test function gradient and basis function term
1620: - g3 - integrand for the test function gradient and basis function gradient term

1622:   Calling sequence of `g0`:
1623: + dim          - the spatial dimension
1624: . Nf           - the number of fields
1625: . NfAux        - the number of auxiliary fields
1626: . uOff         - the offset into u[] and u_t[] for each field
1627: . uOff_x       - the offset into u_x[] for each field
1628: . u            - each field evaluated at the current point
1629: . u_t          - the time derivative of each field evaluated at the current point
1630: . u_x          - the gradient of each field evaluated at the current point
1631: . aOff         - the offset into a[] and a_t[] for each auxiliary field
1632: . aOff_x       - the offset into a_x[] for each auxiliary field
1633: . a            - each auxiliary field evaluated at the current point
1634: . a_t          - the time derivative of each auxiliary field evaluated at the current point
1635: . a_x          - the gradient of auxiliary each field evaluated at the current point
1636: . t            - current time
1637: . u_tShift     - the multiplier a for dF/dU_t
1638: . x            - coordinates of the current point
1639: . numConstants - number of constant parameters
1640: . constants    - constant parameters
1641: - g0           - output values at the current point

1643:   Level: intermediate

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

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

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

1654: .seealso: `PetscDS`, `PetscDSSetJacobian()`
1655: @*/
1656: 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[]))
1657: {
1658:   PetscPointJac *tmp0, *tmp1, *tmp2, *tmp3;
1659:   PetscInt       n0, n1, n2, n3;

1661:   PetscFunctionBegin;
1663:   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);
1664:   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);
1665:   PetscCall(PetscWeakFormGetJacobian(ds->wf, NULL, 0, f, g, 0, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3));
1666:   *g0 = tmp0 ? tmp0[0] : NULL;
1667:   *g1 = tmp1 ? tmp1[0] : NULL;
1668:   *g2 = tmp2 ? tmp2[0] : NULL;
1669:   *g3 = tmp3 ? tmp3[0] : NULL;
1670:   PetscFunctionReturn(PETSC_SUCCESS);
1671: }

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

1676:   Not Collective

1678:   Input Parameters:
1679: + ds - The `PetscDS`
1680: . f  - The test field number
1681: . g  - The field number
1682: . g0 - integrand for the test and basis function term
1683: . g1 - integrand for the test function and basis function gradient term
1684: . g2 - integrand for the test function gradient and basis function term
1685: - g3 - integrand for the test function gradient and basis function gradient term

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

1708:   Level: intermediate

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

1713:   We are using a first order FEM model for the weak form\:
1714:   \int_\Omega \phi g_0(u, u_t, \nabla u, x, t) \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi

1716: .seealso: `PetscDS`, `PetscDSGetJacobian()`
1717: @*/
1718: 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[]))
1719: {
1720:   PetscFunctionBegin;
1726:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
1727:   PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", g);
1728:   PetscCall(PetscWeakFormSetIndexJacobian(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3));
1729:   PetscFunctionReturn(PETSC_SUCCESS);
1730: }

1732: /*@
1733:   PetscDSUseJacobianPreconditioner - Set whether to construct a Jacobian preconditioner

1735:   Not Collective

1737:   Input Parameters:
1738: + prob      - The `PetscDS`
1739: - useJacPre - flag that enables construction of a Jacobian preconditioner

1741:   Level: intermediate

1743:   Developer Note:
1744:   Should be called `PetscDSSetUseJacobianPreconditioner()`

1746: .seealso: `PetscDS`, `PetscDSGetJacobianPreconditioner()`, `PetscDSSetJacobianPreconditioner()`, `PetscDSGetJacobian()`
1747: @*/
1748: PetscErrorCode PetscDSUseJacobianPreconditioner(PetscDS prob, PetscBool useJacPre)
1749: {
1750:   PetscFunctionBegin;
1752:   prob->useJacPre = useJacPre;
1753:   PetscFunctionReturn(PETSC_SUCCESS);
1754: }

1756: /*@
1757:   PetscDSHasJacobianPreconditioner - Checks if a Jacobian preconditioner matrix has been set

1759:   Not Collective

1761:   Input Parameter:
1762: . ds - The `PetscDS`

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

1767:   Level: intermediate

1769: .seealso: `PetscDS`, `PetscDSGetJacobianPreconditioner()`, `PetscDSSetJacobianPreconditioner()`, `PetscDSGetJacobian()`
1770: @*/
1771: PetscErrorCode PetscDSHasJacobianPreconditioner(PetscDS ds, PetscBool *hasJacPre)
1772: {
1773:   PetscFunctionBegin;
1775:   *hasJacPre = PETSC_FALSE;
1776:   if (!ds->useJacPre) PetscFunctionReturn(PETSC_SUCCESS);
1777:   PetscCall(PetscWeakFormHasJacobianPreconditioner(ds->wf, hasJacPre));
1778:   PetscFunctionReturn(PETSC_SUCCESS);
1779: }

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

1785:   Not Collective

1787:   Input Parameters:
1788: + ds - The `PetscDS`
1789: . f  - The test field number
1790: - g  - The field number

1792:   Output Parameters:
1793: + g0 - integrand for the test and basis function term
1794: . g1 - integrand for the test function and basis function gradient term
1795: . g2 - integrand for the test function gradient and basis function term
1796: - g3 - integrand for the test function gradient and basis function gradient term

1798:   Calling sequence of `g0`:
1799: + dim          - the spatial dimension
1800: . Nf           - the number of fields
1801: . NfAux        - the number of auxiliary fields
1802: . uOff         - the offset into u[] and u_t[] for each field
1803: . uOff_x       - the offset into u_x[] for each field
1804: . u            - each field evaluated at the current point
1805: . u_t          - the time derivative of each field evaluated at the current point
1806: . u_x          - the gradient of each field evaluated at the current point
1807: . aOff         - the offset into a[] and a_t[] for each auxiliary field
1808: . aOff_x       - the offset into a_x[] for each auxiliary field
1809: . a            - each auxiliary field evaluated at the current point
1810: . a_t          - the time derivative of each auxiliary field evaluated at the current point
1811: . a_x          - the gradient of auxiliary each field evaluated at the current point
1812: . t            - current time
1813: . u_tShift     - the multiplier a for dF/dU_t
1814: . x            - coordinates of the current point
1815: . numConstants - number of constant parameters
1816: . constants    - constant parameters
1817: - g0           - output values at the current point

1819:   Level: intermediate

1821:   Note:
1822:   `g1`, `g2`, and `g3` have identical calling sequences to `g0` and are omitted for brevity.
1823:   We are using a first order FEM model for the weak form\:
1824:   \int_\Omega \phi g_0(u, u_t, \nabla u, x, t) \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi

1826: .seealso: `PetscDS`, `PetscDSSetJacobianPreconditioner()`, `PetscDSGetJacobian()`
1827: @*/
1828: 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[]))
1829: {
1830:   PetscPointJac *tmp0, *tmp1, *tmp2, *tmp3;
1831:   PetscInt       n0, n1, n2, n3;

1833:   PetscFunctionBegin;
1835:   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);
1836:   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);
1837:   PetscCall(PetscWeakFormGetJacobianPreconditioner(ds->wf, NULL, 0, f, g, 0, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3));
1838:   *g0 = tmp0 ? tmp0[0] : NULL;
1839:   *g1 = tmp1 ? tmp1[0] : NULL;
1840:   *g2 = tmp2 ? tmp2[0] : NULL;
1841:   *g3 = tmp3 ? tmp3[0] : NULL;
1842:   PetscFunctionReturn(PETSC_SUCCESS);
1843: }

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

1849:   Not Collective

1851:   Input Parameters:
1852: + ds - The `PetscDS`
1853: . f  - The test field number
1854: . g  - The field number
1855: . g0 - integrand for the test and basis function term
1856: . g1 - integrand for the test function and basis function gradient term
1857: . g2 - integrand for the test function gradient and basis function term
1858: - g3 - integrand for the test function gradient and basis function gradient term

1860:   Calling sequence of `g0`:
1861: + dim          - the spatial dimension
1862: . Nf           - the number of fields
1863: . NfAux        - the number of auxiliary fields
1864: . uOff         - the offset into u[] and u_t[] for each field
1865: . uOff_x       - the offset into u_x[] for each field
1866: . u            - each field evaluated at the current point
1867: . u_t          - the time derivative of each field evaluated at the current point
1868: . u_x          - the gradient of each field evaluated at the current point
1869: . aOff         - the offset into a[] and a_t[] for each auxiliary field
1870: . aOff_x       - the offset into a_x[] for each auxiliary field
1871: . a            - each auxiliary field evaluated at the current point
1872: . a_t          - the time derivative of each auxiliary field evaluated at the current point
1873: . a_x          - the gradient of auxiliary each field evaluated at the current point
1874: . t            - current time
1875: . u_tShift     - the multiplier a for dF/dU_t
1876: . x            - coordinates of the current point
1877: . numConstants - number of constant parameters
1878: . constants    - constant parameters
1879: - g0           - output values at the current point

1881:   Level: intermediate

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

1886:   We are using a first order FEM model for the weak form\:
1887:   \int_\Omega \phi g_0(u, u_t, \nabla u, x, t) \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi

1889: .seealso: `PetscDS`, `PetscDSGetJacobianPreconditioner()`, `PetscDSSetJacobian()`
1890: @*/
1891: 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[]))
1892: {
1893:   PetscFunctionBegin;
1899:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
1900:   PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", g);
1901:   PetscCall(PetscWeakFormSetIndexJacobianPreconditioner(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3));
1902:   PetscFunctionReturn(PETSC_SUCCESS);
1903: }

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

1908:   Not Collective

1910:   Input Parameter:
1911: . ds - The `PetscDS`

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

1916:   Level: intermediate

1918: .seealso: `PetscDS`, `PetscDSGetDynamicJacobian()`, `PetscDSSetDynamicJacobian()`, `PetscDSGetJacobian()`
1919: @*/
1920: PetscErrorCode PetscDSHasDynamicJacobian(PetscDS ds, PetscBool *hasDynJac)
1921: {
1922:   PetscFunctionBegin;
1924:   PetscCall(PetscWeakFormHasDynamicJacobian(ds->wf, hasDynJac));
1925:   PetscFunctionReturn(PETSC_SUCCESS);
1926: }

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

1931:   Not Collective

1933:   Input Parameters:
1934: + ds - The `PetscDS`
1935: . f  - The test field number
1936: - g  - The field number

1938:   Output Parameters:
1939: + g0 - integrand for the test and basis function term
1940: . g1 - integrand for the test function and basis function gradient term
1941: . g2 - integrand for the test function gradient and basis function term
1942: - g3 - integrand for the test function gradient and basis function gradient term

1944:   Calling sequence of `g0`:
1945: + dim          - the spatial dimension
1946: . Nf           - the number of fields
1947: . NfAux        - the number of auxiliary fields
1948: . uOff         - the offset into u[] and u_t[] for each field
1949: . uOff_x       - the offset into u_x[] for each field
1950: . u            - each field evaluated at the current point
1951: . u_t          - the time derivative of each field evaluated at the current point
1952: . u_x          - the gradient of each field evaluated at the current point
1953: . aOff         - the offset into a[] and a_t[] for each auxiliary field
1954: . aOff_x       - the offset into a_x[] for each auxiliary field
1955: . a            - each auxiliary field evaluated at the current point
1956: . a_t          - the time derivative of each auxiliary field evaluated at the current point
1957: . a_x          - the gradient of auxiliary each field evaluated at the current point
1958: . t            - current time
1959: . u_tShift     - the multiplier a for dF/dU_t
1960: . x            - coordinates of the current point
1961: . numConstants - number of constant parameters
1962: . constants    - constant parameters
1963: - g0           - output values at the current point

1965:   Level: intermediate

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

1970:   We are using a first order FEM model for the weak form\:
1971:   \int_\Omega \phi g_0(u, u_t, \nabla u, x, t) \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi

1973: .seealso: `PetscDS`, `PetscDSSetJacobian()`
1974: @*/
1975: 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[]))
1976: {
1977:   PetscPointJac *tmp0, *tmp1, *tmp2, *tmp3;
1978:   PetscInt       n0, n1, n2, n3;

1980:   PetscFunctionBegin;
1982:   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);
1983:   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);
1984:   PetscCall(PetscWeakFormGetDynamicJacobian(ds->wf, NULL, 0, f, g, 0, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3));
1985:   *g0 = tmp0 ? tmp0[0] : NULL;
1986:   *g1 = tmp1 ? tmp1[0] : NULL;
1987:   *g2 = tmp2 ? tmp2[0] : NULL;
1988:   *g3 = tmp3 ? tmp3[0] : NULL;
1989:   PetscFunctionReturn(PETSC_SUCCESS);
1990: }

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

1995:   Not Collective

1997:   Input Parameters:
1998: + ds - The `PetscDS`
1999: . f  - The test field number
2000: . g  - The field number
2001: . g0 - integrand for the test and basis function term
2002: . g1 - integrand for the test function and basis function gradient term
2003: . g2 - integrand for the test function gradient and basis function term
2004: - g3 - integrand for the test function gradient and basis function gradient term

2006:   Calling sequence of `g0`:
2007: + dim          - the spatial dimension
2008: . Nf           - the number of fields
2009: . NfAux        - the number of auxiliary fields
2010: . uOff         - the offset into u[] and u_t[] for each field
2011: . uOff_x       - the offset into u_x[] for each field
2012: . u            - each field evaluated at the current point
2013: . u_t          - the time derivative of each field evaluated at the current point
2014: . u_x          - the gradient of each field evaluated at the current point
2015: . aOff         - the offset into a[] and a_t[] for each auxiliary field
2016: . aOff_x       - the offset into a_x[] for each auxiliary field
2017: . a            - each auxiliary field evaluated at the current point
2018: . a_t          - the time derivative of each auxiliary field evaluated at the current point
2019: . a_x          - the gradient of auxiliary each field evaluated at the current point
2020: . t            - current time
2021: . u_tShift     - the multiplier a for dF/dU_t
2022: . x            - coordinates of the current point
2023: . numConstants - number of constant parameters
2024: . constants    - constant parameters
2025: - g0           - output values at the current point

2027:   Level: intermediate

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

2032:   We are using a first order FEM model for the weak form\:
2033:   \int_\Omega \phi g_0(u, u_t, \nabla u, x, t) \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi

2035: .seealso: `PetscDS`, `PetscDSGetJacobian()`
2036: @*/
2037: 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[]))
2038: {
2039:   PetscFunctionBegin;
2045:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2046:   PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", g);
2047:   PetscCall(PetscWeakFormSetIndexDynamicJacobian(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3));
2048:   PetscFunctionReturn(PETSC_SUCCESS);
2049: }

2051: /*@C
2052:   PetscDSGetRiemannSolver - Returns the Riemann solver for the given field

2054:   Not Collective

2056:   Input Parameters:
2057: + ds - The `PetscDS` object
2058: - f  - The field number

2060:   Output Parameter:
2061: . r - Riemann solver

2063:   Calling sequence of `r`:
2064: + dim          - The spatial dimension
2065: . Nf           - The number of fields
2066: . x            - The coordinates at a point on the interface
2067: . n            - The normal vector to the interface
2068: . uL           - The state vector to the left of the interface
2069: . uR           - The state vector to the right of the interface
2070: . flux         - output array of flux through the interface
2071: . numConstants - number of constant parameters
2072: . constants    - constant parameters
2073: - ctx          - optional user context

2075:   Level: intermediate

2077: .seealso: `PetscDS`, `PetscDSSetRiemannSolver()`
2078: @*/
2079: 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))
2080: {
2081:   PetscRiemannFunc *tmp;
2082:   PetscInt          n;

2084:   PetscFunctionBegin;
2086:   PetscAssertPointer(r, 3);
2087:   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);
2088:   PetscCall(PetscWeakFormGetRiemannSolver(ds->wf, NULL, 0, f, 0, &n, &tmp));
2089:   *r = tmp ? tmp[0] : NULL;
2090:   PetscFunctionReturn(PETSC_SUCCESS);
2091: }

2093: /*@C
2094:   PetscDSSetRiemannSolver - Sets the Riemann solver for the given field

2096:   Not Collective

2098:   Input Parameters:
2099: + ds - The `PetscDS` object
2100: . f  - The field number
2101: - r  - Riemann solver

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

2115:   Level: intermediate

2117: .seealso: `PetscDS`, `PetscDSGetRiemannSolver()`
2118: @*/
2119: 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))
2120: {
2121:   PetscFunctionBegin;
2124:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2125:   PetscCall(PetscWeakFormSetIndexRiemannSolver(ds->wf, NULL, 0, f, 0, 0, r));
2126:   PetscFunctionReturn(PETSC_SUCCESS);
2127: }

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

2132:   Not Collective

2134:   Input Parameters:
2135: + ds - The `PetscDS`
2136: - f  - The field number

2138:   Output Parameter:
2139: . update - update function

2141:   Calling sequence of `update`:
2142: + dim          - the spatial dimension
2143: . Nf           - the number of fields
2144: . NfAux        - the number of auxiliary fields
2145: . uOff         - the offset into u[] and u_t[] for each field
2146: . uOff_x       - the offset into u_x[] for each field
2147: . u            - each field evaluated at the current point
2148: . u_t          - the time derivative of each field evaluated at the current point
2149: . u_x          - the gradient of each field evaluated at the current point
2150: . aOff         - the offset into a[] and a_t[] for each auxiliary field
2151: . aOff_x       - the offset into a_x[] for each auxiliary field
2152: . a            - each auxiliary field evaluated at the current point
2153: . a_t          - the time derivative of each auxiliary field evaluated at the current point
2154: . a_x          - the gradient of auxiliary each field evaluated at the current point
2155: . t            - current time
2156: . x            - coordinates of the current point
2157: . numConstants - number of constant parameters
2158: . constants    - constant parameters
2159: - uNew         - new value for field at the current point

2161:   Level: intermediate

2163: .seealso: `PetscDS`, `PetscDSSetUpdate()`, `PetscDSSetResidual()`
2164: @*/
2165: 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[]))
2166: {
2167:   PetscFunctionBegin;
2169:   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);
2170:   if (update) {
2171:     PetscAssertPointer(update, 3);
2172:     *update = ds->update[f];
2173:   }
2174:   PetscFunctionReturn(PETSC_SUCCESS);
2175: }

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

2180:   Not Collective

2182:   Input Parameters:
2183: + ds     - The `PetscDS`
2184: . f      - The field number
2185: - update - update function

2187:   Calling sequence of `update`:
2188: + dim          - the spatial dimension
2189: . Nf           - the number of fields
2190: . NfAux        - the number of auxiliary fields
2191: . uOff         - the offset into u[] and u_t[] for each field
2192: . uOff_x       - the offset into u_x[] for each field
2193: . u            - each field evaluated at the current point
2194: . u_t          - the time derivative of each field evaluated at the current point
2195: . u_x          - the gradient of each field evaluated at the current point
2196: . aOff         - the offset into a[] and a_t[] for each auxiliary field
2197: . aOff_x       - the offset into a_x[] for each auxiliary field
2198: . a            - each auxiliary field evaluated at the current point
2199: . a_t          - the time derivative of each auxiliary field evaluated at the current point
2200: . a_x          - the gradient of auxiliary each field evaluated at the current point
2201: . t            - current time
2202: . x            - coordinates of the current point
2203: . numConstants - number of constant parameters
2204: . constants    - constant parameters
2205: - uNew         - new field values at the current point

2207:   Level: intermediate

2209: .seealso: `PetscDS`, `PetscDSGetResidual()`
2210: @*/
2211: 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[]))
2212: {
2213:   PetscFunctionBegin;
2216:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2217:   PetscCall(PetscDSEnlarge_Static(ds, f + 1));
2218:   ds->update[f] = update;
2219:   PetscFunctionReturn(PETSC_SUCCESS);
2220: }

2222: PetscErrorCode PetscDSGetContext(PetscDS ds, PetscInt f, void *ctx)
2223: {
2224:   PetscFunctionBegin;
2226:   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);
2227:   PetscAssertPointer(ctx, 3);
2228:   *(void **)ctx = ds->ctx[f];
2229:   PetscFunctionReturn(PETSC_SUCCESS);
2230: }

2232: PetscErrorCode PetscDSSetContext(PetscDS ds, PetscInt f, void *ctx)
2233: {
2234:   PetscFunctionBegin;
2236:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2237:   PetscCall(PetscDSEnlarge_Static(ds, f + 1));
2238:   ds->ctx[f] = ctx;
2239:   PetscFunctionReturn(PETSC_SUCCESS);
2240: }

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

2245:   Not Collective

2247:   Input Parameters:
2248: + ds - The PetscDS
2249: - f  - The test field number

2251:   Output Parameters:
2252: + f0 - boundary integrand for the test function term
2253: - f1 - boundary integrand for the test function gradient term

2255:   Calling sequence of `f0`:
2256: + dim          - the spatial dimension
2257: . Nf           - the number of fields
2258: . NfAux        - the number of auxiliary fields
2259: . uOff         - the offset into u[] and u_t[] for each field
2260: . uOff_x       - the offset into u_x[] for each field
2261: . u            - each field evaluated at the current point
2262: . u_t          - the time derivative of each field evaluated at the current point
2263: . u_x          - the gradient of each field evaluated at the current point
2264: . aOff         - the offset into a[] and a_t[] for each auxiliary field
2265: . aOff_x       - the offset into a_x[] for each auxiliary field
2266: . a            - each auxiliary field evaluated at the current point
2267: . a_t          - the time derivative of each auxiliary field evaluated at the current point
2268: . a_x          - the gradient of auxiliary each field evaluated at the current point
2269: . t            - current time
2270: . x            - coordinates of the current point
2271: . n            - unit normal at the current point
2272: . numConstants - number of constant parameters
2273: . constants    - constant parameters
2274: - f0           - output values at the current point

2276:   Level: intermediate

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

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

2284: .seealso: `PetscDS`, `PetscDSSetBdResidual()`
2285: @*/
2286: 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[]))
2287: {
2288:   PetscBdPointFunc *tmp0, *tmp1;
2289:   PetscInt          n0, n1;

2291:   PetscFunctionBegin;
2293:   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);
2294:   PetscCall(PetscWeakFormGetBdResidual(ds->wf, NULL, 0, f, 0, &n0, &tmp0, &n1, &tmp1));
2295:   *f0 = tmp0 ? tmp0[0] : NULL;
2296:   *f1 = tmp1 ? tmp1[0] : NULL;
2297:   PetscFunctionReturn(PETSC_SUCCESS);
2298: }

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

2303:   Not Collective

2305:   Input Parameters:
2306: + ds - The `PetscDS`
2307: . f  - The test field number
2308: . f0 - boundary integrand for the test function term
2309: - f1 - boundary integrand for the test function gradient term

2311:   Calling sequence of `f0`:
2312: + dim          - the spatial dimension
2313: . Nf           - the number of fields
2314: . NfAux        - the number of auxiliary fields
2315: . uOff         - the offset into u[] and u_t[] for each field
2316: . uOff_x       - the offset into u_x[] for each field
2317: . u            - each field evaluated at the current point
2318: . u_t          - the time derivative of each field evaluated at the current point
2319: . u_x          - the gradient of each field evaluated at the current point
2320: . aOff         - the offset into a[] and a_t[] for each auxiliary field
2321: . aOff_x       - the offset into a_x[] for each auxiliary field
2322: . a            - each auxiliary field evaluated at the current point
2323: . a_t          - the time derivative of each auxiliary field evaluated at the current point
2324: . a_x          - the gradient of auxiliary each field evaluated at the current point
2325: . t            - current time
2326: . x            - coordinates of the current point
2327: . n            - unit normal at the current point
2328: . numConstants - number of constant parameters
2329: . constants    - constant parameters
2330: - f0           - output values at the current point

2332:   Level: intermediate

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

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

2340: .seealso: `PetscDS`, `PetscDSGetBdResidual()`
2341: @*/
2342: 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[]))
2343: {
2344:   PetscFunctionBegin;
2346:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2347:   PetscCall(PetscWeakFormSetIndexBdResidual(ds->wf, NULL, 0, f, 0, 0, f0, 0, f1));
2348:   PetscFunctionReturn(PETSC_SUCCESS);
2349: }

2351: /*@
2352:   PetscDSHasBdJacobian - Indicates that boundary Jacobian functions have been set

2354:   Not Collective

2356:   Input Parameter:
2357: . ds - The `PetscDS`

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

2362:   Level: intermediate

2364: .seealso: `PetscDS`, `PetscDSHasJacobian()`, `PetscDSSetBdJacobian()`, `PetscDSGetBdJacobian()`
2365: @*/
2366: PetscErrorCode PetscDSHasBdJacobian(PetscDS ds, PetscBool *hasBdJac)
2367: {
2368:   PetscFunctionBegin;
2370:   PetscAssertPointer(hasBdJac, 2);
2371:   PetscCall(PetscWeakFormHasBdJacobian(ds->wf, hasBdJac));
2372:   PetscFunctionReturn(PETSC_SUCCESS);
2373: }

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

2378:   Not Collective

2380:   Input Parameters:
2381: + ds - The `PetscDS`
2382: . f  - The test field number
2383: - g  - The field number

2385:   Output Parameters:
2386: + g0 - integrand for the test and basis function term
2387: . g1 - integrand for the test function and basis function gradient term
2388: . g2 - integrand for the test function gradient and basis function term
2389: - g3 - integrand for the test function gradient and basis function gradient term

2391:   Calling sequence of `g0`:
2392: + dim          - the spatial dimension
2393: . Nf           - the number of fields
2394: . NfAux        - the number of auxiliary fields
2395: . uOff         - the offset into u[] and u_t[] for each field
2396: . uOff_x       - the offset into u_x[] for each field
2397: . u            - each field evaluated at the current point
2398: . u_t          - the time derivative of each field evaluated at the current point
2399: . u_x          - the gradient of each field evaluated at the current point
2400: . aOff         - the offset into a[] and a_t[] for each auxiliary field
2401: . aOff_x       - the offset into a_x[] for each auxiliary field
2402: . a            - each auxiliary field evaluated at the current point
2403: . a_t          - the time derivative of each auxiliary field evaluated at the current point
2404: . a_x          - the gradient of auxiliary each field evaluated at the current point
2405: . t            - current time
2406: . u_tShift     - the multiplier a for dF/dU_t
2407: . x            - coordinates of the current point
2408: . n            - normal at the current point
2409: . numConstants - number of constant parameters
2410: . constants    - constant parameters
2411: - g0           - output values at the current point

2413:   Level: intermediate

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

2418:   We are using a first order FEM model for the weak form\:
2419:   \int_\Gamma \phi {\vec g}_0(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \cdot \hat n \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \hat n \cdot \nabla \psi

2421: .seealso: `PetscDS`, `PetscDSSetBdJacobian()`
2422: @*/
2423: 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[]))
2424: {
2425:   PetscBdPointJac *tmp0, *tmp1, *tmp2, *tmp3;
2426:   PetscInt         n0, n1, n2, n3;

2428:   PetscFunctionBegin;
2430:   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);
2431:   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);
2432:   PetscCall(PetscWeakFormGetBdJacobian(ds->wf, NULL, 0, f, g, 0, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3));
2433:   *g0 = tmp0 ? tmp0[0] : NULL;
2434:   *g1 = tmp1 ? tmp1[0] : NULL;
2435:   *g2 = tmp2 ? tmp2[0] : NULL;
2436:   *g3 = tmp3 ? tmp3[0] : NULL;
2437:   PetscFunctionReturn(PETSC_SUCCESS);
2438: }

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

2443:   Not Collective

2445:   Input Parameters:
2446: + ds - The PetscDS
2447: . f  - The test field number
2448: . g  - The field number
2449: . g0 - integrand for the test and basis function term
2450: . g1 - integrand for the test function and basis function gradient term
2451: . g2 - integrand for the test function gradient and basis function term
2452: - g3 - integrand for the test function gradient and basis function gradient term

2454:   Calling sequence of `g0`:
2455: + dim          - the spatial dimension
2456: . Nf           - the number of fields
2457: . NfAux        - the number of auxiliary fields
2458: . uOff         - the offset into u[] and u_t[] for each field
2459: . uOff_x       - the offset into u_x[] for each field
2460: . u            - each field evaluated at the current point
2461: . u_t          - the time derivative of each field evaluated at the current point
2462: . u_x          - the gradient of each field evaluated at the current point
2463: . aOff         - the offset into a[] and a_t[] for each auxiliary field
2464: . aOff_x       - the offset into a_x[] for each auxiliary field
2465: . a            - each auxiliary field evaluated at the current point
2466: . a_t          - the time derivative of each auxiliary field evaluated at the current point
2467: . a_x          - the gradient of auxiliary each field evaluated at the current point
2468: . t            - current time
2469: . u_tShift     - the multiplier a for dF/dU_t
2470: . x            - coordinates of the current point
2471: . n            - normal at the current point
2472: . numConstants - number of constant parameters
2473: . constants    - constant parameters
2474: - g0           - output values at the current point

2476:   Level: intermediate

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

2481:   We are using a first order FEM model for the weak form\:
2482:   \int_\Gamma \phi {\vec g}_0(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \cdot \hat n \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \hat n \cdot \nabla \psi

2484: .seealso: `PetscDS`, `PetscDSGetBdJacobian()`
2485: @*/
2486: 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[]))
2487: {
2488:   PetscFunctionBegin;
2494:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2495:   PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", g);
2496:   PetscCall(PetscWeakFormSetIndexBdJacobian(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3));
2497:   PetscFunctionReturn(PETSC_SUCCESS);
2498: }

2500: /*@
2501:   PetscDSHasBdJacobianPreconditioner - Signals that boundary Jacobian preconditioner functions have been set

2503:   Not Collective

2505:   Input Parameter:
2506: . ds - The `PetscDS`

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

2511:   Level: intermediate

2513: .seealso: `PetscDS`, `PetscDSHasJacobian()`, `PetscDSSetBdJacobian()`, `PetscDSGetBdJacobian()`
2514: @*/
2515: PetscErrorCode PetscDSHasBdJacobianPreconditioner(PetscDS ds, PetscBool *hasBdJacPre)
2516: {
2517:   PetscFunctionBegin;
2519:   PetscAssertPointer(hasBdJacPre, 2);
2520:   PetscCall(PetscWeakFormHasBdJacobianPreconditioner(ds->wf, hasBdJacPre));
2521:   PetscFunctionReturn(PETSC_SUCCESS);
2522: }

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

2527:   Not Collective; No Fortran Support

2529:   Input Parameters:
2530: + ds - The `PetscDS`
2531: . f  - The test field number
2532: - g  - The field number

2534:   Output Parameters:
2535: + g0 - integrand for the test and basis function term
2536: . g1 - integrand for the test function and basis function gradient term
2537: . g2 - integrand for the test function gradient and basis function term
2538: - g3 - integrand for the test function gradient and basis function gradient term

2540:   Calling sequence of `g0`:
2541: + dim          - the spatial dimension
2542: . Nf           - the number of fields
2543: . NfAux        - the number of auxiliary fields
2544: . uOff         - the offset into u[] and u_t[] for each field
2545: . uOff_x       - the offset into u_x[] for each field
2546: . u            - each field evaluated at the current point
2547: . u_t          - the time derivative of each field evaluated at the current point
2548: . u_x          - the gradient of each field evaluated at the current point
2549: . aOff         - the offset into a[] and a_t[] for each auxiliary field
2550: . aOff_x       - the offset into a_x[] for each auxiliary field
2551: . a            - each auxiliary field evaluated at the current point
2552: . a_t          - the time derivative of each auxiliary field evaluated at the current point
2553: . a_x          - the gradient of auxiliary each field evaluated at the current point
2554: . t            - current time
2555: . u_tShift     - the multiplier a for dF/dU_t
2556: . x            - coordinates of the current point
2557: . n            - normal at the current point
2558: . numConstants - number of constant parameters
2559: . constants    - constant parameters
2560: - g0           - output values at the current point

2562:   Level: intermediate

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

2567:   We are using a first order FEM model for the weak form\:
2568:   \int_\Gamma \phi {\vec g}_0(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \cdot \hat n \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \hat n \cdot \nabla \psi

2570: .seealso: `PetscDS`, `PetscDSSetBdJacobianPreconditioner()`
2571: @*/
2572: 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[]))
2573: {
2574:   PetscBdPointJac *tmp0, *tmp1, *tmp2, *tmp3;
2575:   PetscInt         n0, n1, n2, n3;

2577:   PetscFunctionBegin;
2579:   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);
2580:   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);
2581:   PetscCall(PetscWeakFormGetBdJacobianPreconditioner(ds->wf, NULL, 0, f, g, 0, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3));
2582:   *g0 = tmp0 ? tmp0[0] : NULL;
2583:   *g1 = tmp1 ? tmp1[0] : NULL;
2584:   *g2 = tmp2 ? tmp2[0] : NULL;
2585:   *g3 = tmp3 ? tmp3[0] : NULL;
2586:   PetscFunctionReturn(PETSC_SUCCESS);
2587: }

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

2592:   Not Collective; No Fortran Support

2594:   Input Parameters:
2595: + ds - The `PetscDS`
2596: . f  - The test field number
2597: . g  - The field number
2598: . g0 - integrand for the test and basis function term
2599: . g1 - integrand for the test function and basis function gradient term
2600: . g2 - integrand for the test function gradient and basis function term
2601: - g3 - integrand for the test function gradient and basis function gradient term

2603:   Calling sequence of `g0':
2604: + dim          - the spatial dimension
2605: . Nf           - the number of fields
2606: . NfAux        - the number of auxiliary fields
2607: . uOff         - the offset into u[] and u_t[] for each field
2608: . uOff_x       - the offset into u_x[] for each field
2609: . u            - each field evaluated at the current point
2610: . u_t          - the time derivative of each field evaluated at the current point
2611: . u_x          - the gradient of each field evaluated at the current point
2612: . aOff         - the offset into a[] and a_t[] for each auxiliary field
2613: . aOff_x       - the offset into a_x[] for each auxiliary field
2614: . a            - each auxiliary field evaluated at the current point
2615: . a_t          - the time derivative of each auxiliary field evaluated at the current point
2616: . a_x          - the gradient of auxiliary each field evaluated at the current point
2617: . t            - current time
2618: . u_tShift     - the multiplier a for dF/dU_t
2619: . x            - coordinates of the current point
2620: . n            - normal at the current point
2621: . numConstants - number of constant parameters
2622: . constants    - constant parameters
2623: - g0           - output values at the current point

2625:   Level: intermediate

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

2630:   We are using a first order FEM model for the weak form\:
2631:   \int_\Gamma \phi {\vec g}_0(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \cdot \hat n \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \hat n \cdot \nabla \psi

2633: .seealso: `PetscDS`, `PetscDSGetBdJacobianPreconditioner()`
2634: @*/
2635: 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[]))
2636: {
2637:   PetscFunctionBegin;
2643:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2644:   PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", g);
2645:   PetscCall(PetscWeakFormSetIndexBdJacobianPreconditioner(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3));
2646:   PetscFunctionReturn(PETSC_SUCCESS);
2647: }

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

2652:   Not Collective

2654:   Input Parameters:
2655: + prob - The PetscDS
2656: - f    - The test field number

2658:   Output Parameters:
2659: + sol - exact solution for the test field
2660: - ctx - exact solution context

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

2670:   Level: intermediate

2672: .seealso: `PetscDS`, `PetscDSSetExactSolution()`, `PetscDSGetExactSolutionTimeDerivative()`
2673: @*/
2674: PetscErrorCode PetscDSGetExactSolution(PetscDS prob, PetscInt f, PetscErrorCode (**sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void **ctx)
2675: {
2676:   PetscFunctionBegin;
2678:   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);
2679:   if (sol) {
2680:     PetscAssertPointer(sol, 3);
2681:     *sol = prob->exactSol[f];
2682:   }
2683:   if (ctx) {
2684:     PetscAssertPointer(ctx, 4);
2685:     *ctx = prob->exactCtx[f];
2686:   }
2687:   PetscFunctionReturn(PETSC_SUCCESS);
2688: }

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

2693:   Not Collective

2695:   Input Parameters:
2696: + prob - The `PetscDS`
2697: . f    - The test field number
2698: . sol  - solution function for the test fields
2699: - ctx  - solution context or `NULL`

2701:   Calling sequence of `sol`:
2702: + dim - the spatial dimension
2703: . t   - current time
2704: . x   - coordinates of the current point
2705: . Nc  - the number of field components
2706: . u   - the solution field evaluated at the current point
2707: - ctx - a user context

2709:   Level: intermediate

2711: .seealso: `PetscDS`, `PetscDSGetExactSolution()`
2712: @*/
2713: PetscErrorCode PetscDSSetExactSolution(PetscDS prob, PetscInt f, PetscErrorCode (*sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void *ctx)
2714: {
2715:   PetscFunctionBegin;
2717:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2718:   PetscCall(PetscDSEnlarge_Static(prob, f + 1));
2719:   if (sol) {
2721:     prob->exactSol[f] = sol;
2722:   }
2723:   if (ctx) {
2725:     prob->exactCtx[f] = ctx;
2726:   }
2727:   PetscFunctionReturn(PETSC_SUCCESS);
2728: }

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

2733:   Not Collective

2735:   Input Parameters:
2736: + prob - The `PetscDS`
2737: - f    - The test field number

2739:   Output Parameters:
2740: + sol - time derivative of the exact solution for the test field
2741: - ctx - time derivative of the exact solution context

2743:   Calling sequence of `exactSol`:
2744: + dim - the spatial dimension
2745: . t   - current time
2746: . x   - coordinates of the current point
2747: . Nc  - the number of field components
2748: . u   - the solution field evaluated at the current point
2749: - ctx - a user context

2751:   Level: intermediate

2753: .seealso: `PetscDS`, `PetscDSSetExactSolutionTimeDerivative()`, `PetscDSGetExactSolution()`
2754: @*/
2755: PetscErrorCode PetscDSGetExactSolutionTimeDerivative(PetscDS prob, PetscInt f, PetscErrorCode (**sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void **ctx)
2756: {
2757:   PetscFunctionBegin;
2759:   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);
2760:   if (sol) {
2761:     PetscAssertPointer(sol, 3);
2762:     *sol = prob->exactSol_t[f];
2763:   }
2764:   if (ctx) {
2765:     PetscAssertPointer(ctx, 4);
2766:     *ctx = prob->exactCtx_t[f];
2767:   }
2768:   PetscFunctionReturn(PETSC_SUCCESS);
2769: }

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

2774:   Not Collective

2776:   Input Parameters:
2777: + prob - The `PetscDS`
2778: . f    - The test field number
2779: . sol  - time derivative of the solution function for the test fields
2780: - ctx  - time derivative of the solution context or `NULL`

2782:   Calling sequence of `sol`:
2783: + dim - the spatial dimension
2784: . t   - current time
2785: . x   - coordinates of the current point
2786: . Nc  - the number of field components
2787: . u   - the solution field evaluated at the current point
2788: - ctx - a user context

2790:   Level: intermediate

2792: .seealso: `PetscDS`, `PetscDSGetExactSolutionTimeDerivative()`, `PetscDSSetExactSolution()`
2793: @*/
2794: PetscErrorCode PetscDSSetExactSolutionTimeDerivative(PetscDS prob, PetscInt f, PetscErrorCode (*sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void *ctx)
2795: {
2796:   PetscFunctionBegin;
2798:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2799:   PetscCall(PetscDSEnlarge_Static(prob, f + 1));
2800:   if (sol) {
2802:     prob->exactSol_t[f] = sol;
2803:   }
2804:   if (ctx) {
2806:     prob->exactCtx_t[f] = ctx;
2807:   }
2808:   PetscFunctionReturn(PETSC_SUCCESS);
2809: }

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

2814:   Not Collective

2816:   Input Parameter:
2817: . ds - The `PetscDS` object

2819:   Output Parameters:
2820: + numConstants - The number of constants
2821: - constants    - The array of constants, NULL if there are none

2823:   Level: intermediate

2825: .seealso: `PetscDS`, `PetscDSSetConstants()`, `PetscDSCreate()`
2826: @*/
2827: PetscErrorCode PetscDSGetConstants(PetscDS ds, PetscInt *numConstants, const PetscScalar *constants[])
2828: {
2829:   PetscFunctionBegin;
2831:   if (numConstants) {
2832:     PetscAssertPointer(numConstants, 2);
2833:     *numConstants = ds->numConstants;
2834:   }
2835:   if (constants) {
2836:     PetscAssertPointer(constants, 3);
2837:     *constants = ds->constants;
2838:   }
2839:   PetscFunctionReturn(PETSC_SUCCESS);
2840: }

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

2845:   Not Collective

2847:   Input Parameters:
2848: + ds           - The `PetscDS` object
2849: . numConstants - The number of constants
2850: - constants    - The array of constants, `NULL` if there are none

2852:   Level: intermediate

2854: .seealso: `PetscDS`, `PetscDSGetConstants()`, `PetscDSCreate()`
2855: @*/
2856: PetscErrorCode PetscDSSetConstants(PetscDS ds, PetscInt numConstants, PetscScalar constants[])
2857: {
2858:   PetscFunctionBegin;
2860:   if (numConstants != ds->numConstants) {
2861:     PetscCall(PetscFree(ds->constants));
2862:     ds->numConstants = numConstants;
2863:     PetscCall(PetscMalloc1(ds->numConstants + ds->numFuncConstants, &ds->constants));
2864:   }
2865:   if (ds->numConstants) {
2866:     PetscAssertPointer(constants, 3);
2867:     PetscCall(PetscArraycpy(ds->constants, constants, ds->numConstants));
2868:   }
2869:   PetscFunctionReturn(PETSC_SUCCESS);
2870: }

2872: /*@C
2873:   PetscDSSetIntegrationParameters - Set the parameters for a particular integration

2875:   Not Collective

2877:   Input Parameters:
2878: + ds     - The `PetscDS` object
2879: . fieldI - The test field for a given point function, or PETSC_DETERMINE
2880: - fieldJ - The basis field for a given point function, or PETSC_DETERMINE

2882:   Level: intermediate

2884: .seealso: `PetscDS`, `PetscDSSetConstants()`, `PetscDSGetConstants()`, `PetscDSCreate()`
2885: @*/
2886: PetscErrorCode PetscDSSetIntegrationParameters(PetscDS ds, PetscInt fieldI, PetscInt fieldJ)
2887: {
2888:   PetscFunctionBegin;
2890:   ds->constants[ds->numConstants]     = fieldI;
2891:   ds->constants[ds->numConstants + 1] = fieldJ;
2892:   PetscFunctionReturn(PETSC_SUCCESS);
2893: }

2895: /*@C
2896:   PetscDSSetCellParameters - Set the parameters for a particular cell

2898:   Not Collective

2900:   Input Parameters:
2901: + ds     - The `PetscDS` object
2902: - volume - The cell volume

2904:   Level: intermediate

2906: .seealso: `PetscDS`, `PetscDSSetConstants()`, `PetscDSGetConstants()`, `PetscDSCreate()`
2907: @*/
2908: PetscErrorCode PetscDSSetCellParameters(PetscDS ds, PetscReal volume)
2909: {
2910:   PetscFunctionBegin;
2912:   ds->constants[ds->numConstants + 2] = volume;
2913:   PetscFunctionReturn(PETSC_SUCCESS);
2914: }

2916: /*@
2917:   PetscDSGetFieldIndex - Returns the index of the given field

2919:   Not Collective

2921:   Input Parameters:
2922: + prob - The `PetscDS` object
2923: - disc - The discretization object

2925:   Output Parameter:
2926: . f - The field number

2928:   Level: beginner

2930: .seealso: `PetscDS`, `PetscGetDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
2931: @*/
2932: PetscErrorCode PetscDSGetFieldIndex(PetscDS prob, PetscObject disc, PetscInt *f)
2933: {
2934:   PetscInt g;

2936:   PetscFunctionBegin;
2938:   PetscAssertPointer(f, 3);
2939:   *f = -1;
2940:   for (g = 0; g < prob->Nf; ++g) {
2941:     if (disc == prob->disc[g]) break;
2942:   }
2943:   PetscCheck(g != prob->Nf, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONG, "Field not found in PetscDS.");
2944:   *f = g;
2945:   PetscFunctionReturn(PETSC_SUCCESS);
2946: }

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

2951:   Not Collective

2953:   Input Parameters:
2954: + prob - The `PetscDS` object
2955: - f    - The field number

2957:   Output Parameter:
2958: . size - The size

2960:   Level: beginner

2962: .seealso: `PetscDS`, `PetscDSGetFieldOffset()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
2963: @*/
2964: PetscErrorCode PetscDSGetFieldSize(PetscDS prob, PetscInt f, PetscInt *size)
2965: {
2966:   PetscFunctionBegin;
2968:   PetscAssertPointer(size, 3);
2969:   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);
2970:   PetscCall(PetscDSSetUp(prob));
2971:   *size = prob->Nb[f];
2972:   PetscFunctionReturn(PETSC_SUCCESS);
2973: }

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

2978:   Not Collective

2980:   Input Parameters:
2981: + prob - The `PetscDS` object
2982: - f    - The field number

2984:   Output Parameter:
2985: . off - The offset

2987:   Level: beginner

2989: .seealso: `PetscDS`, `PetscDSGetFieldSize()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
2990: @*/
2991: PetscErrorCode PetscDSGetFieldOffset(PetscDS prob, PetscInt f, PetscInt *off)
2992: {
2993:   PetscInt size, g;

2995:   PetscFunctionBegin;
2997:   PetscAssertPointer(off, 3);
2998:   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);
2999:   *off = 0;
3000:   for (g = 0; g < f; ++g) {
3001:     PetscCall(PetscDSGetFieldSize(prob, g, &size));
3002:     *off += size;
3003:   }
3004:   PetscFunctionReturn(PETSC_SUCCESS);
3005: }

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

3010:   Not Collective

3012:   Input Parameters:
3013: + ds - The `PetscDS` object
3014: - f  - The field number

3016:   Output Parameter:
3017: . off - The offset

3019:   Level: beginner

3021: .seealso: `PetscDS`, `PetscDSGetFieldSize()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3022: @*/
3023: PetscErrorCode PetscDSGetFieldOffsetCohesive(PetscDS ds, PetscInt f, PetscInt *off)
3024: {
3025:   PetscInt size, g;

3027:   PetscFunctionBegin;
3029:   PetscAssertPointer(off, 3);
3030:   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);
3031:   *off = 0;
3032:   for (g = 0; g < f; ++g) {
3033:     PetscBool cohesive;

3035:     PetscCall(PetscDSGetCohesive(ds, g, &cohesive));
3036:     PetscCall(PetscDSGetFieldSize(ds, g, &size));
3037:     *off += cohesive ? size : size * 2;
3038:   }
3039:   PetscFunctionReturn(PETSC_SUCCESS);
3040: }

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

3045:   Not Collective

3047:   Input Parameter:
3048: . prob - The `PetscDS` object

3050:   Output Parameter:
3051: . dimensions - The number of dimensions

3053:   Level: beginner

3055: .seealso: `PetscDS`, `PetscDSGetComponentOffsets()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3056: @*/
3057: PetscErrorCode PetscDSGetDimensions(PetscDS prob, PetscInt *dimensions[])
3058: {
3059:   PetscFunctionBegin;
3061:   PetscCall(PetscDSSetUp(prob));
3062:   PetscAssertPointer(dimensions, 2);
3063:   *dimensions = prob->Nb;
3064:   PetscFunctionReturn(PETSC_SUCCESS);
3065: }

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

3070:   Not Collective

3072:   Input Parameter:
3073: . prob - The `PetscDS` object

3075:   Output Parameter:
3076: . components - The number of components

3078:   Level: beginner

3080: .seealso: `PetscDS`, `PetscDSGetComponentOffsets()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3081: @*/
3082: PetscErrorCode PetscDSGetComponents(PetscDS prob, PetscInt *components[])
3083: {
3084:   PetscFunctionBegin;
3086:   PetscCall(PetscDSSetUp(prob));
3087:   PetscAssertPointer(components, 2);
3088:   *components = prob->Nc;
3089:   PetscFunctionReturn(PETSC_SUCCESS);
3090: }

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

3095:   Not Collective

3097:   Input Parameters:
3098: + prob - The `PetscDS` object
3099: - f    - The field number

3101:   Output Parameter:
3102: . off - The offset

3104:   Level: beginner

3106: .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3107: @*/
3108: PetscErrorCode PetscDSGetComponentOffset(PetscDS prob, PetscInt f, PetscInt *off)
3109: {
3110:   PetscFunctionBegin;
3112:   PetscAssertPointer(off, 3);
3113:   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);
3114:   PetscCall(PetscDSSetUp(prob));
3115:   *off = prob->off[f];
3116:   PetscFunctionReturn(PETSC_SUCCESS);
3117: }

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

3122:   Not Collective

3124:   Input Parameter:
3125: . prob - The `PetscDS` object

3127:   Output Parameter:
3128: . offsets - The offsets

3130:   Level: beginner

3132: .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3133: @*/
3134: PetscErrorCode PetscDSGetComponentOffsets(PetscDS prob, PetscInt *offsets[])
3135: {
3136:   PetscFunctionBegin;
3138:   PetscAssertPointer(offsets, 2);
3139:   PetscCall(PetscDSSetUp(prob));
3140:   *offsets = prob->off;
3141:   PetscFunctionReturn(PETSC_SUCCESS);
3142: }

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

3147:   Not Collective

3149:   Input Parameter:
3150: . prob - The `PetscDS` object

3152:   Output Parameter:
3153: . offsets - The offsets

3155:   Level: beginner

3157: .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3158: @*/
3159: PetscErrorCode PetscDSGetComponentDerivativeOffsets(PetscDS prob, PetscInt *offsets[])
3160: {
3161:   PetscFunctionBegin;
3163:   PetscAssertPointer(offsets, 2);
3164:   PetscCall(PetscDSSetUp(prob));
3165:   *offsets = prob->offDer;
3166:   PetscFunctionReturn(PETSC_SUCCESS);
3167: }

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

3172:   Not Collective

3174:   Input Parameters:
3175: + ds - The `PetscDS` object
3176: - s  - The cohesive side, 0 for negative, 1 for positive, 2 for cohesive

3178:   Output Parameter:
3179: . offsets - The offsets

3181:   Level: beginner

3183: .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3184: @*/
3185: PetscErrorCode PetscDSGetComponentOffsetsCohesive(PetscDS ds, PetscInt s, PetscInt *offsets[])
3186: {
3187:   PetscFunctionBegin;
3189:   PetscAssertPointer(offsets, 3);
3190:   PetscCheck(ds->isCohesive, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cohesive offsets are only valid for a cohesive DS");
3191:   PetscCheck(!(s < 0) && !(s > 2), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cohesive side %" PetscInt_FMT " is not in [0, 2]", s);
3192:   PetscCall(PetscDSSetUp(ds));
3193:   *offsets = ds->offCohesive[s];
3194:   PetscFunctionReturn(PETSC_SUCCESS);
3195: }

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

3200:   Not Collective

3202:   Input Parameters:
3203: + ds - The `PetscDS` object
3204: - s  - The cohesive side, 0 for negative, 1 for positive, 2 for cohesive

3206:   Output Parameter:
3207: . offsets - The offsets

3209:   Level: beginner

3211: .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3212: @*/
3213: PetscErrorCode PetscDSGetComponentDerivativeOffsetsCohesive(PetscDS ds, PetscInt s, PetscInt *offsets[])
3214: {
3215:   PetscFunctionBegin;
3217:   PetscAssertPointer(offsets, 3);
3218:   PetscCheck(ds->isCohesive, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cohesive offsets are only valid for a cohesive DS");
3219:   PetscCheck(!(s < 0) && !(s > 2), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cohesive side %" PetscInt_FMT " is not in [0, 2]", s);
3220:   PetscCall(PetscDSSetUp(ds));
3221:   *offsets = ds->offDerCohesive[s];
3222:   PetscFunctionReturn(PETSC_SUCCESS);
3223: }

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

3228:   Not Collective

3230:   Input Parameter:
3231: . prob - The `PetscDS` object

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

3236:   Level: intermediate

3238: .seealso: `PetscDS`, `PetscTabulation`, `PetscDSCreate()`
3239: @*/
3240: PetscErrorCode PetscDSGetTabulation(PetscDS prob, PetscTabulation *T[])
3241: {
3242:   PetscFunctionBegin;
3244:   PetscAssertPointer(T, 2);
3245:   PetscCall(PetscDSSetUp(prob));
3246:   *T = prob->T;
3247:   PetscFunctionReturn(PETSC_SUCCESS);
3248: }

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

3253:   Not Collective

3255:   Input Parameter:
3256: . prob - The `PetscDS` object

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

3261:   Level: intermediate

3263: .seealso: `PetscTabulation`, `PetscDS`, `PetscDSGetTabulation()`, `PetscDSCreate()`
3264: @*/
3265: PetscErrorCode PetscDSGetFaceTabulation(PetscDS prob, PetscTabulation *Tf[])
3266: {
3267:   PetscFunctionBegin;
3269:   PetscAssertPointer(Tf, 2);
3270:   PetscCall(PetscDSSetUp(prob));
3271:   *Tf = prob->Tf;
3272:   PetscFunctionReturn(PETSC_SUCCESS);
3273: }

3275: PetscErrorCode PetscDSGetEvaluationArrays(PetscDS prob, PetscScalar **u, PetscScalar **u_t, PetscScalar **u_x)
3276: {
3277:   PetscFunctionBegin;
3279:   PetscCall(PetscDSSetUp(prob));
3280:   if (u) {
3281:     PetscAssertPointer(u, 2);
3282:     *u = prob->u;
3283:   }
3284:   if (u_t) {
3285:     PetscAssertPointer(u_t, 3);
3286:     *u_t = prob->u_t;
3287:   }
3288:   if (u_x) {
3289:     PetscAssertPointer(u_x, 4);
3290:     *u_x = prob->u_x;
3291:   }
3292:   PetscFunctionReturn(PETSC_SUCCESS);
3293: }

3295: PetscErrorCode PetscDSGetWeakFormArrays(PetscDS prob, PetscScalar **f0, PetscScalar **f1, PetscScalar **g0, PetscScalar **g1, PetscScalar **g2, PetscScalar **g3)
3296: {
3297:   PetscFunctionBegin;
3299:   PetscCall(PetscDSSetUp(prob));
3300:   if (f0) {
3301:     PetscAssertPointer(f0, 2);
3302:     *f0 = prob->f0;
3303:   }
3304:   if (f1) {
3305:     PetscAssertPointer(f1, 3);
3306:     *f1 = prob->f1;
3307:   }
3308:   if (g0) {
3309:     PetscAssertPointer(g0, 4);
3310:     *g0 = prob->g0;
3311:   }
3312:   if (g1) {
3313:     PetscAssertPointer(g1, 5);
3314:     *g1 = prob->g1;
3315:   }
3316:   if (g2) {
3317:     PetscAssertPointer(g2, 6);
3318:     *g2 = prob->g2;
3319:   }
3320:   if (g3) {
3321:     PetscAssertPointer(g3, 7);
3322:     *g3 = prob->g3;
3323:   }
3324:   PetscFunctionReturn(PETSC_SUCCESS);
3325: }

3327: PetscErrorCode PetscDSGetWorkspace(PetscDS prob, PetscReal **x, PetscScalar **basisReal, PetscScalar **basisDerReal, PetscScalar **testReal, PetscScalar **testDerReal)
3328: {
3329:   PetscFunctionBegin;
3331:   PetscCall(PetscDSSetUp(prob));
3332:   if (x) {
3333:     PetscAssertPointer(x, 2);
3334:     *x = prob->x;
3335:   }
3336:   if (basisReal) {
3337:     PetscAssertPointer(basisReal, 3);
3338:     *basisReal = prob->basisReal;
3339:   }
3340:   if (basisDerReal) {
3341:     PetscAssertPointer(basisDerReal, 4);
3342:     *basisDerReal = prob->basisDerReal;
3343:   }
3344:   if (testReal) {
3345:     PetscAssertPointer(testReal, 5);
3346:     *testReal = prob->testReal;
3347:   }
3348:   if (testDerReal) {
3349:     PetscAssertPointer(testDerReal, 6);
3350:     *testDerReal = prob->testDerReal;
3351:   }
3352:   PetscFunctionReturn(PETSC_SUCCESS);
3353: }

3355: /*@C
3356:   PetscDSAddBoundary - Add a boundary condition to the model.

3358:   Collective

3360:   Input Parameters:
3361: + ds       - The PetscDS object
3362: . type     - The type of condition, e.g. `DM_BC_ESSENTIAL`/`DM_BC_ESSENTIAL_FIELD` (Dirichlet), or `DM_BC_NATURAL` (Neumann)
3363: . name     - The BC name
3364: . label    - The label defining constrained points
3365: . Nv       - The number of `DMLabel` values for constrained points
3366: . values   - An array of label values for constrained points
3367: . field    - The field to constrain
3368: . Nc       - The number of constrained field components (0 will constrain all fields)
3369: . comps    - An array of constrained component numbers
3370: . bcFunc   - A pointwise function giving boundary values
3371: . bcFunc_t - A pointwise function giving the time derivative of the boundary values, or NULL
3372: - ctx      - An optional user context for bcFunc

3374:   Output Parameter:
3375: . bd - The boundary number

3377:   Options Database Keys:
3378: + -bc_<boundary name> <num>      - Overrides the boundary ids
3379: - -bc_<boundary name>_comp <num> - Overrides the boundary components

3381:   Level: developer

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

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

3388:   If the type is `DM_BC_ESSENTIAL_FIELD` or other _FIELD value, then the calling sequence is\:
3389: .vb
3390:   void bcFunc(PetscInt dim, PetscInt Nf, PetscInt NfAux,
3391:               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
3392:               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
3393:               PetscReal time, const PetscReal x[], PetscScalar bcval[])
3394: .ve
3395: + dim - the spatial dimension
3396: . Nf - the number of fields
3397: . uOff - the offset into u[] and u_t[] for each field
3398: . uOff_x - the offset into u_x[] for each field
3399: . u - each field evaluated at the current point
3400: . u_t - the time derivative of each field evaluated at the current point
3401: . u_x - the gradient of each field evaluated at the current point
3402: . aOff - the offset into a[] and a_t[] for each auxiliary field
3403: . aOff_x - the offset into a_x[] for each auxiliary field
3404: . a - each auxiliary field evaluated at the current point
3405: . a_t - the time derivative of each auxiliary field evaluated at the current point
3406: . a_x - the gradient of auxiliary each field evaluated at the current point
3407: . t - current time
3408: . x - coordinates of the current point
3409: . numConstants - number of constant parameters
3410: . constants - constant parameters
3411: - bcval - output values at the current point

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

3419: .seealso: `PetscDS`, `PetscWeakForm`, `DMLabel`, `DMBoundaryConditionType`, `PetscDSAddBoundaryByName()`, `PetscDSGetBoundary()`, `PetscDSSetResidual()`, `PetscDSSetBdResidual()`
3420: @*/
3421: 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)
3422: {
3423:   DSBoundary  head = ds->boundary, b;
3424:   PetscInt    n    = 0;
3425:   const char *lname;

3427:   PetscFunctionBegin;
3430:   PetscAssertPointer(name, 3);
3435:   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);
3436:   if (Nc > 0) {
3437:     PetscInt *fcomps;
3438:     PetscInt  c;

3440:     PetscCall(PetscDSGetComponents(ds, &fcomps));
3441:     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);
3442:     for (c = 0; c < Nc; ++c) {
3443:       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);
3444:     }
3445:   }
3446:   PetscCall(PetscNew(&b));
3447:   PetscCall(PetscStrallocpy(name, (char **)&b->name));
3448:   PetscCall(PetscWeakFormCreate(PETSC_COMM_SELF, &b->wf));
3449:   PetscCall(PetscWeakFormSetNumFields(b->wf, ds->Nf));
3450:   PetscCall(PetscMalloc1(Nv, &b->values));
3451:   if (Nv) PetscCall(PetscArraycpy(b->values, values, Nv));
3452:   PetscCall(PetscMalloc1(Nc, &b->comps));
3453:   if (Nc) PetscCall(PetscArraycpy(b->comps, comps, Nc));
3454:   PetscCall(PetscObjectGetName((PetscObject)label, &lname));
3455:   PetscCall(PetscStrallocpy(lname, (char **)&b->lname));
3456:   b->type   = type;
3457:   b->label  = label;
3458:   b->Nv     = Nv;
3459:   b->field  = field;
3460:   b->Nc     = Nc;
3461:   b->func   = bcFunc;
3462:   b->func_t = bcFunc_t;
3463:   b->ctx    = ctx;
3464:   b->next   = NULL;
3465:   /* Append to linked list so that we can preserve the order */
3466:   if (!head) ds->boundary = b;
3467:   while (head) {
3468:     if (!head->next) {
3469:       head->next = b;
3470:       head       = b;
3471:     }
3472:     head = head->next;
3473:     ++n;
3474:   }
3475:   if (bd) {
3476:     PetscAssertPointer(bd, 13);
3477:     *bd = n;
3478:   }
3479:   PetscFunctionReturn(PETSC_SUCCESS);
3480: }

3482: // PetscClangLinter pragma ignore: -fdoc-section-header-unknown
3483: /*@C
3484:   PetscDSAddBoundaryByName - Add a boundary condition to the model.

3486:   Collective

3488:   Input Parameters:
3489: + ds       - The `PetscDS` object
3490: . type     - The type of condition, e.g. `DM_BC_ESSENTIAL`/`DM_BC_ESSENTIAL_FIELD` (Dirichlet), or `DM_BC_NATURAL` (Neumann)
3491: . name     - The BC name
3492: . lname    - The naem of the label defining constrained points
3493: . Nv       - The number of `DMLabel` values for constrained points
3494: . values   - An array of label values for constrained points
3495: . field    - The field to constrain
3496: . Nc       - The number of constrained field components (0 will constrain all fields)
3497: . comps    - An array of constrained component numbers
3498: . bcFunc   - A pointwise function giving boundary values
3499: . bcFunc_t - A pointwise function giving the time derivative of the boundary values, or NULL
3500: - ctx      - An optional user context for bcFunc

3502:   Output Parameter:
3503: . bd - The boundary number

3505:   Options Database Keys:
3506: + -bc_<boundary name> <num>      - Overrides the boundary ids
3507: - -bc_<boundary name>_comp <num> - Overrides the boundary components

3509:   Calling Sequence of `bcFunc` and `bcFunc_t`:
3510:   If the type is `DM_BC_ESSENTIAL`
3511: .vb
3512:   void bcFunc(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar bcval[])
3513: .ve
3514:   If the type is `DM_BC_ESSENTIAL_FIELD` or other _FIELD value,
3515: .vb
3516:   void bcFunc(PetscInt dim, PetscInt Nf, PetscInt NfAux,
3517:               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
3518:               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
3519:               PetscReal time, const PetscReal x[], PetscScalar bcval[])
3520: .ve
3521: + dim - the spatial dimension
3522: . Nf - the number of fields
3523: . uOff - the offset into u[] and u_t[] for each field
3524: . uOff_x - the offset into u_x[] for each field
3525: . u - each field evaluated at the current point
3526: . u_t - the time derivative of each field evaluated at the current point
3527: . u_x - the gradient of each field evaluated at the current point
3528: . aOff - the offset into a[] and a_t[] for each auxiliary field
3529: . aOff_x - the offset into a_x[] for each auxiliary field
3530: . a - each auxiliary field evaluated at the current point
3531: . a_t - the time derivative of each auxiliary field evaluated at the current point
3532: . a_x - the gradient of auxiliary each field evaluated at the current point
3533: . t - current time
3534: . x - coordinates of the current point
3535: . numConstants - number of constant parameters
3536: . constants - constant parameters
3537: - bcval - output values at the current point

3539:   Level: developer

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

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

3549: .seealso: `PetscDS`, `PetscWeakForm`, `DMLabel`, `DMBoundaryConditionType`, `PetscDSAddBoundary()`, `PetscDSGetBoundary()`, `PetscDSSetResidual()`, `PetscDSSetBdResidual()`
3550: @*/
3551: 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)
3552: {
3553:   DSBoundary head = ds->boundary, b;
3554:   PetscInt   n    = 0;

3556:   PetscFunctionBegin;
3559:   PetscAssertPointer(name, 3);
3560:   PetscAssertPointer(lname, 4);
3564:   PetscCall(PetscNew(&b));
3565:   PetscCall(PetscStrallocpy(name, (char **)&b->name));
3566:   PetscCall(PetscWeakFormCreate(PETSC_COMM_SELF, &b->wf));
3567:   PetscCall(PetscWeakFormSetNumFields(b->wf, ds->Nf));
3568:   PetscCall(PetscMalloc1(Nv, &b->values));
3569:   if (Nv) PetscCall(PetscArraycpy(b->values, values, Nv));
3570:   PetscCall(PetscMalloc1(Nc, &b->comps));
3571:   if (Nc) PetscCall(PetscArraycpy(b->comps, comps, Nc));
3572:   PetscCall(PetscStrallocpy(lname, (char **)&b->lname));
3573:   b->type   = type;
3574:   b->label  = NULL;
3575:   b->Nv     = Nv;
3576:   b->field  = field;
3577:   b->Nc     = Nc;
3578:   b->func   = bcFunc;
3579:   b->func_t = bcFunc_t;
3580:   b->ctx    = ctx;
3581:   b->next   = NULL;
3582:   /* Append to linked list so that we can preserve the order */
3583:   if (!head) ds->boundary = b;
3584:   while (head) {
3585:     if (!head->next) {
3586:       head->next = b;
3587:       head       = b;
3588:     }
3589:     head = head->next;
3590:     ++n;
3591:   }
3592:   if (bd) {
3593:     PetscAssertPointer(bd, 13);
3594:     *bd = n;
3595:   }
3596:   PetscFunctionReturn(PETSC_SUCCESS);
3597: }

3599: /*@C
3600:   PetscDSUpdateBoundary - Change a boundary condition for the model.

3602:   Input Parameters:
3603: + ds       - The `PetscDS` object
3604: . bd       - The boundary condition number
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` ids for constrained points
3609: . values   - An array of ids for constrained points
3610: . field    - The field to constrain
3611: . Nc       - The number of constrained field components
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:   Level: developer

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

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

3628: .seealso: `PetscDS`, `PetscWeakForm`, `DMBoundaryConditionType`, `PetscDSAddBoundary()`, `PetscDSGetBoundary()`, `PetscDSGetNumBoundary()`, `DMLabel`
3629: @*/
3630: 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)
3631: {
3632:   DSBoundary b = ds->boundary;
3633:   PetscInt   n = 0;

3635:   PetscFunctionBegin;
3637:   while (b) {
3638:     if (n == bd) break;
3639:     b = b->next;
3640:     ++n;
3641:   }
3642:   PetscCheck(b, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %" PetscInt_FMT " is not in [0, %" PetscInt_FMT ")", bd, n);
3643:   if (name) {
3644:     PetscCall(PetscFree(b->name));
3645:     PetscCall(PetscStrallocpy(name, (char **)&b->name));
3646:   }
3647:   b->type = type;
3648:   if (label) {
3649:     const char *name;

3651:     b->label = label;
3652:     PetscCall(PetscFree(b->lname));
3653:     PetscCall(PetscObjectGetName((PetscObject)label, &name));
3654:     PetscCall(PetscStrallocpy(name, (char **)&b->lname));
3655:   }
3656:   if (Nv >= 0) {
3657:     b->Nv = Nv;
3658:     PetscCall(PetscFree(b->values));
3659:     PetscCall(PetscMalloc1(Nv, &b->values));
3660:     if (Nv) PetscCall(PetscArraycpy(b->values, values, Nv));
3661:   }
3662:   if (field >= 0) b->field = field;
3663:   if (Nc >= 0) {
3664:     b->Nc = Nc;
3665:     PetscCall(PetscFree(b->comps));
3666:     PetscCall(PetscMalloc1(Nc, &b->comps));
3667:     if (Nc) PetscCall(PetscArraycpy(b->comps, comps, Nc));
3668:   }
3669:   if (bcFunc) b->func = bcFunc;
3670:   if (bcFunc_t) b->func_t = bcFunc_t;
3671:   if (ctx) b->ctx = ctx;
3672:   PetscFunctionReturn(PETSC_SUCCESS);
3673: }

3675: /*@
3676:   PetscDSGetNumBoundary - Get the number of registered BC

3678:   Input Parameter:
3679: . ds - The `PetscDS` object

3681:   Output Parameter:
3682: . numBd - The number of BC

3684:   Level: intermediate

3686: .seealso: `PetscDS`, `PetscDSAddBoundary()`, `PetscDSGetBoundary()`
3687: @*/
3688: PetscErrorCode PetscDSGetNumBoundary(PetscDS ds, PetscInt *numBd)
3689: {
3690:   DSBoundary b = ds->boundary;

3692:   PetscFunctionBegin;
3694:   PetscAssertPointer(numBd, 2);
3695:   *numBd = 0;
3696:   while (b) {
3697:     ++(*numBd);
3698:     b = b->next;
3699:   }
3700:   PetscFunctionReturn(PETSC_SUCCESS);
3701: }

3703: /*@C
3704:   PetscDSGetBoundary - Gets a boundary condition to the model

3706:   Input Parameters:
3707: + ds - The `PetscDS` object
3708: - bd - The BC number

3710:   Output Parameters:
3711: + wf     - The `PetscWeakForm` holding the pointwise functions
3712: . type   - The type of condition, e.g. `DM_BC_ESSENTIAL`/`DM_BC_ESSENTIAL_FIELD` (Dirichlet), or `DM_BC_NATURAL` (Neumann)
3713: . name   - The BC name
3714: . label  - The label defining constrained points
3715: . Nv     - The number of `DMLabel` ids for constrained points
3716: . values - An array of ids for constrained points
3717: . field  - The field to constrain
3718: . Nc     - The number of constrained field components
3719: . comps  - An array of constrained component numbers
3720: . func   - A pointwise function giving boundary values
3721: . func_t - A pointwise function giving the time derivative of the boundary values
3722: - ctx    - An optional user context for bcFunc

3724:   Options Database Keys:
3725: + -bc_<boundary name> <num>      - Overrides the boundary ids
3726: - -bc_<boundary name>_comp <num> - Overrides the boundary components

3728:   Level: developer

3730: .seealso: `PetscDS`, `PetscWeakForm`, `DMBoundaryConditionType`, `PetscDSAddBoundary()`, `DMLabel`
3731: @*/
3732: 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)
3733: {
3734:   DSBoundary b = ds->boundary;
3735:   PetscInt   n = 0;

3737:   PetscFunctionBegin;
3739:   while (b) {
3740:     if (n == bd) break;
3741:     b = b->next;
3742:     ++n;
3743:   }
3744:   PetscCheck(b, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %" PetscInt_FMT " is not in [0, %" PetscInt_FMT ")", bd, n);
3745:   if (wf) {
3746:     PetscAssertPointer(wf, 3);
3747:     *wf = b->wf;
3748:   }
3749:   if (type) {
3750:     PetscAssertPointer(type, 4);
3751:     *type = b->type;
3752:   }
3753:   if (name) {
3754:     PetscAssertPointer(name, 5);
3755:     *name = b->name;
3756:   }
3757:   if (label) {
3758:     PetscAssertPointer(label, 6);
3759:     *label = b->label;
3760:   }
3761:   if (Nv) {
3762:     PetscAssertPointer(Nv, 7);
3763:     *Nv = b->Nv;
3764:   }
3765:   if (values) {
3766:     PetscAssertPointer(values, 8);
3767:     *values = b->values;
3768:   }
3769:   if (field) {
3770:     PetscAssertPointer(field, 9);
3771:     *field = b->field;
3772:   }
3773:   if (Nc) {
3774:     PetscAssertPointer(Nc, 10);
3775:     *Nc = b->Nc;
3776:   }
3777:   if (comps) {
3778:     PetscAssertPointer(comps, 11);
3779:     *comps = b->comps;
3780:   }
3781:   if (func) {
3782:     PetscAssertPointer(func, 12);
3783:     *func = b->func;
3784:   }
3785:   if (func_t) {
3786:     PetscAssertPointer(func_t, 13);
3787:     *func_t = b->func_t;
3788:   }
3789:   if (ctx) {
3790:     PetscAssertPointer(ctx, 14);
3791:     *ctx = b->ctx;
3792:   }
3793:   PetscFunctionReturn(PETSC_SUCCESS);
3794: }

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

3799:   Not Collective

3801:   Input Parameters:
3802: + ds - The source `PetscDS` object
3803: - dm - The `DM` holding labels

3805:   Level: intermediate

3807: .seealso: `PetscDS`, `DMBoundary`, `DM`, `PetscDSCopyBoundary()`, `PetscDSCreate()`, `DMGetLabel()`
3808: @*/
3809: PetscErrorCode PetscDSUpdateBoundaryLabels(PetscDS ds, DM dm)
3810: {
3811:   DSBoundary b;

3813:   PetscFunctionBegin;
3816:   for (b = ds->boundary; b; b = b->next) {
3817:     if (b->lname) PetscCall(DMGetLabel(dm, b->lname, &b->label));
3818:   }
3819:   PetscFunctionReturn(PETSC_SUCCESS);
3820: }

3822: static PetscErrorCode DSBoundaryDuplicate_Internal(DSBoundary b, DSBoundary *bNew)
3823: {
3824:   PetscFunctionBegin;
3825:   PetscCall(PetscNew(bNew));
3826:   PetscCall(PetscWeakFormCreate(PETSC_COMM_SELF, &(*bNew)->wf));
3827:   PetscCall(PetscWeakFormCopy(b->wf, (*bNew)->wf));
3828:   PetscCall(PetscStrallocpy(b->name, (char **)&((*bNew)->name)));
3829:   PetscCall(PetscStrallocpy(b->lname, (char **)&((*bNew)->lname)));
3830:   (*bNew)->type  = b->type;
3831:   (*bNew)->label = b->label;
3832:   (*bNew)->Nv    = b->Nv;
3833:   PetscCall(PetscMalloc1(b->Nv, &(*bNew)->values));
3834:   PetscCall(PetscArraycpy((*bNew)->values, b->values, b->Nv));
3835:   (*bNew)->field = b->field;
3836:   (*bNew)->Nc    = b->Nc;
3837:   PetscCall(PetscMalloc1(b->Nc, &(*bNew)->comps));
3838:   PetscCall(PetscArraycpy((*bNew)->comps, b->comps, b->Nc));
3839:   (*bNew)->func   = b->func;
3840:   (*bNew)->func_t = b->func_t;
3841:   (*bNew)->ctx    = b->ctx;
3842:   PetscFunctionReturn(PETSC_SUCCESS);
3843: }

3845: /*@
3846:   PetscDSCopyBoundary - Copy all boundary condition objects to the new problem

3848:   Not Collective

3850:   Input Parameters:
3851: + ds        - The source `PetscDS` object
3852: . numFields - The number of selected fields, or `PETSC_DEFAULT` for all fields
3853: - fields    - The selected fields, or NULL for all fields

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

3858:   Level: intermediate

3860: .seealso: `PetscDS`, `DMBoundary`, `PetscDSCopyEquations()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
3861: @*/
3862: PetscErrorCode PetscDSCopyBoundary(PetscDS ds, PetscInt numFields, const PetscInt fields[], PetscDS newds)
3863: {
3864:   DSBoundary b, *lastnext;

3866:   PetscFunctionBegin;
3869:   if (ds == newds) PetscFunctionReturn(PETSC_SUCCESS);
3870:   PetscCall(PetscDSDestroyBoundary(newds));
3871:   lastnext = &newds->boundary;
3872:   for (b = ds->boundary; b; b = b->next) {
3873:     DSBoundary bNew;
3874:     PetscInt   fieldNew = -1;

3876:     if (numFields > 0 && fields) {
3877:       PetscInt f;

3879:       for (f = 0; f < numFields; ++f)
3880:         if (b->field == fields[f]) break;
3881:       if (f == numFields) continue;
3882:       fieldNew = f;
3883:     }
3884:     PetscCall(DSBoundaryDuplicate_Internal(b, &bNew));
3885:     bNew->field = fieldNew < 0 ? b->field : fieldNew;
3886:     *lastnext   = bNew;
3887:     lastnext    = &bNew->next;
3888:   }
3889:   PetscFunctionReturn(PETSC_SUCCESS);
3890: }

3892: /*@
3893:   PetscDSDestroyBoundary - Remove all `DMBoundary` objects from the `PetscDS`

3895:   Not Collective

3897:   Input Parameter:
3898: . ds - The `PetscDS` object

3900:   Level: intermediate

3902: .seealso: `PetscDS`, `DMBoundary`, `PetscDSCopyBoundary()`, `PetscDSCopyEquations()`
3903: @*/
3904: PetscErrorCode PetscDSDestroyBoundary(PetscDS ds)
3905: {
3906:   DSBoundary next = ds->boundary;

3908:   PetscFunctionBegin;
3909:   while (next) {
3910:     DSBoundary b = next;

3912:     next = b->next;
3913:     PetscCall(PetscWeakFormDestroy(&b->wf));
3914:     PetscCall(PetscFree(b->name));
3915:     PetscCall(PetscFree(b->lname));
3916:     PetscCall(PetscFree(b->values));
3917:     PetscCall(PetscFree(b->comps));
3918:     PetscCall(PetscFree(b));
3919:   }
3920:   PetscFunctionReturn(PETSC_SUCCESS);
3921: }

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

3926:   Not Collective

3928:   Input Parameters:
3929: + prob      - The `PetscDS` object
3930: . numFields - Number of new fields
3931: . fields    - Old field number for each new field
3932: . minDegree - Minimum degree for a discretization, or `PETSC_DETERMINE` for no limit
3933: - maxDegree - Maximum degree for a discretization, or `PETSC_DETERMINE` for no limit

3935:   Output Parameter:
3936: . newprob - The `PetscDS` copy

3938:   Level: intermediate

3940: .seealso: `PetscDS`, `PetscDSSelectEquations()`, `PetscDSCopyBoundary()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
3941: @*/
3942: PetscErrorCode PetscDSSelectDiscretizations(PetscDS prob, PetscInt numFields, const PetscInt fields[], PetscInt minDegree, PetscInt maxDegree, PetscDS newprob)
3943: {
3944:   PetscInt Nf, Nfn, fn;

3946:   PetscFunctionBegin;
3948:   if (fields) PetscAssertPointer(fields, 3);
3950:   PetscCall(PetscDSGetNumFields(prob, &Nf));
3951:   PetscCall(PetscDSGetNumFields(newprob, &Nfn));
3952:   numFields = numFields < 0 ? Nf : numFields;
3953:   for (fn = 0; fn < numFields; ++fn) {
3954:     const PetscInt f = fields ? fields[fn] : fn;
3955:     PetscObject    disc;
3956:     PetscClassId   id;

3958:     if (f >= Nf) continue;
3959:     PetscCall(PetscDSGetDiscretization(prob, f, &disc));
3960:     PetscCallContinue(PetscObjectGetClassId(disc, &id));
3961:     if (id == PETSCFE_CLASSID) {
3962:       PetscFE fe;

3964:       PetscCall(PetscFELimitDegree((PetscFE)disc, minDegree, maxDegree, &fe));
3965:       PetscCall(PetscDSSetDiscretization(newprob, fn, (PetscObject)fe));
3966:       PetscCall(PetscFEDestroy(&fe));
3967:     } else {
3968:       PetscCall(PetscDSSetDiscretization(newprob, fn, disc));
3969:     }
3970:   }
3971:   PetscFunctionReturn(PETSC_SUCCESS);
3972: }

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

3977:   Not Collective

3979:   Input Parameters:
3980: + prob      - The `PetscDS` object
3981: . numFields - Number of new fields
3982: - fields    - Old field number for each new field

3984:   Output Parameter:
3985: . newprob - The `PetscDS` copy

3987:   Level: intermediate

3989: .seealso: `PetscDS`, `PetscDSSelectDiscretizations()`, `PetscDSCopyBoundary()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
3990: @*/
3991: PetscErrorCode PetscDSSelectEquations(PetscDS prob, PetscInt numFields, const PetscInt fields[], PetscDS newprob)
3992: {
3993:   PetscInt Nf, Nfn, fn, gn;

3995:   PetscFunctionBegin;
3997:   if (fields) PetscAssertPointer(fields, 3);
3999:   PetscCall(PetscDSGetNumFields(prob, &Nf));
4000:   PetscCall(PetscDSGetNumFields(newprob, &Nfn));
4001:   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);
4002:   for (fn = 0; fn < numFields; ++fn) {
4003:     const PetscInt   f = fields ? fields[fn] : fn;
4004:     PetscPointFunc   obj;
4005:     PetscPointFunc   f0, f1;
4006:     PetscBdPointFunc f0Bd, f1Bd;
4007:     PetscRiemannFunc r;

4009:     if (f >= Nf) continue;
4010:     PetscCall(PetscDSGetObjective(prob, f, &obj));
4011:     PetscCall(PetscDSGetResidual(prob, f, &f0, &f1));
4012:     PetscCall(PetscDSGetBdResidual(prob, f, &f0Bd, &f1Bd));
4013:     PetscCall(PetscDSGetRiemannSolver(prob, f, &r));
4014:     PetscCall(PetscDSSetObjective(newprob, fn, obj));
4015:     PetscCall(PetscDSSetResidual(newprob, fn, f0, f1));
4016:     PetscCall(PetscDSSetBdResidual(newprob, fn, f0Bd, f1Bd));
4017:     PetscCall(PetscDSSetRiemannSolver(newprob, fn, r));
4018:     for (gn = 0; gn < numFields; ++gn) {
4019:       const PetscInt  g = fields ? fields[gn] : gn;
4020:       PetscPointJac   g0, g1, g2, g3;
4021:       PetscPointJac   g0p, g1p, g2p, g3p;
4022:       PetscBdPointJac g0Bd, g1Bd, g2Bd, g3Bd;

4024:       if (g >= Nf) continue;
4025:       PetscCall(PetscDSGetJacobian(prob, f, g, &g0, &g1, &g2, &g3));
4026:       PetscCall(PetscDSGetJacobianPreconditioner(prob, f, g, &g0p, &g1p, &g2p, &g3p));
4027:       PetscCall(PetscDSGetBdJacobian(prob, f, g, &g0Bd, &g1Bd, &g2Bd, &g3Bd));
4028:       PetscCall(PetscDSSetJacobian(newprob, fn, gn, g0, g1, g2, g3));
4029:       PetscCall(PetscDSSetJacobianPreconditioner(newprob, fn, gn, g0p, g1p, g2p, g3p));
4030:       PetscCall(PetscDSSetBdJacobian(newprob, fn, gn, g0Bd, g1Bd, g2Bd, g3Bd));
4031:     }
4032:   }
4033:   PetscFunctionReturn(PETSC_SUCCESS);
4034: }

4036: /*@
4037:   PetscDSCopyEquations - Copy all pointwise function pointers to another `PetscDS`

4039:   Not Collective

4041:   Input Parameter:
4042: . prob - The `PetscDS` object

4044:   Output Parameter:
4045: . newprob - The `PetscDS` copy

4047:   Level: intermediate

4049: .seealso: `PetscDS`, `PetscDSCopyBoundary()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
4050: @*/
4051: PetscErrorCode PetscDSCopyEquations(PetscDS prob, PetscDS newprob)
4052: {
4053:   PetscWeakForm wf, newwf;
4054:   PetscInt      Nf, Ng;

4056:   PetscFunctionBegin;
4059:   PetscCall(PetscDSGetNumFields(prob, &Nf));
4060:   PetscCall(PetscDSGetNumFields(newprob, &Ng));
4061:   PetscCheck(Nf == Ng, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_SIZ, "Number of fields must match %" PetscInt_FMT " != %" PetscInt_FMT, Nf, Ng);
4062:   PetscCall(PetscDSGetWeakForm(prob, &wf));
4063:   PetscCall(PetscDSGetWeakForm(newprob, &newwf));
4064:   PetscCall(PetscWeakFormCopy(wf, newwf));
4065:   PetscFunctionReturn(PETSC_SUCCESS);
4066: }

4068: /*@
4069:   PetscDSCopyConstants - Copy all constants to another `PetscDS`

4071:   Not Collective

4073:   Input Parameter:
4074: . prob - The `PetscDS` object

4076:   Output Parameter:
4077: . newprob - The `PetscDS` copy

4079:   Level: intermediate

4081: .seealso: `PetscDS`, `PetscDSCopyBoundary()`, `PetscDSCopyEquations()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
4082: @*/
4083: PetscErrorCode PetscDSCopyConstants(PetscDS prob, PetscDS newprob)
4084: {
4085:   PetscInt           Nc;
4086:   const PetscScalar *constants;

4088:   PetscFunctionBegin;
4091:   PetscCall(PetscDSGetConstants(prob, &Nc, &constants));
4092:   PetscCall(PetscDSSetConstants(newprob, Nc, (PetscScalar *)constants));
4093:   PetscFunctionReturn(PETSC_SUCCESS);
4094: }

4096: /*@
4097:   PetscDSCopyExactSolutions - Copy all exact solutions to another `PetscDS`

4099:   Not Collective

4101:   Input Parameter:
4102: . ds - The `PetscDS` object

4104:   Output Parameter:
4105: . newds - The `PetscDS` copy

4107:   Level: intermediate

4109: .seealso: `PetscDS`, `PetscDSCopyBoundary()`, `PetscDSCopyEquations()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
4110: @*/
4111: PetscErrorCode PetscDSCopyExactSolutions(PetscDS ds, PetscDS newds)
4112: {
4113:   PetscSimplePointFn *sol;
4114:   void               *ctx;
4115:   PetscInt            Nf, f;

4117:   PetscFunctionBegin;
4120:   PetscCall(PetscDSGetNumFields(ds, &Nf));
4121:   for (f = 0; f < Nf; ++f) {
4122:     PetscCall(PetscDSGetExactSolution(ds, f, &sol, &ctx));
4123:     PetscCall(PetscDSSetExactSolution(newds, f, sol, ctx));
4124:     PetscCall(PetscDSGetExactSolutionTimeDerivative(ds, f, &sol, &ctx));
4125:     PetscCall(PetscDSSetExactSolutionTimeDerivative(newds, f, sol, ctx));
4126:   }
4127:   PetscFunctionReturn(PETSC_SUCCESS);
4128: }

4130: PetscErrorCode PetscDSCopy(PetscDS ds, PetscInt minDegree, PetscInt maxDegree, DM dmNew, PetscDS dsNew)
4131: {
4132:   DSBoundary b;
4133:   PetscInt   cdim, Nf, f, d;
4134:   PetscBool  isCohesive;
4135:   void      *ctx;

4137:   PetscFunctionBegin;
4138:   PetscCall(PetscDSCopyConstants(ds, dsNew));
4139:   PetscCall(PetscDSCopyExactSolutions(ds, dsNew));
4140:   PetscCall(PetscDSSelectDiscretizations(ds, PETSC_DETERMINE, NULL, minDegree, maxDegree, dsNew));
4141:   PetscCall(PetscDSCopyEquations(ds, dsNew));
4142:   PetscCall(PetscDSGetNumFields(ds, &Nf));
4143:   for (f = 0; f < Nf; ++f) {
4144:     PetscCall(PetscDSGetContext(ds, f, &ctx));
4145:     PetscCall(PetscDSSetContext(dsNew, f, ctx));
4146:     PetscCall(PetscDSGetCohesive(ds, f, &isCohesive));
4147:     PetscCall(PetscDSSetCohesive(dsNew, f, isCohesive));
4148:     PetscCall(PetscDSGetJetDegree(ds, f, &d));
4149:     PetscCall(PetscDSSetJetDegree(dsNew, f, d));
4150:   }
4151:   if (Nf) {
4152:     PetscCall(PetscDSGetCoordinateDimension(ds, &cdim));
4153:     PetscCall(PetscDSSetCoordinateDimension(dsNew, cdim));
4154:   }
4155:   PetscCall(PetscDSCopyBoundary(ds, PETSC_DETERMINE, NULL, dsNew));
4156:   for (b = dsNew->boundary; b; b = b->next) {
4157:     PetscCall(DMGetLabel(dmNew, b->lname, &b->label));
4158:     /* Do not check if label exists here, since p4est calls this for the reference tree which does not have the labels */
4159:     //PetscCheck(b->label,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Label %s missing in new DM", name);
4160:   }
4161:   PetscFunctionReturn(PETSC_SUCCESS);
4162: }

4164: PetscErrorCode PetscDSGetHeightSubspace(PetscDS prob, PetscInt height, PetscDS *subprob)
4165: {
4166:   PetscInt dim, Nf, f;

4168:   PetscFunctionBegin;
4170:   PetscAssertPointer(subprob, 3);
4171:   if (height == 0) {
4172:     *subprob = prob;
4173:     PetscFunctionReturn(PETSC_SUCCESS);
4174:   }
4175:   PetscCall(PetscDSGetNumFields(prob, &Nf));
4176:   PetscCall(PetscDSGetSpatialDimension(prob, &dim));
4177:   PetscCheck(height <= dim, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_OUTOFRANGE, "DS can only handle height in [0, %" PetscInt_FMT "], not %" PetscInt_FMT, dim, height);
4178:   if (!prob->subprobs) PetscCall(PetscCalloc1(dim, &prob->subprobs));
4179:   if (!prob->subprobs[height - 1]) {
4180:     PetscInt cdim;

4182:     PetscCall(PetscDSCreate(PetscObjectComm((PetscObject)prob), &prob->subprobs[height - 1]));
4183:     PetscCall(PetscDSGetCoordinateDimension(prob, &cdim));
4184:     PetscCall(PetscDSSetCoordinateDimension(prob->subprobs[height - 1], cdim));
4185:     for (f = 0; f < Nf; ++f) {
4186:       PetscFE      subfe;
4187:       PetscObject  obj;
4188:       PetscClassId id;

4190:       PetscCall(PetscDSGetDiscretization(prob, f, &obj));
4191:       PetscCall(PetscObjectGetClassId(obj, &id));
4192:       if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetHeightSubspace((PetscFE)obj, height, &subfe));
4193:       else SETERRQ(PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONG, "Unsupported discretization type for field %" PetscInt_FMT, f);
4194:       PetscCall(PetscDSSetDiscretization(prob->subprobs[height - 1], f, (PetscObject)subfe));
4195:     }
4196:   }
4197:   *subprob = prob->subprobs[height - 1];
4198:   PetscFunctionReturn(PETSC_SUCCESS);
4199: }

4201: PetscErrorCode PetscDSPermuteQuadPoint(PetscDS ds, PetscInt ornt, PetscInt field, PetscInt q, PetscInt *qperm)
4202: {
4203:   IS              permIS;
4204:   PetscQuadrature quad;
4205:   DMPolytopeType  ct;
4206:   const PetscInt *perm;
4207:   PetscInt        Na, Nq;

4209:   PetscFunctionBeginHot;
4210:   PetscCall(PetscFEGetQuadrature((PetscFE)ds->disc[field], &quad));
4211:   PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL));
4212:   PetscCall(PetscQuadratureGetCellType(quad, &ct));
4213:   PetscCheck(q >= 0 && q < Nq, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Quadrature point %" PetscInt_FMT " is not in [0, %" PetscInt_FMT ")", q, Nq);
4214:   Na = DMPolytopeTypeGetNumArrangements(ct) / 2;
4215:   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);
4216:   if (!ds->quadPerm[(PetscInt)ct]) PetscCall(PetscQuadratureComputePermutations(quad, NULL, &ds->quadPerm[(PetscInt)ct]));
4217:   permIS = ds->quadPerm[(PetscInt)ct][ornt + Na];
4218:   PetscCall(ISGetIndices(permIS, &perm));
4219:   *qperm = perm[q];
4220:   PetscCall(ISRestoreIndices(permIS, &perm));
4221:   PetscFunctionReturn(PETSC_SUCCESS);
4222: }

4224: PetscErrorCode PetscDSGetDiscType_Internal(PetscDS ds, PetscInt f, PetscDiscType *disctype)
4225: {
4226:   PetscObject  obj;
4227:   PetscClassId id;
4228:   PetscInt     Nf;

4230:   PetscFunctionBegin;
4232:   PetscAssertPointer(disctype, 3);
4233:   *disctype = PETSC_DISC_NONE;
4234:   PetscCall(PetscDSGetNumFields(ds, &Nf));
4235:   PetscCheck(f < Nf, PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_SIZ, "Field %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, Nf);
4236:   PetscCall(PetscDSGetDiscretization(ds, f, &obj));
4237:   if (obj) {
4238:     PetscCall(PetscObjectGetClassId(obj, &id));
4239:     if (id == PETSCFE_CLASSID) *disctype = PETSC_DISC_FE;
4240:     else *disctype = PETSC_DISC_FV;
4241:   }
4242:   PetscFunctionReturn(PETSC_SUCCESS);
4243: }

4245: static PetscErrorCode PetscDSDestroy_Basic(PetscDS ds)
4246: {
4247:   PetscFunctionBegin;
4248:   PetscCall(PetscFree(ds->data));
4249:   PetscFunctionReturn(PETSC_SUCCESS);
4250: }

4252: static PetscErrorCode PetscDSInitialize_Basic(PetscDS ds)
4253: {
4254:   PetscFunctionBegin;
4255:   ds->ops->setfromoptions = NULL;
4256:   ds->ops->setup          = NULL;
4257:   ds->ops->view           = NULL;
4258:   ds->ops->destroy        = PetscDSDestroy_Basic;
4259:   PetscFunctionReturn(PETSC_SUCCESS);
4260: }

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

4265:   Level: intermediate

4267: .seealso: `PetscDSType`, `PetscDSCreate()`, `PetscDSSetType()`
4268: M*/

4270: PETSC_EXTERN PetscErrorCode PetscDSCreate_Basic(PetscDS ds)
4271: {
4272:   PetscDS_Basic *b;

4274:   PetscFunctionBegin;
4276:   PetscCall(PetscNew(&b));
4277:   ds->data = b;

4279:   PetscCall(PetscDSInitialize_Basic(ds));
4280:   PetscFunctionReturn(PETSC_SUCCESS);
4281: }