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;

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

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

223:   Collective

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

230:   Options Database Key:
231: . -name [viewertype][:...] - option name and values. See `PetscObjectViewFromOptions()` for the possible arguments

233:   Level: intermediate

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

245: /*@
246:   PetscDSView - Views a `PetscDS`

248:   Collective

250:   Input Parameters:
251: + prob - the `PetscDS` object to view
252: - v    - the viewer

254:   Level: developer

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

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

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

275:   Collective

277:   Input Parameter:
278: . prob - the `PetscDS` object to set options for

280:   Options Database Keys:
281: + -petscds_type type            - Set the `PetscDS` type
282: . -petscds_view                 - View the `PetscDS`
283: . -petscds_jac_pre (true|false) - Turn formation of a separate Jacobian preconditioner on or off
284: . -bc_NAME ids                  - comma separated list of label ids for the boundary condition NAME
285: - -bc_NAME_comp comps           - comma separated list of field components to constrain for the boundary condition NAME

287:   Level: intermediate

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

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

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

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

351: /*@
352:   PetscDSSetUp - Construct data structures for the `PetscDS`

354:   Collective

356:   Input Parameter:
357: . prob - the `PetscDS` object to setup

359:   Level: developer

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

614: /*@
615:   PetscDSDestroy - Destroys a `PetscDS` object

617:   Collective

619:   Input Parameter:
620: . ds - the `PetscDS` object to destroy

622:   Level: developer

624: .seealso: `PetscDSView()`
625: @*/
626: PetscErrorCode PetscDSDestroy(PetscDS *ds)
627: {
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;

640:     PetscCall(PetscDSGetSpatialDimension(*ds, &dim));
641:     for (PetscInt d = 0; d < dim; ++d) PetscCall(PetscDSDestroy(&(*ds)->subprobs[d]));
642:   }
643:   PetscCall(PetscFree((*ds)->subprobs));
644:   PetscCall(PetscDSDestroyStructs_Static(*ds));
645:   for (PetscInt 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:   PetscFunctionBegin;
898:   PetscAssertPointer(numCohesive, 2);
899:   *numCohesive = 0;
900:   for (PetscInt f = 0; f < ds->Nf; ++f) *numCohesive += ds->cohesive[f] ? 1 : 0;
901:   PetscFunctionReturn(PETSC_SUCCESS);
902: }

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

907:   Not Collective

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

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

916:   Level: developer

918: .seealso: `PetscDS`, `PetscDSSetCohesive()`, `PetscDSIsCohesive()`, `PetscDSCreate()`
919: @*/
920: PetscErrorCode PetscDSGetCohesive(PetscDS ds, PetscInt f, PetscBool *isCohesive)
921: {
922:   PetscFunctionBegin;
924:   PetscAssertPointer(isCohesive, 3);
925:   PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
926:   *isCohesive = ds->cohesive[f];
927:   PetscFunctionReturn(PETSC_SUCCESS);
928: }

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

933:   Not Collective

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

940:   Level: developer

942: .seealso: `PetscDS`, `PetscDSGetCohesive()`, `PetscDSIsCohesive()`, `PetscDSCreate()`
943: @*/
944: PetscErrorCode PetscDSSetCohesive(PetscDS ds, PetscInt f, PetscBool isCohesive)
945: {
946:   PetscFunctionBegin;
948:   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);
949:   ds->cohesive[f] = isCohesive;
950:   ds->isCohesive  = PETSC_FALSE;
951:   for (PetscInt i = 0; i < ds->Nf; ++i) ds->isCohesive = ds->isCohesive || ds->cohesive[f] ? PETSC_TRUE : PETSC_FALSE;
952:   PetscFunctionReturn(PETSC_SUCCESS);
953: }

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

958:   Not Collective

960:   Input Parameter:
961: . prob - The `PetscDS` object

963:   Output Parameter:
964: . dim - The total problem dimension

966:   Level: beginner

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

980: /*@
981:   PetscDSGetTotalComponents - Returns the total number of components in this system

983:   Not Collective

985:   Input Parameter:
986: . prob - The `PetscDS` object

988:   Output Parameter:
989: . Nc - The total number of components

991:   Level: beginner

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

1005: /*@
1006:   PetscDSGetDiscretization - Returns the discretization object for the given field

1008:   Not Collective

1010:   Input Parameters:
1011: + prob - The `PetscDS` object
1012: - f    - The field number

1014:   Output Parameter:
1015: . disc - The discretization object, this can be a `PetscFE` or a `PetscFV`

1017:   Level: beginner

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

1031: /*@
1032:   PetscDSSetDiscretization - Sets the discretization object for the given field

1034:   Not Collective

1036:   Input Parameters:
1037: + prob - The `PetscDS` object
1038: . f    - The field number
1039: - disc - The discretization object, this can be a `PetscFE` or a `PetscFV`

1041:   Level: beginner

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

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

1069: /*@
1070:   PetscDSGetWeakForm - Returns the weak form object from within the `PetscDS`

1072:   Not Collective

1074:   Input Parameter:
1075: . ds - The `PetscDS` object

1077:   Output Parameter:
1078: . wf - The weak form object

1080:   Level: beginner

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

1093: /*@
1094:   PetscDSSetWeakForm - Sets the weak form object to be used by the `PetscDS`

1096:   Not Collective

1098:   Input Parameters:
1099: + ds - The `PetscDS` object
1100: - wf - The weak form object

1102:   Level: beginner

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

1118: /*@
1119:   PetscDSAddDiscretization - Adds a discretization object

1121:   Not Collective

1123:   Input Parameters:
1124: + prob - The `PetscDS` object
1125: - disc - The discretization object, this can be a `PetscFE` or `PetscFV`

1127:   Level: beginner

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

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

1141:   Not Collective

1143:   Input Parameter:
1144: . prob - The `PetscDS` object

1146:   Output Parameter:
1147: . q - The quadrature object

1149:   Level: intermediate

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

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

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

1172:   Not Collective

1174:   Input Parameters:
1175: + prob - The `PetscDS` object
1176: - f    - The field number

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

1181:   Level: developer

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

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

1198:   Not Collective

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

1205:   Level: developer

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

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

1221:   Not Collective

1223:   Input Parameters:
1224: + ds - The `PetscDS` object
1225: - f  - The field number

1227:   Output Parameter:
1228: . k - The highest derivative we need to tabulate

1230:   Level: developer

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

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

1247:   Not Collective

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

1254:   Level: developer

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

1267: /*@C
1268:   PetscDSGetObjective - Get the pointwise objective function for a given test field that was provided with `PetscDSSetObjective()`

1270:   Not Collective

1272:   Input Parameters:
1273: + ds - The `PetscDS`
1274: - f  - The test field number

1276:   Output Parameter:
1277: . obj - integrand for the test function term, see `PetscPointFn`

1279:   Level: intermediate

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

1284: .seealso: `PetscPointFn`, `PetscDS`, `PetscDSSetObjective()`, `PetscDSGetResidual()`
1285: @*/
1286: PetscErrorCode PetscDSGetObjective(PetscDS ds, PetscInt f, PetscPointFn **obj)
1287: {
1288:   PetscPointFn **tmp;
1289:   PetscInt       n;

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

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

1303:   Not Collective

1305:   Input Parameters:
1306: + ds  - The `PetscDS`
1307: . f   - The test field number
1308: - obj - integrand for the test function term, see `PetscPointFn`

1310:   Level: intermediate

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

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

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

1330:   Not Collective

1332:   Input Parameters:
1333: + ds - The `PetscDS`
1334: - f  - The test field number

1336:   Output Parameters:
1337: + f0 - integrand for the test function term, see `PetscPointFn`
1338: - f1 - integrand for the test function gradient term, see `PetscPointFn`

1340:   Level: intermediate

1342:   Note:
1343:   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)$

1345: .seealso: `PetscPointFn`, `PetscDS`, `PetscDSSetResidual()`
1346: @*/
1347: PetscErrorCode PetscDSGetResidual(PetscDS ds, PetscInt f, PetscPointFn **f0, PetscPointFn **f1)
1348: {
1349:   PetscPointFn **tmp0, **tmp1;
1350:   PetscInt       n0, n1;

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

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

1364:   Not Collective

1366:   Input Parameters:
1367: + ds - The `PetscDS`
1368: . f  - The test field number
1369: . f0 - integrand for the test function term, see `PetscPointFn`
1370: - f1 - integrand for the test function gradient term, see `PetscPointFn`

1372:   Level: intermediate

1374:   Note:
1375:   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)$

1377: .seealso: `PetscPointFn`, `PetscDS`, `PetscDSGetResidual()`
1378: @*/
1379: PetscErrorCode PetscDSSetResidual(PetscDS ds, PetscInt f, PetscPointFn *f0, PetscPointFn *f1)
1380: {
1381:   PetscFunctionBegin;
1385:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
1386:   PetscCall(PetscWeakFormSetIndexResidual(ds->wf, NULL, 0, f, 0, 0, f0, 0, f1));
1387:   PetscFunctionReturn(PETSC_SUCCESS);
1388: }

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

1393:   Not Collective

1395:   Input Parameters:
1396: + ds - The `PetscDS`
1397: - f  - The test field number

1399:   Output Parameters:
1400: + f0 - integrand for the test function term, see `PetscPointFn`
1401: - f1 - integrand for the test function gradient term, see `PetscPointFn`

1403:   Level: intermediate

1405:   Note:
1406:   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)$

1408: .seealso: `PetscPointFn`, `PetscDS`, `PetscDSSetRHSResidual()`
1409: @*/
1410: PetscErrorCode PetscDSGetRHSResidual(PetscDS ds, PetscInt f, PetscPointFn **f0, PetscPointFn **f1)
1411: {
1412:   PetscPointFn **tmp0, **tmp1;
1413:   PetscInt       n0, n1;

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

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

1427:   Not Collective

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

1435:   Level: intermediate

1437:   Note:
1438:   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)$

1440: .seealso: `PetscDS`, `PetscDSGetResidual()`
1441: @*/
1442: PetscErrorCode PetscDSSetRHSResidual(PetscDS ds, PetscInt f, PetscPointFn *f0, PetscPointFn *f1)
1443: {
1444:   PetscFunctionBegin;
1448:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
1449:   PetscCall(PetscWeakFormSetIndexResidual(ds->wf, NULL, 0, f, 100, 0, f0, 0, f1));
1450:   PetscFunctionReturn(PETSC_SUCCESS);
1451: }

1453: /*@
1454:   PetscDSHasJacobian - Checks that the Jacobian functions have been set

1456:   Not Collective

1458:   Input Parameter:
1459: . ds - The `PetscDS`

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

1464:   Level: intermediate

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

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

1479:   Not Collective

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

1486:   Output Parameters:
1487: + g0 - integrand for the test and basis function term, see `PetscPointJacFn`
1488: . g1 - integrand for the test function and basis function gradient term, see `PetscPointJacFn`
1489: . g2 - integrand for the test function gradient and basis function term, see `PetscPointJacFn`
1490: - g3 - integrand for the test function gradient and basis function gradient term, see `PetscPointJacFn`

1492:   Level: intermediate

1494:   Note:
1495:   We are using a first order FEM model for the weak form\:

1497:   $$
1498:   \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
1499:   + \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
1500:   $$

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

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

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

1524:   Not Collective

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

1535:   Level: intermediate

1537:   Note:
1538:   We are using a first order FEM model for the weak form\:

1540:   $$
1541:   \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
1542:   + \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
1543:   $$

1545: .seealso: `PetscDS`, `PetscDSGetJacobian()`, `PetscPointJacFn`
1546: @*/
1547: PetscErrorCode PetscDSSetJacobian(PetscDS ds, PetscInt f, PetscInt g, PetscPointJacFn *g0, PetscPointJacFn *g1, PetscPointJacFn *g2, PetscPointJacFn *g3)
1548: {
1549:   PetscFunctionBegin;
1555:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
1556:   PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", g);
1557:   PetscCall(PetscWeakFormSetIndexJacobian(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3));
1558:   PetscFunctionReturn(PETSC_SUCCESS);
1559: }

1561: /*@
1562:   PetscDSUseJacobianPreconditioner - Set whether to construct a Jacobian preconditioner

1564:   Not Collective

1566:   Input Parameters:
1567: + prob      - The `PetscDS`
1568: - useJacPre - flag that enables construction of a Jacobian preconditioner

1570:   Level: intermediate

1572:   Developer Note:
1573:   Should be called `PetscDSSetUseJacobianPreconditioner()`

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

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

1588:   Not Collective

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

1593:   Output Parameter:
1594: . hasJacPre - the flag

1596:   Level: intermediate

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

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

1614:   Not Collective

1616:   Input Parameters:
1617: + ds - The `PetscDS`
1618: . f  - The test field number
1619: - g  - The field number

1621:   Output Parameters:
1622: + g0 - integrand for the test and basis function term, see `PetscPointJacFn`
1623: . g1 - integrand for the test function and basis function gradient term, see `PetscPointJacFn`
1624: . g2 - integrand for the test function gradient and basis function term, see `PetscPointJacFn`
1625: - g3 - integrand for the test function gradient and basis function gradient term, see `PetscPointJacFn`

1627:   Level: intermediate

1629:   Note:
1630:   We are using a first order FEM model for the weak form\:

1632:   $$
1633:   \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
1634:   + \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
1635:   $$

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

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

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

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

1663:   Not Collective

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

1674:   Level: intermediate

1676:   Note:
1677:   We are using a first order FEM model for the weak form\:

1679:   $$
1680:   \int_\Omega \phi\, g_0(u, u_t, \nabla u, x, t) \psi + \phi\, {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi
1681:   + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi
1682:   $$

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

1687: .seealso: `PetscDS`, `PetscDSGetJacobianPreconditioner()`, `PetscDSSetJacobian()`, `PetscPointJacFn`
1688: @*/
1689: PetscErrorCode PetscDSSetJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g, PetscPointJacFn *g0, PetscPointJacFn *g1, PetscPointJacFn *g2, PetscPointJacFn *g3)
1690: {
1691:   PetscFunctionBegin;
1697:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
1698:   PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", g);
1699:   PetscCall(PetscWeakFormSetIndexJacobianPreconditioner(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3));
1700:   PetscFunctionReturn(PETSC_SUCCESS);
1701: }

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

1706:   Not Collective

1708:   Input Parameter:
1709: . ds - The `PetscDS`

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

1714:   Level: intermediate

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

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

1729:   Not Collective

1731:   Input Parameters:
1732: + ds - The `PetscDS`
1733: . f  - The test field number
1734: - g  - The field number

1736:   Output Parameters:
1737: + g0 - integrand for the test and basis function term, see `PetscPointJacFn`
1738: . g1 - integrand for the test function and basis function gradient term, see `PetscPointJacFn`
1739: . g2 - integrand for the test function gradient and basis function term, see `PetscPointJacFn`
1740: - g3 - integrand for the test function gradient and basis function gradient term, see `PetscPointJacFn`

1742:   Level: intermediate

1744:   Note:
1745:   We are using a first order FEM model for the weak form\:

1747:   $$
1748:   \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
1749:   + \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
1750:   $$

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

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

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

1774:   Not Collective

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

1785:   Level: intermediate

1787:   Note:
1788:   We are using a first order FEM model for the weak form\:

1790:   $$
1791:   \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
1792:   + \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
1793:   $$

1795: .seealso: `PetscDS`, `PetscDSGetDynamicJacobian()`, `PetscDSGetJacobian()`, `PetscPointJacFn`
1796: @*/
1797: PetscErrorCode PetscDSSetDynamicJacobian(PetscDS ds, PetscInt f, PetscInt g, PetscPointJacFn *g0, PetscPointJacFn *g1, PetscPointJacFn *g2, PetscPointJacFn *g3)
1798: {
1799:   PetscFunctionBegin;
1805:   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
1806:   PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", g);
1807:   PetscCall(PetscWeakFormSetIndexDynamicJacobian(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3));
1808:   PetscFunctionReturn(PETSC_SUCCESS);
1809: }

1811: /*@C
1812:   PetscDSGetRiemannSolver - Returns the Riemann solver for the given field

1814:   Not Collective

1816:   Input Parameters:
1817: + ds - The `PetscDS` object
1818: - f  - The field number

1820:   Output Parameter:
1821: . r - Riemann solver, see `PetscRiemannFn`

1823:   Level: intermediate

1825: .seealso: `PetscDS`, `PetscRiemannFn`, `PetscDSSetRiemannSolver()`
1826: @*/
1827: PetscErrorCode PetscDSGetRiemannSolver(PetscDS ds, PetscInt f, PetscRiemannFn **r)
1828: {
1829:   PetscRiemannFn **tmp;
1830:   PetscInt         n;

1832:   PetscFunctionBegin;
1834:   PetscAssertPointer(r, 3);
1835:   PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
1836:   PetscCall(PetscWeakFormGetRiemannSolver(ds->wf, NULL, 0, f, 0, &n, &tmp));
1837:   *r = tmp ? tmp[0] : NULL;
1838:   PetscFunctionReturn(PETSC_SUCCESS);
1839: }

1841: /*@C
1842:   PetscDSSetRiemannSolver - Sets the Riemann solver for the given field

1844:   Not Collective

1846:   Input Parameters:
1847: + ds - The `PetscDS` object
1848: . f  - The field number
1849: - r  - Riemann solver, see `PetscRiemannFn`

1851:   Level: intermediate

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

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

1868:   Not Collective

1870:   Input Parameters:
1871: + ds - The `PetscDS`
1872: - f  - The field number

1874:   Output Parameter:
1875: . update - update function, see `PetscPointFn`

1877:   Level: intermediate

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

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

1896:   Not Collective

1898:   Input Parameters:
1899: + ds     - The `PetscDS`
1900: . f      - The field number
1901: - update - update function, see `PetscPointFn`

1903:   Level: intermediate

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

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

1921:   Not Collective

1923:   Input Parameters:
1924: + ds  - The `PetscDS`
1925: . f   - The field number
1926: - ctx - the context

1928:   Level: intermediate

1930:   Fortran Notes:
1931:   This only works when the context is a Fortran derived type or a `PetscObject`. Define `ctx` with
1932: .vb
1933:   type(tUsertype), pointer :: ctx
1934: .ve

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

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

1951:   Not Collective

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

1958:   Level: intermediate

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

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

1975:   Not Collective

1977:   Input Parameters:
1978: + ds - The PetscDS
1979: - f  - The test field number

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

1985:   Level: intermediate

1987:   Note:
1988:   We are using a first order FEM model for the weak form\:

1990:   $$
1991:   \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
1992:   $$

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

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

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

2013:   Not Collective

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

2021:   Level: intermediate

2023:   Note:
2024:   We are using a first order FEM model for the weak form\:

2026:   $$
2027:   \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
2028:   $$

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

2041: /*@
2042:   PetscDSHasBdJacobian - Indicates that boundary Jacobian functions have been set

2044:   Not Collective

2046:   Input Parameter:
2047: . ds - The `PetscDS`

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

2052:   Level: intermediate

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

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

2068:   Not Collective

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

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

2081:   Level: intermediate

2083:   Note:
2084:   We are using a first order FEM model for the weak form\:

2086:   $$
2087:   \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
2088:   + \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
2089:   $$

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

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

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

2113:   Not Collective

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

2124:   Level: intermediate

2126:   Note:
2127:   We are using a first order FEM model for the weak form\:

2129:   $$
2130:   \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
2131:   + \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
2132:   $$

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

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

2153:   Not Collective

2155:   Input Parameter:
2156: . ds - The `PetscDS`

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

2161:   Level: intermediate

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

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

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

2181:   Not Collective; No Fortran Support

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

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

2194:   Level: intermediate

2196:   Note:
2197:   We are using a first order FEM model for the weak form\:

2199:   $$
2200:   \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
2201:   + \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
2202:   $$

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

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

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

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

2230:   Not Collective; No Fortran Support

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

2241:   Level: intermediate

2243:   Note:
2244:   We are using a first order FEM model for the weak form\:

2246:   $$
2247:   \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
2248:   + \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
2249:   $$

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

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

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

2273:   Not Collective

2275:   Input Parameters:
2276: + prob - The `PetscDS`
2277: - f    - The test field number

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

2283:   Level: intermediate

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

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

2306:   Not Collective

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

2314:   Level: intermediate

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

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

2338:   Not Collective

2340:   Input Parameters:
2341: + prob - The `PetscDS`
2342: - f    - The test field number

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

2348:   Level: intermediate

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

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

2371:   Not Collective

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

2379:   Level: intermediate

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

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

2403:   Not Collective

2405:   Input Parameters:
2406: + ds - The PetscDS
2407: - f  - The field number

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

2413:   Level: intermediate

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

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

2436:   Not Collective

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

2444:   Level: intermediate

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

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

2468:   Not Collective

2470:   Input Parameters:
2471: + ds - The `PetscDS`
2472: - f  - The field number

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

2478:   Level: intermediate

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

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

2501:   Not Collective

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

2509:   Level: intermediate

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

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

2533:   Not Collective

2535:   Input Parameter:
2536: . ds - The `PetscDS` object

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

2542:   Level: intermediate

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

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

2564:   Not Collective

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

2571:   Level: intermediate

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

2591: /*@C
2592:   PetscDSSetIntegrationParameters - Set the parameters for a particular integration

2594:   Not Collective

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

2601:   Level: intermediate

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

2614: /*@C
2615:   PetscDSSetCellParameters - Set the parameters for a particular cell

2617:   Not Collective

2619:   Input Parameters:
2620: + ds     - The `PetscDS` object
2621: - volume - The cell volume

2623:   Level: intermediate

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

2635: /*@
2636:   PetscDSGetFieldIndex - Returns the index of the given field

2638:   Not Collective

2640:   Input Parameters:
2641: + prob - The `PetscDS` object
2642: - disc - The discretization object

2644:   Output Parameter:
2645: . f - The field number

2647:   Level: beginner

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

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

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

2670:   Not Collective

2672:   Input Parameters:
2673: + prob - The `PetscDS` object
2674: - f    - The field number

2676:   Output Parameter:
2677: . size - The size

2679:   Level: beginner

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

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

2697:   Not Collective

2699:   Input Parameters:
2700: + prob - The `PetscDS` object
2701: - f    - The field number

2703:   Output Parameter:
2704: . off - The offset

2706:   Level: beginner

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

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

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

2729:   Not Collective

2731:   Input Parameters:
2732: + ds - The `PetscDS` object
2733: - f  - The field number

2735:   Output Parameter:
2736: . off - The offset

2738:   Level: beginner

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

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

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

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

2764:   Not Collective

2766:   Input Parameter:
2767: . prob - The `PetscDS` object

2769:   Output Parameter:
2770: . dimensions - The number of dimensions

2772:   Level: beginner

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

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

2789:   Not Collective

2791:   Input Parameter:
2792: . prob - The `PetscDS` object

2794:   Output Parameter:
2795: . components - The number of components

2797:   Level: beginner

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

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

2814:   Not Collective

2816:   Input Parameters:
2817: + prob - The `PetscDS` object
2818: - f    - The field number

2820:   Output Parameter:
2821: . off - The offset

2823:   Level: beginner

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

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

2841:   Not Collective

2843:   Input Parameter:
2844: . prob - The `PetscDS` object

2846:   Output Parameter:
2847: . offsets - The offsets

2849:   Level: beginner

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

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

2866:   Not Collective

2868:   Input Parameter:
2869: . prob - The `PetscDS` object

2871:   Output Parameter:
2872: . offsets - The offsets

2874:   Level: beginner

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

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

2891:   Not Collective

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

2897:   Output Parameter:
2898: . offsets - The offsets

2900:   Level: beginner

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

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

2919:   Not Collective

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

2925:   Output Parameter:
2926: . offsets - The offsets

2928:   Level: beginner

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

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

2947:   Not Collective

2949:   Input Parameter:
2950: . prob - The `PetscDS` object

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

2955:   Level: intermediate

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

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

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

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

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

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

2992:   Not Collective

2994:   Input Parameter:
2995: . prob - The `PetscDS` object

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

3000:   Level: intermediate

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

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

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

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

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

3097: /*@C
3098:   PetscDSAddBoundary - Add a boundary condition to the model.

3100:   Collective

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

3116:   Output Parameter:
3117: . bd - The boundary number

3119:   Options Database Keys:
3120: + -bc_NAME values     - comma separated list of values for the boundary condition NAME
3121: - -bc_NAME_comp comps - comma separated list of components for the boundary condition NAME

3123:   Level: developer

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

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

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

3162: .seealso: `PetscDS`, `PetscWeakForm`, `DMLabel`, `DMBoundaryConditionType`, `PetscDSAddBoundaryByName()`, `PetscDSGetBoundary()`, `PetscDSSetResidual()`, `PetscDSSetBdResidual()`
3163: @*/
3164: PetscErrorCode PetscDSAddBoundary(PetscDS ds, DMBoundaryConditionType type, const char name[], DMLabel label, PetscInt Nv, const PetscInt values[], PetscInt field, PetscInt Nc, const PetscInt comps[], PetscVoidFn *bcFunc, PetscVoidFn *bcFunc_t, PetscCtx ctx, PetscInt *bd)
3165: {
3166:   DSBoundary  head = ds->boundary, b;
3167:   PetscInt    n    = 0;
3168:   const char *lname;

3170:   PetscFunctionBegin;
3173:   PetscAssertPointer(name, 3);
3178:   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);
3179:   if (Nc > 0) {
3180:     PetscInt *fcomps;

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

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

3228:   Collective

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

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

3247:   Options Database Keys:
3248: + -bc_NAME values     - comma separated list of values for the boundary condition NAME
3249: - -bc_NAME_comp comps - comma separated list of components for the boundary condition NAME

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

3281:   Level: developer

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

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

3291: .seealso: `PetscDS`, `PetscWeakForm`, `DMLabel`, `DMBoundaryConditionType`, `PetscDSAddBoundary()`, `PetscDSGetBoundary()`, `PetscDSSetResidual()`, `PetscDSSetBdResidual()`
3292: @*/
3293: PetscErrorCode PetscDSAddBoundaryByName(PetscDS ds, DMBoundaryConditionType type, const char name[], const char lname[], PetscInt Nv, const PetscInt values[], PetscInt field, PetscInt Nc, const PetscInt comps[], PetscVoidFn *bcFunc, PetscVoidFn *bcFunc_t, PetscCtx ctx, PetscInt *bd)
3294: {
3295:   DSBoundary head = ds->boundary, b;
3296:   PetscInt   n    = 0;

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

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

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

3359:   Level: developer

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

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

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

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

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

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

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

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

3426:   Level: intermediate

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

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

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

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

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

3466:   Level: developer

3468: .seealso: `PetscDS`, `PetscWeakForm`, `DMBoundaryConditionType`, `PetscDSAddBoundary()`, `DMLabel`
3469: @*/
3470: 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[], PetscVoidFn **func, PetscVoidFn **func_t, void **ctx)
3471: {
3472:   DSBoundary b = ds->boundary;
3473:   PetscInt   n = 0;

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

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

3537:   Not Collective

3539:   Input Parameters:
3540: + ds - The source `PetscDS` object
3541: - dm - The `DM` holding labels

3543:   Level: intermediate

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

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

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

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

3586:   Not Collective

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

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

3596:   Level: intermediate

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

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

3614:     if (numFields > 0 && fields) {
3615:       PetscInt f;

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

3630: /*@
3631:   PetscDSDestroyBoundary - Remove all `DMBoundary` objects from the `PetscDS`

3633:   Not Collective

3635:   Input Parameter:
3636: . ds - The `PetscDS` object

3638:   Level: intermediate

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

3646:   PetscFunctionBegin;
3647:   while (next) {
3648:     DSBoundary b = next;

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

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

3664:   Not Collective

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

3673:   Output Parameter:
3674: . newprob - The `PetscDS` copy

3676:   Level: intermediate

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

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

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

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

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

3715:   Not Collective

3717:   Input Parameters:
3718: + prob      - The `PetscDS` object
3719: . numFields - Number of new fields
3720: - fields    - Old field number for each new field

3722:   Output Parameter:
3723: . newprob - The `PetscDS` copy

3725:   Level: intermediate

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

3733:   PetscFunctionBegin;
3735:   if (fields) PetscAssertPointer(fields, 3);
3737:   PetscCall(PetscDSGetNumFields(prob, &Nf));
3738:   PetscCall(PetscDSGetNumFields(newprob, &Nfn));
3739:   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);
3740:   for (fn = 0; fn < numFields; ++fn) {
3741:     const PetscInt  f = fields ? fields[fn] : fn;
3742:     PetscPointFn   *obj;
3743:     PetscPointFn   *f0, *f1;
3744:     PetscBdPointFn *f0Bd, *f1Bd;
3745:     PetscRiemannFn *r;

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

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

3774: /*@
3775:   PetscDSCopyEquations - Copy all pointwise function pointers to another `PetscDS`

3777:   Not Collective

3779:   Input Parameter:
3780: . prob - The `PetscDS` object

3782:   Output Parameter:
3783: . newprob - The `PetscDS` copy

3785:   Level: intermediate

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

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

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

3809:   Not Collective

3811:   Input Parameter:
3812: . prob - The `PetscDS` object

3814:   Output Parameter:
3815: . newprob - The `PetscDS` copy

3817:   Level: intermediate

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

3826:   PetscFunctionBegin;
3829:   PetscCall(PetscDSGetConstants(prob, &Nc, &constants));
3830:   PetscCall(PetscDSSetConstants(newprob, Nc, (PetscScalar *)constants));
3831:   PetscFunctionReturn(PETSC_SUCCESS);
3832: }

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

3837:   Not Collective

3839:   Input Parameter:
3840: . ds - The `PetscDS` object

3842:   Output Parameter:
3843: . newds - The `PetscDS` copy

3845:   Level: intermediate

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

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

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

3871:   Not Collective

3873:   Input Parameter:
3874: . ds - The `PetscDS` object

3876:   Output Parameter:
3877: . newds - The `PetscDS` copy

3879:   Level: intermediate

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

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

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

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

3937: PetscErrorCode PetscDSGetHeightSubspace(PetscDS prob, PetscInt height, PetscDS *subprob)
3938: {
3939:   PetscInt dim, Nf, f;

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

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

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

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

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

3997: PetscErrorCode PetscDSGetDiscType_Internal(PetscDS ds, PetscInt f, PetscDiscType *disctype)
3998: {
3999:   PetscObject  obj;
4000:   PetscClassId id;
4001:   PetscInt     Nf;

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

4018: static PetscErrorCode PetscDSDestroy_Basic(PetscDS ds)
4019: {
4020:   PetscFunctionBegin;
4021:   PetscCall(PetscFree(ds->data));
4022:   PetscFunctionReturn(PETSC_SUCCESS);
4023: }

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

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

4038:   Level: intermediate

4040: .seealso: `PetscDSType`, `PetscDSCreate()`, `PetscDSSetType()`
4041: M*/

4043: PETSC_EXTERN PetscErrorCode PetscDSCreate_Basic(PetscDS ds)
4044: {
4045:   PetscDS_Basic *b;

4047:   PetscFunctionBegin;
4049:   PetscCall(PetscNew(&b));
4050:   ds->data = b;

4052:   PetscCall(PetscDSInitialize_Basic(ds));
4053:   PetscFunctionReturn(PETSC_SUCCESS);
4054: }