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: /*@C
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: /*@C
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:   *ds = NULL;
685:   PetscCall(PetscDSInitializePackage());

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

689:   p->Nf           = 0;
690:   p->setup        = PETSC_FALSE;
691:   p->numConstants = 0;
692:   p->constants    = NULL;
693:   p->dimEmbed     = -1;
694:   p->useJacPre    = PETSC_TRUE;
695:   p->forceQuad    = PETSC_TRUE;
696:   PetscCall(PetscWeakFormCreate(comm, &p->wf));
697:   PetscCall(PetscArrayzero(p->quadPerm, DM_NUM_POLYTOPES));

699:   *ds = p;
700:   PetscFunctionReturn(PETSC_SUCCESS);
701: }

703: /*@
704:   PetscDSGetNumFields - Returns the number of fields in the `PetscDS`

706:   Not Collective

708:   Input Parameter:
709: . prob - The `PetscDS` object

711:   Output Parameter:
712: . Nf - The number of fields

714:   Level: beginner

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

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

730:   Not Collective

732:   Input Parameter:
733: . prob - The `PetscDS` object

735:   Output Parameter:
736: . dim - The spatial dimension

738:   Level: beginner

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

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

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

766:   Not Collective

768:   Input Parameter:
769: . prob - The `PetscDS` object

771:   Output Parameter:
772: . dimEmbed - The coordinate dimension

774:   Level: beginner

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

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

791:   Logically Collective

793:   Input Parameters:
794: + prob     - The `PetscDS` object
795: - dimEmbed - The coordinate dimension

797:   Level: beginner

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

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

813:   Not collective

815:   Input Parameter:
816: . ds - The `PetscDS` object

818:   Output Parameter:
819: . forceQuad - The flag

821:   Level: intermediate

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

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

837:   Logically collective on ds

839:   Input Parameters:
840: + ds        - The `PetscDS` object
841: - forceQuad - The flag

843:   Level: intermediate

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

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

858:   Not Collective

860:   Input Parameter:
861: . ds - The `PetscDS` object

863:   Output Parameter:
864: . isCohesive - The flag

866:   Level: developer

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

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

882:   Not Collective

884:   Input Parameter:
885: . ds - The `PetscDS` object

887:   Output Parameter:
888: . numCohesive - The number of cohesive fields

890:   Level: developer

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

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

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

909:   Not Collective

911:   Input Parameters:
912: + ds - The `PetscDS` object
913: - f  - The field index

915:   Output Parameter:
916: . isCohesive - The flag

918:   Level: developer

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

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

935:   Not Collective

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

942:   Level: developer

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

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

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

962:   Not Collective

964:   Input Parameter:
965: . prob - The `PetscDS` object

967:   Output Parameter:
968: . dim - The total problem dimension

970:   Level: beginner

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

984: /*@
985:   PetscDSGetTotalComponents - Returns the total number of components in this system

987:   Not Collective

989:   Input Parameter:
990: . prob - The `PetscDS` object

992:   Output Parameter:
993: . Nc - The total number of components

995:   Level: beginner

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

1009: /*@
1010:   PetscDSGetDiscretization - Returns the discretization object for the given field

1012:   Not Collective

1014:   Input Parameters:
1015: + prob - The `PetscDS` object
1016: - f    - The field number

1018:   Output Parameter:
1019: . disc - The discretization object

1021:   Level: beginner

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

1035: /*@
1036:   PetscDSSetDiscretization - Sets the discretization object for the given field

1038:   Not Collective

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

1045:   Level: beginner

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

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

1073: /*@
1074:   PetscDSGetWeakForm - Returns the weak form object

1076:   Not Collective

1078:   Input Parameter:
1079: . ds - The `PetscDS` object

1081:   Output Parameter:
1082: . wf - The weak form object

1084:   Level: beginner

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

1097: /*@
1098:   PetscDSSetWeakForm - Sets the weak form object

1100:   Not Collective

1102:   Input Parameters:
1103: + ds - The `PetscDS` object
1104: - wf - The weak form object

1106:   Level: beginner

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

1122: /*@
1123:   PetscDSAddDiscretization - Adds a discretization object

1125:   Not Collective

1127:   Input Parameters:
1128: + prob - The `PetscDS` object
1129: - disc - The boundary discretization object

1131:   Level: beginner

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

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

1145:   Not Collective

1147:   Input Parameter:
1148: . prob - The `PetscDS` object

1150:   Output Parameter:
1151: . q - The quadrature object

1153:   Level: intermediate

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

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

1173: /*@
1174:   PetscDSGetImplicit - Returns the flag for implicit solve for this field. This is just a guide for `TSIMEX`

1176:   Not Collective

1178:   Input Parameters:
1179: + prob - The `PetscDS` object
1180: - f    - The field number

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

1185:   Level: developer

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

1199: /*@
1200:   PetscDSSetImplicit - Set the flag for implicit solve for this field. This is just a guide for `TSIMEX`

1202:   Not Collective

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

1209:   Level: developer

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

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

1225:   Not Collective

1227:   Input Parameters:
1228: + ds - The `PetscDS` object
1229: - f  - The field number

1231:   Output Parameter:
1232: . k - The highest derivative we need to tabulate

1234:   Level: developer

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

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

1251:   Not Collective

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

1258:   Level: developer

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

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

1274:   Not Collective

1276:   Input Parameters:
1277: + ds - The `PetscDS`
1278: - f  - The test field number

1280:   Output Parameter:
1281: . obj - integrand for the test function term

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

1303:   Level: intermediate

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

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

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

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

1327:   Not Collective

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

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

1354:   Level: intermediate

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

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

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

1374:   Not Collective

1376:   Input Parameters:
1377: + ds - The `PetscDS`
1378: - f  - The test field number

1380:   Output Parameters:
1381: + f0 - integrand for the test function term
1382: - f1 - integrand for the test function gradient term

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

1404:   Level: intermediate

1406:   Note:
1407:   `f1` has an identical form and is omitted for brevity.

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

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

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

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

1430:   Not Collective

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

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

1458:   Level: intermediate

1460:   Note:
1461:   `f1` has an identical form and is omitted for brevity.

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

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

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

1481:   Not Collective

1483:   Input Parameters:
1484: + ds - The `PetscDS`
1485: - f  - The test field number

1487:   Output Parameters:
1488: + f0 - integrand for the test function term
1489: - f1 - integrand for the test function gradient term

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

1511:   Level: intermediate

1513:   Note:
1514:   `f1` has an identical form and is omitted for brevity.

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

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

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

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

1537:   Not Collective

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

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

1565:   Level: intermediate

1567:   Note:
1568:   `f1` has an identical form and is omitted for brevity.

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

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

1585: /*@
1586:   PetscDSHasJacobian - Checks that the Jacobian functions have been set

1588:   Not Collective

1590:   Input Parameter:
1591: . ds - The `PetscDS`

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

1596:   Level: intermediate

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

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

1611:   Not Collective

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

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

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

1645:   Level: intermediate

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

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

1652:   $$
1653:   \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
1654:   $$

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

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

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

1678:   Not Collective

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

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

1710:   Level: intermediate

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

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

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

1734: /*@
1735:   PetscDSUseJacobianPreconditioner - Set whether to construct a Jacobian preconditioner

1737:   Not Collective

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

1743:   Level: intermediate

1745:   Developer Note:
1746:   Should be called `PetscDSSetUseJacobianPreconditioner()`

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

1758: /*@
1759:   PetscDSHasJacobianPreconditioner - Checks if a Jacobian preconditioner matrix has been set

1761:   Not Collective

1763:   Input Parameter:
1764: . ds - The `PetscDS`

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

1769:   Level: intermediate

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

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

1787:   Not Collective

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

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

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

1821:   Level: intermediate

1823:   Note:
1824:   `g1`, `g2`, and `g3` have identical calling sequences to `g0` and are omitted for brevity.
1825:   We are using a first order FEM model for the weak form\:
1826:   \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

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

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

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

1851:   Not Collective

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

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

1883:   Level: intermediate

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

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

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

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

1910:   Not Collective

1912:   Input Parameter:
1913: . ds - The `PetscDS`

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

1918:   Level: intermediate

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

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

1933:   Not Collective

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

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

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

1967:   Level: intermediate

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

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

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

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

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

1997:   Not Collective

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

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

2029:   Level: intermediate

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

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

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

2053: /*@C
2054:   PetscDSGetRiemannSolver - Returns the Riemann solver for the given field

2056:   Not Collective

2058:   Input Parameters:
2059: + ds - The `PetscDS` object
2060: - f  - The field number

2062:   Output Parameter:
2063: . r - Riemann solver

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

2077:   Level: intermediate

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

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

2095: /*@C
2096:   PetscDSSetRiemannSolver - Sets the Riemann solver for the given field

2098:   Not Collective

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

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

2117:   Level: intermediate

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

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

2134:   Not Collective

2136:   Input Parameters:
2137: + ds - The `PetscDS`
2138: - f  - The field number

2140:   Output Parameter:
2141: . update - update function

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

2163:   Level: intermediate

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

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

2182:   Not Collective

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

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

2209:   Level: intermediate

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

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

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

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

2247:   Not Collective

2249:   Input Parameters:
2250: + ds - The PetscDS
2251: - f  - The test field number

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

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

2278:   Level: intermediate

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

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

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

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

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

2305:   Not Collective

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

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

2334:   Level: intermediate

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

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

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

2353: /*@
2354:   PetscDSHasBdJacobian - Indicates that boundary Jacobian functions have been set

2356:   Not Collective

2358:   Input Parameter:
2359: . ds - The `PetscDS`

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

2364:   Level: intermediate

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

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

2380:   Not Collective

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

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

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

2415:   Level: intermediate

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

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

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

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

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

2445:   Not Collective

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

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

2478:   Level: intermediate

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

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

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

2502: /*@
2503:   PetscDSHasBdJacobianPreconditioner - Signals that boundary Jacobian preconditioner functions have been set

2505:   Not Collective

2507:   Input Parameter:
2508: . ds - The `PetscDS`

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

2513:   Level: intermediate

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

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

2529:   Not Collective; No Fortran Support

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

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

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

2564:   Level: intermediate

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

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

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

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

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

2594:   Not Collective; No Fortran Support

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

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

2627:   Level: intermediate

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

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

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

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

2654:   Not Collective

2656:   Input Parameters:
2657: + prob - The PetscDS
2658: - f    - The test field number

2660:   Output Parameters:
2661: + sol - exact solution for the test field
2662: - ctx - exact solution context

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

2672:   Level: intermediate

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

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

2695:   Not Collective

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

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

2711:   Level: intermediate

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

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

2735:   Not Collective

2737:   Input Parameters:
2738: + prob - The `PetscDS`
2739: - f    - The test field number

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

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

2753:   Level: intermediate

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

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

2776:   Not Collective

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

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

2792:   Level: intermediate

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

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

2816:   Not Collective

2818:   Input Parameter:
2819: . prob - The `PetscDS` object

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

2825:   Level: intermediate

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

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

2847:   Not Collective

2849:   Input Parameters:
2850: + prob         - The `PetscDS` object
2851: . numConstants - The number of constants
2852: - constants    - The array of constants, `NULL` if there are none

2854:   Level: intermediate

2856: .seealso: `PetscDS`, `PetscDSGetConstants()`, `PetscDSCreate()`
2857: @*/
2858: PetscErrorCode PetscDSSetConstants(PetscDS prob, PetscInt numConstants, PetscScalar constants[])
2859: {
2860:   PetscFunctionBegin;
2862:   if (numConstants != prob->numConstants) {
2863:     PetscCall(PetscFree(prob->constants));
2864:     prob->numConstants = numConstants;
2865:     if (prob->numConstants) {
2866:       PetscCall(PetscMalloc1(prob->numConstants, &prob->constants));
2867:     } else {
2868:       prob->constants = NULL;
2869:     }
2870:   }
2871:   if (prob->numConstants) {
2872:     PetscAssertPointer(constants, 3);
2873:     PetscCall(PetscArraycpy(prob->constants, constants, prob->numConstants));
2874:   }
2875:   PetscFunctionReturn(PETSC_SUCCESS);
2876: }

2878: /*@
2879:   PetscDSGetFieldIndex - Returns the index of the given field

2881:   Not Collective

2883:   Input Parameters:
2884: + prob - The `PetscDS` object
2885: - disc - The discretization object

2887:   Output Parameter:
2888: . f - The field number

2890:   Level: beginner

2892: .seealso: `PetscDS`, `PetscGetDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
2893: @*/
2894: PetscErrorCode PetscDSGetFieldIndex(PetscDS prob, PetscObject disc, PetscInt *f)
2895: {
2896:   PetscInt g;

2898:   PetscFunctionBegin;
2900:   PetscAssertPointer(f, 3);
2901:   *f = -1;
2902:   for (g = 0; g < prob->Nf; ++g) {
2903:     if (disc == prob->disc[g]) break;
2904:   }
2905:   PetscCheck(g != prob->Nf, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONG, "Field not found in PetscDS.");
2906:   *f = g;
2907:   PetscFunctionReturn(PETSC_SUCCESS);
2908: }

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

2913:   Not Collective

2915:   Input Parameters:
2916: + prob - The `PetscDS` object
2917: - f    - The field number

2919:   Output Parameter:
2920: . size - The size

2922:   Level: beginner

2924: .seealso: `PetscDS`, `PetscDSGetFieldOffset()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
2925: @*/
2926: PetscErrorCode PetscDSGetFieldSize(PetscDS prob, PetscInt f, PetscInt *size)
2927: {
2928:   PetscFunctionBegin;
2930:   PetscAssertPointer(size, 3);
2931:   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);
2932:   PetscCall(PetscDSSetUp(prob));
2933:   *size = prob->Nb[f];
2934:   PetscFunctionReturn(PETSC_SUCCESS);
2935: }

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

2940:   Not Collective

2942:   Input Parameters:
2943: + prob - The `PetscDS` object
2944: - f    - The field number

2946:   Output Parameter:
2947: . off - The offset

2949:   Level: beginner

2951: .seealso: `PetscDS`, `PetscDSGetFieldSize()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
2952: @*/
2953: PetscErrorCode PetscDSGetFieldOffset(PetscDS prob, PetscInt f, PetscInt *off)
2954: {
2955:   PetscInt size, g;

2957:   PetscFunctionBegin;
2959:   PetscAssertPointer(off, 3);
2960:   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);
2961:   *off = 0;
2962:   for (g = 0; g < f; ++g) {
2963:     PetscCall(PetscDSGetFieldSize(prob, g, &size));
2964:     *off += size;
2965:   }
2966:   PetscFunctionReturn(PETSC_SUCCESS);
2967: }

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

2972:   Not Collective

2974:   Input Parameters:
2975: + ds - The `PetscDS` object
2976: - f  - The field number

2978:   Output Parameter:
2979: . off - The offset

2981:   Level: beginner

2983: .seealso: `PetscDS`, `PetscDSGetFieldSize()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
2984: @*/
2985: PetscErrorCode PetscDSGetFieldOffsetCohesive(PetscDS ds, PetscInt f, PetscInt *off)
2986: {
2987:   PetscInt size, g;

2989:   PetscFunctionBegin;
2991:   PetscAssertPointer(off, 3);
2992:   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);
2993:   *off = 0;
2994:   for (g = 0; g < f; ++g) {
2995:     PetscBool cohesive;

2997:     PetscCall(PetscDSGetCohesive(ds, g, &cohesive));
2998:     PetscCall(PetscDSGetFieldSize(ds, g, &size));
2999:     *off += cohesive ? size : size * 2;
3000:   }
3001:   PetscFunctionReturn(PETSC_SUCCESS);
3002: }

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

3007:   Not Collective

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

3012:   Output Parameter:
3013: . dimensions - The number of dimensions

3015:   Level: beginner

3017: .seealso: `PetscDS`, `PetscDSGetComponentOffsets()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3018: @*/
3019: PetscErrorCode PetscDSGetDimensions(PetscDS prob, PetscInt *dimensions[])
3020: {
3021:   PetscFunctionBegin;
3023:   PetscCall(PetscDSSetUp(prob));
3024:   PetscAssertPointer(dimensions, 2);
3025:   *dimensions = prob->Nb;
3026:   PetscFunctionReturn(PETSC_SUCCESS);
3027: }

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

3032:   Not Collective

3034:   Input Parameter:
3035: . prob - The `PetscDS` object

3037:   Output Parameter:
3038: . components - The number of components

3040:   Level: beginner

3042: .seealso: `PetscDS`, `PetscDSGetComponentOffsets()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3043: @*/
3044: PetscErrorCode PetscDSGetComponents(PetscDS prob, PetscInt *components[])
3045: {
3046:   PetscFunctionBegin;
3048:   PetscCall(PetscDSSetUp(prob));
3049:   PetscAssertPointer(components, 2);
3050:   *components = prob->Nc;
3051:   PetscFunctionReturn(PETSC_SUCCESS);
3052: }

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

3057:   Not Collective

3059:   Input Parameters:
3060: + prob - The `PetscDS` object
3061: - f    - The field number

3063:   Output Parameter:
3064: . off - The offset

3066:   Level: beginner

3068: .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3069: @*/
3070: PetscErrorCode PetscDSGetComponentOffset(PetscDS prob, PetscInt f, PetscInt *off)
3071: {
3072:   PetscFunctionBegin;
3074:   PetscAssertPointer(off, 3);
3075:   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);
3076:   PetscCall(PetscDSSetUp(prob));
3077:   *off = prob->off[f];
3078:   PetscFunctionReturn(PETSC_SUCCESS);
3079: }

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

3084:   Not Collective

3086:   Input Parameter:
3087: . prob - The `PetscDS` object

3089:   Output Parameter:
3090: . offsets - The offsets

3092:   Level: beginner

3094: .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3095: @*/
3096: PetscErrorCode PetscDSGetComponentOffsets(PetscDS prob, PetscInt *offsets[])
3097: {
3098:   PetscFunctionBegin;
3100:   PetscAssertPointer(offsets, 2);
3101:   PetscCall(PetscDSSetUp(prob));
3102:   *offsets = prob->off;
3103:   PetscFunctionReturn(PETSC_SUCCESS);
3104: }

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

3109:   Not Collective

3111:   Input Parameter:
3112: . prob - The `PetscDS` object

3114:   Output Parameter:
3115: . offsets - The offsets

3117:   Level: beginner

3119: .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3120: @*/
3121: PetscErrorCode PetscDSGetComponentDerivativeOffsets(PetscDS prob, PetscInt *offsets[])
3122: {
3123:   PetscFunctionBegin;
3125:   PetscAssertPointer(offsets, 2);
3126:   PetscCall(PetscDSSetUp(prob));
3127:   *offsets = prob->offDer;
3128:   PetscFunctionReturn(PETSC_SUCCESS);
3129: }

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

3134:   Not Collective

3136:   Input Parameters:
3137: + ds - The `PetscDS` object
3138: - s  - The cohesive side, 0 for negative, 1 for positive, 2 for cohesive

3140:   Output Parameter:
3141: . offsets - The offsets

3143:   Level: beginner

3145: .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3146: @*/
3147: PetscErrorCode PetscDSGetComponentOffsetsCohesive(PetscDS ds, PetscInt s, PetscInt *offsets[])
3148: {
3149:   PetscFunctionBegin;
3151:   PetscAssertPointer(offsets, 3);
3152:   PetscCheck(ds->isCohesive, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cohesive offsets are only valid for a cohesive DS");
3153:   PetscCheck(!(s < 0) && !(s > 2), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cohesive side %" PetscInt_FMT " is not in [0, 2]", s);
3154:   PetscCall(PetscDSSetUp(ds));
3155:   *offsets = ds->offCohesive[s];
3156:   PetscFunctionReturn(PETSC_SUCCESS);
3157: }

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

3162:   Not Collective

3164:   Input Parameters:
3165: + ds - The `PetscDS` object
3166: - s  - The cohesive side, 0 for negative, 1 for positive, 2 for cohesive

3168:   Output Parameter:
3169: . offsets - The offsets

3171:   Level: beginner

3173: .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3174: @*/
3175: PetscErrorCode PetscDSGetComponentDerivativeOffsetsCohesive(PetscDS ds, PetscInt s, PetscInt *offsets[])
3176: {
3177:   PetscFunctionBegin;
3179:   PetscAssertPointer(offsets, 3);
3180:   PetscCheck(ds->isCohesive, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cohesive offsets are only valid for a cohesive DS");
3181:   PetscCheck(!(s < 0) && !(s > 2), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cohesive side %" PetscInt_FMT " is not in [0, 2]", s);
3182:   PetscCall(PetscDSSetUp(ds));
3183:   *offsets = ds->offDerCohesive[s];
3184:   PetscFunctionReturn(PETSC_SUCCESS);
3185: }

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

3190:   Not Collective

3192:   Input Parameter:
3193: . prob - The `PetscDS` object

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

3198:   Level: intermediate

3200: .seealso: `PetscDS`, `PetscTabulation`, `PetscDSCreate()`
3201: @*/
3202: PetscErrorCode PetscDSGetTabulation(PetscDS prob, PetscTabulation *T[])
3203: {
3204:   PetscFunctionBegin;
3206:   PetscAssertPointer(T, 2);
3207:   PetscCall(PetscDSSetUp(prob));
3208:   *T = prob->T;
3209:   PetscFunctionReturn(PETSC_SUCCESS);
3210: }

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

3215:   Not Collective

3217:   Input Parameter:
3218: . prob - The `PetscDS` object

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

3223:   Level: intermediate

3225: .seealso: `PetscTabulation`, `PetscDS`, `PetscDSGetTabulation()`, `PetscDSCreate()`
3226: @*/
3227: PetscErrorCode PetscDSGetFaceTabulation(PetscDS prob, PetscTabulation *Tf[])
3228: {
3229:   PetscFunctionBegin;
3231:   PetscAssertPointer(Tf, 2);
3232:   PetscCall(PetscDSSetUp(prob));
3233:   *Tf = prob->Tf;
3234:   PetscFunctionReturn(PETSC_SUCCESS);
3235: }

3237: PetscErrorCode PetscDSGetEvaluationArrays(PetscDS prob, PetscScalar **u, PetscScalar **u_t, PetscScalar **u_x)
3238: {
3239:   PetscFunctionBegin;
3241:   PetscCall(PetscDSSetUp(prob));
3242:   if (u) {
3243:     PetscAssertPointer(u, 2);
3244:     *u = prob->u;
3245:   }
3246:   if (u_t) {
3247:     PetscAssertPointer(u_t, 3);
3248:     *u_t = prob->u_t;
3249:   }
3250:   if (u_x) {
3251:     PetscAssertPointer(u_x, 4);
3252:     *u_x = prob->u_x;
3253:   }
3254:   PetscFunctionReturn(PETSC_SUCCESS);
3255: }

3257: PetscErrorCode PetscDSGetWeakFormArrays(PetscDS prob, PetscScalar **f0, PetscScalar **f1, PetscScalar **g0, PetscScalar **g1, PetscScalar **g2, PetscScalar **g3)
3258: {
3259:   PetscFunctionBegin;
3261:   PetscCall(PetscDSSetUp(prob));
3262:   if (f0) {
3263:     PetscAssertPointer(f0, 2);
3264:     *f0 = prob->f0;
3265:   }
3266:   if (f1) {
3267:     PetscAssertPointer(f1, 3);
3268:     *f1 = prob->f1;
3269:   }
3270:   if (g0) {
3271:     PetscAssertPointer(g0, 4);
3272:     *g0 = prob->g0;
3273:   }
3274:   if (g1) {
3275:     PetscAssertPointer(g1, 5);
3276:     *g1 = prob->g1;
3277:   }
3278:   if (g2) {
3279:     PetscAssertPointer(g2, 6);
3280:     *g2 = prob->g2;
3281:   }
3282:   if (g3) {
3283:     PetscAssertPointer(g3, 7);
3284:     *g3 = prob->g3;
3285:   }
3286:   PetscFunctionReturn(PETSC_SUCCESS);
3287: }

3289: PetscErrorCode PetscDSGetWorkspace(PetscDS prob, PetscReal **x, PetscScalar **basisReal, PetscScalar **basisDerReal, PetscScalar **testReal, PetscScalar **testDerReal)
3290: {
3291:   PetscFunctionBegin;
3293:   PetscCall(PetscDSSetUp(prob));
3294:   if (x) {
3295:     PetscAssertPointer(x, 2);
3296:     *x = prob->x;
3297:   }
3298:   if (basisReal) {
3299:     PetscAssertPointer(basisReal, 3);
3300:     *basisReal = prob->basisReal;
3301:   }
3302:   if (basisDerReal) {
3303:     PetscAssertPointer(basisDerReal, 4);
3304:     *basisDerReal = prob->basisDerReal;
3305:   }
3306:   if (testReal) {
3307:     PetscAssertPointer(testReal, 5);
3308:     *testReal = prob->testReal;
3309:   }
3310:   if (testDerReal) {
3311:     PetscAssertPointer(testDerReal, 6);
3312:     *testDerReal = prob->testDerReal;
3313:   }
3314:   PetscFunctionReturn(PETSC_SUCCESS);
3315: }

3317: /*@C
3318:   PetscDSAddBoundary - Add a boundary condition to the model.

3320:   Collective

3322:   Input Parameters:
3323: + ds       - The PetscDS object
3324: . type     - The type of condition, e.g. `DM_BC_ESSENTIAL`/`DM_BC_ESSENTIAL_FIELD` (Dirichlet), or `DM_BC_NATURAL` (Neumann)
3325: . name     - The BC name
3326: . label    - The label defining constrained points
3327: . Nv       - The number of `DMLabel` values for constrained points
3328: . values   - An array of label values for constrained points
3329: . field    - The field to constrain
3330: . Nc       - The number of constrained field components (0 will constrain all fields)
3331: . comps    - An array of constrained component numbers
3332: . bcFunc   - A pointwise function giving boundary values
3333: . bcFunc_t - A pointwise function giving the time derivative of the boundary values, or NULL
3334: - ctx      - An optional user context for bcFunc

3336:   Output Parameter:
3337: . bd - The boundary number

3339:   Options Database Keys:
3340: + -bc_<boundary name> <num>      - Overrides the boundary ids
3341: - -bc_<boundary name>_comp <num> - Overrides the boundary components

3343:   Level: developer

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

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

3350:   If the type is `DM_BC_ESSENTIAL_FIELD` or other _FIELD value, then the calling sequence is\:
3351: .vb
3352:   void bcFunc(PetscInt dim, PetscInt Nf, PetscInt NfAux,
3353:               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
3354:               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
3355:               PetscReal time, const PetscReal x[], PetscScalar bcval[])
3356: .ve
3357: + dim - the spatial dimension
3358: . Nf - the number of fields
3359: . uOff - the offset into u[] and u_t[] for each field
3360: . uOff_x - the offset into u_x[] for each field
3361: . u - each field evaluated at the current point
3362: . u_t - the time derivative of each field evaluated at the current point
3363: . u_x - the gradient of each field evaluated at the current point
3364: . aOff - the offset into a[] and a_t[] for each auxiliary field
3365: . aOff_x - the offset into a_x[] for each auxiliary field
3366: . a - each auxiliary field evaluated at the current point
3367: . a_t - the time derivative of each auxiliary field evaluated at the current point
3368: . a_x - the gradient of auxiliary each field evaluated at the current point
3369: . t - current time
3370: . x - coordinates of the current point
3371: . numConstants - number of constant parameters
3372: . constants - constant parameters
3373: - bcval - output values at the current point

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

3381: .seealso: `PetscDS`, `PetscWeakForm`, `DMLabel`, `DMBoundaryConditionType`, `PetscDSAddBoundaryByName()`, `PetscDSGetBoundary()`, `PetscDSSetResidual()`, `PetscDSSetBdResidual()`
3382: @*/
3383: 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)
3384: {
3385:   DSBoundary  head = ds->boundary, b;
3386:   PetscInt    n    = 0;
3387:   const char *lname;

3389:   PetscFunctionBegin;
3392:   PetscAssertPointer(name, 3);
3397:   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);
3398:   if (Nc > 0) {
3399:     PetscInt *fcomps;
3400:     PetscInt  c;

3402:     PetscCall(PetscDSGetComponents(ds, &fcomps));
3403:     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);
3404:     for (c = 0; c < Nc; ++c) {
3405:       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);
3406:     }
3407:   }
3408:   PetscCall(PetscNew(&b));
3409:   PetscCall(PetscStrallocpy(name, (char **)&b->name));
3410:   PetscCall(PetscWeakFormCreate(PETSC_COMM_SELF, &b->wf));
3411:   PetscCall(PetscWeakFormSetNumFields(b->wf, ds->Nf));
3412:   PetscCall(PetscMalloc1(Nv, &b->values));
3413:   if (Nv) PetscCall(PetscArraycpy(b->values, values, Nv));
3414:   PetscCall(PetscMalloc1(Nc, &b->comps));
3415:   if (Nc) PetscCall(PetscArraycpy(b->comps, comps, Nc));
3416:   PetscCall(PetscObjectGetName((PetscObject)label, &lname));
3417:   PetscCall(PetscStrallocpy(lname, (char **)&b->lname));
3418:   b->type   = type;
3419:   b->label  = label;
3420:   b->Nv     = Nv;
3421:   b->field  = field;
3422:   b->Nc     = Nc;
3423:   b->func   = bcFunc;
3424:   b->func_t = bcFunc_t;
3425:   b->ctx    = ctx;
3426:   b->next   = NULL;
3427:   /* Append to linked list so that we can preserve the order */
3428:   if (!head) ds->boundary = b;
3429:   while (head) {
3430:     if (!head->next) {
3431:       head->next = b;
3432:       head       = b;
3433:     }
3434:     head = head->next;
3435:     ++n;
3436:   }
3437:   if (bd) {
3438:     PetscAssertPointer(bd, 13);
3439:     *bd = n;
3440:   }
3441:   PetscFunctionReturn(PETSC_SUCCESS);
3442: }

3444: // PetscClangLinter pragma ignore: -fdoc-section-header-unknown
3445: /*@C
3446:   PetscDSAddBoundaryByName - Add a boundary condition to the model.

3448:   Collective

3450:   Input Parameters:
3451: + ds       - The `PetscDS` object
3452: . type     - The type of condition, e.g. `DM_BC_ESSENTIAL`/`DM_BC_ESSENTIAL_FIELD` (Dirichlet), or `DM_BC_NATURAL` (Neumann)
3453: . name     - The BC name
3454: . lname    - The naem of the label defining constrained points
3455: . Nv       - The number of `DMLabel` values for constrained points
3456: . values   - An array of label values for constrained points
3457: . field    - The field to constrain
3458: . Nc       - The number of constrained field components (0 will constrain all fields)
3459: . comps    - An array of constrained component numbers
3460: . bcFunc   - A pointwise function giving boundary values
3461: . bcFunc_t - A pointwise function giving the time derivative of the boundary values, or NULL
3462: - ctx      - An optional user context for bcFunc

3464:   Output Parameter:
3465: . bd - The boundary number

3467:   Options Database Keys:
3468: + -bc_<boundary name> <num>      - Overrides the boundary ids
3469: - -bc_<boundary name>_comp <num> - Overrides the boundary components

3471:   Calling Sequence of `bcFunc` and `bcFunc_t`:
3472:   If the type is `DM_BC_ESSENTIAL`
3473: .vb
3474:   void bcFunc(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar bcval[])
3475: .ve
3476:   If the type is `DM_BC_ESSENTIAL_FIELD` or other _FIELD value,
3477: .vb
3478:   void bcFunc(PetscInt dim, PetscInt Nf, PetscInt NfAux,
3479:               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
3480:               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
3481:               PetscReal time, const PetscReal x[], PetscScalar bcval[])
3482: .ve
3483: + dim - the spatial dimension
3484: . Nf - the number of fields
3485: . uOff - the offset into u[] and u_t[] for each field
3486: . uOff_x - the offset into u_x[] for each field
3487: . u - each field evaluated at the current point
3488: . u_t - the time derivative of each field evaluated at the current point
3489: . u_x - the gradient of each field evaluated at the current point
3490: . aOff - the offset into a[] and a_t[] for each auxiliary field
3491: . aOff_x - the offset into a_x[] for each auxiliary field
3492: . a - each auxiliary field evaluated at the current point
3493: . a_t - the time derivative of each auxiliary field evaluated at the current point
3494: . a_x - the gradient of auxiliary each field evaluated at the current point
3495: . t - current time
3496: . x - coordinates of the current point
3497: . numConstants - number of constant parameters
3498: . constants - constant parameters
3499: - bcval - output values at the current point

3501:   Level: developer

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

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

3511: .seealso: `PetscDS`, `PetscWeakForm`, `DMLabel`, `DMBoundaryConditionType`, `PetscDSAddBoundary()`, `PetscDSGetBoundary()`, `PetscDSSetResidual()`, `PetscDSSetBdResidual()`
3512: @*/
3513: 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)
3514: {
3515:   DSBoundary head = ds->boundary, b;
3516:   PetscInt   n    = 0;

3518:   PetscFunctionBegin;
3521:   PetscAssertPointer(name, 3);
3522:   PetscAssertPointer(lname, 4);
3526:   PetscCall(PetscNew(&b));
3527:   PetscCall(PetscStrallocpy(name, (char **)&b->name));
3528:   PetscCall(PetscWeakFormCreate(PETSC_COMM_SELF, &b->wf));
3529:   PetscCall(PetscWeakFormSetNumFields(b->wf, ds->Nf));
3530:   PetscCall(PetscMalloc1(Nv, &b->values));
3531:   if (Nv) PetscCall(PetscArraycpy(b->values, values, Nv));
3532:   PetscCall(PetscMalloc1(Nc, &b->comps));
3533:   if (Nc) PetscCall(PetscArraycpy(b->comps, comps, Nc));
3534:   PetscCall(PetscStrallocpy(lname, (char **)&b->lname));
3535:   b->type   = type;
3536:   b->label  = NULL;
3537:   b->Nv     = Nv;
3538:   b->field  = field;
3539:   b->Nc     = Nc;
3540:   b->func   = bcFunc;
3541:   b->func_t = bcFunc_t;
3542:   b->ctx    = ctx;
3543:   b->next   = NULL;
3544:   /* Append to linked list so that we can preserve the order */
3545:   if (!head) ds->boundary = b;
3546:   while (head) {
3547:     if (!head->next) {
3548:       head->next = b;
3549:       head       = b;
3550:     }
3551:     head = head->next;
3552:     ++n;
3553:   }
3554:   if (bd) {
3555:     PetscAssertPointer(bd, 13);
3556:     *bd = n;
3557:   }
3558:   PetscFunctionReturn(PETSC_SUCCESS);
3559: }

3561: /*@C
3562:   PetscDSUpdateBoundary - Change a boundary condition for the model.

3564:   Input Parameters:
3565: + ds       - The `PetscDS` object
3566: . bd       - The boundary condition number
3567: . type     - The type of condition, e.g. `DM_BC_ESSENTIAL`/`DM_BC_ESSENTIAL_FIELD` (Dirichlet), or `DM_BC_NATURAL` (Neumann)
3568: . name     - The BC name
3569: . label    - The label defining constrained points
3570: . Nv       - The number of `DMLabel` ids for constrained points
3571: . values   - An array of ids for constrained points
3572: . field    - The field to constrain
3573: . Nc       - The number of constrained field components
3574: . comps    - An array of constrained component numbers
3575: . bcFunc   - A pointwise function giving boundary values
3576: . bcFunc_t - A pointwise function giving the time derivative of the boundary values, or NULL
3577: - ctx      - An optional user context for bcFunc

3579:   Level: developer

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

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

3590: .seealso: `PetscDS`, `PetscWeakForm`, `DMBoundaryConditionType`, `PetscDSAddBoundary()`, `PetscDSGetBoundary()`, `PetscDSGetNumBoundary()`, `DMLabel`
3591: @*/
3592: 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)
3593: {
3594:   DSBoundary b = ds->boundary;
3595:   PetscInt   n = 0;

3597:   PetscFunctionBegin;
3599:   while (b) {
3600:     if (n == bd) break;
3601:     b = b->next;
3602:     ++n;
3603:   }
3604:   PetscCheck(b, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %" PetscInt_FMT " is not in [0, %" PetscInt_FMT ")", bd, n);
3605:   if (name) {
3606:     PetscCall(PetscFree(b->name));
3607:     PetscCall(PetscStrallocpy(name, (char **)&b->name));
3608:   }
3609:   b->type = type;
3610:   if (label) {
3611:     const char *name;

3613:     b->label = label;
3614:     PetscCall(PetscFree(b->lname));
3615:     PetscCall(PetscObjectGetName((PetscObject)label, &name));
3616:     PetscCall(PetscStrallocpy(name, (char **)&b->lname));
3617:   }
3618:   if (Nv >= 0) {
3619:     b->Nv = Nv;
3620:     PetscCall(PetscFree(b->values));
3621:     PetscCall(PetscMalloc1(Nv, &b->values));
3622:     if (Nv) PetscCall(PetscArraycpy(b->values, values, Nv));
3623:   }
3624:   if (field >= 0) b->field = field;
3625:   if (Nc >= 0) {
3626:     b->Nc = Nc;
3627:     PetscCall(PetscFree(b->comps));
3628:     PetscCall(PetscMalloc1(Nc, &b->comps));
3629:     if (Nc) PetscCall(PetscArraycpy(b->comps, comps, Nc));
3630:   }
3631:   if (bcFunc) b->func = bcFunc;
3632:   if (bcFunc_t) b->func_t = bcFunc_t;
3633:   if (ctx) b->ctx = ctx;
3634:   PetscFunctionReturn(PETSC_SUCCESS);
3635: }

3637: /*@
3638:   PetscDSGetNumBoundary - Get the number of registered BC

3640:   Input Parameter:
3641: . ds - The `PetscDS` object

3643:   Output Parameter:
3644: . numBd - The number of BC

3646:   Level: intermediate

3648: .seealso: `PetscDS`, `PetscDSAddBoundary()`, `PetscDSGetBoundary()`
3649: @*/
3650: PetscErrorCode PetscDSGetNumBoundary(PetscDS ds, PetscInt *numBd)
3651: {
3652:   DSBoundary b = ds->boundary;

3654:   PetscFunctionBegin;
3656:   PetscAssertPointer(numBd, 2);
3657:   *numBd = 0;
3658:   while (b) {
3659:     ++(*numBd);
3660:     b = b->next;
3661:   }
3662:   PetscFunctionReturn(PETSC_SUCCESS);
3663: }

3665: /*@C
3666:   PetscDSGetBoundary - Gets a boundary condition to the model

3668:   Input Parameters:
3669: + ds - The `PetscDS` object
3670: - bd - The BC number

3672:   Output Parameters:
3673: + wf     - The `PetscWeakForm` holding the pointwise functions
3674: . type   - The type of condition, e.g. `DM_BC_ESSENTIAL`/`DM_BC_ESSENTIAL_FIELD` (Dirichlet), or `DM_BC_NATURAL` (Neumann)
3675: . name   - The BC name
3676: . label  - The label defining constrained points
3677: . Nv     - The number of `DMLabel` ids for constrained points
3678: . values - An array of ids for constrained points
3679: . field  - The field to constrain
3680: . Nc     - The number of constrained field components
3681: . comps  - An array of constrained component numbers
3682: . func   - A pointwise function giving boundary values
3683: . func_t - A pointwise function giving the time derivative of the boundary values
3684: - ctx    - An optional user context for bcFunc

3686:   Options Database Keys:
3687: + -bc_<boundary name> <num>      - Overrides the boundary ids
3688: - -bc_<boundary name>_comp <num> - Overrides the boundary components

3690:   Level: developer

3692: .seealso: `PetscDS`, `PetscWeakForm`, `DMBoundaryConditionType`, `PetscDSAddBoundary()`, `DMLabel`
3693: @*/
3694: 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)
3695: {
3696:   DSBoundary b = ds->boundary;
3697:   PetscInt   n = 0;

3699:   PetscFunctionBegin;
3701:   while (b) {
3702:     if (n == bd) break;
3703:     b = b->next;
3704:     ++n;
3705:   }
3706:   PetscCheck(b, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %" PetscInt_FMT " is not in [0, %" PetscInt_FMT ")", bd, n);
3707:   if (wf) {
3708:     PetscAssertPointer(wf, 3);
3709:     *wf = b->wf;
3710:   }
3711:   if (type) {
3712:     PetscAssertPointer(type, 4);
3713:     *type = b->type;
3714:   }
3715:   if (name) {
3716:     PetscAssertPointer(name, 5);
3717:     *name = b->name;
3718:   }
3719:   if (label) {
3720:     PetscAssertPointer(label, 6);
3721:     *label = b->label;
3722:   }
3723:   if (Nv) {
3724:     PetscAssertPointer(Nv, 7);
3725:     *Nv = b->Nv;
3726:   }
3727:   if (values) {
3728:     PetscAssertPointer(values, 8);
3729:     *values = b->values;
3730:   }
3731:   if (field) {
3732:     PetscAssertPointer(field, 9);
3733:     *field = b->field;
3734:   }
3735:   if (Nc) {
3736:     PetscAssertPointer(Nc, 10);
3737:     *Nc = b->Nc;
3738:   }
3739:   if (comps) {
3740:     PetscAssertPointer(comps, 11);
3741:     *comps = b->comps;
3742:   }
3743:   if (func) {
3744:     PetscAssertPointer(func, 12);
3745:     *func = b->func;
3746:   }
3747:   if (func_t) {
3748:     PetscAssertPointer(func_t, 13);
3749:     *func_t = b->func_t;
3750:   }
3751:   if (ctx) {
3752:     PetscAssertPointer(ctx, 14);
3753:     *ctx = b->ctx;
3754:   }
3755:   PetscFunctionReturn(PETSC_SUCCESS);
3756: }

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

3761:   Not Collective

3763:   Input Parameters:
3764: + ds - The source `PetscDS` object
3765: - dm - The `DM` holding labels

3767:   Level: intermediate

3769: .seealso: `PetscDS`, `DMBoundary`, `DM`, `PetscDSCopyBoundary()`, `PetscDSCreate()`, `DMGetLabel()`
3770: @*/
3771: PetscErrorCode PetscDSUpdateBoundaryLabels(PetscDS ds, DM dm)
3772: {
3773:   DSBoundary b;

3775:   PetscFunctionBegin;
3778:   for (b = ds->boundary; b; b = b->next) {
3779:     if (b->lname) PetscCall(DMGetLabel(dm, b->lname, &b->label));
3780:   }
3781:   PetscFunctionReturn(PETSC_SUCCESS);
3782: }

3784: static PetscErrorCode DSBoundaryDuplicate_Internal(DSBoundary b, DSBoundary *bNew)
3785: {
3786:   PetscFunctionBegin;
3787:   PetscCall(PetscNew(bNew));
3788:   PetscCall(PetscWeakFormCreate(PETSC_COMM_SELF, &(*bNew)->wf));
3789:   PetscCall(PetscWeakFormCopy(b->wf, (*bNew)->wf));
3790:   PetscCall(PetscStrallocpy(b->name, (char **)&((*bNew)->name)));
3791:   PetscCall(PetscStrallocpy(b->lname, (char **)&((*bNew)->lname)));
3792:   (*bNew)->type  = b->type;
3793:   (*bNew)->label = b->label;
3794:   (*bNew)->Nv    = b->Nv;
3795:   PetscCall(PetscMalloc1(b->Nv, &(*bNew)->values));
3796:   PetscCall(PetscArraycpy((*bNew)->values, b->values, b->Nv));
3797:   (*bNew)->field = b->field;
3798:   (*bNew)->Nc    = b->Nc;
3799:   PetscCall(PetscMalloc1(b->Nc, &(*bNew)->comps));
3800:   PetscCall(PetscArraycpy((*bNew)->comps, b->comps, b->Nc));
3801:   (*bNew)->func   = b->func;
3802:   (*bNew)->func_t = b->func_t;
3803:   (*bNew)->ctx    = b->ctx;
3804:   PetscFunctionReturn(PETSC_SUCCESS);
3805: }

3807: /*@
3808:   PetscDSCopyBoundary - Copy all boundary condition objects to the new problem

3810:   Not Collective

3812:   Input Parameters:
3813: + ds        - The source `PetscDS` object
3814: . numFields - The number of selected fields, or `PETSC_DEFAULT` for all fields
3815: - fields    - The selected fields, or NULL for all fields

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

3820:   Level: intermediate

3822: .seealso: `PetscDS`, `DMBoundary`, `PetscDSCopyEquations()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
3823: @*/
3824: PetscErrorCode PetscDSCopyBoundary(PetscDS ds, PetscInt numFields, const PetscInt fields[], PetscDS newds)
3825: {
3826:   DSBoundary b, *lastnext;

3828:   PetscFunctionBegin;
3831:   if (ds == newds) PetscFunctionReturn(PETSC_SUCCESS);
3832:   PetscCall(PetscDSDestroyBoundary(newds));
3833:   lastnext = &newds->boundary;
3834:   for (b = ds->boundary; b; b = b->next) {
3835:     DSBoundary bNew;
3836:     PetscInt   fieldNew = -1;

3838:     if (numFields > 0 && fields) {
3839:       PetscInt f;

3841:       for (f = 0; f < numFields; ++f)
3842:         if (b->field == fields[f]) break;
3843:       if (f == numFields) continue;
3844:       fieldNew = f;
3845:     }
3846:     PetscCall(DSBoundaryDuplicate_Internal(b, &bNew));
3847:     bNew->field = fieldNew < 0 ? b->field : fieldNew;
3848:     *lastnext   = bNew;
3849:     lastnext    = &bNew->next;
3850:   }
3851:   PetscFunctionReturn(PETSC_SUCCESS);
3852: }

3854: /*@
3855:   PetscDSDestroyBoundary - Remove all `DMBoundary` objects from the `PetscDS`

3857:   Not Collective

3859:   Input Parameter:
3860: . ds - The `PetscDS` object

3862:   Level: intermediate

3864: .seealso: `PetscDS`, `DMBoundary`, `PetscDSCopyBoundary()`, `PetscDSCopyEquations()`
3865: @*/
3866: PetscErrorCode PetscDSDestroyBoundary(PetscDS ds)
3867: {
3868:   DSBoundary next = ds->boundary;

3870:   PetscFunctionBegin;
3871:   while (next) {
3872:     DSBoundary b = next;

3874:     next = b->next;
3875:     PetscCall(PetscWeakFormDestroy(&b->wf));
3876:     PetscCall(PetscFree(b->name));
3877:     PetscCall(PetscFree(b->lname));
3878:     PetscCall(PetscFree(b->values));
3879:     PetscCall(PetscFree(b->comps));
3880:     PetscCall(PetscFree(b));
3881:   }
3882:   PetscFunctionReturn(PETSC_SUCCESS);
3883: }

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

3888:   Not Collective

3890:   Input Parameters:
3891: + prob      - The `PetscDS` object
3892: . numFields - Number of new fields
3893: - fields    - Old field number for each new field

3895:   Output Parameter:
3896: . newprob - The `PetscDS` copy

3898:   Level: intermediate

3900: .seealso: `PetscDS`, `PetscDSSelectEquations()`, `PetscDSCopyBoundary()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
3901: @*/
3902: PetscErrorCode PetscDSSelectDiscretizations(PetscDS prob, PetscInt numFields, const PetscInt fields[], PetscDS newprob)
3903: {
3904:   PetscInt Nf, Nfn, fn;

3906:   PetscFunctionBegin;
3908:   if (fields) PetscAssertPointer(fields, 3);
3910:   PetscCall(PetscDSGetNumFields(prob, &Nf));
3911:   PetscCall(PetscDSGetNumFields(newprob, &Nfn));
3912:   numFields = numFields < 0 ? Nf : numFields;
3913:   for (fn = 0; fn < numFields; ++fn) {
3914:     const PetscInt f = fields ? fields[fn] : fn;
3915:     PetscObject    disc;

3917:     if (f >= Nf) continue;
3918:     PetscCall(PetscDSGetDiscretization(prob, f, &disc));
3919:     PetscCall(PetscDSSetDiscretization(newprob, fn, disc));
3920:   }
3921:   PetscFunctionReturn(PETSC_SUCCESS);
3922: }

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

3927:   Not Collective

3929:   Input Parameters:
3930: + prob      - The `PetscDS` object
3931: . numFields - Number of new fields
3932: - fields    - Old field number for each new field

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

3937:   Level: intermediate

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

3945:   PetscFunctionBegin;
3947:   if (fields) PetscAssertPointer(fields, 3);
3949:   PetscCall(PetscDSGetNumFields(prob, &Nf));
3950:   PetscCall(PetscDSGetNumFields(newprob, &Nfn));
3951:   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);
3952:   for (fn = 0; fn < numFields; ++fn) {
3953:     const PetscInt   f = fields ? fields[fn] : fn;
3954:     PetscPointFunc   obj;
3955:     PetscPointFunc   f0, f1;
3956:     PetscBdPointFunc f0Bd, f1Bd;
3957:     PetscRiemannFunc r;

3959:     if (f >= Nf) continue;
3960:     PetscCall(PetscDSGetObjective(prob, f, &obj));
3961:     PetscCall(PetscDSGetResidual(prob, f, &f0, &f1));
3962:     PetscCall(PetscDSGetBdResidual(prob, f, &f0Bd, &f1Bd));
3963:     PetscCall(PetscDSGetRiemannSolver(prob, f, &r));
3964:     PetscCall(PetscDSSetObjective(newprob, fn, obj));
3965:     PetscCall(PetscDSSetResidual(newprob, fn, f0, f1));
3966:     PetscCall(PetscDSSetBdResidual(newprob, fn, f0Bd, f1Bd));
3967:     PetscCall(PetscDSSetRiemannSolver(newprob, fn, r));
3968:     for (gn = 0; gn < numFields; ++gn) {
3969:       const PetscInt  g = fields ? fields[gn] : gn;
3970:       PetscPointJac   g0, g1, g2, g3;
3971:       PetscPointJac   g0p, g1p, g2p, g3p;
3972:       PetscBdPointJac g0Bd, g1Bd, g2Bd, g3Bd;

3974:       if (g >= Nf) continue;
3975:       PetscCall(PetscDSGetJacobian(prob, f, g, &g0, &g1, &g2, &g3));
3976:       PetscCall(PetscDSGetJacobianPreconditioner(prob, f, g, &g0p, &g1p, &g2p, &g3p));
3977:       PetscCall(PetscDSGetBdJacobian(prob, f, g, &g0Bd, &g1Bd, &g2Bd, &g3Bd));
3978:       PetscCall(PetscDSSetJacobian(newprob, fn, gn, g0, g1, g2, g3));
3979:       PetscCall(PetscDSSetJacobianPreconditioner(newprob, fn, gn, g0p, g1p, g2p, g3p));
3980:       PetscCall(PetscDSSetBdJacobian(newprob, fn, gn, g0Bd, g1Bd, g2Bd, g3Bd));
3981:     }
3982:   }
3983:   PetscFunctionReturn(PETSC_SUCCESS);
3984: }

3986: /*@
3987:   PetscDSCopyEquations - Copy all pointwise function pointers to another `PetscDS`

3989:   Not Collective

3991:   Input Parameter:
3992: . prob - The `PetscDS` object

3994:   Output Parameter:
3995: . newprob - The `PetscDS` copy

3997:   Level: intermediate

3999: .seealso: `PetscDS`, `PetscDSCopyBoundary()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
4000: @*/
4001: PetscErrorCode PetscDSCopyEquations(PetscDS prob, PetscDS newprob)
4002: {
4003:   PetscWeakForm wf, newwf;
4004:   PetscInt      Nf, Ng;

4006:   PetscFunctionBegin;
4009:   PetscCall(PetscDSGetNumFields(prob, &Nf));
4010:   PetscCall(PetscDSGetNumFields(newprob, &Ng));
4011:   PetscCheck(Nf == Ng, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_SIZ, "Number of fields must match %" PetscInt_FMT " != %" PetscInt_FMT, Nf, Ng);
4012:   PetscCall(PetscDSGetWeakForm(prob, &wf));
4013:   PetscCall(PetscDSGetWeakForm(newprob, &newwf));
4014:   PetscCall(PetscWeakFormCopy(wf, newwf));
4015:   PetscFunctionReturn(PETSC_SUCCESS);
4016: }

4018: /*@
4019:   PetscDSCopyConstants - Copy all constants to another `PetscDS`

4021:   Not Collective

4023:   Input Parameter:
4024: . prob - The `PetscDS` object

4026:   Output Parameter:
4027: . newprob - The `PetscDS` copy

4029:   Level: intermediate

4031: .seealso: `PetscDS`, `PetscDSCopyBoundary()`, `PetscDSCopyEquations()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
4032: @*/
4033: PetscErrorCode PetscDSCopyConstants(PetscDS prob, PetscDS newprob)
4034: {
4035:   PetscInt           Nc;
4036:   const PetscScalar *constants;

4038:   PetscFunctionBegin;
4041:   PetscCall(PetscDSGetConstants(prob, &Nc, &constants));
4042:   PetscCall(PetscDSSetConstants(newprob, Nc, (PetscScalar *)constants));
4043:   PetscFunctionReturn(PETSC_SUCCESS);
4044: }

4046: /*@
4047:   PetscDSCopyExactSolutions - Copy all exact solutions to another `PetscDS`

4049:   Not Collective

4051:   Input Parameter:
4052: . ds - The `PetscDS` object

4054:   Output Parameter:
4055: . newds - The `PetscDS` copy

4057:   Level: intermediate

4059: .seealso: `PetscDS`, `PetscDSCopyBoundary()`, `PetscDSCopyEquations()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
4060: @*/
4061: PetscErrorCode PetscDSCopyExactSolutions(PetscDS ds, PetscDS newds)
4062: {
4063:   PetscSimplePointFn *sol;
4064:   void               *ctx;
4065:   PetscInt            Nf, f;

4067:   PetscFunctionBegin;
4070:   PetscCall(PetscDSGetNumFields(ds, &Nf));
4071:   for (f = 0; f < Nf; ++f) {
4072:     PetscCall(PetscDSGetExactSolution(ds, f, &sol, &ctx));
4073:     PetscCall(PetscDSSetExactSolution(newds, f, sol, ctx));
4074:     PetscCall(PetscDSGetExactSolutionTimeDerivative(ds, f, &sol, &ctx));
4075:     PetscCall(PetscDSSetExactSolutionTimeDerivative(newds, f, sol, ctx));
4076:   }
4077:   PetscFunctionReturn(PETSC_SUCCESS);
4078: }

4080: PetscErrorCode PetscDSCopy(PetscDS ds, DM dmNew, PetscDS dsNew)
4081: {
4082:   DSBoundary b;
4083:   PetscInt   cdim, Nf, f, d;
4084:   PetscBool  isCohesive;
4085:   void      *ctx;

4087:   PetscFunctionBegin;
4088:   PetscCall(PetscDSCopyConstants(ds, dsNew));
4089:   PetscCall(PetscDSCopyExactSolutions(ds, dsNew));
4090:   PetscCall(PetscDSSelectDiscretizations(ds, PETSC_DETERMINE, NULL, dsNew));
4091:   PetscCall(PetscDSCopyEquations(ds, dsNew));
4092:   PetscCall(PetscDSGetNumFields(ds, &Nf));
4093:   for (f = 0; f < Nf; ++f) {
4094:     PetscCall(PetscDSGetContext(ds, f, &ctx));
4095:     PetscCall(PetscDSSetContext(dsNew, f, ctx));
4096:     PetscCall(PetscDSGetCohesive(ds, f, &isCohesive));
4097:     PetscCall(PetscDSSetCohesive(dsNew, f, isCohesive));
4098:     PetscCall(PetscDSGetJetDegree(ds, f, &d));
4099:     PetscCall(PetscDSSetJetDegree(dsNew, f, d));
4100:   }
4101:   if (Nf) {
4102:     PetscCall(PetscDSGetCoordinateDimension(ds, &cdim));
4103:     PetscCall(PetscDSSetCoordinateDimension(dsNew, cdim));
4104:   }
4105:   PetscCall(PetscDSCopyBoundary(ds, PETSC_DETERMINE, NULL, dsNew));
4106:   for (b = dsNew->boundary; b; b = b->next) {
4107:     PetscCall(DMGetLabel(dmNew, b->lname, &b->label));
4108:     /* Do not check if label exists here, since p4est calls this for the reference tree which does not have the labels */
4109:     //PetscCheck(b->label,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Label %s missing in new DM", name);
4110:   }
4111:   PetscFunctionReturn(PETSC_SUCCESS);
4112: }

4114: PetscErrorCode PetscDSGetHeightSubspace(PetscDS prob, PetscInt height, PetscDS *subprob)
4115: {
4116:   PetscInt dim, Nf, f;

4118:   PetscFunctionBegin;
4120:   PetscAssertPointer(subprob, 3);
4121:   if (height == 0) {
4122:     *subprob = prob;
4123:     PetscFunctionReturn(PETSC_SUCCESS);
4124:   }
4125:   PetscCall(PetscDSGetNumFields(prob, &Nf));
4126:   PetscCall(PetscDSGetSpatialDimension(prob, &dim));
4127:   PetscCheck(height <= dim, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_OUTOFRANGE, "DS can only handle height in [0, %" PetscInt_FMT "], not %" PetscInt_FMT, dim, height);
4128:   if (!prob->subprobs) PetscCall(PetscCalloc1(dim, &prob->subprobs));
4129:   if (!prob->subprobs[height - 1]) {
4130:     PetscInt cdim;

4132:     PetscCall(PetscDSCreate(PetscObjectComm((PetscObject)prob), &prob->subprobs[height - 1]));
4133:     PetscCall(PetscDSGetCoordinateDimension(prob, &cdim));
4134:     PetscCall(PetscDSSetCoordinateDimension(prob->subprobs[height - 1], cdim));
4135:     for (f = 0; f < Nf; ++f) {
4136:       PetscFE      subfe;
4137:       PetscObject  obj;
4138:       PetscClassId id;

4140:       PetscCall(PetscDSGetDiscretization(prob, f, &obj));
4141:       PetscCall(PetscObjectGetClassId(obj, &id));
4142:       if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetHeightSubspace((PetscFE)obj, height, &subfe));
4143:       else SETERRQ(PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONG, "Unsupported discretization type for field %" PetscInt_FMT, f);
4144:       PetscCall(PetscDSSetDiscretization(prob->subprobs[height - 1], f, (PetscObject)subfe));
4145:     }
4146:   }
4147:   *subprob = prob->subprobs[height - 1];
4148:   PetscFunctionReturn(PETSC_SUCCESS);
4149: }

4151: PetscErrorCode PetscDSPermuteQuadPoint(PetscDS ds, PetscInt ornt, PetscInt field, PetscInt q, PetscInt *qperm)
4152: {
4153:   IS              permIS;
4154:   PetscQuadrature quad;
4155:   DMPolytopeType  ct;
4156:   const PetscInt *perm;
4157:   PetscInt        Na, Nq;

4159:   PetscFunctionBeginHot;
4160:   PetscCall(PetscFEGetQuadrature((PetscFE)ds->disc[field], &quad));
4161:   PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL));
4162:   PetscCall(PetscQuadratureGetCellType(quad, &ct));
4163:   PetscCheck(q >= 0 && q < Nq, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Quadrature point %" PetscInt_FMT " is not in [0, %" PetscInt_FMT ")", q, Nq);
4164:   Na = DMPolytopeTypeGetNumArrangements(ct) / 2;
4165:   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);
4166:   if (!ds->quadPerm[(PetscInt)ct]) PetscCall(PetscQuadratureComputePermutations(quad, NULL, &ds->quadPerm[(PetscInt)ct]));
4167:   permIS = ds->quadPerm[(PetscInt)ct][ornt + Na];
4168:   PetscCall(ISGetIndices(permIS, &perm));
4169:   *qperm = perm[q];
4170:   PetscCall(ISRestoreIndices(permIS, &perm));
4171:   PetscFunctionReturn(PETSC_SUCCESS);
4172: }

4174: PetscErrorCode PetscDSGetDiscType_Internal(PetscDS ds, PetscInt f, PetscDiscType *disctype)
4175: {
4176:   PetscObject  obj;
4177:   PetscClassId id;
4178:   PetscInt     Nf;

4180:   PetscFunctionBegin;
4182:   PetscAssertPointer(disctype, 3);
4183:   *disctype = PETSC_DISC_NONE;
4184:   PetscCall(PetscDSGetNumFields(ds, &Nf));
4185:   PetscCheck(f < Nf, PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_SIZ, "Field %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, Nf);
4186:   PetscCall(PetscDSGetDiscretization(ds, f, &obj));
4187:   if (obj) {
4188:     PetscCall(PetscObjectGetClassId(obj, &id));
4189:     if (id == PETSCFE_CLASSID) *disctype = PETSC_DISC_FE;
4190:     else *disctype = PETSC_DISC_FV;
4191:   }
4192:   PetscFunctionReturn(PETSC_SUCCESS);
4193: }

4195: static PetscErrorCode PetscDSDestroy_Basic(PetscDS ds)
4196: {
4197:   PetscFunctionBegin;
4198:   PetscCall(PetscFree(ds->data));
4199:   PetscFunctionReturn(PETSC_SUCCESS);
4200: }

4202: static PetscErrorCode PetscDSInitialize_Basic(PetscDS ds)
4203: {
4204:   PetscFunctionBegin;
4205:   ds->ops->setfromoptions = NULL;
4206:   ds->ops->setup          = NULL;
4207:   ds->ops->view           = NULL;
4208:   ds->ops->destroy        = PetscDSDestroy_Basic;
4209:   PetscFunctionReturn(PETSC_SUCCESS);
4210: }

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

4215:   Level: intermediate

4217: .seealso: `PetscDSType`, `PetscDSCreate()`, `PetscDSSetType()`
4218: M*/

4220: PETSC_EXTERN PetscErrorCode PetscDSCreate_Basic(PetscDS ds)
4221: {
4222:   PetscDS_Basic *b;

4224:   PetscFunctionBegin;
4226:   PetscCall(PetscNew(&b));
4227:   ds->data = b;

4229:   PetscCall(PetscDSInitialize_Basic(ds));
4230:   PetscFunctionReturn(PETSC_SUCCESS);
4231: }