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: /*@C
  9:   PetscDSRegister - Adds a new `PetscDS` implementation

 11:   Not Collective; No Fortran Support

 13:   Input Parameters:
 14: + sname    - The name of a new user-defined creation routine
 15: - function - The creation routine itself

 17:   Example Usage:
 18: .vb
 19:     PetscDSRegister("my_ds", MyPetscDSCreate);
 20: .ve

 22:   Then, your PetscDS type can be chosen with the procedural interface via
 23: .vb
 24:     PetscDSCreate(MPI_Comm, PetscDS *);
 25:     PetscDSSetType(PetscDS, "my_ds");
 26: .ve
 27:   or at runtime via the option
 28: .vb
 29:     -petscds_type my_ds
 30: .ve

 32:   Level: advanced

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

 37: .seealso: `PetscDSType`, `PetscDS`, `PetscDSRegisterAll()`, `PetscDSRegisterDestroy()`
 38: @*/
 39: PetscErrorCode PetscDSRegister(const char sname[], PetscErrorCode (*function)(PetscDS))
 40: {
 41:   PetscFunctionBegin;
 42:   PetscCall(PetscFunctionListAdd(&PetscDSList, sname, function));
 43:   PetscFunctionReturn(PETSC_SUCCESS);
 44: }

 46: /*@
 47:   PetscDSSetType - Builds a particular `PetscDS`

 49:   Collective; No Fortran Support

 51:   Input Parameters:
 52: + prob - The `PetscDS` object
 53: - name - The `PetscDSType`

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

 58:   Level: intermediate

 60: .seealso: `PetscDSType`, `PetscDS`, `PetscDSGetType()`, `PetscDSCreate()`
 61: @*/
 62: PetscErrorCode PetscDSSetType(PetscDS prob, PetscDSType name)
 63: {
 64:   PetscErrorCode (*r)(PetscDS);
 65:   PetscBool match;

 67:   PetscFunctionBegin;
 69:   PetscCall(PetscObjectTypeCompare((PetscObject)prob, name, &match));
 70:   if (match) PetscFunctionReturn(PETSC_SUCCESS);

 72:   PetscCall(PetscDSRegisterAll());
 73:   PetscCall(PetscFunctionListFind(PetscDSList, name, &r));
 74:   PetscCheck(r, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscDS type: %s", name);

 76:   PetscTryTypeMethod(prob, destroy);
 77:   prob->ops->destroy = NULL;

 79:   PetscCall((*r)(prob));
 80:   PetscCall(PetscObjectChangeTypeName((PetscObject)prob, name));
 81:   PetscFunctionReturn(PETSC_SUCCESS);
 82: }

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

 87:   Not Collective; No Fortran Support

 89:   Input Parameter:
 90: . prob - The `PetscDS`

 92:   Output Parameter:
 93: . name - The `PetscDSType` name

 95:   Level: intermediate

 97: .seealso: `PetscDSType`, `PetscDS`, `PetscDSSetType()`, `PetscDSCreate()`
 98: @*/
 99: PetscErrorCode PetscDSGetType(PetscDS prob, PetscDSType *name)
100: {
101:   PetscFunctionBegin;
103:   PetscAssertPointer(name, 2);
104:   PetscCall(PetscDSRegisterAll());
105:   *name = ((PetscObject)prob)->type_name;
106:   PetscFunctionReturn(PETSC_SUCCESS);
107: }

109: static PetscErrorCode PetscDSView_Ascii(PetscDS ds, PetscViewer viewer)
110: {
111:   PetscViewerFormat  format;
112:   const PetscScalar *constants;
113:   PetscInt           Nf, numConstants, f;

115:   PetscFunctionBegin;
116:   PetscCall(PetscDSGetNumFields(ds, &Nf));
117:   PetscCall(PetscViewerGetFormat(viewer, &format));
118:   PetscCall(PetscViewerASCIIPrintf(viewer, "Discrete System with %" PetscInt_FMT " fields\n", Nf));
119:   PetscCall(PetscViewerASCIIPushTab(viewer));
120:   PetscCall(PetscViewerASCIIPrintf(viewer, "  cell total dim %" PetscInt_FMT " total comp %" PetscInt_FMT "\n", ds->totDim, ds->totComp));
121:   if (ds->isCohesive) PetscCall(PetscViewerASCIIPrintf(viewer, "  cohesive cell\n"));
122:   for (f = 0; f < Nf; ++f) {
123:     DSBoundary      b;
124:     PetscObject     obj;
125:     PetscClassId    id;
126:     PetscQuadrature q;
127:     const char     *name;
128:     PetscInt        Nc, Nq, Nqc;

130:     PetscCall(PetscDSGetDiscretization(ds, f, &obj));
131:     PetscCall(PetscObjectGetClassId(obj, &id));
132:     PetscCall(PetscObjectGetName(obj, &name));
133:     PetscCall(PetscViewerASCIIPrintf(viewer, "Field %s", name ? name : "<unknown>"));
134:     PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
135:     if (id == PETSCFE_CLASSID) {
136:       PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
137:       PetscCall(PetscFEGetQuadrature((PetscFE)obj, &q));
138:       PetscCall(PetscViewerASCIIPrintf(viewer, " FEM"));
139:     } else if (id == PETSCFV_CLASSID) {
140:       PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc));
141:       PetscCall(PetscFVGetQuadrature((PetscFV)obj, &q));
142:       PetscCall(PetscViewerASCIIPrintf(viewer, " FVM"));
143:     } else SETERRQ(PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %" PetscInt_FMT, f);
144:     if (Nc > 1) PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT " components", Nc));
145:     else PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT " component ", Nc));
146:     if (ds->implicit[f]) PetscCall(PetscViewerASCIIPrintf(viewer, " (implicit)"));
147:     else PetscCall(PetscViewerASCIIPrintf(viewer, " (explicit)"));
148:     if (q) {
149:       PetscCall(PetscQuadratureGetData(q, NULL, &Nqc, &Nq, NULL, NULL));
150:       PetscCall(PetscViewerASCIIPrintf(viewer, " (Nq %" PetscInt_FMT " Nqc %" PetscInt_FMT ")", Nq, Nqc));
151:     }
152:     PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT "-jet", ds->jetDegree[f]));
153:     PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
154:     PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
155:     PetscCall(PetscViewerASCIIPushTab(viewer));
156:     if (id == PETSCFE_CLASSID) PetscCall(PetscFEView((PetscFE)obj, viewer));
157:     else if (id == PETSCFV_CLASSID) PetscCall(PetscFVView((PetscFV)obj, viewer));
158:     PetscCall(PetscViewerASCIIPopTab(viewer));

160:     for (b = ds->boundary; b; b = b->next) {
161:       char    *name;
162:       PetscInt c, i;

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

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

224:   Collective

226:   Input Parameters:
227: + A    - the `PetscDS` object
228: . obj  - Optional object that provides the options prefix used in the search of the options database
229: - name - command line option

231:   Level: intermediate

233: .seealso: `PetscDSType`, `PetscDS`, `PetscDSView()`, `PetscObjectViewFromOptions()`, `PetscDSCreate()`
234: @*/
235: PetscErrorCode PetscDSViewFromOptions(PetscDS A, PetscObject obj, const char name[])
236: {
237:   PetscFunctionBegin;
239:   PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
240:   PetscFunctionReturn(PETSC_SUCCESS);
241: }

243: /*@
244:   PetscDSView - Views a `PetscDS`

246:   Collective

248:   Input Parameters:
249: + prob - the `PetscDS` object to view
250: - v    - the viewer

252:   Level: developer

254: .seealso: `PetscDSType`, `PetscDS`, `PetscViewer`, `PetscDSDestroy()`, `PetscDSViewFromOptions()`
255: @*/
256: PetscErrorCode PetscDSView(PetscDS prob, PetscViewer v)
257: {
258:   PetscBool isascii;

260:   PetscFunctionBegin;
262:   if (!v) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)prob), &v));
264:   PetscCall(PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &isascii));
265:   if (isascii) PetscCall(PetscDSView_Ascii(prob, v));
266:   PetscTryTypeMethod(prob, view, v);
267:   PetscFunctionReturn(PETSC_SUCCESS);
268: }

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

273:   Collective

275:   Input Parameter:
276: . prob - the `PetscDS` object to set options for

278:   Options Database Keys:
279: + -petscds_type <type>     - Set the `PetscDS` type
280: . -petscds_view <view opt> - View the `PetscDS`
281: . -petscds_jac_pre         - Turn formation of a separate Jacobian preconditioner on or off
282: . -bc_<name> <ids>         - Specify a list of label ids for a boundary condition
283: - -bc_<name>_comp <comps>  - Specify a list of field components to constrain for a boundary condition

285:   Level: intermediate

287: .seealso: `PetscDS`, `PetscDSView()`
288: @*/
289: PetscErrorCode PetscDSSetFromOptions(PetscDS prob)
290: {
291:   DSBoundary  b;
292:   const char *defaultType;
293:   char        name[256];
294:   PetscBool   flg;

296:   PetscFunctionBegin;
298:   if (!((PetscObject)prob)->type_name) {
299:     defaultType = PETSCDSBASIC;
300:   } else {
301:     defaultType = ((PetscObject)prob)->type_name;
302:   }
303:   PetscCall(PetscDSRegisterAll());

305:   PetscObjectOptionsBegin((PetscObject)prob);
306:   for (b = prob->boundary; b; b = b->next) {
307:     char      optname[1024];
308:     PetscInt  ids[1024], len = 1024;
309:     PetscBool flg;

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

349: /*@
350:   PetscDSSetUp - Construct data structures for the `PetscDS`

352:   Collective

354:   Input Parameter:
355: . prob - the `PetscDS` object to setup

357:   Level: developer

359: .seealso: `PetscDS`, `PetscDSView()`, `PetscDSDestroy()`
360: @*/
361: PetscErrorCode PetscDSSetUp(PetscDS prob)
362: {
363:   const PetscInt Nf          = prob->Nf;
364:   PetscBool      hasH        = PETSC_FALSE;
365:   PetscInt       maxOrder[4] = {-2, -2, -2, -2};
366:   PetscInt       dim, dimEmbed, NbMax = 0, NcMax = 0, NqMax = 0, NsMax = 1, f;

368:   PetscFunctionBegin;
370:   if (prob->setup) PetscFunctionReturn(PETSC_SUCCESS);
371:   /* Calculate sizes */
372:   PetscCall(PetscDSGetSpatialDimension(prob, &dim));
373:   PetscCall(PetscDSGetCoordinateDimension(prob, &dimEmbed));
374:   prob->totDim = prob->totComp = 0;
375:   PetscCall(PetscMalloc2(Nf, &prob->Nc, Nf, &prob->Nb));
376:   PetscCall(PetscCalloc2(Nf + 1, &prob->off, Nf + 1, &prob->offDer));
377:   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]));
378:   PetscCall(PetscMalloc2(Nf, &prob->T, Nf, &prob->Tf));
379:   if (prob->forceQuad) {
380:     // Note: This assumes we have one kind of cell at each dimension.
381:     //       We can fix this by having quadrature hold the celltype
382:     PetscQuadrature maxQuad[4] = {NULL, NULL, NULL, NULL};

384:     for (f = 0; f < Nf; ++f) {
385:       PetscObject     obj;
386:       PetscClassId    id;
387:       PetscQuadrature q = NULL, fq = NULL;
388:       PetscInt        dim = -1, order = -1, forder = -1;

390:       PetscCall(PetscDSGetDiscretization(prob, f, &obj));
391:       if (!obj) continue;
392:       PetscCall(PetscObjectGetClassId(obj, &id));
393:       if (id == PETSCFE_CLASSID) {
394:         PetscFE fe = (PetscFE)obj;

396:         PetscCall(PetscFEGetQuadrature(fe, &q));
397:         PetscCall(PetscFEGetFaceQuadrature(fe, &fq));
398:       } else if (id == PETSCFV_CLASSID) {
399:         PetscFV fv = (PetscFV)obj;

401:         PetscCall(PetscFVGetQuadrature(fv, &q));
402:       }
403:       if (q) {
404:         PetscCall(PetscQuadratureGetData(q, &dim, NULL, NULL, NULL, NULL));
405:         PetscCall(PetscQuadratureGetOrder(q, &order));
406:         if (order > maxOrder[dim]) {
407:           maxOrder[dim] = order;
408:           maxQuad[dim]  = q;
409:         }
410:       }
411:       if (fq) {
412:         PetscCall(PetscQuadratureGetData(fq, &dim, NULL, NULL, NULL, NULL));
413:         PetscCall(PetscQuadratureGetOrder(fq, &forder));
414:         if (forder > maxOrder[dim]) {
415:           maxOrder[dim] = forder;
416:           maxQuad[dim]  = fq;
417:         }
418:       }
419:     }
420:     for (f = 0; f < Nf; ++f) {
421:       PetscObject     obj;
422:       PetscClassId    id;
423:       PetscQuadrature q;
424:       PetscInt        dim;

426:       PetscCall(PetscDSGetDiscretization(prob, f, &obj));
427:       if (!obj) continue;
428:       PetscCall(PetscObjectGetClassId(obj, &id));
429:       if (id == PETSCFE_CLASSID) {
430:         PetscFE fe = (PetscFE)obj;

432:         PetscCall(PetscFEGetQuadrature(fe, &q));
433:         PetscCall(PetscQuadratureGetData(q, &dim, NULL, NULL, NULL, NULL));
434:         PetscCall(PetscFESetQuadrature(fe, maxQuad[dim]));
435:         PetscCall(PetscFESetFaceQuadrature(fe, dim ? maxQuad[dim - 1] : NULL));
436:       } else if (id == PETSCFV_CLASSID) {
437:         PetscFV fv = (PetscFV)obj;

439:         PetscCall(PetscFVGetQuadrature(fv, &q));
440:         PetscCall(PetscQuadratureGetData(q, &dim, NULL, NULL, NULL, NULL));
441:         PetscCall(PetscFVSetQuadrature(fv, maxQuad[dim]));
442:       }
443:     }
444:   }
445:   for (f = 0; f < Nf; ++f) {
446:     PetscObject     obj;
447:     PetscClassId    id;
448:     PetscQuadrature q  = NULL;
449:     PetscInt        Nq = 0, Nb, Nc;

451:     PetscCall(PetscDSGetDiscretization(prob, f, &obj));
452:     if (prob->jetDegree[f] > 1) hasH = PETSC_TRUE;
453:     if (!obj) {
454:       /* Empty mesh */
455:       Nb = Nc    = 0;
456:       prob->T[f] = prob->Tf[f] = NULL;
457:     } else {
458:       PetscCall(PetscObjectGetClassId(obj, &id));
459:       if (id == PETSCFE_CLASSID) {
460:         PetscFE fe = (PetscFE)obj;

462:         PetscCall(PetscFEGetQuadrature(fe, &q));
463:         {
464:           PetscQuadrature fq;
465:           PetscInt        dim, order;

467:           PetscCall(PetscQuadratureGetData(q, &dim, NULL, NULL, NULL, NULL));
468:           PetscCall(PetscQuadratureGetOrder(q, &order));
469:           if (maxOrder[dim] < 0) maxOrder[dim] = order;
470:           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]);
471:           PetscCall(PetscFEGetFaceQuadrature(fe, &fq));
472:           if (fq) {
473:             PetscCall(PetscQuadratureGetData(fq, &dim, NULL, NULL, NULL, NULL));
474:             PetscCall(PetscQuadratureGetOrder(fq, &order));
475:             if (maxOrder[dim] < 0) maxOrder[dim] = order;
476:             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]);
477:           }
478:         }
479:         PetscCall(PetscFEGetDimension(fe, &Nb));
480:         PetscCall(PetscFEGetNumComponents(fe, &Nc));
481:         PetscCall(PetscFEGetCellTabulation(fe, prob->jetDegree[f], &prob->T[f]));
482:         PetscCall(PetscFEGetFaceTabulation(fe, prob->jetDegree[f], &prob->Tf[f]));
483:       } else if (id == PETSCFV_CLASSID) {
484:         PetscFV fv = (PetscFV)obj;

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

525: static PetscErrorCode PetscDSDestroyStructs_Static(PetscDS prob)
526: {
527:   PetscFunctionBegin;
528:   PetscCall(PetscFree2(prob->Nc, prob->Nb));
529:   PetscCall(PetscFree2(prob->off, prob->offDer));
530:   PetscCall(PetscFree6(prob->offCohesive[0], prob->offCohesive[1], prob->offCohesive[2], prob->offDerCohesive[0], prob->offDerCohesive[1], prob->offDerCohesive[2]));
531:   PetscCall(PetscFree2(prob->T, prob->Tf));
532:   PetscCall(PetscFree3(prob->u, prob->u_t, prob->u_x));
533:   PetscCall(PetscFree5(prob->x, prob->basisReal, prob->basisDerReal, prob->testReal, prob->testDerReal));
534:   PetscCall(PetscFree6(prob->f0, prob->f1, prob->g0, prob->g1, prob->g2, prob->g3));
535:   PetscFunctionReturn(PETSC_SUCCESS);
536: }

538: static PetscErrorCode PetscDSEnlarge_Static(PetscDS prob, PetscInt NfNew)
539: {
540:   PetscObject         *tmpd;
541:   PetscBool           *tmpi;
542:   PetscInt            *tmpk;
543:   PetscBool           *tmpc;
544:   PetscPointFn       **tmpup;
545:   PetscSimplePointFn **tmpexactSol, **tmpexactSol_t, **tmplowerBound, **tmpupperBound;
546:   void               **tmpexactCtx, **tmpexactCtx_t, **tmplowerCtx, **tmpupperCtx;
547:   void               **tmpctx;
548:   PetscInt             Nf = prob->Nf, f;

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

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

615:   Collective

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

620:   Level: developer

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

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

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

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

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

668:   Collective

670:   Input Parameter:
671: . comm - The communicator for the `PetscDS` object

673:   Output Parameter:
674: . ds - The `PetscDS` object

676:   Level: beginner

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

684:   PetscFunctionBegin;
685:   PetscAssertPointer(ds, 2);
686:   PetscCall(PetscDSInitializePackage());

688:   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->numFuncConstants = 3; // Row and col fields, cell size
693:   p->dimEmbed         = -1;
694:   p->useJacPre        = PETSC_TRUE;
695:   p->forceQuad        = PETSC_TRUE;
696:   PetscCall(PetscMalloc1(p->numConstants + p->numFuncConstants, &p->constants));
697:   PetscCall(PetscWeakFormCreate(comm, &p->wf));
698:   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, this can be a `PetscFE` or a `PetscFV`

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, this can be a `PetscFE` or a `PetscFV`

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 from within the `PetscDS`

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 to be used by the `PetscDS`

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 discretization object, this can be a `PetscFE` or `PetscFV`

1131:   Level: beginner

1133: .seealso: `PetscWeakForm`, `PetscFE`, `PetscFV`, `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 `TSARKIMEX`

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: `TSARKIMEX`, `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 `TSARKIMEX`

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: `TSARKIMEX`, `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 that was provided with `PetscDSSetObjective()`

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, see `PetscPointFn`

1283:   Level: intermediate

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

1288: .seealso: `PetscPointFn`, `PetscDS`, `PetscDSSetObjective()`, `PetscDSGetResidual()`
1289: @*/
1290: PetscErrorCode PetscDSGetObjective(PetscDS ds, PetscInt f, PetscPointFn **obj)
1291: {
1292:   PetscPointFn **tmp;
1293:   PetscInt       n;

1295:   PetscFunctionBegin;
1297:   PetscAssertPointer(obj, 3);
1298:   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);
1299:   PetscCall(PetscWeakFormGetObjective(ds->wf, NULL, 0, f, 0, &n, &tmp));
1300:   *obj = tmp ? tmp[0] : NULL;
1301:   PetscFunctionReturn(PETSC_SUCCESS);
1302: }

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

1307:   Not Collective

1309:   Input Parameters:
1310: + ds  - The `PetscDS`
1311: . f   - The test field number
1312: - obj - integrand for the test function term, see `PetscPointFn`

1314:   Level: intermediate

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

1319: .seealso: `PetscPointFn`, `PetscDS`, `PetscDSGetObjective()`, `PetscDSSetResidual()`
1320: @*/
1321: PetscErrorCode PetscDSSetObjective(PetscDS ds, PetscInt f, PetscPointFn *obj)
1322: {
1323:   PetscFunctionBegin;
1326:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
1327:   PetscCall(PetscWeakFormSetIndexObjective(ds->wf, NULL, 0, f, 0, 0, obj));
1328:   PetscFunctionReturn(PETSC_SUCCESS);
1329: }

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

1334:   Not Collective

1336:   Input Parameters:
1337: + ds - The `PetscDS`
1338: - f  - The test field number

1340:   Output Parameters:
1341: + f0 - integrand for the test function term, see `PetscPointFn`
1342: - f1 - integrand for the test function gradient term, see `PetscPointFn`

1344:   Level: intermediate

1346:   Note:
1347:   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)$

1349: .seealso: `PetscPointFn`, `PetscDS`, `PetscDSSetResidual()`
1350: @*/
1351: PetscErrorCode PetscDSGetResidual(PetscDS ds, PetscInt f, PetscPointFn **f0, PetscPointFn **f1)
1352: {
1353:   PetscPointFn **tmp0, **tmp1;
1354:   PetscInt       n0, n1;

1356:   PetscFunctionBegin;
1358:   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);
1359:   PetscCall(PetscWeakFormGetResidual(ds->wf, NULL, 0, f, 0, &n0, &tmp0, &n1, &tmp1));
1360:   *f0 = tmp0 ? tmp0[0] : NULL;
1361:   *f1 = tmp1 ? tmp1[0] : NULL;
1362:   PetscFunctionReturn(PETSC_SUCCESS);
1363: }

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

1368:   Not Collective

1370:   Input Parameters:
1371: + ds - The `PetscDS`
1372: . f  - The test field number
1373: . f0 - integrand for the test function term, see `PetscPointFn`
1374: - f1 - integrand for the test function gradient term, see `PetscPointFn`

1376:   Level: intermediate

1378:   Note:
1379:   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)$

1381: .seealso: `PetscPointFn`, `PetscDS`, `PetscDSGetResidual()`
1382: @*/
1383: PetscErrorCode PetscDSSetResidual(PetscDS ds, PetscInt f, PetscPointFn *f0, PetscPointFn *f1)
1384: {
1385:   PetscFunctionBegin;
1389:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
1390:   PetscCall(PetscWeakFormSetIndexResidual(ds->wf, NULL, 0, f, 0, 0, f0, 0, f1));
1391:   PetscFunctionReturn(PETSC_SUCCESS);
1392: }

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

1397:   Not Collective

1399:   Input Parameters:
1400: + ds - The `PetscDS`
1401: - f  - The test field number

1403:   Output Parameters:
1404: + f0 - integrand for the test function term, see `PetscPointFn`
1405: - f1 - integrand for the test function gradient term, see `PetscPointFn`

1407:   Level: intermediate

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

1412: .seealso: `PetscPointFn`, `PetscDS`, `PetscDSSetRHSResidual()`
1413: @*/
1414: PetscErrorCode PetscDSGetRHSResidual(PetscDS ds, PetscInt f, PetscPointFn **f0, PetscPointFn **f1)
1415: {
1416:   PetscPointFn **tmp0, **tmp1;
1417:   PetscInt       n0, n1;

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

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

1431:   Not Collective

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

1439:   Level: intermediate

1441:   Note:
1442:   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)$

1444: .seealso: `PetscDS`, `PetscDSGetResidual()`
1445: @*/
1446: PetscErrorCode PetscDSSetRHSResidual(PetscDS ds, PetscInt f, PetscPointFn *f0, PetscPointFn *f1)
1447: {
1448:   PetscFunctionBegin;
1452:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
1453:   PetscCall(PetscWeakFormSetIndexResidual(ds->wf, NULL, 0, f, 100, 0, f0, 0, f1));
1454:   PetscFunctionReturn(PETSC_SUCCESS);
1455: }

1457: /*@
1458:   PetscDSHasJacobian - Checks that the Jacobian functions have been set

1460:   Not Collective

1462:   Input Parameter:
1463: . ds - The `PetscDS`

1465:   Output Parameter:
1466: . hasJac - flag that indicates the pointwise function for the Jacobian has been set

1468:   Level: intermediate

1470: .seealso: `PetscDS`, `PetscDSGetJacobianPreconditioner()`, `PetscDSSetJacobianPreconditioner()`, `PetscDSGetJacobian()`
1471: @*/
1472: PetscErrorCode PetscDSHasJacobian(PetscDS ds, PetscBool *hasJac)
1473: {
1474:   PetscFunctionBegin;
1476:   PetscCall(PetscWeakFormHasJacobian(ds->wf, hasJac));
1477:   PetscFunctionReturn(PETSC_SUCCESS);
1478: }

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

1483:   Not Collective

1485:   Input Parameters:
1486: + ds - The `PetscDS`
1487: . f  - The test field number
1488: - g  - The field number

1490:   Output Parameters:
1491: + g0 - integrand for the test and basis function term, see `PetscPointJacFn`
1492: . g1 - integrand for the test function and basis function gradient term, see `PetscPointJacFn`
1493: . g2 - integrand for the test function gradient and basis function term, see `PetscPointJacFn`
1494: - g3 - integrand for the test function gradient and basis function gradient term, see `PetscPointJacFn`

1496:   Level: intermediate

1498:   Note:
1499:   We are using a first order FEM model for the weak form\:

1501:   $$
1502:   \int_\Omega \phi\, g_0(u, u_t, \nabla u, x, t) \psi + \phi\, {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi
1503:   + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi
1504:   $$

1506: .seealso: `PetscDS`, `PetscDSSetJacobian()`, `PetscPointJacFn`
1507: @*/
1508: PetscErrorCode PetscDSGetJacobian(PetscDS ds, PetscInt f, PetscInt g, PetscPointJacFn **g0, PetscPointJacFn **g1, PetscPointJacFn **g2, PetscPointJacFn **g3)
1509: {
1510:   PetscPointJacFn **tmp0, **tmp1, **tmp2, **tmp3;
1511:   PetscInt          n0, n1, n2, n3;

1513:   PetscFunctionBegin;
1515:   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);
1516:   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);
1517:   PetscCall(PetscWeakFormGetJacobian(ds->wf, NULL, 0, f, g, 0, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3));
1518:   *g0 = tmp0 ? tmp0[0] : NULL;
1519:   *g1 = tmp1 ? tmp1[0] : NULL;
1520:   *g2 = tmp2 ? tmp2[0] : NULL;
1521:   *g3 = tmp3 ? tmp3[0] : NULL;
1522:   PetscFunctionReturn(PETSC_SUCCESS);
1523: }

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

1528:   Not Collective

1530:   Input Parameters:
1531: + ds - The `PetscDS`
1532: . f  - The test field number
1533: . g  - The field number
1534: . g0 - integrand for the test and basis function term, see `PetscPointJacFn`
1535: . g1 - integrand for the test function and basis function gradient term, see `PetscPointJacFn`
1536: . g2 - integrand for the test function gradient and basis function term, see `PetscPointJacFn`
1537: - g3 - integrand for the test function gradient and basis function gradient term, see `PetscPointJacFn`

1539:   Level: intermediate

1541:   Note:
1542:   We are using a first order FEM model for the weak form\:

1544:   $$
1545:   \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
1546:   + \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
1547:   $$

1549: .seealso: `PetscDS`, `PetscDSGetJacobian()`, `PetscPointJacFn`
1550: @*/
1551: PetscErrorCode PetscDSSetJacobian(PetscDS ds, PetscInt f, PetscInt g, PetscPointJacFn *g0, PetscPointJacFn *g1, PetscPointJacFn *g2, PetscPointJacFn *g3)
1552: {
1553:   PetscFunctionBegin;
1559:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
1560:   PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", g);
1561:   PetscCall(PetscWeakFormSetIndexJacobian(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3));
1562:   PetscFunctionReturn(PETSC_SUCCESS);
1563: }

1565: /*@
1566:   PetscDSUseJacobianPreconditioner - Set whether to construct a Jacobian preconditioner

1568:   Not Collective

1570:   Input Parameters:
1571: + prob      - The `PetscDS`
1572: - useJacPre - flag that enables construction of a Jacobian preconditioner

1574:   Level: intermediate

1576:   Developer Note:
1577:   Should be called `PetscDSSetUseJacobianPreconditioner()`

1579: .seealso: `PetscDS`, `PetscDSGetJacobianPreconditioner()`, `PetscDSSetJacobianPreconditioner()`, `PetscDSGetJacobian()`
1580: @*/
1581: PetscErrorCode PetscDSUseJacobianPreconditioner(PetscDS prob, PetscBool useJacPre)
1582: {
1583:   PetscFunctionBegin;
1585:   prob->useJacPre = useJacPre;
1586:   PetscFunctionReturn(PETSC_SUCCESS);
1587: }

1589: /*@
1590:   PetscDSHasJacobianPreconditioner - Checks if a Jacobian matrix for constructing a preconditioner has been set

1592:   Not Collective

1594:   Input Parameter:
1595: . ds - The `PetscDS`

1597:   Output Parameter:
1598: . hasJacPre - the flag

1600:   Level: intermediate

1602: .seealso: `PetscDS`, `PetscDSGetJacobianPreconditioner()`, `PetscDSSetJacobianPreconditioner()`, `PetscDSGetJacobian()`
1603: @*/
1604: PetscErrorCode PetscDSHasJacobianPreconditioner(PetscDS ds, PetscBool *hasJacPre)
1605: {
1606:   PetscFunctionBegin;
1608:   *hasJacPre = PETSC_FALSE;
1609:   if (!ds->useJacPre) PetscFunctionReturn(PETSC_SUCCESS);
1610:   PetscCall(PetscWeakFormHasJacobianPreconditioner(ds->wf, hasJacPre));
1611:   PetscFunctionReturn(PETSC_SUCCESS);
1612: }

1614: /*@C
1615:   PetscDSGetJacobianPreconditioner - Get the pointwise Jacobian function for given test and basis field that constructs the matrix used
1616:   to compute the preconditioner. If this is missing, the system matrix is used to build the preconditioner.

1618:   Not Collective

1620:   Input Parameters:
1621: + ds - The `PetscDS`
1622: . f  - The test field number
1623: - g  - The field number

1625:   Output Parameters:
1626: + g0 - integrand for the test and basis function term, see `PetscPointJacFn`
1627: . g1 - integrand for the test function and basis function gradient term, see `PetscPointJacFn`
1628: . g2 - integrand for the test function gradient and basis function term, see `PetscPointJacFn`
1629: - g3 - integrand for the test function gradient and basis function gradient term, see `PetscPointJacFn`

1631:   Level: intermediate

1633:   Note:
1634:   We are using a first order FEM model for the weak form\:

1636:   $$
1637:   \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
1638:   + \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
1639:   $$

1641:   Developer Note:
1642:   The name is confusing since the function computes a matrix used to construct the preconditioner, not a preconditioner.

1644: .seealso: `PetscDS`, `PetscDSSetJacobianPreconditioner()`, `PetscDSGetJacobian()`, `PetscPointJacFn`
1645: @*/
1646: PetscErrorCode PetscDSGetJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g, PetscPointJacFn **g0, PetscPointJacFn **g1, PetscPointJacFn **g2, PetscPointJacFn **g3)
1647: {
1648:   PetscPointJacFn **tmp0, **tmp1, **tmp2, **tmp3;
1649:   PetscInt          n0, n1, n2, n3;

1651:   PetscFunctionBegin;
1653:   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);
1654:   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);
1655:   PetscCall(PetscWeakFormGetJacobianPreconditioner(ds->wf, NULL, 0, f, g, 0, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3));
1656:   *g0 = tmp0 ? tmp0[0] : NULL;
1657:   *g1 = tmp1 ? tmp1[0] : NULL;
1658:   *g2 = tmp2 ? tmp2[0] : NULL;
1659:   *g3 = tmp3 ? tmp3[0] : NULL;
1660:   PetscFunctionReturn(PETSC_SUCCESS);
1661: }

1663: /*@C
1664:   PetscDSSetJacobianPreconditioner - Set the pointwise Jacobian function for given test and basis fields that constructs the matrix used
1665:   to compute the preconditioner. If this is missing, the system matrix is used to build the preconditioner.

1667:   Not Collective

1669:   Input Parameters:
1670: + ds - The `PetscDS`
1671: . f  - The test field number
1672: . g  - The field number
1673: . g0 - integrand for the test and basis function term, see `PetscPointJacFn`
1674: . g1 - integrand for the test function and basis function gradient term, see `PetscPointJacFn`
1675: . g2 - integrand for the test function gradient and basis function term, see `PetscPointJacFn`
1676: - g3 - integrand for the test function gradient and basis function gradient term, see `PetscPointJacFn`

1678:   Level: intermediate

1680:   Note:
1681:   We are using a first order FEM model for the weak form\:

1683:   $$
1684:   \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
1685:   + \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
1686:   $$

1688:   Developer Note:
1689:   The name is confusing since the function computes a matrix used to construct the preconditioner, not a preconditioner.

1691: .seealso: `PetscDS`, `PetscDSGetJacobianPreconditioner()`, `PetscDSSetJacobian()`, `PetscPointJacFn`
1692: @*/
1693: PetscErrorCode PetscDSSetJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g, PetscPointJacFn *g0, PetscPointJacFn *g1, PetscPointJacFn *g2, PetscPointJacFn *g3)
1694: {
1695:   PetscFunctionBegin;
1701:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
1702:   PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", g);
1703:   PetscCall(PetscWeakFormSetIndexJacobianPreconditioner(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3));
1704:   PetscFunctionReturn(PETSC_SUCCESS);
1705: }

1707: /*@
1708:   PetscDSHasDynamicJacobian - Signals that a dynamic Jacobian, $dF/du_t$, has been set

1710:   Not Collective

1712:   Input Parameter:
1713: . ds - The `PetscDS`

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

1718:   Level: intermediate

1720: .seealso: `PetscDS`, `PetscDSGetDynamicJacobian()`, `PetscDSSetDynamicJacobian()`, `PetscDSGetJacobian()`
1721: @*/
1722: PetscErrorCode PetscDSHasDynamicJacobian(PetscDS ds, PetscBool *hasDynJac)
1723: {
1724:   PetscFunctionBegin;
1726:   PetscCall(PetscWeakFormHasDynamicJacobian(ds->wf, hasDynJac));
1727:   PetscFunctionReturn(PETSC_SUCCESS);
1728: }

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

1733:   Not Collective

1735:   Input Parameters:
1736: + ds - The `PetscDS`
1737: . f  - The test field number
1738: - g  - The field number

1740:   Output Parameters:
1741: + g0 - integrand for the test and basis function term, see `PetscPointJacFn`
1742: . g1 - integrand for the test function and basis function gradient term, see `PetscPointJacFn`
1743: . g2 - integrand for the test function gradient and basis function term, see `PetscPointJacFn`
1744: - g3 - integrand for the test function gradient and basis function gradient term, see `PetscPointJacFn`

1746:   Level: intermediate

1748:   Note:
1749:   We are using a first order FEM model for the weak form\:

1751:   $$
1752:   \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
1753:   + \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
1754:   $$

1756: .seealso: `PetscDS`, `PetscDSSetJacobian()`, `PetscDSSetDynamicJacobian()`, `PetscPointJacFn`
1757: @*/
1758: PetscErrorCode PetscDSGetDynamicJacobian(PetscDS ds, PetscInt f, PetscInt g, PetscPointJacFn **g0, PetscPointJacFn **g1, PetscPointJacFn **g2, PetscPointJacFn **g3)
1759: {
1760:   PetscPointJacFn **tmp0, **tmp1, **tmp2, **tmp3;
1761:   PetscInt          n0, n1, n2, n3;

1763:   PetscFunctionBegin;
1765:   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);
1766:   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);
1767:   PetscCall(PetscWeakFormGetDynamicJacobian(ds->wf, NULL, 0, f, g, 0, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3));
1768:   *g0 = tmp0 ? tmp0[0] : NULL;
1769:   *g1 = tmp1 ? tmp1[0] : NULL;
1770:   *g2 = tmp2 ? tmp2[0] : NULL;
1771:   *g3 = tmp3 ? tmp3[0] : NULL;
1772:   PetscFunctionReturn(PETSC_SUCCESS);
1773: }

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

1778:   Not Collective

1780:   Input Parameters:
1781: + ds - The `PetscDS`
1782: . f  - The test field number
1783: . g  - The field number
1784: . g0 - integrand for the test and basis function term, see `PetscPointJacFn`
1785: . g1 - integrand for the test function and basis function gradient term, see `PetscPointJacFn`
1786: . g2 - integrand for the test function gradient and basis function term, see `PetscPointJacFn`
1787: - g3 - integrand for the test function gradient and basis function gradient term, see `PetscPointJacFn`

1789:   Level: intermediate

1791:   Note:
1792:   We are using a first order FEM model for the weak form\:

1794:   $$
1795:   \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
1796:   + \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
1797:   $$

1799: .seealso: `PetscDS`, `PetscDSGetDynamicJacobian()`, `PetscDSGetJacobian()`, `PetscPointJacFn`
1800: @*/
1801: PetscErrorCode PetscDSSetDynamicJacobian(PetscDS ds, PetscInt f, PetscInt g, PetscPointJacFn *g0, PetscPointJacFn *g1, PetscPointJacFn *g2, PetscPointJacFn *g3)
1802: {
1803:   PetscFunctionBegin;
1809:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
1810:   PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", g);
1811:   PetscCall(PetscWeakFormSetIndexDynamicJacobian(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3));
1812:   PetscFunctionReturn(PETSC_SUCCESS);
1813: }

1815: /*@C
1816:   PetscDSGetRiemannSolver - Returns the Riemann solver for the given field

1818:   Not Collective

1820:   Input Parameters:
1821: + ds - The `PetscDS` object
1822: - f  - The field number

1824:   Output Parameter:
1825: . r - Riemann solver, see `PetscRiemannFn`

1827:   Level: intermediate

1829: .seealso: `PetscDS`, `PetscRiemannFn`, `PetscDSSetRiemannSolver()`
1830: @*/
1831: PetscErrorCode PetscDSGetRiemannSolver(PetscDS ds, PetscInt f, PetscRiemannFn **r)
1832: {
1833:   PetscRiemannFn **tmp;
1834:   PetscInt         n;

1836:   PetscFunctionBegin;
1838:   PetscAssertPointer(r, 3);
1839:   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);
1840:   PetscCall(PetscWeakFormGetRiemannSolver(ds->wf, NULL, 0, f, 0, &n, &tmp));
1841:   *r = tmp ? tmp[0] : NULL;
1842:   PetscFunctionReturn(PETSC_SUCCESS);
1843: }

1845: /*@C
1846:   PetscDSSetRiemannSolver - Sets the Riemann solver for the given field

1848:   Not Collective

1850:   Input Parameters:
1851: + ds - The `PetscDS` object
1852: . f  - The field number
1853: - r  - Riemann solver, see `PetscRiemannFn`

1855:   Level: intermediate

1857: .seealso: `PetscDS`, `PetscRiemannFn`, `PetscDSGetRiemannSolver()`
1858: @*/
1859: PetscErrorCode PetscDSSetRiemannSolver(PetscDS ds, PetscInt f, PetscRiemannFn *r)
1860: {
1861:   PetscFunctionBegin;
1864:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
1865:   PetscCall(PetscWeakFormSetIndexRiemannSolver(ds->wf, NULL, 0, f, 0, 0, r));
1866:   PetscFunctionReturn(PETSC_SUCCESS);
1867: }

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

1872:   Not Collective

1874:   Input Parameters:
1875: + ds - The `PetscDS`
1876: - f  - The field number

1878:   Output Parameter:
1879: . update - update function, see `PetscPointFn`

1881:   Level: intermediate

1883: .seealso: `PetscDS`, `PetscPointFn`, `PetscDSSetUpdate()`, `PetscDSSetResidual()`
1884: @*/
1885: PetscErrorCode PetscDSGetUpdate(PetscDS ds, PetscInt f, PetscPointFn **update)
1886: {
1887:   PetscFunctionBegin;
1889:   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);
1890:   if (update) {
1891:     PetscAssertPointer(update, 3);
1892:     *update = ds->update[f];
1893:   }
1894:   PetscFunctionReturn(PETSC_SUCCESS);
1895: }

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

1900:   Not Collective

1902:   Input Parameters:
1903: + ds     - The `PetscDS`
1904: . f      - The field number
1905: - update - update function, see `PetscPointFn`

1907:   Level: intermediate

1909: .seealso: `PetscDS`, `PetscPointFn`, `PetscDSGetResidual()`
1910: @*/
1911: PetscErrorCode PetscDSSetUpdate(PetscDS ds, PetscInt f, PetscPointFn *update)
1912: {
1913:   PetscFunctionBegin;
1916:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
1917:   PetscCall(PetscDSEnlarge_Static(ds, f + 1));
1918:   ds->update[f] = update;
1919:   PetscFunctionReturn(PETSC_SUCCESS);
1920: }

1922: /*@C
1923:   PetscDSGetContext - Returns the context that was passed by `PetscDSSetContext()`

1925:   Not Collective

1927:   Input Parameters:
1928: + ds  - The `PetscDS`
1929: . f   - The field number
1930: - ctx - the context

1932:   Level: intermediate

1934: .seealso: `PetscDS`, `PetscPointFn`, `PetscDSSetContext()`
1935: @*/
1936: PetscErrorCode PetscDSGetContext(PetscDS ds, PetscInt f, void *ctx)
1937: {
1938:   PetscFunctionBegin;
1940:   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);
1941:   PetscAssertPointer(ctx, 3);
1942:   *(void **)ctx = ds->ctx[f];
1943:   PetscFunctionReturn(PETSC_SUCCESS);
1944: }

1946: /*@C
1947:   PetscDSSetContext - Sets the context that is passed back to some of the pointwise function callbacks used by this `PetscDS`

1949:   Not Collective

1951:   Input Parameters:
1952: + ds  - The `PetscDS`
1953: . f   - The field number
1954: - ctx - the context

1956:   Level: intermediate

1958: .seealso: `PetscDS`, `PetscPointFn`, `PetscDSGetContext()`
1959: @*/
1960: PetscErrorCode PetscDSSetContext(PetscDS ds, PetscInt f, void *ctx)
1961: {
1962:   PetscFunctionBegin;
1964:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
1965:   PetscCall(PetscDSEnlarge_Static(ds, f + 1));
1966:   ds->ctx[f] = ctx;
1967:   PetscFunctionReturn(PETSC_SUCCESS);
1968: }

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

1973:   Not Collective

1975:   Input Parameters:
1976: + ds - The PetscDS
1977: - f  - The test field number

1979:   Output Parameters:
1980: + f0 - boundary integrand for the test function term, see `PetscBdPointFn`
1981: - f1 - boundary integrand for the test function gradient term, see `PetscBdPointFn`

1983:   Level: intermediate

1985:   Note:
1986:   We are using a first order FEM model for the weak form\:

1988:   $$
1989:   \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
1990:   $$

1992: .seealso: `PetscDS`, `PetscBdPointFn`, `PetscDSSetBdResidual()`
1993: @*/
1994: PetscErrorCode PetscDSGetBdResidual(PetscDS ds, PetscInt f, PetscBdPointFn **f0, PetscBdPointFn **f1)
1995: {
1996:   PetscBdPointFn **tmp0, **tmp1;
1997:   PetscInt         n0, n1;

1999:   PetscFunctionBegin;
2001:   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);
2002:   PetscCall(PetscWeakFormGetBdResidual(ds->wf, NULL, 0, f, 0, &n0, &tmp0, &n1, &tmp1));
2003:   *f0 = tmp0 ? tmp0[0] : NULL;
2004:   *f1 = tmp1 ? tmp1[0] : NULL;
2005:   PetscFunctionReturn(PETSC_SUCCESS);
2006: }

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

2011:   Not Collective

2013:   Input Parameters:
2014: + ds - The `PetscDS`
2015: . f  - The test field number
2016: . f0 - boundary integrand for the test function term, see `PetscBdPointFn`
2017: - f1 - boundary integrand for the test function gradient term, see `PetscBdPointFn`

2019:   Level: intermediate

2021:   Note:
2022:   We are using a first order FEM model for the weak form\:

2024:   $$
2025:   \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
2026:   $$

2028: .seealso: `PetscDS`, `PetscBdPointFn`, `PetscDSGetBdResidual()`
2029: @*/
2030: PetscErrorCode PetscDSSetBdResidual(PetscDS ds, PetscInt f, PetscBdPointFn *f0, PetscBdPointFn *f1)
2031: {
2032:   PetscFunctionBegin;
2034:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2035:   PetscCall(PetscWeakFormSetIndexBdResidual(ds->wf, NULL, 0, f, 0, 0, f0, 0, f1));
2036:   PetscFunctionReturn(PETSC_SUCCESS);
2037: }

2039: /*@
2040:   PetscDSHasBdJacobian - Indicates that boundary Jacobian functions have been set

2042:   Not Collective

2044:   Input Parameter:
2045: . ds - The `PetscDS`

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

2050:   Level: intermediate

2052: .seealso: `PetscDS`, `PetscDSHasJacobian()`, `PetscDSSetBdJacobian()`, `PetscDSGetBdJacobian()`
2053: @*/
2054: PetscErrorCode PetscDSHasBdJacobian(PetscDS ds, PetscBool *hasBdJac)
2055: {
2056:   PetscFunctionBegin;
2058:   PetscAssertPointer(hasBdJac, 2);
2059:   PetscCall(PetscWeakFormHasBdJacobian(ds->wf, hasBdJac));
2060:   PetscFunctionReturn(PETSC_SUCCESS);
2061: }

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

2066:   Not Collective

2068:   Input Parameters:
2069: + ds - The `PetscDS`
2070: . f  - The test field number
2071: - g  - The field number

2073:   Output Parameters:
2074: + g0 - integrand for the test and basis function term, see `PetscBdPointJacFn`
2075: . g1 - integrand for the test function and basis function gradient term, see `PetscBdPointJacFn`
2076: . g2 - integrand for the test function gradient and basis function term, see `PetscBdPointJacFn`
2077: - g3 - integrand for the test function gradient and basis function gradient term, see `PetscBdPointJacFn`

2079:   Level: intermediate

2081:   Note:
2082:   We are using a first order FEM model for the weak form\:

2084:   $$
2085:   \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
2086:   + \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
2087:   $$

2089: .seealso: `PetscDS`, `PetscBdPointJacFn`, `PetscDSSetBdJacobian()`
2090: @*/
2091: PetscErrorCode PetscDSGetBdJacobian(PetscDS ds, PetscInt f, PetscInt g, PetscBdPointJacFn **g0, PetscBdPointJacFn **g1, PetscBdPointJacFn **g2, PetscBdPointJacFn **g3)
2092: {
2093:   PetscBdPointJacFn **tmp0, **tmp1, **tmp2, **tmp3;
2094:   PetscInt            n0, n1, n2, n3;

2096:   PetscFunctionBegin;
2098:   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);
2099:   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);
2100:   PetscCall(PetscWeakFormGetBdJacobian(ds->wf, NULL, 0, f, g, 0, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3));
2101:   *g0 = tmp0 ? tmp0[0] : NULL;
2102:   *g1 = tmp1 ? tmp1[0] : NULL;
2103:   *g2 = tmp2 ? tmp2[0] : NULL;
2104:   *g3 = tmp3 ? tmp3[0] : NULL;
2105:   PetscFunctionReturn(PETSC_SUCCESS);
2106: }

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

2111:   Not Collective

2113:   Input Parameters:
2114: + ds - The PetscDS
2115: . f  - The test field number
2116: . g  - The field number
2117: . g0 - integrand for the test and basis function term, see `PetscBdPointJacFn`
2118: . g1 - integrand for the test function and basis function gradient term, see `PetscBdPointJacFn`
2119: . g2 - integrand for the test function gradient and basis function term, see `PetscBdPointJacFn`
2120: - g3 - integrand for the test function gradient and basis function gradient term, see `PetscBdPointJacFn`

2122:   Level: intermediate

2124:   Note:
2125:   We are using a first order FEM model for the weak form\:

2127:   $$
2128:   \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
2129:   + \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
2130:   $$

2132: .seealso: `PetscDS`, `PetscBdPointJacFn`, `PetscDSGetBdJacobian()`
2133: @*/
2134: PetscErrorCode PetscDSSetBdJacobian(PetscDS ds, PetscInt f, PetscInt g, PetscBdPointJacFn *g0, PetscBdPointJacFn *g1, PetscBdPointJacFn *g2, PetscBdPointJacFn *g3)
2135: {
2136:   PetscFunctionBegin;
2142:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2143:   PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", g);
2144:   PetscCall(PetscWeakFormSetIndexBdJacobian(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3));
2145:   PetscFunctionReturn(PETSC_SUCCESS);
2146: }

2148: /*@
2149:   PetscDSHasBdJacobianPreconditioner - Signals that boundary Jacobian preconditioner functions have been set with `PetscDSSetBdJacobianPreconditioner()`

2151:   Not Collective

2153:   Input Parameter:
2154: . ds - The `PetscDS`

2156:   Output Parameter:
2157: . hasBdJacPre - flag that pointwise function for the boundary Jacobian matrix to construct the preconditioner has been set

2159:   Level: intermediate

2161:   Developer Note:
2162:   The name is confusing since the function computes a matrix used to construct the preconditioner, not a preconditioner.

2164: .seealso: `PetscDS`, `PetscDSHasJacobian()`, `PetscDSSetBdJacobian()`, `PetscDSGetBdJacobian()`
2165: @*/
2166: PetscErrorCode PetscDSHasBdJacobianPreconditioner(PetscDS ds, PetscBool *hasBdJacPre)
2167: {
2168:   PetscFunctionBegin;
2170:   PetscAssertPointer(hasBdJacPre, 2);
2171:   PetscCall(PetscWeakFormHasBdJacobianPreconditioner(ds->wf, hasBdJacPre));
2172:   PetscFunctionReturn(PETSC_SUCCESS);
2173: }

2175: /*@C
2176:   PetscDSGetBdJacobianPreconditioner - Get the pointwise boundary Jacobian function for given test and basis field that constructs the
2177:   matrix used to construct the preconditioner

2179:   Not Collective; No Fortran Support

2181:   Input Parameters:
2182: + ds - The `PetscDS`
2183: . f  - The test field number
2184: - g  - The field number

2186:   Output Parameters:
2187: + g0 - integrand for the test and basis function term, see `PetscBdPointJacFn`
2188: . g1 - integrand for the test function and basis function gradient term, see `PetscBdPointJacFn`
2189: . g2 - integrand for the test function gradient and basis function term, see `PetscBdPointJacFn`
2190: - g3 - integrand for the test function gradient and basis function gradient term, see `PetscBdPointJacFn`

2192:   Level: intermediate

2194:   Note:
2195:   We are using a first order FEM model for the weak form\:

2197:   $$
2198:   \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
2199:   + \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
2200:   $$

2202:   Developer Note:
2203:   The name is confusing since the function computes a matrix used to construct the preconditioner, not a preconditioner.

2205: .seealso: `PetscDS`, `PetscBdPointJacFn`, `PetscDSSetBdJacobianPreconditioner()`
2206: @*/
2207: PetscErrorCode PetscDSGetBdJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g, PetscBdPointJacFn **g0, PetscBdPointJacFn **g1, PetscBdPointJacFn **g2, PetscBdPointJacFn **g3)
2208: {
2209:   PetscBdPointJacFn **tmp0, **tmp1, **tmp2, **tmp3;
2210:   PetscInt            n0, n1, n2, n3;

2212:   PetscFunctionBegin;
2214:   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);
2215:   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);
2216:   PetscCall(PetscWeakFormGetBdJacobianPreconditioner(ds->wf, NULL, 0, f, g, 0, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3));
2217:   *g0 = tmp0 ? tmp0[0] : NULL;
2218:   *g1 = tmp1 ? tmp1[0] : NULL;
2219:   *g2 = tmp2 ? tmp2[0] : NULL;
2220:   *g3 = tmp3 ? tmp3[0] : NULL;
2221:   PetscFunctionReturn(PETSC_SUCCESS);
2222: }

2224: /*@C
2225:   PetscDSSetBdJacobianPreconditioner - Set the pointwise boundary Jacobian preconditioner function for given test and basis field that constructs the
2226:   matrix used to construct the preconditioner

2228:   Not Collective; No Fortran Support

2230:   Input Parameters:
2231: + ds - The `PetscDS`
2232: . f  - The test field number
2233: . g  - The field number
2234: . g0 - integrand for the test and basis function term, see `PetscBdPointJacFn`
2235: . g1 - integrand for the test function and basis function gradient term, see `PetscBdPointJacFn`
2236: . g2 - integrand for the test function gradient and basis function term, see `PetscBdPointJacFn`
2237: - g3 - integrand for the test function gradient and basis function gradient term, see `PetscBdPointJacFn`

2239:   Level: intermediate

2241:   Note:
2242:   We are using a first order FEM model for the weak form\:

2244:   $$
2245:   \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
2246:   + \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
2247:   $$

2249:   Developer Note:
2250:   The name is confusing since the function computes a matrix used to construct the preconditioner, not a preconditioner.

2252: .seealso: `PetscDS`, `PetscBdPointJacFn`, `PetscDSGetBdJacobianPreconditioner()`
2253: @*/
2254: PetscErrorCode PetscDSSetBdJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g, PetscBdPointJacFn *g0, PetscBdPointJacFn *g1, PetscBdPointJacFn *g2, PetscBdPointJacFn *g3)
2255: {
2256:   PetscFunctionBegin;
2262:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2263:   PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", g);
2264:   PetscCall(PetscWeakFormSetIndexBdJacobianPreconditioner(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3));
2265:   PetscFunctionReturn(PETSC_SUCCESS);
2266: }

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

2271:   Not Collective

2273:   Input Parameters:
2274: + prob - The `PetscDS`
2275: - f    - The test field number

2277:   Output Parameters:
2278: + sol - exact solution function for the test field, see `PetscPointExactSolutionFn`
2279: - ctx - exact solution context

2281:   Level: intermediate

2283: .seealso: `PetscDS`, `PetscPointExactSolutionFn`, `PetscDSSetExactSolution()`, `PetscDSGetExactSolutionTimeDerivative()`
2284: @*/
2285: PetscErrorCode PetscDSGetExactSolution(PetscDS prob, PetscInt f, PetscPointExactSolutionFn **sol, void **ctx)
2286: {
2287:   PetscFunctionBegin;
2289:   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);
2290:   if (sol) {
2291:     PetscAssertPointer(sol, 3);
2292:     *sol = prob->exactSol[f];
2293:   }
2294:   if (ctx) {
2295:     PetscAssertPointer(ctx, 4);
2296:     *ctx = prob->exactCtx[f];
2297:   }
2298:   PetscFunctionReturn(PETSC_SUCCESS);
2299: }

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

2304:   Not Collective

2306:   Input Parameters:
2307: + prob - The `PetscDS`
2308: . f    - The test field number
2309: . sol  - solution function for the test fields, see `PetscPointExactSolutionFn`
2310: - ctx  - solution context or `NULL`

2312:   Level: intermediate

2314: .seealso: `PetscDS`, `PetscPointExactSolutionFn`, `PetscDSGetExactSolution()`
2315: @*/
2316: PetscErrorCode PetscDSSetExactSolution(PetscDS prob, PetscInt f, PetscPointExactSolutionFn *sol, void *ctx)
2317: {
2318:   PetscFunctionBegin;
2320:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2321:   PetscCall(PetscDSEnlarge_Static(prob, f + 1));
2322:   if (sol) {
2324:     prob->exactSol[f] = sol;
2325:   }
2326:   if (ctx) {
2328:     prob->exactCtx[f] = ctx;
2329:   }
2330:   PetscFunctionReturn(PETSC_SUCCESS);
2331: }

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

2336:   Not Collective

2338:   Input Parameters:
2339: + prob - The `PetscDS`
2340: - f    - The test field number

2342:   Output Parameters:
2343: + sol - time derivative of the exact solution for the test field, see `PetscPointExactSolutionFn`
2344: - ctx - the exact solution context

2346:   Level: intermediate

2348: .seealso: `PetscDS`, `PetscPointExactSolutionFn`, `PetscDSSetExactSolutionTimeDerivative()`, `PetscDSGetExactSolution()`
2349: @*/
2350: PetscErrorCode PetscDSGetExactSolutionTimeDerivative(PetscDS prob, PetscInt f, PetscPointExactSolutionFn **sol, void **ctx)
2351: {
2352:   PetscFunctionBegin;
2354:   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);
2355:   if (sol) {
2356:     PetscAssertPointer(sol, 3);
2357:     *sol = prob->exactSol_t[f];
2358:   }
2359:   if (ctx) {
2360:     PetscAssertPointer(ctx, 4);
2361:     *ctx = prob->exactCtx_t[f];
2362:   }
2363:   PetscFunctionReturn(PETSC_SUCCESS);
2364: }

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

2369:   Not Collective

2371:   Input Parameters:
2372: + prob - The `PetscDS`
2373: . f    - The test field number
2374: . sol  - time derivative of the solution function for the test fields, see `PetscPointExactSolutionFn`
2375: - ctx  - the solution context or `NULL`

2377:   Level: intermediate

2379: .seealso: `PetscDS`, `PetscPointExactSolutionFn`, `PetscDSGetExactSolutionTimeDerivative()`, `PetscDSSetExactSolution()`
2380: @*/
2381: PetscErrorCode PetscDSSetExactSolutionTimeDerivative(PetscDS prob, PetscInt f, PetscPointExactSolutionFn *sol, void *ctx)
2382: {
2383:   PetscFunctionBegin;
2385:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2386:   PetscCall(PetscDSEnlarge_Static(prob, f + 1));
2387:   if (sol) {
2389:     prob->exactSol_t[f] = sol;
2390:   }
2391:   if (ctx) {
2393:     prob->exactCtx_t[f] = ctx;
2394:   }
2395:   PetscFunctionReturn(PETSC_SUCCESS);
2396: }

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

2401:   Not Collective

2403:   Input Parameters:
2404: + ds - The PetscDS
2405: - f  - The field number

2407:   Output Parameters:
2408: + lb  - lower bound function for the field, see `PetscPointBoundFn`
2409: - ctx - lower bound context that was set with `PetscDSSetLowerBound()`

2411:   Level: intermediate

2413: .seealso: `PetscDS`, `PetscPointBoundFn`, `PetscDSSetLowerBound()`, `PetscDSGetUpperBound()`, `PetscDSGetExactSolution()`
2414: @*/
2415: PetscErrorCode PetscDSGetLowerBound(PetscDS ds, PetscInt f, PetscPointBoundFn **lb, void **ctx)
2416: {
2417:   PetscFunctionBegin;
2419:   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);
2420:   if (lb) {
2421:     PetscAssertPointer(lb, 3);
2422:     *lb = ds->lowerBound[f];
2423:   }
2424:   if (ctx) {
2425:     PetscAssertPointer(ctx, 4);
2426:     *ctx = ds->lowerCtx[f];
2427:   }
2428:   PetscFunctionReturn(PETSC_SUCCESS);
2429: }

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

2434:   Not Collective

2436:   Input Parameters:
2437: + ds  - The `PetscDS`
2438: . f   - The field number
2439: . lb  - lower bound function for the test fields, see `PetscPointBoundFn`
2440: - ctx - lower bound context or `NULL` which will be passed to `lb`

2442:   Level: intermediate

2444: .seealso: `PetscDS`, `PetscPointBoundFn`, `PetscDSGetLowerBound()`, `PetscDSGetUpperBound()`, `PetscDSGetExactSolution()`
2445: @*/
2446: PetscErrorCode PetscDSSetLowerBound(PetscDS ds, PetscInt f, PetscPointBoundFn *lb, void *ctx)
2447: {
2448:   PetscFunctionBegin;
2450:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2451:   PetscCall(PetscDSEnlarge_Static(ds, f + 1));
2452:   if (lb) {
2454:     ds->lowerBound[f] = lb;
2455:   }
2456:   if (ctx) {
2458:     ds->lowerCtx[f] = ctx;
2459:   }
2460:   PetscFunctionReturn(PETSC_SUCCESS);
2461: }

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

2466:   Not Collective

2468:   Input Parameters:
2469: + ds - The `PetscDS`
2470: - f  - The field number

2472:   Output Parameters:
2473: + ub  - upper bound function for the field, see `PetscPointBoundFn`
2474: - ctx - upper bound context that was set with `PetscDSSetUpperBound()`

2476:   Level: intermediate

2478: .seealso: `PetscDS`, `PetscPointBoundFn`, `PetscDSSetUpperBound()`, `PetscDSGetLowerBound()`, `PetscDSGetExactSolution()`
2479: @*/
2480: PetscErrorCode PetscDSGetUpperBound(PetscDS ds, PetscInt f, PetscPointBoundFn **ub, void **ctx)
2481: {
2482:   PetscFunctionBegin;
2484:   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);
2485:   if (ub) {
2486:     PetscAssertPointer(ub, 3);
2487:     *ub = ds->upperBound[f];
2488:   }
2489:   if (ctx) {
2490:     PetscAssertPointer(ctx, 4);
2491:     *ctx = ds->upperCtx[f];
2492:   }
2493:   PetscFunctionReturn(PETSC_SUCCESS);
2494: }

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

2499:   Not Collective

2501:   Input Parameters:
2502: + ds  - The `PetscDS`
2503: . f   - The field number
2504: . ub  - upper bound function for the test fields, see `PetscPointBoundFn`
2505: - ctx - context or `NULL` that will be passed to `ub`

2507:   Level: intermediate

2509: .seealso: `PetscDS`, `PetscPointBoundFn`, `PetscDSGetUpperBound()`, `PetscDSGetLowerBound()`, `PetscDSGetExactSolution()`
2510: @*/
2511: PetscErrorCode PetscDSSetUpperBound(PetscDS ds, PetscInt f, PetscPointBoundFn *ub, void *ctx)
2512: {
2513:   PetscFunctionBegin;
2515:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2516:   PetscCall(PetscDSEnlarge_Static(ds, f + 1));
2517:   if (ub) {
2519:     ds->upperBound[f] = ub;
2520:   }
2521:   if (ctx) {
2523:     ds->upperCtx[f] = ctx;
2524:   }
2525:   PetscFunctionReturn(PETSC_SUCCESS);
2526: }

2528: /*@C
2529:   PetscDSGetConstants - Returns the array of constants passed to point functions from a `PetscDS` object

2531:   Not Collective

2533:   Input Parameter:
2534: . ds - The `PetscDS` object

2536:   Output Parameters:
2537: + numConstants - The number of constants, or pass in `NULL` if not required
2538: - constants    - The array of constants, `NULL` if there are none

2540:   Level: intermediate

2542: .seealso: `PetscDS`, `PetscDSSetConstants()`, `PetscDSCreate()`
2543: @*/
2544: PetscErrorCode PetscDSGetConstants(PetscDS ds, PeOp PetscInt *numConstants, PeOp const PetscScalar *constants[])
2545: {
2546:   PetscFunctionBegin;
2548:   if (numConstants) {
2549:     PetscAssertPointer(numConstants, 2);
2550:     *numConstants = ds->numConstants;
2551:   }
2552:   if (constants) {
2553:     PetscAssertPointer(constants, 3);
2554:     *constants = ds->constants;
2555:   }
2556:   PetscFunctionReturn(PETSC_SUCCESS);
2557: }

2559: /*@C
2560:   PetscDSSetConstants - Set the array of constants passed to point functions from a `PetscDS`

2562:   Not Collective

2564:   Input Parameters:
2565: + ds           - The `PetscDS` object
2566: . numConstants - The number of constants
2567: - constants    - The array of constants, `NULL` if there are none

2569:   Level: intermediate

2571: .seealso: `PetscDS`, `PetscDSGetConstants()`, `PetscDSCreate()`
2572: @*/
2573: PetscErrorCode PetscDSSetConstants(PetscDS ds, PetscInt numConstants, PetscScalar constants[])
2574: {
2575:   PetscFunctionBegin;
2577:   if (numConstants != ds->numConstants) {
2578:     PetscCall(PetscFree(ds->constants));
2579:     ds->numConstants = numConstants;
2580:     PetscCall(PetscMalloc1(ds->numConstants + ds->numFuncConstants, &ds->constants));
2581:   }
2582:   if (ds->numConstants) {
2583:     PetscAssertPointer(constants, 3);
2584:     PetscCall(PetscArraycpy(ds->constants, constants, ds->numConstants));
2585:   }
2586:   PetscFunctionReturn(PETSC_SUCCESS);
2587: }

2589: /*@C
2590:   PetscDSSetIntegrationParameters - Set the parameters for a particular integration

2592:   Not Collective

2594:   Input Parameters:
2595: + ds     - The `PetscDS` object
2596: . fieldI - The test field for a given point function, or `PETSC_DETERMINE`
2597: - fieldJ - The basis field for a given point function, or `PETSC_DETERMINE`

2599:   Level: intermediate

2601: .seealso: `PetscDS`, `PetscDSSetConstants()`, `PetscDSGetConstants()`, `PetscDSCreate()`
2602: @*/
2603: PetscErrorCode PetscDSSetIntegrationParameters(PetscDS ds, PetscInt fieldI, PetscInt fieldJ)
2604: {
2605:   PetscFunctionBegin;
2607:   ds->constants[ds->numConstants]     = fieldI;
2608:   ds->constants[ds->numConstants + 1] = fieldJ;
2609:   PetscFunctionReturn(PETSC_SUCCESS);
2610: }

2612: /*@C
2613:   PetscDSSetCellParameters - Set the parameters for a particular cell

2615:   Not Collective

2617:   Input Parameters:
2618: + ds     - The `PetscDS` object
2619: - volume - The cell volume

2621:   Level: intermediate

2623: .seealso: `PetscDS`, `PetscDSSetConstants()`, `PetscDSGetConstants()`, `PetscDSCreate()`
2624: @*/
2625: PetscErrorCode PetscDSSetCellParameters(PetscDS ds, PetscReal volume)
2626: {
2627:   PetscFunctionBegin;
2629:   ds->constants[ds->numConstants + 2] = volume;
2630:   PetscFunctionReturn(PETSC_SUCCESS);
2631: }

2633: /*@
2634:   PetscDSGetFieldIndex - Returns the index of the given field

2636:   Not Collective

2638:   Input Parameters:
2639: + prob - The `PetscDS` object
2640: - disc - The discretization object

2642:   Output Parameter:
2643: . f - The field number

2645:   Level: beginner

2647: .seealso: `PetscDS`, `PetscGetDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
2648: @*/
2649: PetscErrorCode PetscDSGetFieldIndex(PetscDS prob, PetscObject disc, PetscInt *f)
2650: {
2651:   PetscInt g;

2653:   PetscFunctionBegin;
2655:   PetscAssertPointer(f, 3);
2656:   *f = -1;
2657:   for (g = 0; g < prob->Nf; ++g) {
2658:     if (disc == prob->disc[g]) break;
2659:   }
2660:   PetscCheck(g != prob->Nf, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONG, "Field not found in PetscDS.");
2661:   *f = g;
2662:   PetscFunctionReturn(PETSC_SUCCESS);
2663: }

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

2668:   Not Collective

2670:   Input Parameters:
2671: + prob - The `PetscDS` object
2672: - f    - The field number

2674:   Output Parameter:
2675: . size - The size

2677:   Level: beginner

2679: .seealso: `PetscDS`, `PetscDSGetFieldOffset()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
2680: @*/
2681: PetscErrorCode PetscDSGetFieldSize(PetscDS prob, PetscInt f, PetscInt *size)
2682: {
2683:   PetscFunctionBegin;
2685:   PetscAssertPointer(size, 3);
2686:   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);
2687:   PetscCall(PetscDSSetUp(prob));
2688:   *size = prob->Nb[f];
2689:   PetscFunctionReturn(PETSC_SUCCESS);
2690: }

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

2695:   Not Collective

2697:   Input Parameters:
2698: + prob - The `PetscDS` object
2699: - f    - The field number

2701:   Output Parameter:
2702: . off - The offset

2704:   Level: beginner

2706: .seealso: `PetscDS`, `PetscDSGetFieldSize()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
2707: @*/
2708: PetscErrorCode PetscDSGetFieldOffset(PetscDS prob, PetscInt f, PetscInt *off)
2709: {
2710:   PetscInt size, g;

2712:   PetscFunctionBegin;
2714:   PetscAssertPointer(off, 3);
2715:   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);
2716:   *off = 0;
2717:   for (g = 0; g < f; ++g) {
2718:     PetscCall(PetscDSGetFieldSize(prob, g, &size));
2719:     *off += size;
2720:   }
2721:   PetscFunctionReturn(PETSC_SUCCESS);
2722: }

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

2727:   Not Collective

2729:   Input Parameters:
2730: + ds - The `PetscDS` object
2731: - f  - The field number

2733:   Output Parameter:
2734: . off - The offset

2736:   Level: beginner

2738: .seealso: `PetscDS`, `PetscDSGetFieldSize()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
2739: @*/
2740: PetscErrorCode PetscDSGetFieldOffsetCohesive(PetscDS ds, PetscInt f, PetscInt *off)
2741: {
2742:   PetscInt size, g;

2744:   PetscFunctionBegin;
2746:   PetscAssertPointer(off, 3);
2747:   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);
2748:   *off = 0;
2749:   for (g = 0; g < f; ++g) {
2750:     PetscBool cohesive;

2752:     PetscCall(PetscDSGetCohesive(ds, g, &cohesive));
2753:     PetscCall(PetscDSGetFieldSize(ds, g, &size));
2754:     *off += cohesive ? size : size * 2;
2755:   }
2756:   PetscFunctionReturn(PETSC_SUCCESS);
2757: }

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

2762:   Not Collective

2764:   Input Parameter:
2765: . prob - The `PetscDS` object

2767:   Output Parameter:
2768: . dimensions - The number of dimensions

2770:   Level: beginner

2772: .seealso: `PetscDS`, `PetscDSGetComponentOffsets()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
2773: @*/
2774: PetscErrorCode PetscDSGetDimensions(PetscDS prob, PetscInt *dimensions[])
2775: {
2776:   PetscFunctionBegin;
2778:   PetscCall(PetscDSSetUp(prob));
2779:   PetscAssertPointer(dimensions, 2);
2780:   *dimensions = prob->Nb;
2781:   PetscFunctionReturn(PETSC_SUCCESS);
2782: }

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

2787:   Not Collective

2789:   Input Parameter:
2790: . prob - The `PetscDS` object

2792:   Output Parameter:
2793: . components - The number of components

2795:   Level: beginner

2797: .seealso: `PetscDS`, `PetscDSGetComponentOffsets()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
2798: @*/
2799: PetscErrorCode PetscDSGetComponents(PetscDS prob, PetscInt *components[])
2800: {
2801:   PetscFunctionBegin;
2803:   PetscCall(PetscDSSetUp(prob));
2804:   PetscAssertPointer(components, 2);
2805:   *components = prob->Nc;
2806:   PetscFunctionReturn(PETSC_SUCCESS);
2807: }

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

2812:   Not Collective

2814:   Input Parameters:
2815: + prob - The `PetscDS` object
2816: - f    - The field number

2818:   Output Parameter:
2819: . off - The offset

2821:   Level: beginner

2823: .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
2824: @*/
2825: PetscErrorCode PetscDSGetComponentOffset(PetscDS prob, PetscInt f, PetscInt *off)
2826: {
2827:   PetscFunctionBegin;
2829:   PetscAssertPointer(off, 3);
2830:   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);
2831:   PetscCall(PetscDSSetUp(prob));
2832:   *off = prob->off[f];
2833:   PetscFunctionReturn(PETSC_SUCCESS);
2834: }

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

2839:   Not Collective

2841:   Input Parameter:
2842: . prob - The `PetscDS` object

2844:   Output Parameter:
2845: . offsets - The offsets

2847:   Level: beginner

2849: .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
2850: @*/
2851: PetscErrorCode PetscDSGetComponentOffsets(PetscDS prob, PetscInt *offsets[])
2852: {
2853:   PetscFunctionBegin;
2855:   PetscAssertPointer(offsets, 2);
2856:   PetscCall(PetscDSSetUp(prob));
2857:   *offsets = prob->off;
2858:   PetscFunctionReturn(PETSC_SUCCESS);
2859: }

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

2864:   Not Collective

2866:   Input Parameter:
2867: . prob - The `PetscDS` object

2869:   Output Parameter:
2870: . offsets - The offsets

2872:   Level: beginner

2874: .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
2875: @*/
2876: PetscErrorCode PetscDSGetComponentDerivativeOffsets(PetscDS prob, PetscInt *offsets[])
2877: {
2878:   PetscFunctionBegin;
2880:   PetscAssertPointer(offsets, 2);
2881:   PetscCall(PetscDSSetUp(prob));
2882:   *offsets = prob->offDer;
2883:   PetscFunctionReturn(PETSC_SUCCESS);
2884: }

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

2889:   Not Collective

2891:   Input Parameters:
2892: + ds - The `PetscDS` object
2893: - s  - The cohesive side, 0 for negative, 1 for positive, 2 for cohesive

2895:   Output Parameter:
2896: . offsets - The offsets

2898:   Level: beginner

2900: .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
2901: @*/
2902: PetscErrorCode PetscDSGetComponentOffsetsCohesive(PetscDS ds, PetscInt s, PetscInt *offsets[])
2903: {
2904:   PetscFunctionBegin;
2906:   PetscAssertPointer(offsets, 3);
2907:   PetscCheck(ds->isCohesive, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cohesive offsets are only valid for a cohesive DS");
2908:   PetscCheck(!(s < 0) && !(s > 2), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cohesive side %" PetscInt_FMT " is not in [0, 2]", s);
2909:   PetscCall(PetscDSSetUp(ds));
2910:   *offsets = ds->offCohesive[s];
2911:   PetscFunctionReturn(PETSC_SUCCESS);
2912: }

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

2917:   Not Collective

2919:   Input Parameters:
2920: + ds - The `PetscDS` object
2921: - s  - The cohesive side, 0 for negative, 1 for positive, 2 for cohesive

2923:   Output Parameter:
2924: . offsets - The offsets

2926:   Level: beginner

2928: .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
2929: @*/
2930: PetscErrorCode PetscDSGetComponentDerivativeOffsetsCohesive(PetscDS ds, PetscInt s, PetscInt *offsets[])
2931: {
2932:   PetscFunctionBegin;
2934:   PetscAssertPointer(offsets, 3);
2935:   PetscCheck(ds->isCohesive, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cohesive offsets are only valid for a cohesive DS");
2936:   PetscCheck(!(s < 0) && !(s > 2), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cohesive side %" PetscInt_FMT " is not in [0, 2]", s);
2937:   PetscCall(PetscDSSetUp(ds));
2938:   *offsets = ds->offDerCohesive[s];
2939:   PetscFunctionReturn(PETSC_SUCCESS);
2940: }

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

2945:   Not Collective

2947:   Input Parameter:
2948: . prob - The `PetscDS` object

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

2953:   Level: intermediate

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

2958:   Fortran Note:
2959:   Use the declaration
2960: .vb
2961:   PetscTabulation, pointer :: tab(:)
2962: .ve
2963:   and access the values using, for example,
2964: .vb
2965:   tab(i)%ptr%K
2966:   tab(i)%ptr%T(j)%ptr
2967: .ve
2968:   where $ i = 1, 2, ..., Nf $ and $ j = 1, 2, ..., tab(i)%ptr%K+1 $.

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

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

2975: .seealso: `PetscDS`, `PetscTabulation`, `PetscDSCreate()`
2976: @*/
2977: PetscErrorCode PetscDSGetTabulation(PetscDS prob, PetscTabulation *T[]) PeNS
2978: {
2979:   PetscFunctionBegin;
2981:   PetscAssertPointer(T, 2);
2982:   PetscCall(PetscDSSetUp(prob));
2983:   *T = prob->T;
2984:   PetscFunctionReturn(PETSC_SUCCESS);
2985: }

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

2990:   Not Collective

2992:   Input Parameter:
2993: . prob - The `PetscDS` object

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

2998:   Level: intermediate

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

3003: .seealso: `PetscTabulation`, `PetscDS`, `PetscDSGetTabulation()`, `PetscDSCreate()`
3004: @*/
3005: PetscErrorCode PetscDSGetFaceTabulation(PetscDS prob, PetscTabulation *Tf[])
3006: {
3007:   PetscFunctionBegin;
3009:   PetscAssertPointer(Tf, 2);
3010:   PetscCall(PetscDSSetUp(prob));
3011:   *Tf = prob->Tf;
3012:   PetscFunctionReturn(PETSC_SUCCESS);
3013: }

3015: PetscErrorCode PetscDSGetEvaluationArrays(PetscDS prob, PetscScalar *u[], PetscScalar *u_t[], PetscScalar *u_x[])
3016: {
3017:   PetscFunctionBegin;
3019:   PetscCall(PetscDSSetUp(prob));
3020:   if (u) {
3021:     PetscAssertPointer(u, 2);
3022:     *u = prob->u;
3023:   }
3024:   if (u_t) {
3025:     PetscAssertPointer(u_t, 3);
3026:     *u_t = prob->u_t;
3027:   }
3028:   if (u_x) {
3029:     PetscAssertPointer(u_x, 4);
3030:     *u_x = prob->u_x;
3031:   }
3032:   PetscFunctionReturn(PETSC_SUCCESS);
3033: }

3035: PetscErrorCode PetscDSGetWeakFormArrays(PetscDS prob, PetscScalar *f0[], PetscScalar *f1[], PetscScalar *g0[], PetscScalar *g1[], PetscScalar *g2[], PetscScalar *g3[])
3036: {
3037:   PetscFunctionBegin;
3039:   PetscCall(PetscDSSetUp(prob));
3040:   if (f0) {
3041:     PetscAssertPointer(f0, 2);
3042:     *f0 = prob->f0;
3043:   }
3044:   if (f1) {
3045:     PetscAssertPointer(f1, 3);
3046:     *f1 = prob->f1;
3047:   }
3048:   if (g0) {
3049:     PetscAssertPointer(g0, 4);
3050:     *g0 = prob->g0;
3051:   }
3052:   if (g1) {
3053:     PetscAssertPointer(g1, 5);
3054:     *g1 = prob->g1;
3055:   }
3056:   if (g2) {
3057:     PetscAssertPointer(g2, 6);
3058:     *g2 = prob->g2;
3059:   }
3060:   if (g3) {
3061:     PetscAssertPointer(g3, 7);
3062:     *g3 = prob->g3;
3063:   }
3064:   PetscFunctionReturn(PETSC_SUCCESS);
3065: }

3067: PetscErrorCode PetscDSGetWorkspace(PetscDS prob, PetscReal **x, PetscScalar **basisReal, PetscScalar **basisDerReal, PetscScalar **testReal, PetscScalar **testDerReal)
3068: {
3069:   PetscFunctionBegin;
3071:   PetscCall(PetscDSSetUp(prob));
3072:   if (x) {
3073:     PetscAssertPointer(x, 2);
3074:     *x = prob->x;
3075:   }
3076:   if (basisReal) {
3077:     PetscAssertPointer(basisReal, 3);
3078:     *basisReal = prob->basisReal;
3079:   }
3080:   if (basisDerReal) {
3081:     PetscAssertPointer(basisDerReal, 4);
3082:     *basisDerReal = prob->basisDerReal;
3083:   }
3084:   if (testReal) {
3085:     PetscAssertPointer(testReal, 5);
3086:     *testReal = prob->testReal;
3087:   }
3088:   if (testDerReal) {
3089:     PetscAssertPointer(testDerReal, 6);
3090:     *testDerReal = prob->testDerReal;
3091:   }
3092:   PetscFunctionReturn(PETSC_SUCCESS);
3093: }

3095: /*@C
3096:   PetscDSAddBoundary - Add a boundary condition to the model.

3098:   Collective

3100:   Input Parameters:
3101: + ds       - The `PetscDS` object
3102: . type     - The type of condition, e.g. `DM_BC_ESSENTIAL`/`DM_BC_ESSENTIAL_FIELD` (Dirichlet), or `DM_BC_NATURAL` (Neumann)
3103: . name     - The name for the boundary condition
3104: . label    - The label defining constrained points
3105: . Nv       - The number of `DMLabel` values for constrained points
3106: . values   - An array of label values for constrained points
3107: . field    - The field to constrain
3108: . Nc       - The number of constrained field components (0 will constrain all fields)
3109: . comps    - An array of constrained component numbers
3110: . bcFunc   - A pointwise function giving boundary values
3111: . bcFunc_t - A pointwise function giving the time derivative of the boundary values, or `NULL`
3112: - ctx      - An optional user context for `bcFunc`

3114:   Output Parameter:
3115: . bd - The boundary number

3117:   Options Database Keys:
3118: + -bc_<boundary name> <num>      - Overrides the boundary ids
3119: - -bc_<boundary name>_comp <num> - Overrides the boundary components

3121:   Level: developer

3123:   Note:
3124:   Both `bcFunc` and `bcFunc_t` will depend on the boundary condition type. If the type if `DM_BC_ESSENTIAL`, then the calling sequence is\:
3125: .vb
3126:   void bcFunc(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar bcval[])
3127: .ve

3129:   If the type is `DM_BC_ESSENTIAL_FIELD` or other _FIELD value, then the calling sequence is\:
3130: .vb
3131:   void bcFunc(PetscInt dim, PetscInt Nf, PetscInt NfAux,
3132:               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
3133:               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
3134:               PetscReal time, const PetscReal x[], PetscScalar bcval[])
3135: .ve
3136: + dim          - the coordinate dimension
3137: . Nf           - the number of fields
3138: . uOff         - the offset into u[] and u_t[] for each field
3139: . uOff_x       - the offset into u_x[] for each field
3140: . u            - each field evaluated at the current point
3141: . u_t          - the time derivative of each field evaluated at the current point
3142: . u_x          - the gradient of each field evaluated at the current point
3143: . aOff         - the offset into a[] and a_t[] for each auxiliary field
3144: . aOff_x       - the offset into a_x[] for each auxiliary field
3145: . a            - each auxiliary field evaluated at the current point
3146: . a_t          - the time derivative of each auxiliary field evaluated at the current point
3147: . a_x          - the gradient of auxiliary each field evaluated at the current point
3148: . t            - current time
3149: . x            - coordinates of the current point
3150: . numConstants - number of constant parameters
3151: . constants    - constant parameters
3152: - bcval        - output values at the current point

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

3160: .seealso: `PetscDS`, `PetscWeakForm`, `DMLabel`, `DMBoundaryConditionType`, `PetscDSAddBoundaryByName()`, `PetscDSGetBoundary()`, `PetscDSSetResidual()`, `PetscDSSetBdResidual()`
3161: @*/
3162: 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)
3163: {
3164:   DSBoundary  head = ds->boundary, b;
3165:   PetscInt    n    = 0;
3166:   const char *lname;

3168:   PetscFunctionBegin;
3171:   PetscAssertPointer(name, 3);
3176:   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);
3177:   if (Nc > 0) {
3178:     PetscInt *fcomps;
3179:     PetscInt  c;

3181:     PetscCall(PetscDSGetComponents(ds, &fcomps));
3182:     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);
3183:     for (c = 0; c < Nc; ++c) {
3184:       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);
3185:     }
3186:   }
3187:   PetscCall(PetscNew(&b));
3188:   PetscCall(PetscStrallocpy(name, (char **)&b->name));
3189:   PetscCall(PetscWeakFormCreate(PETSC_COMM_SELF, &b->wf));
3190:   PetscCall(PetscWeakFormSetNumFields(b->wf, ds->Nf));
3191:   PetscCall(PetscMalloc1(Nv, &b->values));
3192:   if (Nv) PetscCall(PetscArraycpy(b->values, values, Nv));
3193:   PetscCall(PetscMalloc1(Nc, &b->comps));
3194:   if (Nc) PetscCall(PetscArraycpy(b->comps, comps, Nc));
3195:   PetscCall(PetscObjectGetName((PetscObject)label, &lname));
3196:   PetscCall(PetscStrallocpy(lname, (char **)&b->lname));
3197:   b->type   = type;
3198:   b->label  = label;
3199:   b->Nv     = Nv;
3200:   b->field  = field;
3201:   b->Nc     = Nc;
3202:   b->func   = bcFunc;
3203:   b->func_t = bcFunc_t;
3204:   b->ctx    = ctx;
3205:   b->next   = NULL;
3206:   /* Append to linked list so that we can preserve the order */
3207:   if (!head) ds->boundary = b;
3208:   while (head) {
3209:     if (!head->next) {
3210:       head->next = b;
3211:       head       = b;
3212:     }
3213:     head = head->next;
3214:     ++n;
3215:   }
3216:   if (bd) {
3217:     PetscAssertPointer(bd, 13);
3218:     *bd = n;
3219:   }
3220:   PetscFunctionReturn(PETSC_SUCCESS);
3221: }

3223: // PetscClangLinter pragma ignore: -fdoc-section-header-unknown
3224: /*@C
3225:   PetscDSAddBoundaryByName - Add a boundary condition to the model.

3227:   Collective

3229:   Input Parameters:
3230: + ds       - The `PetscDS` object
3231: . type     - The type of condition, e.g. `DM_BC_ESSENTIAL`/`DM_BC_ESSENTIAL_FIELD` (Dirichlet), or `DM_BC_NATURAL` (Neumann)
3232: . name     - The boundary condition name
3233: . lname    - The name of the label defining constrained points
3234: . Nv       - The number of `DMLabel` values for constrained points
3235: . values   - An array of label values for constrained points
3236: . field    - The field to constrain
3237: . Nc       - The number of constrained field components (0 will constrain all fields)
3238: . comps    - An array of constrained component numbers
3239: . bcFunc   - A pointwise function giving boundary values
3240: . bcFunc_t - A pointwise function giving the time derivative of the boundary values, or `NULL`
3241: - ctx      - An optional user context for `bcFunc`

3243:   Output Parameter:
3244: . bd - The boundary number

3246:   Options Database Keys:
3247: + -bc_<boundary name> <num>      - Overrides the boundary ids
3248: - -bc_<boundary name>_comp <num> - Overrides the boundary components

3250:   Calling Sequence of `bcFunc` and `bcFunc_t`:
3251:   If the type is `DM_BC_ESSENTIAL`
3252: .vb
3253:   void bcFunc(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar bcval[])
3254: .ve
3255:   If the type is `DM_BC_ESSENTIAL_FIELD` or other _FIELD value,
3256: .vb
3257:   void bcFunc(PetscInt dim, PetscInt Nf, PetscInt NfAux,
3258:               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
3259:               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
3260:               PetscReal time, const PetscReal x[], PetscScalar bcval[])
3261: .ve
3262: + dim          - the coordinate dimension
3263: . Nf           - the number of fields
3264: . uOff         - the offset into `u`[] and `u_t`[] for each field
3265: . uOff_x       - the offset into `u_x`[] for each field
3266: . u            - each field evaluated at the current point
3267: . u_t          - the time derivative of each field evaluated at the current point
3268: . u_x          - the gradient of each field evaluated at the current point
3269: . aOff         - the offset into `a`[] and `a_t`[] for each auxiliary field
3270: . aOff_x       - the offset into `a_x`[] for each auxiliary field
3271: . a            - each auxiliary field evaluated at the current point
3272: . a_t          - the time derivative of each auxiliary field evaluated at the current point
3273: . a_x          - the gradient of auxiliary each field evaluated at the current point
3274: . t            - current time
3275: . x            - coordinates of the current point
3276: . numConstants - number of constant parameters
3277: . constants    - constant parameters
3278: - bcval        - output values at the current point

3280:   Level: developer

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

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

3290: .seealso: `PetscDS`, `PetscWeakForm`, `DMLabel`, `DMBoundaryConditionType`, `PetscDSAddBoundary()`, `PetscDSGetBoundary()`, `PetscDSSetResidual()`, `PetscDSSetBdResidual()`
3291: @*/
3292: 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)
3293: {
3294:   DSBoundary head = ds->boundary, b;
3295:   PetscInt   n    = 0;

3297:   PetscFunctionBegin;
3300:   PetscAssertPointer(name, 3);
3301:   PetscAssertPointer(lname, 4);
3305:   PetscCall(PetscNew(&b));
3306:   PetscCall(PetscStrallocpy(name, (char **)&b->name));
3307:   PetscCall(PetscWeakFormCreate(PETSC_COMM_SELF, &b->wf));
3308:   PetscCall(PetscWeakFormSetNumFields(b->wf, ds->Nf));
3309:   PetscCall(PetscMalloc1(Nv, &b->values));
3310:   if (Nv) PetscCall(PetscArraycpy(b->values, values, Nv));
3311:   PetscCall(PetscMalloc1(Nc, &b->comps));
3312:   if (Nc) PetscCall(PetscArraycpy(b->comps, comps, Nc));
3313:   PetscCall(PetscStrallocpy(lname, (char **)&b->lname));
3314:   b->type   = type;
3315:   b->label  = NULL;
3316:   b->Nv     = Nv;
3317:   b->field  = field;
3318:   b->Nc     = Nc;
3319:   b->func   = bcFunc;
3320:   b->func_t = bcFunc_t;
3321:   b->ctx    = ctx;
3322:   b->next   = NULL;
3323:   /* Append to linked list so that we can preserve the order */
3324:   if (!head) ds->boundary = b;
3325:   while (head) {
3326:     if (!head->next) {
3327:       head->next = b;
3328:       head       = b;
3329:     }
3330:     head = head->next;
3331:     ++n;
3332:   }
3333:   if (bd) {
3334:     PetscAssertPointer(bd, 13);
3335:     *bd = n;
3336:   }
3337:   PetscFunctionReturn(PETSC_SUCCESS);
3338: }

3340: /*@C
3341:   PetscDSUpdateBoundary - Change a boundary condition for the model.

3343:   Input Parameters:
3344: + ds       - The `PetscDS` object
3345: . bd       - The boundary condition number
3346: . type     - The type of condition, e.g. `DM_BC_ESSENTIAL`/`DM_BC_ESSENTIAL_FIELD` (Dirichlet), or `DM_BC_NATURAL` (Neumann)
3347: . name     - The boundary condition name
3348: . label    - The label defining constrained points
3349: . Nv       - The number of `DMLabel` ids for constrained points
3350: . values   - An array of ids for constrained points
3351: . field    - The field to constrain
3352: . Nc       - The number of constrained field components
3353: . comps    - An array of constrained component numbers
3354: . bcFunc   - A pointwise function giving boundary values
3355: . bcFunc_t - A pointwise function giving the time derivative of the boundary values, or `NULL`
3356: - ctx      - An optional user context for `bcFunc`

3358:   Level: developer

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

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

3369: .seealso: `PetscDS`, `PetscWeakForm`, `DMBoundaryConditionType`, `PetscDSAddBoundary()`, `PetscDSGetBoundary()`, `PetscDSGetNumBoundary()`, `DMLabel`
3370: @*/
3371: 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)
3372: {
3373:   DSBoundary b = ds->boundary;
3374:   PetscInt   n = 0;

3376:   PetscFunctionBegin;
3378:   while (b) {
3379:     if (n == bd) break;
3380:     b = b->next;
3381:     ++n;
3382:   }
3383:   PetscCheck(b, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %" PetscInt_FMT " is not in [0, %" PetscInt_FMT ")", bd, n);
3384:   if (name) {
3385:     PetscCall(PetscFree(b->name));
3386:     PetscCall(PetscStrallocpy(name, (char **)&b->name));
3387:   }
3388:   b->type = type;
3389:   if (label) {
3390:     const char *name;

3392:     b->label = label;
3393:     PetscCall(PetscFree(b->lname));
3394:     PetscCall(PetscObjectGetName((PetscObject)label, &name));
3395:     PetscCall(PetscStrallocpy(name, (char **)&b->lname));
3396:   }
3397:   if (Nv >= 0) {
3398:     b->Nv = Nv;
3399:     PetscCall(PetscFree(b->values));
3400:     PetscCall(PetscMalloc1(Nv, &b->values));
3401:     if (Nv) PetscCall(PetscArraycpy(b->values, values, Nv));
3402:   }
3403:   if (field >= 0) b->field = field;
3404:   if (Nc >= 0) {
3405:     b->Nc = Nc;
3406:     PetscCall(PetscFree(b->comps));
3407:     PetscCall(PetscMalloc1(Nc, &b->comps));
3408:     if (Nc) PetscCall(PetscArraycpy(b->comps, comps, Nc));
3409:   }
3410:   if (bcFunc) b->func = bcFunc;
3411:   if (bcFunc_t) b->func_t = bcFunc_t;
3412:   if (ctx) b->ctx = ctx;
3413:   PetscFunctionReturn(PETSC_SUCCESS);
3414: }

3416: /*@
3417:   PetscDSGetNumBoundary - Get the number of registered boundary conditions

3419:   Input Parameter:
3420: . ds - The `PetscDS` object

3422:   Output Parameter:
3423: . numBd - The number of boundary conditions

3425:   Level: intermediate

3427: .seealso: `PetscDS`, `PetscDSAddBoundary()`, `PetscDSGetBoundary()`
3428: @*/
3429: PetscErrorCode PetscDSGetNumBoundary(PetscDS ds, PetscInt *numBd)
3430: {
3431:   DSBoundary b = ds->boundary;

3433:   PetscFunctionBegin;
3435:   PetscAssertPointer(numBd, 2);
3436:   *numBd = 0;
3437:   while (b) {
3438:     ++(*numBd);
3439:     b = b->next;
3440:   }
3441:   PetscFunctionReturn(PETSC_SUCCESS);
3442: }

3444: /*@C
3445:   PetscDSGetBoundary - Gets a boundary condition from the model

3447:   Input Parameters:
3448: + ds - The `PetscDS` object
3449: - bd - The boundary condition number

3451:   Output Parameters:
3452: + wf     - The `PetscWeakForm` holding the pointwise functions
3453: . type   - The type of condition, e.g. `DM_BC_ESSENTIAL`/`DM_BC_ESSENTIAL_FIELD` (Dirichlet), or `DM_BC_NATURAL` (Neumann)
3454: . name   - The boundary condition name
3455: . label  - The label defining constrained points
3456: . Nv     - The number of `DMLabel` ids for constrained points
3457: . values - An array of ids for constrained points
3458: . field  - The field to constrain
3459: . Nc     - The number of constrained field components
3460: . comps  - An array of constrained component numbers
3461: . func   - A pointwise function giving boundary values
3462: . func_t - A pointwise function giving the time derivative of the boundary values
3463: - ctx    - An optional user context for `bcFunc`

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

3469:   Level: developer

3471: .seealso: `PetscDS`, `PetscWeakForm`, `DMBoundaryConditionType`, `PetscDSAddBoundary()`, `DMLabel`
3472: @*/
3473: 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)
3474: {
3475:   DSBoundary b = ds->boundary;
3476:   PetscInt   n = 0;

3478:   PetscFunctionBegin;
3480:   while (b) {
3481:     if (n == bd) break;
3482:     b = b->next;
3483:     ++n;
3484:   }
3485:   PetscCheck(b, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %" PetscInt_FMT " is not in [0, %" PetscInt_FMT ")", bd, n);
3486:   if (wf) {
3487:     PetscAssertPointer(wf, 3);
3488:     *wf = b->wf;
3489:   }
3490:   if (type) {
3491:     PetscAssertPointer(type, 4);
3492:     *type = b->type;
3493:   }
3494:   if (name) {
3495:     PetscAssertPointer(name, 5);
3496:     *name = b->name;
3497:   }
3498:   if (label) {
3499:     PetscAssertPointer(label, 6);
3500:     *label = b->label;
3501:   }
3502:   if (Nv) {
3503:     PetscAssertPointer(Nv, 7);
3504:     *Nv = b->Nv;
3505:   }
3506:   if (values) {
3507:     PetscAssertPointer(values, 8);
3508:     *values = b->values;
3509:   }
3510:   if (field) {
3511:     PetscAssertPointer(field, 9);
3512:     *field = b->field;
3513:   }
3514:   if (Nc) {
3515:     PetscAssertPointer(Nc, 10);
3516:     *Nc = b->Nc;
3517:   }
3518:   if (comps) {
3519:     PetscAssertPointer(comps, 11);
3520:     *comps = b->comps;
3521:   }
3522:   if (func) {
3523:     PetscAssertPointer(func, 12);
3524:     *func = b->func;
3525:   }
3526:   if (func_t) {
3527:     PetscAssertPointer(func_t, 13);
3528:     *func_t = b->func_t;
3529:   }
3530:   if (ctx) {
3531:     PetscAssertPointer(ctx, 14);
3532:     *ctx = b->ctx;
3533:   }
3534:   PetscFunctionReturn(PETSC_SUCCESS);
3535: }

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

3540:   Not Collective

3542:   Input Parameters:
3543: + ds - The source `PetscDS` object
3544: - dm - The `DM` holding labels

3546:   Level: intermediate

3548: .seealso: `PetscDS`, `DMBoundary`, `DM`, `PetscDSCopyBoundary()`, `PetscDSCreate()`, `DMGetLabel()`
3549: @*/
3550: PetscErrorCode PetscDSUpdateBoundaryLabels(PetscDS ds, DM dm)
3551: {
3552:   DSBoundary b;

3554:   PetscFunctionBegin;
3557:   for (b = ds->boundary; b; b = b->next) {
3558:     if (b->lname) PetscCall(DMGetLabel(dm, b->lname, &b->label));
3559:   }
3560:   PetscFunctionReturn(PETSC_SUCCESS);
3561: }

3563: static PetscErrorCode DSBoundaryDuplicate_Internal(DSBoundary b, DSBoundary *bNew)
3564: {
3565:   PetscFunctionBegin;
3566:   PetscCall(PetscNew(bNew));
3567:   PetscCall(PetscWeakFormCreate(PETSC_COMM_SELF, &(*bNew)->wf));
3568:   PetscCall(PetscWeakFormCopy(b->wf, (*bNew)->wf));
3569:   PetscCall(PetscStrallocpy(b->name, (char **)&((*bNew)->name)));
3570:   PetscCall(PetscStrallocpy(b->lname, (char **)&((*bNew)->lname)));
3571:   (*bNew)->type  = b->type;
3572:   (*bNew)->label = b->label;
3573:   (*bNew)->Nv    = b->Nv;
3574:   PetscCall(PetscMalloc1(b->Nv, &(*bNew)->values));
3575:   PetscCall(PetscArraycpy((*bNew)->values, b->values, b->Nv));
3576:   (*bNew)->field = b->field;
3577:   (*bNew)->Nc    = b->Nc;
3578:   PetscCall(PetscMalloc1(b->Nc, &(*bNew)->comps));
3579:   PetscCall(PetscArraycpy((*bNew)->comps, b->comps, b->Nc));
3580:   (*bNew)->func   = b->func;
3581:   (*bNew)->func_t = b->func_t;
3582:   (*bNew)->ctx    = b->ctx;
3583:   PetscFunctionReturn(PETSC_SUCCESS);
3584: }

3586: /*@
3587:   PetscDSCopyBoundary - Copy all boundary condition objects to the new `PetscDS`

3589:   Not Collective

3591:   Input Parameters:
3592: + ds        - The source `PetscDS` object
3593: . numFields - The number of selected fields, or `PETSC_DEFAULT` for all fields
3594: - fields    - The selected fields, or `NULL` for all fields

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

3599:   Level: intermediate

3601: .seealso: `PetscDS`, `DMBoundary`, `PetscDSCopyEquations()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
3602: @*/
3603: PetscErrorCode PetscDSCopyBoundary(PetscDS ds, PetscInt numFields, const PetscInt fields[], PetscDS newds)
3604: {
3605:   DSBoundary b, *lastnext;

3607:   PetscFunctionBegin;
3610:   if (ds == newds) PetscFunctionReturn(PETSC_SUCCESS);
3611:   PetscCall(PetscDSDestroyBoundary(newds));
3612:   lastnext = &newds->boundary;
3613:   for (b = ds->boundary; b; b = b->next) {
3614:     DSBoundary bNew;
3615:     PetscInt   fieldNew = -1;

3617:     if (numFields > 0 && fields) {
3618:       PetscInt f;

3620:       for (f = 0; f < numFields; ++f)
3621:         if (b->field == fields[f]) break;
3622:       if (f == numFields) continue;
3623:       fieldNew = f;
3624:     }
3625:     PetscCall(DSBoundaryDuplicate_Internal(b, &bNew));
3626:     bNew->field = fieldNew < 0 ? b->field : fieldNew;
3627:     *lastnext   = bNew;
3628:     lastnext    = &bNew->next;
3629:   }
3630:   PetscFunctionReturn(PETSC_SUCCESS);
3631: }

3633: /*@
3634:   PetscDSDestroyBoundary - Remove all `DMBoundary` objects from the `PetscDS`

3636:   Not Collective

3638:   Input Parameter:
3639: . ds - The `PetscDS` object

3641:   Level: intermediate

3643: .seealso: `PetscDS`, `DMBoundary`, `PetscDSCopyBoundary()`, `PetscDSCopyEquations()`
3644: @*/
3645: PetscErrorCode PetscDSDestroyBoundary(PetscDS ds)
3646: {
3647:   DSBoundary next = ds->boundary;

3649:   PetscFunctionBegin;
3650:   while (next) {
3651:     DSBoundary b = next;

3653:     next = b->next;
3654:     PetscCall(PetscWeakFormDestroy(&b->wf));
3655:     PetscCall(PetscFree(b->name));
3656:     PetscCall(PetscFree(b->lname));
3657:     PetscCall(PetscFree(b->values));
3658:     PetscCall(PetscFree(b->comps));
3659:     PetscCall(PetscFree(b));
3660:   }
3661:   PetscFunctionReturn(PETSC_SUCCESS);
3662: }

3664: /*@
3665:   PetscDSSelectDiscretizations - Copy discretizations to the new `PetscDS` with different field layout

3667:   Not Collective

3669:   Input Parameters:
3670: + prob      - The `PetscDS` object
3671: . numFields - Number of new fields
3672: . fields    - Old field number for each new field
3673: . minDegree - Minimum degree for a discretization, or `PETSC_DETERMINE` for no limit
3674: - maxDegree - Maximum degree for a discretization, or `PETSC_DETERMINE` for no limit

3676:   Output Parameter:
3677: . newprob - The `PetscDS` copy

3679:   Level: intermediate

3681: .seealso: `PetscDS`, `PetscDSSelectEquations()`, `PetscDSCopyBoundary()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
3682: @*/
3683: PetscErrorCode PetscDSSelectDiscretizations(PetscDS prob, PetscInt numFields, const PetscInt fields[], PetscInt minDegree, PetscInt maxDegree, PetscDS newprob)
3684: {
3685:   PetscInt Nf, Nfn, fn;

3687:   PetscFunctionBegin;
3689:   if (fields) PetscAssertPointer(fields, 3);
3691:   PetscCall(PetscDSGetNumFields(prob, &Nf));
3692:   PetscCall(PetscDSGetNumFields(newprob, &Nfn));
3693:   numFields = numFields < 0 ? Nf : numFields;
3694:   for (fn = 0; fn < numFields; ++fn) {
3695:     const PetscInt f = fields ? fields[fn] : fn;
3696:     PetscObject    disc;
3697:     PetscClassId   id;

3699:     if (f >= Nf) continue;
3700:     PetscCall(PetscDSGetDiscretization(prob, f, &disc));
3701:     PetscCallContinue(PetscObjectGetClassId(disc, &id));
3702:     if (id == PETSCFE_CLASSID) {
3703:       PetscFE fe;

3705:       PetscCall(PetscFELimitDegree((PetscFE)disc, minDegree, maxDegree, &fe));
3706:       PetscCall(PetscDSSetDiscretization(newprob, fn, (PetscObject)fe));
3707:       PetscCall(PetscFEDestroy(&fe));
3708:     } else {
3709:       PetscCall(PetscDSSetDiscretization(newprob, fn, disc));
3710:     }
3711:   }
3712:   PetscFunctionReturn(PETSC_SUCCESS);
3713: }

3715: /*@
3716:   PetscDSSelectEquations - Copy pointwise function pointers to the new `PetscDS` with different field layout

3718:   Not Collective

3720:   Input Parameters:
3721: + prob      - The `PetscDS` object
3722: . numFields - Number of new fields
3723: - fields    - Old field number for each new field

3725:   Output Parameter:
3726: . newprob - The `PetscDS` copy

3728:   Level: intermediate

3730: .seealso: `PetscDS`, `PetscDSSelectDiscretizations()`, `PetscDSCopyBoundary()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
3731: @*/
3732: PetscErrorCode PetscDSSelectEquations(PetscDS prob, PetscInt numFields, const PetscInt fields[], PetscDS newprob)
3733: {
3734:   PetscInt Nf, Nfn, fn, gn;

3736:   PetscFunctionBegin;
3738:   if (fields) PetscAssertPointer(fields, 3);
3740:   PetscCall(PetscDSGetNumFields(prob, &Nf));
3741:   PetscCall(PetscDSGetNumFields(newprob, &Nfn));
3742:   PetscCheck(numFields <= Nfn, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_SIZ, "Number of fields %" PetscInt_FMT " to transfer must not be greater than the total number of fields %" PetscInt_FMT, numFields, Nfn);
3743:   for (fn = 0; fn < numFields; ++fn) {
3744:     const PetscInt  f = fields ? fields[fn] : fn;
3745:     PetscPointFn   *obj;
3746:     PetscPointFn   *f0, *f1;
3747:     PetscBdPointFn *f0Bd, *f1Bd;
3748:     PetscRiemannFn *r;

3750:     if (f >= Nf) continue;
3751:     PetscCall(PetscDSGetObjective(prob, f, &obj));
3752:     PetscCall(PetscDSGetResidual(prob, f, &f0, &f1));
3753:     PetscCall(PetscDSGetBdResidual(prob, f, &f0Bd, &f1Bd));
3754:     PetscCall(PetscDSGetRiemannSolver(prob, f, &r));
3755:     PetscCall(PetscDSSetObjective(newprob, fn, obj));
3756:     PetscCall(PetscDSSetResidual(newprob, fn, f0, f1));
3757:     PetscCall(PetscDSSetBdResidual(newprob, fn, f0Bd, f1Bd));
3758:     PetscCall(PetscDSSetRiemannSolver(newprob, fn, r));
3759:     for (gn = 0; gn < numFields; ++gn) {
3760:       const PetscInt     g = fields ? fields[gn] : gn;
3761:       PetscPointJacFn   *g0, *g1, *g2, *g3;
3762:       PetscPointJacFn   *g0p, *g1p, *g2p, *g3p;
3763:       PetscBdPointJacFn *g0Bd, *g1Bd, *g2Bd, *g3Bd;

3765:       if (g >= Nf) continue;
3766:       PetscCall(PetscDSGetJacobian(prob, f, g, &g0, &g1, &g2, &g3));
3767:       PetscCall(PetscDSGetJacobianPreconditioner(prob, f, g, &g0p, &g1p, &g2p, &g3p));
3768:       PetscCall(PetscDSGetBdJacobian(prob, f, g, &g0Bd, &g1Bd, &g2Bd, &g3Bd));
3769:       PetscCall(PetscDSSetJacobian(newprob, fn, gn, g0, g1, g2, g3));
3770:       PetscCall(PetscDSSetJacobianPreconditioner(newprob, fn, gn, g0p, g1p, g2p, g3p));
3771:       PetscCall(PetscDSSetBdJacobian(newprob, fn, gn, g0Bd, g1Bd, g2Bd, g3Bd));
3772:     }
3773:   }
3774:   PetscFunctionReturn(PETSC_SUCCESS);
3775: }

3777: /*@
3778:   PetscDSCopyEquations - Copy all pointwise function pointers to another `PetscDS`

3780:   Not Collective

3782:   Input Parameter:
3783: . prob - The `PetscDS` object

3785:   Output Parameter:
3786: . newprob - The `PetscDS` copy

3788:   Level: intermediate

3790: .seealso: `PetscDS`, `PetscDSCopyBoundary()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
3791: @*/
3792: PetscErrorCode PetscDSCopyEquations(PetscDS prob, PetscDS newprob)
3793: {
3794:   PetscWeakForm wf, newwf;
3795:   PetscInt      Nf, Ng;

3797:   PetscFunctionBegin;
3800:   PetscCall(PetscDSGetNumFields(prob, &Nf));
3801:   PetscCall(PetscDSGetNumFields(newprob, &Ng));
3802:   PetscCheck(Nf == Ng, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_SIZ, "Number of fields must match %" PetscInt_FMT " != %" PetscInt_FMT, Nf, Ng);
3803:   PetscCall(PetscDSGetWeakForm(prob, &wf));
3804:   PetscCall(PetscDSGetWeakForm(newprob, &newwf));
3805:   PetscCall(PetscWeakFormCopy(wf, newwf));
3806:   PetscFunctionReturn(PETSC_SUCCESS);
3807: }

3809: /*@
3810:   PetscDSCopyConstants - Copy all constants set with `PetscDSSetConstants()` to another `PetscDS`

3812:   Not Collective

3814:   Input Parameter:
3815: . prob - The `PetscDS` object

3817:   Output Parameter:
3818: . newprob - The `PetscDS` copy

3820:   Level: intermediate

3822: .seealso: `PetscDS`, `PetscDSCopyBoundary()`, `PetscDSCopyEquations()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
3823: @*/
3824: PetscErrorCode PetscDSCopyConstants(PetscDS prob, PetscDS newprob)
3825: {
3826:   PetscInt           Nc;
3827:   const PetscScalar *constants;

3829:   PetscFunctionBegin;
3832:   PetscCall(PetscDSGetConstants(prob, &Nc, &constants));
3833:   PetscCall(PetscDSSetConstants(newprob, Nc, (PetscScalar *)constants));
3834:   PetscFunctionReturn(PETSC_SUCCESS);
3835: }

3837: /*@
3838:   PetscDSCopyExactSolutions - Copy all exact solutions set with `PetscDSSetExactSolution()` and `PetscDSSetExactSolutionTimeDerivative()` to another `PetscDS`

3840:   Not Collective

3842:   Input Parameter:
3843: . ds - The `PetscDS` object

3845:   Output Parameter:
3846: . newds - The `PetscDS` copy

3848:   Level: intermediate

3850: .seealso: `PetscDS`, `PetscDSCopyBoundary()`, `PetscDSCopyEquations()`, `PetscDSCopyBounds()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
3851: @*/
3852: PetscErrorCode PetscDSCopyExactSolutions(PetscDS ds, PetscDS newds)
3853: {
3854:   PetscSimplePointFn *sol;
3855:   void               *ctx;
3856:   PetscInt            Nf, f;

3858:   PetscFunctionBegin;
3861:   PetscCall(PetscDSGetNumFields(ds, &Nf));
3862:   for (f = 0; f < Nf; ++f) {
3863:     PetscCall(PetscDSGetExactSolution(ds, f, &sol, &ctx));
3864:     PetscCall(PetscDSSetExactSolution(newds, f, sol, ctx));
3865:     PetscCall(PetscDSGetExactSolutionTimeDerivative(ds, f, &sol, &ctx));
3866:     PetscCall(PetscDSSetExactSolutionTimeDerivative(newds, f, sol, ctx));
3867:   }
3868:   PetscFunctionReturn(PETSC_SUCCESS);
3869: }

3871: /*@
3872:   PetscDSCopyBounds - Copy lower and upper solution bounds set with `PetscDSSetLowerBound()` and `PetscDSSetLowerBound()` to another `PetscDS`

3874:   Not Collective

3876:   Input Parameter:
3877: . ds - The `PetscDS` object

3879:   Output Parameter:
3880: . newds - The `PetscDS` copy

3882:   Level: intermediate

3884: .seealso: `PetscDS`, `PetscDSCopyBoundary()`, `PetscDSCopyEquations()`, `PetscDSCopyExactSolutions()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
3885: @*/
3886: PetscErrorCode PetscDSCopyBounds(PetscDS ds, PetscDS newds)
3887: {
3888:   PetscSimplePointFn *bound;
3889:   void               *ctx;
3890:   PetscInt            Nf, f;

3892:   PetscFunctionBegin;
3895:   PetscCall(PetscDSGetNumFields(ds, &Nf));
3896:   for (f = 0; f < Nf; ++f) {
3897:     PetscCall(PetscDSGetLowerBound(ds, f, &bound, &ctx));
3898:     PetscCall(PetscDSSetLowerBound(newds, f, bound, ctx));
3899:     PetscCall(PetscDSGetUpperBound(ds, f, &bound, &ctx));
3900:     PetscCall(PetscDSSetUpperBound(newds, f, bound, ctx));
3901:   }
3902:   PetscFunctionReturn(PETSC_SUCCESS);
3903: }

3905: PetscErrorCode PetscDSCopy(PetscDS ds, PetscInt minDegree, PetscInt maxDegree, DM dmNew, PetscDS dsNew)
3906: {
3907:   DSBoundary b;
3908:   PetscInt   cdim, Nf, f, d;
3909:   PetscBool  isCohesive;
3910:   void      *ctx;

3912:   PetscFunctionBegin;
3913:   PetscCall(PetscDSCopyConstants(ds, dsNew));
3914:   PetscCall(PetscDSCopyExactSolutions(ds, dsNew));
3915:   PetscCall(PetscDSCopyBounds(ds, dsNew));
3916:   PetscCall(PetscDSSelectDiscretizations(ds, PETSC_DETERMINE, NULL, minDegree, maxDegree, dsNew));
3917:   PetscCall(PetscDSCopyEquations(ds, dsNew));
3918:   PetscCall(PetscDSGetNumFields(ds, &Nf));
3919:   for (f = 0; f < Nf; ++f) {
3920:     PetscCall(PetscDSGetContext(ds, f, &ctx));
3921:     PetscCall(PetscDSSetContext(dsNew, f, ctx));
3922:     PetscCall(PetscDSGetCohesive(ds, f, &isCohesive));
3923:     PetscCall(PetscDSSetCohesive(dsNew, f, isCohesive));
3924:     PetscCall(PetscDSGetJetDegree(ds, f, &d));
3925:     PetscCall(PetscDSSetJetDegree(dsNew, f, d));
3926:   }
3927:   if (Nf) {
3928:     PetscCall(PetscDSGetCoordinateDimension(ds, &cdim));
3929:     PetscCall(PetscDSSetCoordinateDimension(dsNew, cdim));
3930:   }
3931:   PetscCall(PetscDSCopyBoundary(ds, PETSC_DETERMINE, NULL, dsNew));
3932:   for (b = dsNew->boundary; b; b = b->next) {
3933:     PetscCall(DMGetLabel(dmNew, b->lname, &b->label));
3934:     /* Do not check if label exists here, since p4est calls this for the reference tree which does not have the labels */
3935:     //PetscCheck(b->label,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Label %s missing in new DM", name);
3936:   }
3937:   PetscFunctionReturn(PETSC_SUCCESS);
3938: }

3940: PetscErrorCode PetscDSGetHeightSubspace(PetscDS prob, PetscInt height, PetscDS *subprob)
3941: {
3942:   PetscInt dim, Nf, f;

3944:   PetscFunctionBegin;
3946:   PetscAssertPointer(subprob, 3);
3947:   if (height == 0) {
3948:     *subprob = prob;
3949:     PetscFunctionReturn(PETSC_SUCCESS);
3950:   }
3951:   PetscCall(PetscDSGetNumFields(prob, &Nf));
3952:   PetscCall(PetscDSGetSpatialDimension(prob, &dim));
3953:   PetscCheck(height <= dim, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_OUTOFRANGE, "DS can only handle height in [0, %" PetscInt_FMT "], not %" PetscInt_FMT, dim, height);
3954:   if (!prob->subprobs) PetscCall(PetscCalloc1(dim, &prob->subprobs));
3955:   if (!prob->subprobs[height - 1]) {
3956:     PetscInt cdim;

3958:     PetscCall(PetscDSCreate(PetscObjectComm((PetscObject)prob), &prob->subprobs[height - 1]));
3959:     PetscCall(PetscDSGetCoordinateDimension(prob, &cdim));
3960:     PetscCall(PetscDSSetCoordinateDimension(prob->subprobs[height - 1], cdim));
3961:     for (f = 0; f < Nf; ++f) {
3962:       PetscFE      subfe;
3963:       PetscObject  obj;
3964:       PetscClassId id;

3966:       PetscCall(PetscDSGetDiscretization(prob, f, &obj));
3967:       PetscCall(PetscObjectGetClassId(obj, &id));
3968:       PetscCheck(id == PETSCFE_CLASSID, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONG, "Unsupported discretization type for field %" PetscInt_FMT, f);
3969:       PetscCall(PetscFEGetHeightSubspace((PetscFE)obj, height, &subfe));
3970:       PetscCall(PetscDSSetDiscretization(prob->subprobs[height - 1], f, (PetscObject)subfe));
3971:     }
3972:   }
3973:   *subprob = prob->subprobs[height - 1];
3974:   PetscFunctionReturn(PETSC_SUCCESS);
3975: }

3977: PetscErrorCode PetscDSPermuteQuadPoint(PetscDS ds, PetscInt ornt, PetscInt field, PetscInt q, PetscInt *qperm)
3978: {
3979:   IS              permIS;
3980:   PetscQuadrature quad;
3981:   DMPolytopeType  ct;
3982:   const PetscInt *perm;
3983:   PetscInt        Na, Nq;

3985:   PetscFunctionBeginHot;
3986:   PetscCall(PetscFEGetQuadrature((PetscFE)ds->disc[field], &quad));
3987:   PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL));
3988:   PetscCall(PetscQuadratureGetCellType(quad, &ct));
3989:   PetscCheck(q >= 0 && q < Nq, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Quadrature point %" PetscInt_FMT " is not in [0, %" PetscInt_FMT ")", q, Nq);
3990:   Na = DMPolytopeTypeGetNumArrangements(ct) / 2;
3991:   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);
3992:   if (!ds->quadPerm[(PetscInt)ct]) PetscCall(PetscQuadratureComputePermutations(quad, NULL, &ds->quadPerm[(PetscInt)ct]));
3993:   permIS = ds->quadPerm[(PetscInt)ct][ornt + Na];
3994:   PetscCall(ISGetIndices(permIS, &perm));
3995:   *qperm = perm[q];
3996:   PetscCall(ISRestoreIndices(permIS, &perm));
3997:   PetscFunctionReturn(PETSC_SUCCESS);
3998: }

4000: PetscErrorCode PetscDSGetDiscType_Internal(PetscDS ds, PetscInt f, PetscDiscType *disctype)
4001: {
4002:   PetscObject  obj;
4003:   PetscClassId id;
4004:   PetscInt     Nf;

4006:   PetscFunctionBegin;
4008:   PetscAssertPointer(disctype, 3);
4009:   *disctype = PETSC_DISC_NONE;
4010:   PetscCall(PetscDSGetNumFields(ds, &Nf));
4011:   PetscCheck(f < Nf, PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_SIZ, "Field %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, Nf);
4012:   PetscCall(PetscDSGetDiscretization(ds, f, &obj));
4013:   if (obj) {
4014:     PetscCall(PetscObjectGetClassId(obj, &id));
4015:     if (id == PETSCFE_CLASSID) *disctype = PETSC_DISC_FE;
4016:     else *disctype = PETSC_DISC_FV;
4017:   }
4018:   PetscFunctionReturn(PETSC_SUCCESS);
4019: }

4021: static PetscErrorCode PetscDSDestroy_Basic(PetscDS ds)
4022: {
4023:   PetscFunctionBegin;
4024:   PetscCall(PetscFree(ds->data));
4025:   PetscFunctionReturn(PETSC_SUCCESS);
4026: }

4028: static PetscErrorCode PetscDSInitialize_Basic(PetscDS ds)
4029: {
4030:   PetscFunctionBegin;
4031:   ds->ops->setfromoptions = NULL;
4032:   ds->ops->setup          = NULL;
4033:   ds->ops->view           = NULL;
4034:   ds->ops->destroy        = PetscDSDestroy_Basic;
4035:   PetscFunctionReturn(PETSC_SUCCESS);
4036: }

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

4041:   Level: intermediate

4043: .seealso: `PetscDSType`, `PetscDSCreate()`, `PetscDSSetType()`
4044: M*/

4046: PETSC_EXTERN PetscErrorCode PetscDSCreate_Basic(PetscDS ds)
4047: {
4048:   PetscDS_Basic *b;

4050:   PetscFunctionBegin;
4052:   PetscCall(PetscNew(&b));
4053:   ds->data = b;

4055:   PetscCall(PetscDSInitialize_Basic(ds));
4056:   PetscFunctionReturn(PETSC_SUCCESS);
4057: }