xref: /petsc/src/dm/dt/interface/dtds.c (revision c3b34e8a033ac2eae760bfa25d5dd353c7b9f646)
1 #include <petsc/private/petscdsimpl.h> /*I "petscds.h" I*/
2 
3 PetscClassId PETSCDS_CLASSID = 0;
4 
5 PetscFunctionList PetscDSList              = NULL;
6 PetscBool         PetscDSRegisterAllCalled = PETSC_FALSE;
7 
8 /* A PetscDS (Discrete System) encodes a set of equations posed in a discrete space, which represents a set of
9    nonlinear continuum equations. The equations can have multiple fields, each field having a different
10    discretization. In addition, different pieces of the domain can have different field combinations and equations.
11 
12    The DS provides the user a description of the approximation space on any given cell. It also gives pointwise
13    functions representing the equations.
14 
15    Each field is associated with a label, marking the cells on which it is supported. Note that a field can be
16    supported on the closure of a cell not in the label due to overlap of the boundary of neighboring cells. The DM
17    then creates a DS for each set of cells with identical approximation spaces. When assembling, the user asks for
18    the space associated with a given cell. DMPlex uses the labels associated with each DS in the default integration loop.
19 */
20 
21 /*@C
22   PetscDSRegister - Adds a new `PetscDS` implementation
23 
24   Not Collective; No Fortran Support
25 
26   Input Parameters:
27 + sname    - The name of a new user-defined creation routine
28 - function - The creation routine itself
29 
30   Example Usage:
31 .vb
32     PetscDSRegister("my_ds", MyPetscDSCreate);
33 .ve
34 
35   Then, your PetscDS type can be chosen with the procedural interface via
36 .vb
37     PetscDSCreate(MPI_Comm, PetscDS *);
38     PetscDSSetType(PetscDS, "my_ds");
39 .ve
40   or at runtime via the option
41 .vb
42     -petscds_type my_ds
43 .ve
44 
45   Level: advanced
46 
47   Note:
48   `PetscDSRegister()` may be called multiple times to add several user-defined `PetscDSs`
49 
50 .seealso: `PetscDSType`, `PetscDS`, `PetscDSRegisterAll()`, `PetscDSRegisterDestroy()`
51 @*/
52 PetscErrorCode PetscDSRegister(const char sname[], PetscErrorCode (*function)(PetscDS))
53 {
54   PetscFunctionBegin;
55   PetscCall(PetscFunctionListAdd(&PetscDSList, sname, function));
56   PetscFunctionReturn(PETSC_SUCCESS);
57 }
58 
59 /*@C
60   PetscDSSetType - Builds a particular `PetscDS`
61 
62   Collective; No Fortran Support
63 
64   Input Parameters:
65 + prob - The `PetscDS` object
66 - name - The `PetscDSType`
67 
68   Options Database Key:
69 . -petscds_type <type> - Sets the PetscDS type; use -help for a list of available types
70 
71   Level: intermediate
72 
73 .seealso: `PetscDSType`, `PetscDS`, `PetscDSGetType()`, `PetscDSCreate()`
74 @*/
75 PetscErrorCode PetscDSSetType(PetscDS prob, PetscDSType name)
76 {
77   PetscErrorCode (*r)(PetscDS);
78   PetscBool match;
79 
80   PetscFunctionBegin;
81   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
82   PetscCall(PetscObjectTypeCompare((PetscObject)prob, name, &match));
83   if (match) PetscFunctionReturn(PETSC_SUCCESS);
84 
85   PetscCall(PetscDSRegisterAll());
86   PetscCall(PetscFunctionListFind(PetscDSList, name, &r));
87   PetscCheck(r, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscDS type: %s", name);
88 
89   PetscTryTypeMethod(prob, destroy);
90   prob->ops->destroy = NULL;
91 
92   PetscCall((*r)(prob));
93   PetscCall(PetscObjectChangeTypeName((PetscObject)prob, name));
94   PetscFunctionReturn(PETSC_SUCCESS);
95 }
96 
97 /*@C
98   PetscDSGetType - Gets the `PetscDSType` name (as a string) from the `PetscDS`
99 
100   Not Collective; No Fortran Support
101 
102   Input Parameter:
103 . prob - The `PetscDS`
104 
105   Output Parameter:
106 . name - The `PetscDSType` name
107 
108   Level: intermediate
109 
110 .seealso: `PetscDSType`, `PetscDS`, `PetscDSSetType()`, `PetscDSCreate()`
111 @*/
112 PetscErrorCode PetscDSGetType(PetscDS prob, PetscDSType *name)
113 {
114   PetscFunctionBegin;
115   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
116   PetscAssertPointer(name, 2);
117   PetscCall(PetscDSRegisterAll());
118   *name = ((PetscObject)prob)->type_name;
119   PetscFunctionReturn(PETSC_SUCCESS);
120 }
121 
122 static PetscErrorCode PetscDSView_Ascii(PetscDS ds, PetscViewer viewer)
123 {
124   PetscViewerFormat  format;
125   const PetscScalar *constants;
126   PetscInt           Nf, numConstants, f;
127 
128   PetscFunctionBegin;
129   PetscCall(PetscDSGetNumFields(ds, &Nf));
130   PetscCall(PetscViewerGetFormat(viewer, &format));
131   PetscCall(PetscViewerASCIIPrintf(viewer, "Discrete System with %" PetscInt_FMT " fields\n", Nf));
132   PetscCall(PetscViewerASCIIPushTab(viewer));
133   PetscCall(PetscViewerASCIIPrintf(viewer, "  cell total dim %" PetscInt_FMT " total comp %" PetscInt_FMT "\n", ds->totDim, ds->totComp));
134   if (ds->isCohesive) PetscCall(PetscViewerASCIIPrintf(viewer, "  cohesive cell\n"));
135   for (f = 0; f < Nf; ++f) {
136     DSBoundary      b;
137     PetscObject     obj;
138     PetscClassId    id;
139     PetscQuadrature q;
140     const char     *name;
141     PetscInt        Nc, Nq, Nqc;
142 
143     PetscCall(PetscDSGetDiscretization(ds, f, &obj));
144     PetscCall(PetscObjectGetClassId(obj, &id));
145     PetscCall(PetscObjectGetName(obj, &name));
146     PetscCall(PetscViewerASCIIPrintf(viewer, "Field %s", name ? name : "<unknown>"));
147     PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
148     if (id == PETSCFE_CLASSID) {
149       PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
150       PetscCall(PetscFEGetQuadrature((PetscFE)obj, &q));
151       PetscCall(PetscViewerASCIIPrintf(viewer, " FEM"));
152     } else if (id == PETSCFV_CLASSID) {
153       PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc));
154       PetscCall(PetscFVGetQuadrature((PetscFV)obj, &q));
155       PetscCall(PetscViewerASCIIPrintf(viewer, " FVM"));
156     } else SETERRQ(PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %" PetscInt_FMT, f);
157     if (Nc > 1) PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT " components", Nc));
158     else PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT " component ", Nc));
159     if (ds->implicit[f]) PetscCall(PetscViewerASCIIPrintf(viewer, " (implicit)"));
160     else PetscCall(PetscViewerASCIIPrintf(viewer, " (explicit)"));
161     if (q) {
162       PetscCall(PetscQuadratureGetData(q, NULL, &Nqc, &Nq, NULL, NULL));
163       PetscCall(PetscViewerASCIIPrintf(viewer, " (Nq %" PetscInt_FMT " Nqc %" PetscInt_FMT ")", Nq, Nqc));
164     }
165     PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT "-jet", ds->jetDegree[f]));
166     PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
167     PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
168     PetscCall(PetscViewerASCIIPushTab(viewer));
169     if (id == PETSCFE_CLASSID) PetscCall(PetscFEView((PetscFE)obj, viewer));
170     else if (id == PETSCFV_CLASSID) PetscCall(PetscFVView((PetscFV)obj, viewer));
171     PetscCall(PetscViewerASCIIPopTab(viewer));
172 
173     for (b = ds->boundary; b; b = b->next) {
174       char    *name;
175       PetscInt c, i;
176 
177       if (b->field != f) continue;
178       PetscCall(PetscViewerASCIIPushTab(viewer));
179       PetscCall(PetscViewerASCIIPrintf(viewer, "Boundary %s (%s) %s\n", b->name, b->lname, DMBoundaryConditionTypes[b->type]));
180       if (!b->Nc) {
181         PetscCall(PetscViewerASCIIPrintf(viewer, "  all components\n"));
182       } else {
183         PetscCall(PetscViewerASCIIPrintf(viewer, "  components: "));
184         PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
185         for (c = 0; c < b->Nc; ++c) {
186           if (c > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ", "));
187           PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT, b->comps[c]));
188         }
189         PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
190         PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
191       }
192       PetscCall(PetscViewerASCIIPrintf(viewer, "  values: "));
193       PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
194       for (i = 0; i < b->Nv; ++i) {
195         if (i > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ", "));
196         PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT, b->values[i]));
197       }
198       PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
199       PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
200 #if defined(__clang__)
201       PETSC_PRAGMA_DIAGNOSTIC_IGNORED_BEGIN("-Wformat-pedantic")
202 #elif defined(__GNUC__) || defined(__GNUG__)
203       PETSC_PRAGMA_DIAGNOSTIC_IGNORED_BEGIN("-Wformat")
204 #endif
205       if (b->func) {
206         PetscCall(PetscDLAddr(b->func, &name));
207         if (name) PetscCall(PetscViewerASCIIPrintf(viewer, "  func: %s\n", name));
208         else PetscCall(PetscViewerASCIIPrintf(viewer, "  func: %p\n", b->func));
209         PetscCall(PetscFree(name));
210       }
211       if (b->func_t) {
212         PetscCall(PetscDLAddr(b->func_t, &name));
213         if (name) PetscCall(PetscViewerASCIIPrintf(viewer, "  func_t: %s\n", name));
214         else PetscCall(PetscViewerASCIIPrintf(viewer, "  func_t: %p\n", b->func_t));
215         PetscCall(PetscFree(name));
216       }
217       PETSC_PRAGMA_DIAGNOSTIC_IGNORED_END()
218       PetscCall(PetscWeakFormView(b->wf, viewer));
219       PetscCall(PetscViewerASCIIPopTab(viewer));
220     }
221   }
222   PetscCall(PetscDSGetConstants(ds, &numConstants, &constants));
223   if (numConstants) {
224     PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " constants\n", numConstants));
225     PetscCall(PetscViewerASCIIPushTab(viewer));
226     for (f = 0; f < numConstants; ++f) PetscCall(PetscViewerASCIIPrintf(viewer, "%g\n", (double)PetscRealPart(constants[f])));
227     PetscCall(PetscViewerASCIIPopTab(viewer));
228   }
229   PetscCall(PetscWeakFormView(ds->wf, viewer));
230   PetscCall(PetscViewerASCIIPopTab(viewer));
231   PetscFunctionReturn(PETSC_SUCCESS);
232 }
233 
234 /*@C
235   PetscDSViewFromOptions - View a `PetscDS` based on values in the options database
236 
237   Collective
238 
239   Input Parameters:
240 + A    - the `PetscDS` object
241 . obj  - Optional object that provides the options prefix used in the search
242 - name - command line option
243 
244   Level: intermediate
245 
246 .seealso: `PetscDSType`, `PetscDS`, `PetscDSView()`, `PetscObjectViewFromOptions()`, `PetscDSCreate()`
247 @*/
248 PetscErrorCode PetscDSViewFromOptions(PetscDS A, PetscObject obj, const char name[])
249 {
250   PetscFunctionBegin;
251   PetscValidHeaderSpecific(A, PETSCDS_CLASSID, 1);
252   PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
253   PetscFunctionReturn(PETSC_SUCCESS);
254 }
255 
256 /*@C
257   PetscDSView - Views a `PetscDS`
258 
259   Collective
260 
261   Input Parameters:
262 + prob - the `PetscDS` object to view
263 - v    - the viewer
264 
265   Level: developer
266 
267 .seealso: `PetscDSType`, `PetscDS`, `PetscViewer`, `PetscDSDestroy()`, `PetscDSViewFromOptions()`
268 @*/
269 PetscErrorCode PetscDSView(PetscDS prob, PetscViewer v)
270 {
271   PetscBool iascii;
272 
273   PetscFunctionBegin;
274   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
275   if (!v) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)prob), &v));
276   else PetscValidHeaderSpecific(v, PETSC_VIEWER_CLASSID, 2);
277   PetscCall(PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &iascii));
278   if (iascii) PetscCall(PetscDSView_Ascii(prob, v));
279   PetscTryTypeMethod(prob, view, v);
280   PetscFunctionReturn(PETSC_SUCCESS);
281 }
282 
283 /*@
284   PetscDSSetFromOptions - sets parameters in a `PetscDS` from the options database
285 
286   Collective
287 
288   Input Parameter:
289 . prob - the `PetscDS` object to set options for
290 
291   Options Database Keys:
292 + -petscds_type <type>     - Set the `PetscDS` type
293 . -petscds_view <view opt> - View the `PetscDS`
294 . -petscds_jac_pre         - Turn formation of a separate Jacobian preconditioner on or off
295 . -bc_<name> <ids>         - Specify a list of label ids for a boundary condition
296 - -bc_<name>_comp <comps>  - Specify a list of field components to constrain for a boundary condition
297 
298   Level: intermediate
299 
300 .seealso: `PetscDS`, `PetscDSView()`
301 @*/
302 PetscErrorCode PetscDSSetFromOptions(PetscDS prob)
303 {
304   DSBoundary  b;
305   const char *defaultType;
306   char        name[256];
307   PetscBool   flg;
308 
309   PetscFunctionBegin;
310   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
311   if (!((PetscObject)prob)->type_name) {
312     defaultType = PETSCDSBASIC;
313   } else {
314     defaultType = ((PetscObject)prob)->type_name;
315   }
316   PetscCall(PetscDSRegisterAll());
317 
318   PetscObjectOptionsBegin((PetscObject)prob);
319   for (b = prob->boundary; b; b = b->next) {
320     char      optname[1024];
321     PetscInt  ids[1024], len = 1024;
322     PetscBool flg;
323 
324     PetscCall(PetscSNPrintf(optname, sizeof(optname), "-bc_%s", b->name));
325     PetscCall(PetscMemzero(ids, sizeof(ids)));
326     PetscCall(PetscOptionsIntArray(optname, "List of boundary IDs", "", ids, &len, &flg));
327     if (flg) {
328       b->Nv = len;
329       PetscCall(PetscFree(b->values));
330       PetscCall(PetscMalloc1(len, &b->values));
331       PetscCall(PetscArraycpy(b->values, ids, len));
332       PetscCall(PetscWeakFormRewriteKeys(b->wf, b->label, len, b->values));
333     }
334     len = 1024;
335     PetscCall(PetscSNPrintf(optname, sizeof(optname), "-bc_%s_comp", b->name));
336     PetscCall(PetscMemzero(ids, sizeof(ids)));
337     PetscCall(PetscOptionsIntArray(optname, "List of boundary field components", "", ids, &len, &flg));
338     if (flg) {
339       b->Nc = len;
340       PetscCall(PetscFree(b->comps));
341       PetscCall(PetscMalloc1(len, &b->comps));
342       PetscCall(PetscArraycpy(b->comps, ids, len));
343     }
344   }
345   PetscCall(PetscOptionsFList("-petscds_type", "Discrete System", "PetscDSSetType", PetscDSList, defaultType, name, 256, &flg));
346   if (flg) {
347     PetscCall(PetscDSSetType(prob, name));
348   } else if (!((PetscObject)prob)->type_name) {
349     PetscCall(PetscDSSetType(prob, defaultType));
350   }
351   PetscCall(PetscOptionsBool("-petscds_jac_pre", "Discrete System", "PetscDSUseJacobianPreconditioner", prob->useJacPre, &prob->useJacPre, &flg));
352   PetscCall(PetscOptionsBool("-petscds_force_quad", "Discrete System", "PetscDSSetForceQuad", prob->forceQuad, &prob->forceQuad, &flg));
353   PetscTryTypeMethod(prob, setfromoptions);
354   /* process any options handlers added with PetscObjectAddOptionsHandler() */
355   PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)prob, PetscOptionsObject));
356   PetscOptionsEnd();
357   if (prob->Nf) PetscCall(PetscDSViewFromOptions(prob, NULL, "-petscds_view"));
358   PetscFunctionReturn(PETSC_SUCCESS);
359 }
360 
361 /*@C
362   PetscDSSetUp - Construct data structures for the `PetscDS`
363 
364   Collective
365 
366   Input Parameter:
367 . prob - the `PetscDS` object to setup
368 
369   Level: developer
370 
371 .seealso: `PetscDS`, `PetscDSView()`, `PetscDSDestroy()`
372 @*/
373 PetscErrorCode PetscDSSetUp(PetscDS prob)
374 {
375   const PetscInt Nf          = prob->Nf;
376   PetscBool      hasH        = PETSC_FALSE;
377   PetscInt       maxOrder[4] = {-1, -1, -1, -1};
378   PetscInt       dim, dimEmbed, NbMax = 0, NcMax = 0, NqMax = 0, NsMax = 1, f;
379 
380   PetscFunctionBegin;
381   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
382   if (prob->setup) PetscFunctionReturn(PETSC_SUCCESS);
383   /* Calculate sizes */
384   PetscCall(PetscDSGetSpatialDimension(prob, &dim));
385   PetscCall(PetscDSGetCoordinateDimension(prob, &dimEmbed));
386   prob->totDim = prob->totComp = 0;
387   PetscCall(PetscMalloc2(Nf, &prob->Nc, Nf, &prob->Nb));
388   PetscCall(PetscCalloc2(Nf + 1, &prob->off, Nf + 1, &prob->offDer));
389   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]));
390   PetscCall(PetscMalloc2(Nf, &prob->T, Nf, &prob->Tf));
391   if (prob->forceQuad) {
392     // Note: This assumes we have one kind of cell at each dimension.
393     //       We can fix this by having quadrature hold the celltype
394     PetscQuadrature maxQuad[4] = {NULL, NULL, NULL, NULL};
395 
396     for (f = 0; f < Nf; ++f) {
397       PetscObject     obj;
398       PetscClassId    id;
399       PetscQuadrature q = NULL, fq = NULL;
400       PetscInt        dim = -1, order = -1, forder = -1;
401 
402       PetscCall(PetscDSGetDiscretization(prob, f, &obj));
403       if (!obj) continue;
404       PetscCall(PetscObjectGetClassId(obj, &id));
405       if (id == PETSCFE_CLASSID) {
406         PetscFE fe = (PetscFE)obj;
407 
408         PetscCall(PetscFEGetQuadrature(fe, &q));
409         PetscCall(PetscFEGetFaceQuadrature(fe, &fq));
410       } else if (id == PETSCFV_CLASSID) {
411         PetscFV fv = (PetscFV)obj;
412 
413         PetscCall(PetscFVGetQuadrature(fv, &q));
414       }
415       if (q) {
416         PetscCall(PetscQuadratureGetData(q, &dim, NULL, NULL, NULL, NULL));
417         PetscCall(PetscQuadratureGetOrder(q, &order));
418         if (order > maxOrder[dim]) {
419           maxOrder[dim] = order;
420           maxQuad[dim]  = q;
421         }
422       }
423       if (fq) {
424         PetscCall(PetscQuadratureGetData(fq, &dim, NULL, NULL, NULL, NULL));
425         PetscCall(PetscQuadratureGetOrder(fq, &forder));
426         if (forder > maxOrder[dim]) {
427           maxOrder[dim] = forder;
428           maxQuad[dim]  = fq;
429         }
430       }
431     }
432     for (f = 0; f < Nf; ++f) {
433       PetscObject     obj;
434       PetscClassId    id;
435       PetscQuadrature q;
436       PetscInt        dim;
437 
438       PetscCall(PetscDSGetDiscretization(prob, f, &obj));
439       if (!obj) continue;
440       PetscCall(PetscObjectGetClassId(obj, &id));
441       if (id == PETSCFE_CLASSID) {
442         PetscFE fe = (PetscFE)obj;
443 
444         PetscCall(PetscFEGetQuadrature(fe, &q));
445         PetscCall(PetscQuadratureGetData(q, &dim, NULL, NULL, NULL, NULL));
446         PetscCall(PetscFESetQuadrature(fe, maxQuad[dim]));
447         PetscCall(PetscFESetFaceQuadrature(fe, dim ? maxQuad[dim - 1] : NULL));
448       } else if (id == PETSCFV_CLASSID) {
449         PetscFV fv = (PetscFV)obj;
450 
451         PetscCall(PetscFVGetQuadrature(fv, &q));
452         PetscCall(PetscQuadratureGetData(q, &dim, NULL, NULL, NULL, NULL));
453         PetscCall(PetscFVSetQuadrature(fv, maxQuad[dim]));
454       }
455     }
456   }
457   for (f = 0; f < Nf; ++f) {
458     PetscObject     obj;
459     PetscClassId    id;
460     PetscQuadrature q  = NULL;
461     PetscInt        Nq = 0, Nb, Nc;
462 
463     PetscCall(PetscDSGetDiscretization(prob, f, &obj));
464     if (prob->jetDegree[f] > 1) hasH = PETSC_TRUE;
465     if (!obj) {
466       /* Empty mesh */
467       Nb = Nc    = 0;
468       prob->T[f] = prob->Tf[f] = NULL;
469     } else {
470       PetscCall(PetscObjectGetClassId(obj, &id));
471       if (id == PETSCFE_CLASSID) {
472         PetscFE fe = (PetscFE)obj;
473 
474         PetscCall(PetscFEGetQuadrature(fe, &q));
475         {
476           PetscQuadrature fq;
477           PetscInt        dim, order;
478 
479           PetscCall(PetscQuadratureGetData(q, &dim, NULL, NULL, NULL, NULL));
480           PetscCall(PetscQuadratureGetOrder(q, &order));
481           if (maxOrder[dim] < 0) maxOrder[dim] = order;
482           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]);
483           PetscCall(PetscFEGetFaceQuadrature(fe, &fq));
484           if (fq) {
485             PetscCall(PetscQuadratureGetData(fq, &dim, NULL, NULL, NULL, NULL));
486             PetscCall(PetscQuadratureGetOrder(fq, &order));
487             if (maxOrder[dim] < 0) maxOrder[dim] = order;
488             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]);
489           }
490         }
491         PetscCall(PetscFEGetDimension(fe, &Nb));
492         PetscCall(PetscFEGetNumComponents(fe, &Nc));
493         PetscCall(PetscFEGetCellTabulation(fe, prob->jetDegree[f], &prob->T[f]));
494         PetscCall(PetscFEGetFaceTabulation(fe, prob->jetDegree[f], &prob->Tf[f]));
495       } else if (id == PETSCFV_CLASSID) {
496         PetscFV fv = (PetscFV)obj;
497 
498         PetscCall(PetscFVGetQuadrature(fv, &q));
499         PetscCall(PetscFVGetNumComponents(fv, &Nc));
500         Nb = Nc;
501         PetscCall(PetscFVGetCellTabulation(fv, &prob->T[f]));
502         /* TODO: should PetscFV also have face tabulation? Otherwise there will be a null pointer in prob->basisFace */
503       } else SETERRQ(PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %" PetscInt_FMT, f);
504     }
505     prob->Nc[f]                    = Nc;
506     prob->Nb[f]                    = Nb;
507     prob->off[f + 1]               = Nc + prob->off[f];
508     prob->offDer[f + 1]            = Nc * dim + prob->offDer[f];
509     prob->offCohesive[0][f + 1]    = (prob->cohesive[f] ? Nc : Nc * 2) + prob->offCohesive[0][f];
510     prob->offDerCohesive[0][f + 1] = (prob->cohesive[f] ? Nc : Nc * 2) * dimEmbed + prob->offDerCohesive[0][f];
511     prob->offCohesive[1][f]        = (prob->cohesive[f] ? 0 : Nc) + prob->offCohesive[0][f];
512     prob->offDerCohesive[1][f]     = (prob->cohesive[f] ? 0 : Nc) * dimEmbed + prob->offDerCohesive[0][f];
513     prob->offCohesive[2][f + 1]    = (prob->cohesive[f] ? Nc : Nc * 2) + prob->offCohesive[2][f];
514     prob->offDerCohesive[2][f + 1] = (prob->cohesive[f] ? Nc : Nc * 2) * dimEmbed + prob->offDerCohesive[2][f];
515     if (q) PetscCall(PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL));
516     NqMax = PetscMax(NqMax, Nq);
517     NbMax = PetscMax(NbMax, Nb);
518     NcMax = PetscMax(NcMax, Nc);
519     prob->totDim += Nb;
520     prob->totComp += Nc;
521     /* There are two faces for all fields on a cohesive cell, except for cohesive fields */
522     if (prob->isCohesive && !prob->cohesive[f]) prob->totDim += Nb;
523   }
524   prob->offCohesive[1][Nf]    = prob->offCohesive[0][Nf];
525   prob->offDerCohesive[1][Nf] = prob->offDerCohesive[0][Nf];
526   /* Allocate works space */
527   NsMax = 2; /* A non-cohesive discretizations can be used on a cohesive cell, so we need this extra workspace for all DS */
528   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));
529   PetscCall(PetscMalloc5(dimEmbed, &prob->x, NbMax * NcMax, &prob->basisReal, NbMax * NcMax * dimEmbed, &prob->basisDerReal, NbMax * NcMax, &prob->testReal, NbMax * NcMax * dimEmbed, &prob->testDerReal));
530   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,
531                          &prob->g2, NsMax * NsMax * NqMax * NcMax * NcMax * dimEmbed * dimEmbed, &prob->g3));
532   PetscTryTypeMethod(prob, setup);
533   prob->setup = PETSC_TRUE;
534   PetscFunctionReturn(PETSC_SUCCESS);
535 }
536 
537 static PetscErrorCode PetscDSDestroyStructs_Static(PetscDS prob)
538 {
539   PetscFunctionBegin;
540   PetscCall(PetscFree2(prob->Nc, prob->Nb));
541   PetscCall(PetscFree2(prob->off, prob->offDer));
542   PetscCall(PetscFree6(prob->offCohesive[0], prob->offCohesive[1], prob->offCohesive[2], prob->offDerCohesive[0], prob->offDerCohesive[1], prob->offDerCohesive[2]));
543   PetscCall(PetscFree2(prob->T, prob->Tf));
544   PetscCall(PetscFree3(prob->u, prob->u_t, prob->u_x));
545   PetscCall(PetscFree5(prob->x, prob->basisReal, prob->basisDerReal, prob->testReal, prob->testDerReal));
546   PetscCall(PetscFree6(prob->f0, prob->f1, prob->g0, prob->g1, prob->g2, prob->g3));
547   PetscFunctionReturn(PETSC_SUCCESS);
548 }
549 
550 static PetscErrorCode PetscDSEnlarge_Static(PetscDS prob, PetscInt NfNew)
551 {
552   PetscObject          *tmpd;
553   PetscBool            *tmpi;
554   PetscInt             *tmpk;
555   PetscBool            *tmpc;
556   PetscPointFunc       *tmpup;
557   PetscSimplePointFunc *tmpexactSol, *tmpexactSol_t;
558   void                **tmpexactCtx, **tmpexactCtx_t;
559   void                **tmpctx;
560   PetscInt              Nf = prob->Nf, f;
561 
562   PetscFunctionBegin;
563   if (Nf >= NfNew) PetscFunctionReturn(PETSC_SUCCESS);
564   prob->setup = PETSC_FALSE;
565   PetscCall(PetscDSDestroyStructs_Static(prob));
566   PetscCall(PetscMalloc4(NfNew, &tmpd, NfNew, &tmpi, NfNew, &tmpc, NfNew, &tmpk));
567   for (f = 0; f < Nf; ++f) {
568     tmpd[f] = prob->disc[f];
569     tmpi[f] = prob->implicit[f];
570     tmpc[f] = prob->cohesive[f];
571     tmpk[f] = prob->jetDegree[f];
572   }
573   for (f = Nf; f < NfNew; ++f) {
574     tmpd[f] = NULL;
575     tmpi[f] = PETSC_TRUE, tmpc[f] = PETSC_FALSE;
576     tmpk[f] = 1;
577   }
578   PetscCall(PetscFree4(prob->disc, prob->implicit, prob->cohesive, prob->jetDegree));
579   PetscCall(PetscWeakFormSetNumFields(prob->wf, NfNew));
580   prob->Nf        = NfNew;
581   prob->disc      = tmpd;
582   prob->implicit  = tmpi;
583   prob->cohesive  = tmpc;
584   prob->jetDegree = tmpk;
585   PetscCall(PetscCalloc2(NfNew, &tmpup, NfNew, &tmpctx));
586   for (f = 0; f < Nf; ++f) tmpup[f] = prob->update[f];
587   for (f = 0; f < Nf; ++f) tmpctx[f] = prob->ctx[f];
588   for (f = Nf; f < NfNew; ++f) tmpup[f] = NULL;
589   for (f = Nf; f < NfNew; ++f) tmpctx[f] = NULL;
590   PetscCall(PetscFree2(prob->update, prob->ctx));
591   prob->update = tmpup;
592   prob->ctx    = tmpctx;
593   PetscCall(PetscCalloc4(NfNew, &tmpexactSol, NfNew, &tmpexactCtx, NfNew, &tmpexactSol_t, NfNew, &tmpexactCtx_t));
594   for (f = 0; f < Nf; ++f) tmpexactSol[f] = prob->exactSol[f];
595   for (f = 0; f < Nf; ++f) tmpexactCtx[f] = prob->exactCtx[f];
596   for (f = 0; f < Nf; ++f) tmpexactSol_t[f] = prob->exactSol_t[f];
597   for (f = 0; f < Nf; ++f) tmpexactCtx_t[f] = prob->exactCtx_t[f];
598   for (f = Nf; f < NfNew; ++f) tmpexactSol[f] = NULL;
599   for (f = Nf; f < NfNew; ++f) tmpexactCtx[f] = NULL;
600   for (f = Nf; f < NfNew; ++f) tmpexactSol_t[f] = NULL;
601   for (f = Nf; f < NfNew; ++f) tmpexactCtx_t[f] = NULL;
602   PetscCall(PetscFree4(prob->exactSol, prob->exactCtx, prob->exactSol_t, prob->exactCtx_t));
603   prob->exactSol   = tmpexactSol;
604   prob->exactCtx   = tmpexactCtx;
605   prob->exactSol_t = tmpexactSol_t;
606   prob->exactCtx_t = tmpexactCtx_t;
607   PetscFunctionReturn(PETSC_SUCCESS);
608 }
609 
610 /*@
611   PetscDSDestroy - Destroys a `PetscDS` object
612 
613   Collective
614 
615   Input Parameter:
616 . ds - the `PetscDS` object to destroy
617 
618   Level: developer
619 
620 .seealso: `PetscDSView()`
621 @*/
622 PetscErrorCode PetscDSDestroy(PetscDS *ds)
623 {
624   PetscInt f;
625 
626   PetscFunctionBegin;
627   if (!*ds) PetscFunctionReturn(PETSC_SUCCESS);
628   PetscValidHeaderSpecific((*ds), PETSCDS_CLASSID, 1);
629 
630   if (--((PetscObject)(*ds))->refct > 0) {
631     *ds = NULL;
632     PetscFunctionReturn(PETSC_SUCCESS);
633   }
634   ((PetscObject)(*ds))->refct = 0;
635   if ((*ds)->subprobs) {
636     PetscInt dim, d;
637 
638     PetscCall(PetscDSGetSpatialDimension(*ds, &dim));
639     for (d = 0; d < dim; ++d) PetscCall(PetscDSDestroy(&(*ds)->subprobs[d]));
640   }
641   PetscCall(PetscFree((*ds)->subprobs));
642   PetscCall(PetscDSDestroyStructs_Static(*ds));
643   for (f = 0; f < (*ds)->Nf; ++f) PetscCall(PetscObjectDereference((*ds)->disc[f]));
644   PetscCall(PetscFree4((*ds)->disc, (*ds)->implicit, (*ds)->cohesive, (*ds)->jetDegree));
645   PetscCall(PetscWeakFormDestroy(&(*ds)->wf));
646   PetscCall(PetscFree2((*ds)->update, (*ds)->ctx));
647   PetscCall(PetscFree4((*ds)->exactSol, (*ds)->exactCtx, (*ds)->exactSol_t, (*ds)->exactCtx_t));
648   PetscTryTypeMethod((*ds), destroy);
649   PetscCall(PetscDSDestroyBoundary(*ds));
650   PetscCall(PetscFree((*ds)->constants));
651   for (PetscInt c = 0; c < DM_NUM_POLYTOPES; ++c) {
652     const PetscInt Na = DMPolytopeTypeGetNumArrangments((DMPolytopeType)c);
653     if ((*ds)->quadPerm[c])
654       for (PetscInt o = 0; o < Na; ++o) PetscCall(ISDestroy(&(*ds)->quadPerm[c][o]));
655     PetscCall(PetscFree((*ds)->quadPerm[c]));
656     (*ds)->quadPerm[c] = NULL;
657   }
658   PetscCall(PetscHeaderDestroy(ds));
659   PetscFunctionReturn(PETSC_SUCCESS);
660 }
661 
662 /*@
663   PetscDSCreate - Creates an empty `PetscDS` object. The type can then be set with `PetscDSSetType()`.
664 
665   Collective
666 
667   Input Parameter:
668 . comm - The communicator for the `PetscDS` object
669 
670   Output Parameter:
671 . ds - The `PetscDS` object
672 
673   Level: beginner
674 
675 .seealso: `PetscDS`, `PetscDSSetType()`, `PETSCDSBASIC`, `PetscDSType`
676 @*/
677 PetscErrorCode PetscDSCreate(MPI_Comm comm, PetscDS *ds)
678 {
679   PetscDS p;
680 
681   PetscFunctionBegin;
682   PetscAssertPointer(ds, 2);
683   *ds = NULL;
684   PetscCall(PetscDSInitializePackage());
685 
686   PetscCall(PetscHeaderCreate(p, PETSCDS_CLASSID, "PetscDS", "Discrete System", "PetscDS", comm, PetscDSDestroy, PetscDSView));
687 
688   p->Nf           = 0;
689   p->setup        = PETSC_FALSE;
690   p->numConstants = 0;
691   p->constants    = NULL;
692   p->dimEmbed     = -1;
693   p->useJacPre    = PETSC_TRUE;
694   p->forceQuad    = PETSC_TRUE;
695   PetscCall(PetscWeakFormCreate(comm, &p->wf));
696   PetscCall(PetscArrayzero(p->quadPerm, DM_NUM_POLYTOPES));
697 
698   *ds = p;
699   PetscFunctionReturn(PETSC_SUCCESS);
700 }
701 
702 /*@
703   PetscDSGetNumFields - Returns the number of fields in the `PetscDS`
704 
705   Not Collective
706 
707   Input Parameter:
708 . prob - The `PetscDS` object
709 
710   Output Parameter:
711 . Nf - The number of fields
712 
713   Level: beginner
714 
715 .seealso: `PetscDS`, `PetscDSGetSpatialDimension()`, `PetscDSCreate()`
716 @*/
717 PetscErrorCode PetscDSGetNumFields(PetscDS prob, PetscInt *Nf)
718 {
719   PetscFunctionBegin;
720   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
721   PetscAssertPointer(Nf, 2);
722   *Nf = prob->Nf;
723   PetscFunctionReturn(PETSC_SUCCESS);
724 }
725 
726 /*@
727   PetscDSGetSpatialDimension - Returns the spatial dimension of the `PetscDS`, meaning the topological dimension of the discretizations
728 
729   Not Collective
730 
731   Input Parameter:
732 . prob - The `PetscDS` object
733 
734   Output Parameter:
735 . dim - The spatial dimension
736 
737   Level: beginner
738 
739 .seealso: `PetscDS`, `PetscDSGetCoordinateDimension()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
740 @*/
741 PetscErrorCode PetscDSGetSpatialDimension(PetscDS prob, PetscInt *dim)
742 {
743   PetscFunctionBegin;
744   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
745   PetscAssertPointer(dim, 2);
746   *dim = 0;
747   if (prob->Nf) {
748     PetscObject  obj;
749     PetscClassId id;
750 
751     PetscCall(PetscDSGetDiscretization(prob, 0, &obj));
752     if (obj) {
753       PetscCall(PetscObjectGetClassId(obj, &id));
754       if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetSpatialDimension((PetscFE)obj, dim));
755       else if (id == PETSCFV_CLASSID) PetscCall(PetscFVGetSpatialDimension((PetscFV)obj, dim));
756       else SETERRQ(PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", 0);
757     }
758   }
759   PetscFunctionReturn(PETSC_SUCCESS);
760 }
761 
762 /*@
763   PetscDSGetCoordinateDimension - Returns the coordinate dimension of the `PetscDS`, meaning the dimension of the space into which the discretiaztions are embedded
764 
765   Not Collective
766 
767   Input Parameter:
768 . prob - The `PetscDS` object
769 
770   Output Parameter:
771 . dimEmbed - The coordinate dimension
772 
773   Level: beginner
774 
775 .seealso: `PetscDS`, `PetscDSSetCoordinateDimension()`, `PetscDSGetSpatialDimension()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
776 @*/
777 PetscErrorCode PetscDSGetCoordinateDimension(PetscDS prob, PetscInt *dimEmbed)
778 {
779   PetscFunctionBegin;
780   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
781   PetscAssertPointer(dimEmbed, 2);
782   PetscCheck(prob->dimEmbed >= 0, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONGSTATE, "No coordinate dimension set for this DS");
783   *dimEmbed = prob->dimEmbed;
784   PetscFunctionReturn(PETSC_SUCCESS);
785 }
786 
787 /*@
788   PetscDSSetCoordinateDimension - Set the coordinate dimension of the `PetscDS`, meaning the dimension of the space into which the discretiaztions are embedded
789 
790   Logically Collective
791 
792   Input Parameters:
793 + prob     - The `PetscDS` object
794 - dimEmbed - The coordinate dimension
795 
796   Level: beginner
797 
798 .seealso: `PetscDS`, `PetscDSGetCoordinateDimension()`, `PetscDSGetSpatialDimension()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
799 @*/
800 PetscErrorCode PetscDSSetCoordinateDimension(PetscDS prob, PetscInt dimEmbed)
801 {
802   PetscFunctionBegin;
803   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
804   PetscCheck(dimEmbed >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Coordinate dimension must be non-negative, not %" PetscInt_FMT, dimEmbed);
805   prob->dimEmbed = dimEmbed;
806   PetscFunctionReturn(PETSC_SUCCESS);
807 }
808 
809 /*@
810   PetscDSGetForceQuad - Returns the flag to force matching quadratures among the field discretizations
811 
812   Not collective
813 
814   Input Parameter:
815 . ds - The `PetscDS` object
816 
817   Output Parameter:
818 . forceQuad - The flag
819 
820   Level: intermediate
821 
822 .seealso: `PetscDS`, `PetscDSSetForceQuad()`, `PetscDSGetDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
823 @*/
824 PetscErrorCode PetscDSGetForceQuad(PetscDS ds, PetscBool *forceQuad)
825 {
826   PetscFunctionBegin;
827   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
828   PetscAssertPointer(forceQuad, 2);
829   *forceQuad = ds->forceQuad;
830   PetscFunctionReturn(PETSC_SUCCESS);
831 }
832 
833 /*@
834   PetscDSSetForceQuad - Set the flag to force matching quadratures among the field discretizations
835 
836   Logically collective on ds
837 
838   Input Parameters:
839 + ds        - The `PetscDS` object
840 - forceQuad - The flag
841 
842   Level: intermediate
843 
844 .seealso: `PetscDS`, `PetscDSGetForceQuad()`, `PetscDSGetDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
845 @*/
846 PetscErrorCode PetscDSSetForceQuad(PetscDS ds, PetscBool forceQuad)
847 {
848   PetscFunctionBegin;
849   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
850   ds->forceQuad = forceQuad;
851   PetscFunctionReturn(PETSC_SUCCESS);
852 }
853 
854 /*@
855   PetscDSIsCohesive - Returns the flag indicating that this `PetscDS` is for a cohesive cell
856 
857   Not Collective
858 
859   Input Parameter:
860 . ds - The `PetscDS` object
861 
862   Output Parameter:
863 . isCohesive - The flag
864 
865   Level: developer
866 
867 .seealso: `PetscDS`, `PetscDSGetNumCohesive()`, `PetscDSGetCohesive()`, `PetscDSSetCohesive()`, `PetscDSCreate()`
868 @*/
869 PetscErrorCode PetscDSIsCohesive(PetscDS ds, PetscBool *isCohesive)
870 {
871   PetscFunctionBegin;
872   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
873   PetscAssertPointer(isCohesive, 2);
874   *isCohesive = ds->isCohesive;
875   PetscFunctionReturn(PETSC_SUCCESS);
876 }
877 
878 /*@
879   PetscDSGetNumCohesive - Returns the number of cohesive fields, meaning those defined on the interior of a cohesive cell
880 
881   Not Collective
882 
883   Input Parameter:
884 . ds - The `PetscDS` object
885 
886   Output Parameter:
887 . numCohesive - The number of cohesive fields
888 
889   Level: developer
890 
891 .seealso: `PetscDS`, `PetscDSSetCohesive()`, `PetscDSCreate()`
892 @*/
893 PetscErrorCode PetscDSGetNumCohesive(PetscDS ds, PetscInt *numCohesive)
894 {
895   PetscInt f;
896 
897   PetscFunctionBegin;
898   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
899   PetscAssertPointer(numCohesive, 2);
900   *numCohesive = 0;
901   for (f = 0; f < ds->Nf; ++f) *numCohesive += ds->cohesive[f] ? 1 : 0;
902   PetscFunctionReturn(PETSC_SUCCESS);
903 }
904 
905 /*@
906   PetscDSGetCohesive - Returns the flag indicating that a field is cohesive, meaning it is defined on the interior of a cohesive cell
907 
908   Not Collective
909 
910   Input Parameters:
911 + ds - The `PetscDS` object
912 - f  - The field index
913 
914   Output Parameter:
915 . isCohesive - The flag
916 
917   Level: developer
918 
919 .seealso: `PetscDS`, `PetscDSSetCohesive()`, `PetscDSIsCohesive()`, `PetscDSCreate()`
920 @*/
921 PetscErrorCode PetscDSGetCohesive(PetscDS ds, PetscInt f, PetscBool *isCohesive)
922 {
923   PetscFunctionBegin;
924   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
925   PetscAssertPointer(isCohesive, 3);
926   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);
927   *isCohesive = ds->cohesive[f];
928   PetscFunctionReturn(PETSC_SUCCESS);
929 }
930 
931 /*@
932   PetscDSSetCohesive - Set the flag indicating that a field is cohesive, meaning it is defined on the interior of a cohesive cell
933 
934   Not Collective
935 
936   Input Parameters:
937 + ds         - The `PetscDS` object
938 . f          - The field index
939 - isCohesive - The flag for a cohesive field
940 
941   Level: developer
942 
943 .seealso: `PetscDS`, `PetscDSGetCohesive()`, `PetscDSIsCohesive()`, `PetscDSCreate()`
944 @*/
945 PetscErrorCode PetscDSSetCohesive(PetscDS ds, PetscInt f, PetscBool isCohesive)
946 {
947   PetscInt i;
948 
949   PetscFunctionBegin;
950   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
951   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);
952   ds->cohesive[f] = isCohesive;
953   ds->isCohesive  = PETSC_FALSE;
954   for (i = 0; i < ds->Nf; ++i) ds->isCohesive = ds->isCohesive || ds->cohesive[f] ? PETSC_TRUE : PETSC_FALSE;
955   PetscFunctionReturn(PETSC_SUCCESS);
956 }
957 
958 /*@
959   PetscDSGetTotalDimension - Returns the total size of the approximation space for this system
960 
961   Not Collective
962 
963   Input Parameter:
964 . prob - The `PetscDS` object
965 
966   Output Parameter:
967 . dim - The total problem dimension
968 
969   Level: beginner
970 
971 .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
972 @*/
973 PetscErrorCode PetscDSGetTotalDimension(PetscDS prob, PetscInt *dim)
974 {
975   PetscFunctionBegin;
976   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
977   PetscCall(PetscDSSetUp(prob));
978   PetscAssertPointer(dim, 2);
979   *dim = prob->totDim;
980   PetscFunctionReturn(PETSC_SUCCESS);
981 }
982 
983 /*@
984   PetscDSGetTotalComponents - Returns the total number of components in this system
985 
986   Not Collective
987 
988   Input Parameter:
989 . prob - The `PetscDS` object
990 
991   Output Parameter:
992 . Nc - The total number of components
993 
994   Level: beginner
995 
996 .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
997 @*/
998 PetscErrorCode PetscDSGetTotalComponents(PetscDS prob, PetscInt *Nc)
999 {
1000   PetscFunctionBegin;
1001   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
1002   PetscCall(PetscDSSetUp(prob));
1003   PetscAssertPointer(Nc, 2);
1004   *Nc = prob->totComp;
1005   PetscFunctionReturn(PETSC_SUCCESS);
1006 }
1007 
1008 /*@
1009   PetscDSGetDiscretization - Returns the discretization object for the given field
1010 
1011   Not Collective
1012 
1013   Input Parameters:
1014 + prob - The `PetscDS` object
1015 - f    - The field number
1016 
1017   Output Parameter:
1018 . disc - The discretization object
1019 
1020   Level: beginner
1021 
1022 .seealso: `PetscDS`, `PetscFE`, `PetscFV`, `PetscDSSetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1023 @*/
1024 PetscErrorCode PetscDSGetDiscretization(PetscDS prob, PetscInt f, PetscObject *disc)
1025 {
1026   PetscFunctionBeginHot;
1027   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
1028   PetscAssertPointer(disc, 3);
1029   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);
1030   *disc = prob->disc[f];
1031   PetscFunctionReturn(PETSC_SUCCESS);
1032 }
1033 
1034 /*@
1035   PetscDSSetDiscretization - Sets the discretization object for the given field
1036 
1037   Not Collective
1038 
1039   Input Parameters:
1040 + prob - The `PetscDS` object
1041 . f    - The field number
1042 - disc - The discretization object
1043 
1044   Level: beginner
1045 
1046 .seealso: `PetscDS`, `PetscFE`, `PetscFV`, `PetscDSGetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1047 @*/
1048 PetscErrorCode PetscDSSetDiscretization(PetscDS prob, PetscInt f, PetscObject disc)
1049 {
1050   PetscFunctionBegin;
1051   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
1052   if (disc) PetscAssertPointer(disc, 3);
1053   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
1054   PetscCall(PetscDSEnlarge_Static(prob, f + 1));
1055   PetscCall(PetscObjectDereference(prob->disc[f]));
1056   prob->disc[f] = disc;
1057   PetscCall(PetscObjectReference(disc));
1058   if (disc) {
1059     PetscClassId id;
1060 
1061     PetscCall(PetscObjectGetClassId(disc, &id));
1062     if (id == PETSCFE_CLASSID) {
1063       PetscCall(PetscDSSetImplicit(prob, f, PETSC_TRUE));
1064     } else if (id == PETSCFV_CLASSID) {
1065       PetscCall(PetscDSSetImplicit(prob, f, PETSC_FALSE));
1066     }
1067     PetscCall(PetscDSSetJetDegree(prob, f, 1));
1068   }
1069   PetscFunctionReturn(PETSC_SUCCESS);
1070 }
1071 
1072 /*@
1073   PetscDSGetWeakForm - Returns the weak form object
1074 
1075   Not Collective
1076 
1077   Input Parameter:
1078 . ds - The `PetscDS` object
1079 
1080   Output Parameter:
1081 . wf - The weak form object
1082 
1083   Level: beginner
1084 
1085 .seealso: `PetscWeakForm`, `PetscDSSetWeakForm()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1086 @*/
1087 PetscErrorCode PetscDSGetWeakForm(PetscDS ds, PetscWeakForm *wf)
1088 {
1089   PetscFunctionBegin;
1090   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
1091   PetscAssertPointer(wf, 2);
1092   *wf = ds->wf;
1093   PetscFunctionReturn(PETSC_SUCCESS);
1094 }
1095 
1096 /*@
1097   PetscDSSetWeakForm - Sets the weak form object
1098 
1099   Not Collective
1100 
1101   Input Parameters:
1102 + ds - The `PetscDS` object
1103 - wf - The weak form object
1104 
1105   Level: beginner
1106 
1107 .seealso: `PetscWeakForm`, `PetscDSGetWeakForm()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1108 @*/
1109 PetscErrorCode PetscDSSetWeakForm(PetscDS ds, PetscWeakForm wf)
1110 {
1111   PetscFunctionBegin;
1112   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
1113   PetscValidHeaderSpecific(wf, PETSCWEAKFORM_CLASSID, 2);
1114   PetscCall(PetscObjectDereference((PetscObject)ds->wf));
1115   ds->wf = wf;
1116   PetscCall(PetscObjectReference((PetscObject)wf));
1117   PetscCall(PetscWeakFormSetNumFields(wf, ds->Nf));
1118   PetscFunctionReturn(PETSC_SUCCESS);
1119 }
1120 
1121 /*@
1122   PetscDSAddDiscretization - Adds a discretization object
1123 
1124   Not Collective
1125 
1126   Input Parameters:
1127 + prob - The `PetscDS` object
1128 - disc - The boundary discretization object
1129 
1130   Level: beginner
1131 
1132 .seealso: `PetscWeakForm`, `PetscDSGetDiscretization()`, `PetscDSSetDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1133 @*/
1134 PetscErrorCode PetscDSAddDiscretization(PetscDS prob, PetscObject disc)
1135 {
1136   PetscFunctionBegin;
1137   PetscCall(PetscDSSetDiscretization(prob, prob->Nf, disc));
1138   PetscFunctionReturn(PETSC_SUCCESS);
1139 }
1140 
1141 /*@
1142   PetscDSGetQuadrature - Returns the quadrature, which must agree for all fields in the `PetscDS`
1143 
1144   Not Collective
1145 
1146   Input Parameter:
1147 . prob - The `PetscDS` object
1148 
1149   Output Parameter:
1150 . q - The quadrature object
1151 
1152   Level: intermediate
1153 
1154 .seealso: `PetscDS`, `PetscQuadrature`, `PetscDSSetImplicit()`, `PetscDSSetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1155 @*/
1156 PetscErrorCode PetscDSGetQuadrature(PetscDS prob, PetscQuadrature *q)
1157 {
1158   PetscObject  obj;
1159   PetscClassId id;
1160 
1161   PetscFunctionBegin;
1162   *q = NULL;
1163   if (!prob->Nf) PetscFunctionReturn(PETSC_SUCCESS);
1164   PetscCall(PetscDSGetDiscretization(prob, 0, &obj));
1165   PetscCall(PetscObjectGetClassId(obj, &id));
1166   if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetQuadrature((PetscFE)obj, q));
1167   else if (id == PETSCFV_CLASSID) PetscCall(PetscFVGetQuadrature((PetscFV)obj, q));
1168   else SETERRQ(PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", 0);
1169   PetscFunctionReturn(PETSC_SUCCESS);
1170 }
1171 
1172 /*@
1173   PetscDSGetImplicit - Returns the flag for implicit solve for this field. This is just a guide for `TSIMEX`
1174 
1175   Not Collective
1176 
1177   Input Parameters:
1178 + prob - The `PetscDS` object
1179 - f    - The field number
1180 
1181   Output Parameter:
1182 . implicit - The flag indicating what kind of solve to use for this field
1183 
1184   Level: developer
1185 
1186 .seealso: `TSIMEX`, `PetscDS`, `PetscDSSetImplicit()`, `PetscDSSetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1187 @*/
1188 PetscErrorCode PetscDSGetImplicit(PetscDS prob, PetscInt f, PetscBool *implicit)
1189 {
1190   PetscFunctionBegin;
1191   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
1192   PetscAssertPointer(implicit, 3);
1193   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);
1194   *implicit = prob->implicit[f];
1195   PetscFunctionReturn(PETSC_SUCCESS);
1196 }
1197 
1198 /*@
1199   PetscDSSetImplicit - Set the flag for implicit solve for this field. This is just a guide for `TSIMEX`
1200 
1201   Not Collective
1202 
1203   Input Parameters:
1204 + prob     - The `PetscDS` object
1205 . f        - The field number
1206 - implicit - The flag indicating what kind of solve to use for this field
1207 
1208   Level: developer
1209 
1210 .seealso: `TSIMEX`, `PetscDSGetImplicit()`, `PetscDSSetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1211 @*/
1212 PetscErrorCode PetscDSSetImplicit(PetscDS prob, PetscInt f, PetscBool implicit)
1213 {
1214   PetscFunctionBegin;
1215   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
1216   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);
1217   prob->implicit[f] = implicit;
1218   PetscFunctionReturn(PETSC_SUCCESS);
1219 }
1220 
1221 /*@
1222   PetscDSGetJetDegree - Returns the highest derivative for this field equation, or the k-jet that the discretization needs to tabulate.
1223 
1224   Not Collective
1225 
1226   Input Parameters:
1227 + ds - The `PetscDS` object
1228 - f  - The field number
1229 
1230   Output Parameter:
1231 . k - The highest derivative we need to tabulate
1232 
1233   Level: developer
1234 
1235 .seealso: `PetscDS`, `PetscDSSetJetDegree()`, `PetscDSSetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1236 @*/
1237 PetscErrorCode PetscDSGetJetDegree(PetscDS ds, PetscInt f, PetscInt *k)
1238 {
1239   PetscFunctionBegin;
1240   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
1241   PetscAssertPointer(k, 3);
1242   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);
1243   *k = ds->jetDegree[f];
1244   PetscFunctionReturn(PETSC_SUCCESS);
1245 }
1246 
1247 /*@
1248   PetscDSSetJetDegree - Set the highest derivative for this field equation, or the k-jet that the discretization needs to tabulate.
1249 
1250   Not Collective
1251 
1252   Input Parameters:
1253 + ds - The `PetscDS` object
1254 . f  - The field number
1255 - k  - The highest derivative we need to tabulate
1256 
1257   Level: developer
1258 
1259 .seealso: ``PetscDS`, `PetscDSGetJetDegree()`, `PetscDSSetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1260 @*/
1261 PetscErrorCode PetscDSSetJetDegree(PetscDS ds, PetscInt f, PetscInt k)
1262 {
1263   PetscFunctionBegin;
1264   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
1265   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);
1266   ds->jetDegree[f] = k;
1267   PetscFunctionReturn(PETSC_SUCCESS);
1268 }
1269 
1270 /*@C
1271   PetscDSGetObjective - Get the pointwise objective function for a given test field
1272 
1273   Not Collective
1274 
1275   Input Parameters:
1276 + ds - The `PetscDS`
1277 - f  - The test field number
1278 
1279   Output Parameter:
1280 . obj - integrand for the test function term
1281 
1282   Calling sequence of `obj`:
1283 + dim          - the spatial dimension
1284 . Nf           - the number of fields
1285 . NfAux        - the number of auxiliary fields
1286 . uOff         - the offset into u[] and u_t[] for each field
1287 . uOff_x       - the offset into u_x[] for each field
1288 . u            - each field evaluated at the current point
1289 . u_t          - the time derivative of each field evaluated at the current point
1290 . u_x          - the gradient of each field evaluated at the current point
1291 . aOff         - the offset into a[] and a_t[] for each auxiliary field
1292 . aOff_x       - the offset into a_x[] for each auxiliary field
1293 . a            - each auxiliary field evaluated at the current point
1294 . a_t          - the time derivative of each auxiliary field evaluated at the current point
1295 . a_x          - the gradient of auxiliary each field evaluated at the current point
1296 . t            - current time
1297 . x            - coordinates of the current point
1298 . numConstants - number of constant parameters
1299 . constants    - constant parameters
1300 - obj          - output values at the current point
1301 
1302   Level: intermediate
1303 
1304   Note:
1305   We are using a first order FEM model for the weak form\: $  \int_\Omega \phi obj(u, u_t, \nabla u, x, t)$
1306 
1307 .seealso: `PetscDS`, `PetscDSSetObjective()`, `PetscDSGetResidual()`
1308 @*/
1309 PetscErrorCode PetscDSGetObjective(PetscDS ds, PetscInt f, void (**obj)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[]))
1310 {
1311   PetscPointFunc *tmp;
1312   PetscInt        n;
1313 
1314   PetscFunctionBegin;
1315   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
1316   PetscAssertPointer(obj, 3);
1317   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);
1318   PetscCall(PetscWeakFormGetObjective(ds->wf, NULL, 0, f, 0, &n, &tmp));
1319   *obj = tmp ? tmp[0] : NULL;
1320   PetscFunctionReturn(PETSC_SUCCESS);
1321 }
1322 
1323 /*@C
1324   PetscDSSetObjective - Set the pointwise objective function for a given test field
1325 
1326   Not Collective
1327 
1328   Input Parameters:
1329 + ds  - The `PetscDS`
1330 . f   - The test field number
1331 - obj - integrand for the test function term
1332 
1333   Calling sequence of `obj`:
1334 + dim          - the spatial dimension
1335 . Nf           - the number of fields
1336 . NfAux        - the number of auxiliary fields
1337 . uOff         - the offset into u[] and u_t[] for each field
1338 . uOff_x       - the offset into u_x[] for each field
1339 . u            - each field evaluated at the current point
1340 . u_t          - the time derivative of each field evaluated at the current point
1341 . u_x          - the gradient of each field evaluated at the current point
1342 . aOff         - the offset into a[] and a_t[] for each auxiliary field
1343 . aOff_x       - the offset into a_x[] for each auxiliary field
1344 . a            - each auxiliary field evaluated at the current point
1345 . a_t          - the time derivative of each auxiliary field evaluated at the current point
1346 . a_x          - the gradient of auxiliary each field evaluated at the current point
1347 . t            - current time
1348 . x            - coordinates of the current point
1349 . numConstants - number of constant parameters
1350 . constants    - constant parameters
1351 - obj          - output values at the current point
1352 
1353   Level: intermediate
1354 
1355   Note:
1356   We are using a first order FEM model for the weak form\: $  \int_\Omega \phi obj(u, u_t, \nabla u, x, t)$
1357 
1358 .seealso: `PetscDS`, `PetscDSGetObjective()`, `PetscDSSetResidual()`
1359 @*/
1360 PetscErrorCode PetscDSSetObjective(PetscDS ds, PetscInt f, void (*obj)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[]))
1361 {
1362   PetscFunctionBegin;
1363   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
1364   if (obj) PetscValidFunction(obj, 3);
1365   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
1366   PetscCall(PetscWeakFormSetIndexObjective(ds->wf, NULL, 0, f, 0, 0, obj));
1367   PetscFunctionReturn(PETSC_SUCCESS);
1368 }
1369 
1370 /*@C
1371   PetscDSGetResidual - Get the pointwise residual function for a given test field
1372 
1373   Not Collective
1374 
1375   Input Parameters:
1376 + ds - The `PetscDS`
1377 - f  - The test field number
1378 
1379   Output Parameters:
1380 + f0 - integrand for the test function term
1381 - f1 - integrand for the test function gradient term
1382 
1383   Calling sequence of `f0`:
1384 + dim          - the spatial dimension
1385 . Nf           - the number of fields
1386 . NfAux        - the number of auxiliary fields
1387 . uOff         - the offset into u[] and u_t[] for each field
1388 . uOff_x       - the offset into u_x[] for each field
1389 . u            - each field evaluated at the current point
1390 . u_t          - the time derivative of each field evaluated at the current point
1391 . u_x          - the gradient of each field evaluated at the current point
1392 . aOff         - the offset into a[] and a_t[] for each auxiliary field
1393 . aOff_x       - the offset into a_x[] for each auxiliary field
1394 . a            - each auxiliary field evaluated at the current point
1395 . a_t          - the time derivative of each auxiliary field evaluated at the current point
1396 . a_x          - the gradient of auxiliary each field evaluated at the current point
1397 . t            - current time
1398 . x            - coordinates of the current point
1399 . numConstants - number of constant parameters
1400 . constants    - constant parameters
1401 - f0           - output values at the current point
1402 
1403   Level: intermediate
1404 
1405   Note:
1406   `f1` has an identical form and is omitted for brevity.
1407 
1408   We are using a first order FEM model for the weak form\: $  \int_\Omega \phi f_0(u, u_t, \nabla u, x, t) + \nabla\phi \cdot {\vec f}_1(u, u_t, \nabla u, x, t)$
1409 
1410 .seealso: `PetscDS`, `PetscDSSetResidual()`
1411 @*/
1412 PetscErrorCode PetscDSGetResidual(PetscDS ds, PetscInt f, void (**f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]), void (**f1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))
1413 {
1414   PetscPointFunc *tmp0, *tmp1;
1415   PetscInt        n0, n1;
1416 
1417   PetscFunctionBegin;
1418   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
1419   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);
1420   PetscCall(PetscWeakFormGetResidual(ds->wf, NULL, 0, f, 0, &n0, &tmp0, &n1, &tmp1));
1421   *f0 = tmp0 ? tmp0[0] : NULL;
1422   *f1 = tmp1 ? tmp1[0] : NULL;
1423   PetscFunctionReturn(PETSC_SUCCESS);
1424 }
1425 
1426 /*@C
1427   PetscDSSetResidual - Set the pointwise residual function for a given test field
1428 
1429   Not Collective
1430 
1431   Input Parameters:
1432 + ds - The `PetscDS`
1433 . f  - The test field number
1434 . f0 - integrand for the test function term
1435 - f1 - integrand for the test function gradient term
1436 
1437   Calling sequence of `f0`:
1438 + dim          - the spatial dimension
1439 . Nf           - the number of fields
1440 . NfAux        - the number of auxiliary fields
1441 . uOff         - the offset into u[] and u_t[] for each field
1442 . uOff_x       - the offset into u_x[] for each field
1443 . u            - each field evaluated at the current point
1444 . u_t          - the time derivative of each field evaluated at the current point
1445 . u_x          - the gradient of each field evaluated at the current point
1446 . aOff         - the offset into a[] and a_t[] for each auxiliary field
1447 . aOff_x       - the offset into a_x[] for each auxiliary field
1448 . a            - each auxiliary field evaluated at the current point
1449 . a_t          - the time derivative of each auxiliary field evaluated at the current point
1450 . a_x          - the gradient of auxiliary each field evaluated at the current point
1451 . t            - current time
1452 . x            - coordinates of the current point
1453 . numConstants - number of constant parameters
1454 . constants    - constant parameters
1455 - f0           - output values at the current point
1456 
1457   Level: intermediate
1458 
1459   Note:
1460   `f1` has an identical form and is omitted for brevity.
1461 
1462   We are using a first order FEM model for the weak form\: $  \int_\Omega \phi f_0(u, u_t, \nabla u, x, t) + \nabla\phi \cdot {\vec f}_1(u, u_t, \nabla u, x, t)$
1463 
1464 .seealso: `PetscDS`, `PetscDSGetResidual()`
1465 @*/
1466 PetscErrorCode PetscDSSetResidual(PetscDS ds, PetscInt f, void (*f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]), void (*f1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))
1467 {
1468   PetscFunctionBegin;
1469   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
1470   if (f0) PetscValidFunction(f0, 3);
1471   if (f1) PetscValidFunction(f1, 4);
1472   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
1473   PetscCall(PetscWeakFormSetIndexResidual(ds->wf, NULL, 0, f, 0, 0, f0, 0, f1));
1474   PetscFunctionReturn(PETSC_SUCCESS);
1475 }
1476 
1477 /*@C
1478   PetscDSGetRHSResidual - Get the pointwise RHS residual function for explicit timestepping for a given test field
1479 
1480   Not Collective
1481 
1482   Input Parameters:
1483 + ds - The `PetscDS`
1484 - f  - The test field number
1485 
1486   Output Parameters:
1487 + f0 - integrand for the test function term
1488 - f1 - integrand for the test function gradient term
1489 
1490   Calling sequence of `f0`:
1491 + dim          - the spatial dimension
1492 . Nf           - the number of fields
1493 . NfAux        - the number of auxiliary fields
1494 . uOff         - the offset into u[] and u_t[] for each field
1495 . uOff_x       - the offset into u_x[] for each field
1496 . u            - each field evaluated at the current point
1497 . u_t          - the time derivative of each field evaluated at the current point
1498 . u_x          - the gradient of each field evaluated at the current point
1499 . aOff         - the offset into a[] and a_t[] for each auxiliary field
1500 . aOff_x       - the offset into a_x[] for each auxiliary field
1501 . a            - each auxiliary field evaluated at the current point
1502 . a_t          - the time derivative of each auxiliary field evaluated at the current point
1503 . a_x          - the gradient of auxiliary each field evaluated at the current point
1504 . t            - current time
1505 . x            - coordinates of the current point
1506 . numConstants - number of constant parameters
1507 . constants    - constant parameters
1508 - f0           - output values at the current point
1509 
1510   Level: intermediate
1511 
1512   Note:
1513   `f1` has an identical form and is omitted for brevity.
1514 
1515   We are using a first order FEM model for the weak form\: $ \int_\Omega \phi f_0(u, u_t, \nabla u, x, t) + \nabla\phi \cdot {\vec f}_1(u, u_t, \nabla u, x, t)$
1516 
1517 .seealso: `PetscDS`, `PetscDSSetRHSResidual()`
1518 @*/
1519 PetscErrorCode PetscDSGetRHSResidual(PetscDS ds, PetscInt f, void (**f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]), void (**f1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))
1520 {
1521   PetscPointFunc *tmp0, *tmp1;
1522   PetscInt        n0, n1;
1523 
1524   PetscFunctionBegin;
1525   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
1526   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);
1527   PetscCall(PetscWeakFormGetResidual(ds->wf, NULL, 0, f, 100, &n0, &tmp0, &n1, &tmp1));
1528   *f0 = tmp0 ? tmp0[0] : NULL;
1529   *f1 = tmp1 ? tmp1[0] : NULL;
1530   PetscFunctionReturn(PETSC_SUCCESS);
1531 }
1532 
1533 /*@C
1534   PetscDSSetRHSResidual - Set the pointwise residual function for explicit timestepping for a given test field
1535 
1536   Not Collective
1537 
1538   Input Parameters:
1539 + ds - The `PetscDS`
1540 . f  - The test field number
1541 . f0 - integrand for the test function term
1542 - f1 - integrand for the test function gradient term
1543 
1544   Calling sequence for the callbacks `f0`:
1545 + dim          - the spatial dimension
1546 . Nf           - the number of fields
1547 . NfAux        - the number of auxiliary fields
1548 . uOff         - the offset into u[] and u_t[] for each field
1549 . uOff_x       - the offset into u_x[] for each field
1550 . u            - each field evaluated at the current point
1551 . u_t          - the time derivative of each field evaluated at the current point
1552 . u_x          - the gradient of each field evaluated at the current point
1553 . aOff         - the offset into a[] and a_t[] for each auxiliary field
1554 . aOff_x       - the offset into a_x[] for each auxiliary field
1555 . a            - each auxiliary field evaluated at the current point
1556 . a_t          - the time derivative of each auxiliary field evaluated at the current point
1557 . a_x          - the gradient of auxiliary each field evaluated at the current point
1558 . t            - current time
1559 . x            - coordinates of the current point
1560 . numConstants - number of constant parameters
1561 . constants    - constant parameters
1562 - f0           - output values at the current point
1563 
1564   Level: intermediate
1565 
1566   Note:
1567   `f1` has an identical form and is omitted for brevity.
1568 
1569   We are using a first order FEM model for the weak form\: $ \int_\Omega \phi f_0(u, u_t, \nabla u, x, t) + \nabla\phi \cdot {\vec f}_1(u, u_t, \nabla u, x, t)$
1570 
1571 .seealso: `PetscDS`, `PetscDSGetResidual()`
1572 @*/
1573 PetscErrorCode PetscDSSetRHSResidual(PetscDS ds, PetscInt f, void (*f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]), void (*f1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))
1574 {
1575   PetscFunctionBegin;
1576   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
1577   if (f0) PetscValidFunction(f0, 3);
1578   if (f1) PetscValidFunction(f1, 4);
1579   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
1580   PetscCall(PetscWeakFormSetIndexResidual(ds->wf, NULL, 0, f, 100, 0, f0, 0, f1));
1581   PetscFunctionReturn(PETSC_SUCCESS);
1582 }
1583 
1584 /*@C
1585   PetscDSHasJacobian - Checks that the Jacobian functions have been set
1586 
1587   Not Collective
1588 
1589   Input Parameter:
1590 . ds - The `PetscDS`
1591 
1592   Output Parameter:
1593 . hasJac - flag that pointwise function for the Jacobian has been set
1594 
1595   Level: intermediate
1596 
1597 .seealso: `PetscDS`, `PetscDSGetJacobianPreconditioner()`, `PetscDSSetJacobianPreconditioner()`, `PetscDSGetJacobian()`
1598 @*/
1599 PetscErrorCode PetscDSHasJacobian(PetscDS ds, PetscBool *hasJac)
1600 {
1601   PetscFunctionBegin;
1602   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
1603   PetscCall(PetscWeakFormHasJacobian(ds->wf, hasJac));
1604   PetscFunctionReturn(PETSC_SUCCESS);
1605 }
1606 
1607 /*@C
1608   PetscDSGetJacobian - Get the pointwise Jacobian function for given test and basis field
1609 
1610   Not Collective
1611 
1612   Input Parameters:
1613 + ds - The `PetscDS`
1614 . f  - The test field number
1615 - g  - The field number
1616 
1617   Output Parameters:
1618 + g0 - integrand for the test and basis function term
1619 . g1 - integrand for the test function and basis function gradient term
1620 . g2 - integrand for the test function gradient and basis function term
1621 - g3 - integrand for the test function gradient and basis function gradient term
1622 
1623   Calling sequence of `g0`:
1624 + dim          - the spatial dimension
1625 . Nf           - the number of fields
1626 . NfAux        - the number of auxiliary fields
1627 . uOff         - the offset into u[] and u_t[] for each field
1628 . uOff_x       - the offset into u_x[] for each field
1629 . u            - each field evaluated at the current point
1630 . u_t          - the time derivative of each field evaluated at the current point
1631 . u_x          - the gradient of each field evaluated at the current point
1632 . aOff         - the offset into a[] and a_t[] for each auxiliary field
1633 . aOff_x       - the offset into a_x[] for each auxiliary field
1634 . a            - each auxiliary field evaluated at the current point
1635 . a_t          - the time derivative of each auxiliary field evaluated at the current point
1636 . a_x          - the gradient of auxiliary each field evaluated at the current point
1637 . t            - current time
1638 . u_tShift     - the multiplier a for dF/dU_t
1639 . x            - coordinates of the current point
1640 . numConstants - number of constant parameters
1641 . constants    - constant parameters
1642 - g0           - output values at the current point
1643 
1644   Level: intermediate
1645 
1646   Note:
1647   `g1`, `g2`, and `g3` have identical calling sequences to `g0` and are omitted for brevity.
1648 
1649   We are using a first order FEM model for the weak form\:
1650 
1651   $$
1652   \int_\Omega \phi g_0(u, u_t, \nabla u, x, t) \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi
1653   $$
1654 
1655 .seealso: `PetscDS`, `PetscDSSetJacobian()`
1656 @*/
1657 PetscErrorCode PetscDSGetJacobian(PetscDS ds, PetscInt f, PetscInt g, void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]), void (**g1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (**g2)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (**g3)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))
1658 {
1659   PetscPointJac *tmp0, *tmp1, *tmp2, *tmp3;
1660   PetscInt       n0, n1, n2, n3;
1661 
1662   PetscFunctionBegin;
1663   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
1664   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);
1665   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);
1666   PetscCall(PetscWeakFormGetJacobian(ds->wf, NULL, 0, f, g, 0, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3));
1667   *g0 = tmp0 ? tmp0[0] : NULL;
1668   *g1 = tmp1 ? tmp1[0] : NULL;
1669   *g2 = tmp2 ? tmp2[0] : NULL;
1670   *g3 = tmp3 ? tmp3[0] : NULL;
1671   PetscFunctionReturn(PETSC_SUCCESS);
1672 }
1673 
1674 /*@C
1675   PetscDSSetJacobian - Set the pointwise Jacobian function for given test and basis fields
1676 
1677   Not Collective
1678 
1679   Input Parameters:
1680 + ds - The `PetscDS`
1681 . f  - The test field number
1682 . g  - The field number
1683 . g0 - integrand for the test and basis function term
1684 . g1 - integrand for the test function and basis function gradient term
1685 . g2 - integrand for the test function gradient and basis function term
1686 - g3 - integrand for the test function gradient and basis function gradient term
1687 
1688   Calling sequence of `g0`:
1689 + dim          - the spatial dimension
1690 . Nf           - the number of fields
1691 . NfAux        - the number of auxiliary fields
1692 . uOff         - the offset into u[] and u_t[] for each field
1693 . uOff_x       - the offset into u_x[] for each field
1694 . u            - each field evaluated at the current point
1695 . u_t          - the time derivative of each field evaluated at the current point
1696 . u_x          - the gradient of each field evaluated at the current point
1697 . aOff         - the offset into a[] and a_t[] for each auxiliary field
1698 . aOff_x       - the offset into a_x[] for each auxiliary field
1699 . a            - each auxiliary field evaluated at the current point
1700 . a_t          - the time derivative of each auxiliary field evaluated at the current point
1701 . a_x          - the gradient of auxiliary each field evaluated at the current point
1702 . t            - current time
1703 . u_tShift     - the multiplier a for dF/dU_t
1704 . x            - coordinates of the current point
1705 . numConstants - number of constant parameters
1706 . constants    - constant parameters
1707 - g0           - output values at the current point
1708 
1709   Level: intermediate
1710 
1711   Note:
1712   `g1`, `g2`, and `g3` have identical calling sequences to `g0` and are omitted for brevity.
1713 
1714   We are using a first order FEM model for the weak form\:
1715   \int_\Omega \phi g_0(u, u_t, \nabla u, x, t) \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi
1716 
1717 .seealso: `PetscDS`, `PetscDSGetJacobian()`
1718 @*/
1719 PetscErrorCode PetscDSSetJacobian(PetscDS ds, PetscInt f, PetscInt g, void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]), void (*g1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (*g2)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (*g3)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))
1720 {
1721   PetscFunctionBegin;
1722   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
1723   if (g0) PetscValidFunction(g0, 4);
1724   if (g1) PetscValidFunction(g1, 5);
1725   if (g2) PetscValidFunction(g2, 6);
1726   if (g3) PetscValidFunction(g3, 7);
1727   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
1728   PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", g);
1729   PetscCall(PetscWeakFormSetIndexJacobian(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3));
1730   PetscFunctionReturn(PETSC_SUCCESS);
1731 }
1732 
1733 /*@C
1734   PetscDSUseJacobianPreconditioner - Set whether to construct a Jacobian preconditioner
1735 
1736   Not Collective
1737 
1738   Input Parameters:
1739 + prob      - The `PetscDS`
1740 - useJacPre - flag that enables construction of a Jacobian preconditioner
1741 
1742   Level: intermediate
1743 
1744 .seealso: `PetscDS`, `PetscDSGetJacobianPreconditioner()`, `PetscDSSetJacobianPreconditioner()`, `PetscDSGetJacobian()`
1745 @*/
1746 PetscErrorCode PetscDSUseJacobianPreconditioner(PetscDS prob, PetscBool useJacPre)
1747 {
1748   PetscFunctionBegin;
1749   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
1750   prob->useJacPre = useJacPre;
1751   PetscFunctionReturn(PETSC_SUCCESS);
1752 }
1753 
1754 /*@C
1755   PetscDSHasJacobianPreconditioner - Checks if a Jacobian preconditioner matrix has been set
1756 
1757   Not Collective
1758 
1759   Input Parameter:
1760 . ds - The `PetscDS`
1761 
1762   Output Parameter:
1763 . hasJacPre - flag that pointwise function for Jacobian preconditioner matrix has been set
1764 
1765   Level: intermediate
1766 
1767 .seealso: `PetscDS`, `PetscDSGetJacobianPreconditioner()`, `PetscDSSetJacobianPreconditioner()`, `PetscDSGetJacobian()`
1768 @*/
1769 PetscErrorCode PetscDSHasJacobianPreconditioner(PetscDS ds, PetscBool *hasJacPre)
1770 {
1771   PetscFunctionBegin;
1772   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
1773   *hasJacPre = PETSC_FALSE;
1774   if (!ds->useJacPre) PetscFunctionReturn(PETSC_SUCCESS);
1775   PetscCall(PetscWeakFormHasJacobianPreconditioner(ds->wf, hasJacPre));
1776   PetscFunctionReturn(PETSC_SUCCESS);
1777 }
1778 
1779 /*@C
1780   PetscDSGetJacobianPreconditioner - Get the pointwise Jacobian preconditioner function for given test and basis field. If this is missing,
1781   the system matrix is used to build the preconditioner.
1782 
1783   Not Collective
1784 
1785   Input Parameters:
1786 + ds - The `PetscDS`
1787 . f  - The test field number
1788 - g  - The field number
1789 
1790   Output Parameters:
1791 + g0 - integrand for the test and basis function term
1792 . g1 - integrand for the test function and basis function gradient term
1793 . g2 - integrand for the test function gradient and basis function term
1794 - g3 - integrand for the test function gradient and basis function gradient term
1795 
1796   Calling sequence of `g0`:
1797 + dim          - the spatial dimension
1798 . Nf           - the number of fields
1799 . NfAux        - the number of auxiliary fields
1800 . uOff         - the offset into u[] and u_t[] for each field
1801 . uOff_x       - the offset into u_x[] for each field
1802 . u            - each field evaluated at the current point
1803 . u_t          - the time derivative of each field evaluated at the current point
1804 . u_x          - the gradient of each field evaluated at the current point
1805 . aOff         - the offset into a[] and a_t[] for each auxiliary field
1806 . aOff_x       - the offset into a_x[] for each auxiliary field
1807 . a            - each auxiliary field evaluated at the current point
1808 . a_t          - the time derivative of each auxiliary field evaluated at the current point
1809 . a_x          - the gradient of auxiliary each field evaluated at the current point
1810 . t            - current time
1811 . u_tShift     - the multiplier a for dF/dU_t
1812 . x            - coordinates of the current point
1813 . numConstants - number of constant parameters
1814 . constants    - constant parameters
1815 - g0           - output values at the current point
1816 
1817   Level: intermediate
1818 
1819   Note:
1820   `g1`, `g2`, and `g3` have identical calling sequences to `g0` and are omitted for brevity.
1821   We are using a first order FEM model for the weak form\:
1822   \int_\Omega \phi g_0(u, u_t, \nabla u, x, t) \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi
1823 
1824 .seealso: `PetscDS`, `PetscDSSetJacobianPreconditioner()`, `PetscDSGetJacobian()`
1825 @*/
1826 PetscErrorCode PetscDSGetJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g, void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]), void (**g1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (**g2)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (**g3)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))
1827 {
1828   PetscPointJac *tmp0, *tmp1, *tmp2, *tmp3;
1829   PetscInt       n0, n1, n2, n3;
1830 
1831   PetscFunctionBegin;
1832   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
1833   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);
1834   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);
1835   PetscCall(PetscWeakFormGetJacobianPreconditioner(ds->wf, NULL, 0, f, g, 0, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3));
1836   *g0 = tmp0 ? tmp0[0] : NULL;
1837   *g1 = tmp1 ? tmp1[0] : NULL;
1838   *g2 = tmp2 ? tmp2[0] : NULL;
1839   *g3 = tmp3 ? tmp3[0] : NULL;
1840   PetscFunctionReturn(PETSC_SUCCESS);
1841 }
1842 
1843 /*@C
1844   PetscDSSetJacobianPreconditioner - Set the pointwise Jacobian preconditioner function for given test and basis fields.
1845   If this is missing, the system matrix is used to build the preconditioner.
1846 
1847   Not Collective
1848 
1849   Input Parameters:
1850 + ds - The `PetscDS`
1851 . f  - The test field number
1852 . g  - The field number
1853 . g0 - integrand for the test and basis function term
1854 . g1 - integrand for the test function and basis function gradient term
1855 . g2 - integrand for the test function gradient and basis function term
1856 - g3 - integrand for the test function gradient and basis function gradient term
1857 
1858   Calling sequence of `g0`:
1859 + dim          - the spatial dimension
1860 . Nf           - the number of fields
1861 . NfAux        - the number of auxiliary fields
1862 . uOff         - the offset into u[] and u_t[] for each field
1863 . uOff_x       - the offset into u_x[] for each field
1864 . u            - each field evaluated at the current point
1865 . u_t          - the time derivative of each field evaluated at the current point
1866 . u_x          - the gradient of each field evaluated at the current point
1867 . aOff         - the offset into a[] and a_t[] for each auxiliary field
1868 . aOff_x       - the offset into a_x[] for each auxiliary field
1869 . a            - each auxiliary field evaluated at the current point
1870 . a_t          - the time derivative of each auxiliary field evaluated at the current point
1871 . a_x          - the gradient of auxiliary each field evaluated at the current point
1872 . t            - current time
1873 . u_tShift     - the multiplier a for dF/dU_t
1874 . x            - coordinates of the current point
1875 . numConstants - number of constant parameters
1876 . constants    - constant parameters
1877 - g0           - output values at the current point
1878 
1879   Level: intermediate
1880 
1881   Note:
1882   `g1`, `g2`, and `g3` have identical calling sequences to `g0` and are omitted for brevity.
1883 
1884   We are using a first order FEM model for the weak form\:
1885   \int_\Omega \phi g_0(u, u_t, \nabla u, x, t) \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi
1886 
1887 .seealso: `PetscDS`, `PetscDSGetJacobianPreconditioner()`, `PetscDSSetJacobian()`
1888 @*/
1889 PetscErrorCode PetscDSSetJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g, void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]), void (*g1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (*g2)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (*g3)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))
1890 {
1891   PetscFunctionBegin;
1892   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
1893   if (g0) PetscValidFunction(g0, 4);
1894   if (g1) PetscValidFunction(g1, 5);
1895   if (g2) PetscValidFunction(g2, 6);
1896   if (g3) PetscValidFunction(g3, 7);
1897   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
1898   PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", g);
1899   PetscCall(PetscWeakFormSetIndexJacobianPreconditioner(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3));
1900   PetscFunctionReturn(PETSC_SUCCESS);
1901 }
1902 
1903 /*@C
1904   PetscDSHasDynamicJacobian - Signals that a dynamic Jacobian, dF/du_t, has been set
1905 
1906   Not Collective
1907 
1908   Input Parameter:
1909 . ds - The `PetscDS`
1910 
1911   Output Parameter:
1912 . hasDynJac - flag that pointwise function for dynamic Jacobian has been set
1913 
1914   Level: intermediate
1915 
1916 .seealso: `PetscDS`, `PetscDSGetDynamicJacobian()`, `PetscDSSetDynamicJacobian()`, `PetscDSGetJacobian()`
1917 @*/
1918 PetscErrorCode PetscDSHasDynamicJacobian(PetscDS ds, PetscBool *hasDynJac)
1919 {
1920   PetscFunctionBegin;
1921   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
1922   PetscCall(PetscWeakFormHasDynamicJacobian(ds->wf, hasDynJac));
1923   PetscFunctionReturn(PETSC_SUCCESS);
1924 }
1925 
1926 /*@C
1927   PetscDSGetDynamicJacobian - Get the pointwise dynamic Jacobian, dF/du_t, function for given test and basis field
1928 
1929   Not Collective
1930 
1931   Input Parameters:
1932 + ds - The `PetscDS`
1933 . f  - The test field number
1934 - g  - The field number
1935 
1936   Output Parameters:
1937 + g0 - integrand for the test and basis function term
1938 . g1 - integrand for the test function and basis function gradient term
1939 . g2 - integrand for the test function gradient and basis function term
1940 - g3 - integrand for the test function gradient and basis function gradient term
1941 
1942   Calling sequence of `g0`:
1943 + dim          - the spatial dimension
1944 . Nf           - the number of fields
1945 . NfAux        - the number of auxiliary fields
1946 . uOff         - the offset into u[] and u_t[] for each field
1947 . uOff_x       - the offset into u_x[] for each field
1948 . u            - each field evaluated at the current point
1949 . u_t          - the time derivative of each field evaluated at the current point
1950 . u_x          - the gradient of each field evaluated at the current point
1951 . aOff         - the offset into a[] and a_t[] for each auxiliary field
1952 . aOff_x       - the offset into a_x[] for each auxiliary field
1953 . a            - each auxiliary field evaluated at the current point
1954 . a_t          - the time derivative of each auxiliary field evaluated at the current point
1955 . a_x          - the gradient of auxiliary each field evaluated at the current point
1956 . t            - current time
1957 . u_tShift     - the multiplier a for dF/dU_t
1958 . x            - coordinates of the current point
1959 . numConstants - number of constant parameters
1960 . constants    - constant parameters
1961 - g0           - output values at the current point
1962 
1963   Level: intermediate
1964 
1965   Note:
1966   `g1`, `g2`, and `g3` have identical calling sequences to `g0` and are omitted for brevity.
1967 
1968   We are using a first order FEM model for the weak form\:
1969   \int_\Omega \phi g_0(u, u_t, \nabla u, x, t) \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi
1970 
1971 .seealso: `PetscDS`, `PetscDSSetJacobian()`
1972 @*/
1973 PetscErrorCode PetscDSGetDynamicJacobian(PetscDS ds, PetscInt f, PetscInt g, void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]), void (**g1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (**g2)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (**g3)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))
1974 {
1975   PetscPointJac *tmp0, *tmp1, *tmp2, *tmp3;
1976   PetscInt       n0, n1, n2, n3;
1977 
1978   PetscFunctionBegin;
1979   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
1980   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);
1981   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);
1982   PetscCall(PetscWeakFormGetDynamicJacobian(ds->wf, NULL, 0, f, g, 0, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3));
1983   *g0 = tmp0 ? tmp0[0] : NULL;
1984   *g1 = tmp1 ? tmp1[0] : NULL;
1985   *g2 = tmp2 ? tmp2[0] : NULL;
1986   *g3 = tmp3 ? tmp3[0] : NULL;
1987   PetscFunctionReturn(PETSC_SUCCESS);
1988 }
1989 
1990 /*@C
1991   PetscDSSetDynamicJacobian - Set the pointwise dynamic Jacobian, dF/du_t, function for given test and basis fields
1992 
1993   Not Collective
1994 
1995   Input Parameters:
1996 + ds - The `PetscDS`
1997 . f  - The test field number
1998 . g  - The field number
1999 . g0 - integrand for the test and basis function term
2000 . g1 - integrand for the test function and basis function gradient term
2001 . g2 - integrand for the test function gradient and basis function term
2002 - g3 - integrand for the test function gradient and basis function gradient term
2003 
2004   Calling sequence of `g0`:
2005 + dim          - the spatial dimension
2006 . Nf           - the number of fields
2007 . NfAux        - the number of auxiliary fields
2008 . uOff         - the offset into u[] and u_t[] for each field
2009 . uOff_x       - the offset into u_x[] for each field
2010 . u            - each field evaluated at the current point
2011 . u_t          - the time derivative of each field evaluated at the current point
2012 . u_x          - the gradient of each field evaluated at the current point
2013 . aOff         - the offset into a[] and a_t[] for each auxiliary field
2014 . aOff_x       - the offset into a_x[] for each auxiliary field
2015 . a            - each auxiliary field evaluated at the current point
2016 . a_t          - the time derivative of each auxiliary field evaluated at the current point
2017 . a_x          - the gradient of auxiliary each field evaluated at the current point
2018 . t            - current time
2019 . u_tShift     - the multiplier a for dF/dU_t
2020 . x            - coordinates of the current point
2021 . numConstants - number of constant parameters
2022 . constants    - constant parameters
2023 - g0           - output values at the current point
2024 
2025   Level: intermediate
2026 
2027   Note:
2028   `g1`, `g2`, and `g3` have identical calling sequences to `g0` and are omitted for brevity.
2029 
2030   We are using a first order FEM model for the weak form\:
2031   \int_\Omega \phi g_0(u, u_t, \nabla u, x, t) \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi
2032 
2033 .seealso: `PetscDS`, `PetscDSGetJacobian()`
2034 @*/
2035 PetscErrorCode PetscDSSetDynamicJacobian(PetscDS ds, PetscInt f, PetscInt g, void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]), void (*g1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (*g2)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (*g3)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))
2036 {
2037   PetscFunctionBegin;
2038   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
2039   if (g0) PetscValidFunction(g0, 4);
2040   if (g1) PetscValidFunction(g1, 5);
2041   if (g2) PetscValidFunction(g2, 6);
2042   if (g3) PetscValidFunction(g3, 7);
2043   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2044   PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", g);
2045   PetscCall(PetscWeakFormSetIndexDynamicJacobian(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3));
2046   PetscFunctionReturn(PETSC_SUCCESS);
2047 }
2048 
2049 /*@C
2050   PetscDSGetRiemannSolver - Returns the Riemann solver for the given field
2051 
2052   Not Collective
2053 
2054   Input Parameters:
2055 + ds - The `PetscDS` object
2056 - f  - The field number
2057 
2058   Output Parameter:
2059 . r - Riemann solver
2060 
2061   Calling sequence of `r`:
2062 + dim          - The spatial dimension
2063 . Nf           - The number of fields
2064 . x            - The coordinates at a point on the interface
2065 . n            - The normal vector to the interface
2066 . uL           - The state vector to the left of the interface
2067 . uR           - The state vector to the right of the interface
2068 . flux         - output array of flux through the interface
2069 . numConstants - number of constant parameters
2070 . constants    - constant parameters
2071 - ctx          - optional user context
2072 
2073   Level: intermediate
2074 
2075 .seealso: `PetscDS`, `PetscDSSetRiemannSolver()`
2076 @*/
2077 PetscErrorCode PetscDSGetRiemannSolver(PetscDS ds, PetscInt f, void (**r)(PetscInt dim, PetscInt Nf, const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscInt numConstants, const PetscScalar constants[], PetscScalar flux[], void *ctx))
2078 {
2079   PetscRiemannFunc *tmp;
2080   PetscInt          n;
2081 
2082   PetscFunctionBegin;
2083   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
2084   PetscAssertPointer(r, 3);
2085   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);
2086   PetscCall(PetscWeakFormGetRiemannSolver(ds->wf, NULL, 0, f, 0, &n, &tmp));
2087   *r = tmp ? tmp[0] : NULL;
2088   PetscFunctionReturn(PETSC_SUCCESS);
2089 }
2090 
2091 /*@C
2092   PetscDSSetRiemannSolver - Sets the Riemann solver for the given field
2093 
2094   Not Collective
2095 
2096   Input Parameters:
2097 + ds - The `PetscDS` object
2098 . f  - The field number
2099 - r  - Riemann solver
2100 
2101   Calling sequence of `r`:
2102 + dim          - The spatial dimension
2103 . Nf           - The number of fields
2104 . x            - The coordinates at a point on the interface
2105 . n            - The normal vector to the interface
2106 . uL           - The state vector to the left of the interface
2107 . uR           - The state vector to the right of the interface
2108 . flux         - output array of flux through the interface
2109 . numConstants - number of constant parameters
2110 . constants    - constant parameters
2111 - ctx          - optional user context
2112 
2113   Level: intermediate
2114 
2115 .seealso: `PetscDS`, `PetscDSGetRiemannSolver()`
2116 @*/
2117 PetscErrorCode PetscDSSetRiemannSolver(PetscDS ds, PetscInt f, void (*r)(PetscInt dim, PetscInt Nf, const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscInt numConstants, const PetscScalar constants[], PetscScalar flux[], void *ctx))
2118 {
2119   PetscFunctionBegin;
2120   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
2121   if (r) PetscValidFunction(r, 3);
2122   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2123   PetscCall(PetscWeakFormSetIndexRiemannSolver(ds->wf, NULL, 0, f, 0, 0, r));
2124   PetscFunctionReturn(PETSC_SUCCESS);
2125 }
2126 
2127 /*@C
2128   PetscDSGetUpdate - Get the pointwise update function for a given field
2129 
2130   Not Collective
2131 
2132   Input Parameters:
2133 + ds - The `PetscDS`
2134 - f  - The field number
2135 
2136   Output Parameter:
2137 . update - update function
2138 
2139   Calling sequence of `update`:
2140 + dim          - the spatial dimension
2141 . Nf           - the number of fields
2142 . NfAux        - the number of auxiliary fields
2143 . uOff         - the offset into u[] and u_t[] for each field
2144 . uOff_x       - the offset into u_x[] for each field
2145 . u            - each field evaluated at the current point
2146 . u_t          - the time derivative of each field evaluated at the current point
2147 . u_x          - the gradient of each field evaluated at the current point
2148 . aOff         - the offset into a[] and a_t[] for each auxiliary field
2149 . aOff_x       - the offset into a_x[] for each auxiliary field
2150 . a            - each auxiliary field evaluated at the current point
2151 . a_t          - the time derivative of each auxiliary field evaluated at the current point
2152 . a_x          - the gradient of auxiliary each field evaluated at the current point
2153 . t            - current time
2154 . x            - coordinates of the current point
2155 . numConstants - number of constant parameters
2156 . constants    - constant parameters
2157 - uNew         - new value for field at the current point
2158 
2159   Level: intermediate
2160 
2161 .seealso: `PetscDS`, `PetscDSSetUpdate()`, `PetscDSSetResidual()`
2162 @*/
2163 PetscErrorCode PetscDSGetUpdate(PetscDS ds, PetscInt f, void (**update)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar uNew[]))
2164 {
2165   PetscFunctionBegin;
2166   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
2167   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);
2168   if (update) {
2169     PetscAssertPointer(update, 3);
2170     *update = ds->update[f];
2171   }
2172   PetscFunctionReturn(PETSC_SUCCESS);
2173 }
2174 
2175 /*@C
2176   PetscDSSetUpdate - Set the pointwise update function for a given field
2177 
2178   Not Collective
2179 
2180   Input Parameters:
2181 + ds     - The `PetscDS`
2182 . f      - The field number
2183 - update - update function
2184 
2185   Calling sequence of `update`:
2186 + dim          - the spatial dimension
2187 . Nf           - the number of fields
2188 . NfAux        - the number of auxiliary fields
2189 . uOff         - the offset into u[] and u_t[] for each field
2190 . uOff_x       - the offset into u_x[] for each field
2191 . u            - each field evaluated at the current point
2192 . u_t          - the time derivative of each field evaluated at the current point
2193 . u_x          - the gradient of each field evaluated at the current point
2194 . aOff         - the offset into a[] and a_t[] for each auxiliary field
2195 . aOff_x       - the offset into a_x[] for each auxiliary field
2196 . a            - each auxiliary field evaluated at the current point
2197 . a_t          - the time derivative of each auxiliary field evaluated at the current point
2198 . a_x          - the gradient of auxiliary each field evaluated at the current point
2199 . t            - current time
2200 . x            - coordinates of the current point
2201 . numConstants - number of constant parameters
2202 . constants    - constant parameters
2203 - uNew         - new field values at the current point
2204 
2205   Level: intermediate
2206 
2207 .seealso: `PetscDS`, `PetscDSGetResidual()`
2208 @*/
2209 PetscErrorCode PetscDSSetUpdate(PetscDS ds, PetscInt f, void (*update)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar uNew[]))
2210 {
2211   PetscFunctionBegin;
2212   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
2213   if (update) PetscValidFunction(update, 3);
2214   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2215   PetscCall(PetscDSEnlarge_Static(ds, f + 1));
2216   ds->update[f] = update;
2217   PetscFunctionReturn(PETSC_SUCCESS);
2218 }
2219 
2220 PetscErrorCode PetscDSGetContext(PetscDS ds, PetscInt f, void *ctx)
2221 {
2222   PetscFunctionBegin;
2223   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
2224   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);
2225   PetscAssertPointer(ctx, 3);
2226   *(void **)ctx = ds->ctx[f];
2227   PetscFunctionReturn(PETSC_SUCCESS);
2228 }
2229 
2230 PetscErrorCode PetscDSSetContext(PetscDS ds, PetscInt f, void *ctx)
2231 {
2232   PetscFunctionBegin;
2233   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
2234   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2235   PetscCall(PetscDSEnlarge_Static(ds, f + 1));
2236   ds->ctx[f] = ctx;
2237   PetscFunctionReturn(PETSC_SUCCESS);
2238 }
2239 
2240 /*@C
2241   PetscDSGetBdResidual - Get the pointwise boundary residual function for a given test field
2242 
2243   Not Collective
2244 
2245   Input Parameters:
2246 + ds - The PetscDS
2247 - f  - The test field number
2248 
2249   Output Parameters:
2250 + f0 - boundary integrand for the test function term
2251 - f1 - boundary integrand for the test function gradient term
2252 
2253   Calling sequence of `f0`:
2254 + dim          - the spatial dimension
2255 . Nf           - the number of fields
2256 . NfAux        - the number of auxiliary fields
2257 . uOff         - the offset into u[] and u_t[] for each field
2258 . uOff_x       - the offset into u_x[] for each field
2259 . u            - each field evaluated at the current point
2260 . u_t          - the time derivative of each field evaluated at the current point
2261 . u_x          - the gradient of each field evaluated at the current point
2262 . aOff         - the offset into a[] and a_t[] for each auxiliary field
2263 . aOff_x       - the offset into a_x[] for each auxiliary field
2264 . a            - each auxiliary field evaluated at the current point
2265 . a_t          - the time derivative of each auxiliary field evaluated at the current point
2266 . a_x          - the gradient of auxiliary each field evaluated at the current point
2267 . t            - current time
2268 . x            - coordinates of the current point
2269 . n            - unit normal at the current point
2270 . numConstants - number of constant parameters
2271 . constants    - constant parameters
2272 - f0           - output values at the current point
2273 
2274   Level: intermediate
2275 
2276   Note:
2277   The calling sequence of `f1` is identical, and therefore omitted for brevity.
2278 
2279   We are using a first order FEM model for the weak form\:
2280   \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
2281 
2282 .seealso: `PetscDS`, `PetscDSSetBdResidual()`
2283 @*/
2284 PetscErrorCode PetscDSGetBdResidual(PetscDS ds, PetscInt f, void (**f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]), void (**f1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))
2285 {
2286   PetscBdPointFunc *tmp0, *tmp1;
2287   PetscInt          n0, n1;
2288 
2289   PetscFunctionBegin;
2290   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
2291   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);
2292   PetscCall(PetscWeakFormGetBdResidual(ds->wf, NULL, 0, f, 0, &n0, &tmp0, &n1, &tmp1));
2293   *f0 = tmp0 ? tmp0[0] : NULL;
2294   *f1 = tmp1 ? tmp1[0] : NULL;
2295   PetscFunctionReturn(PETSC_SUCCESS);
2296 }
2297 
2298 /*@C
2299   PetscDSSetBdResidual - Get the pointwise boundary residual function for a given test field
2300 
2301   Not Collective
2302 
2303   Input Parameters:
2304 + ds - The `PetscDS`
2305 . f  - The test field number
2306 . f0 - boundary integrand for the test function term
2307 - f1 - boundary integrand for the test function gradient term
2308 
2309   Calling sequence of `f0`:
2310 + dim          - the spatial dimension
2311 . Nf           - the number of fields
2312 . NfAux        - the number of auxiliary fields
2313 . uOff         - the offset into u[] and u_t[] for each field
2314 . uOff_x       - the offset into u_x[] for each field
2315 . u            - each field evaluated at the current point
2316 . u_t          - the time derivative of each field evaluated at the current point
2317 . u_x          - the gradient of each field evaluated at the current point
2318 . aOff         - the offset into a[] and a_t[] for each auxiliary field
2319 . aOff_x       - the offset into a_x[] for each auxiliary field
2320 . a            - each auxiliary field evaluated at the current point
2321 . a_t          - the time derivative of each auxiliary field evaluated at the current point
2322 . a_x          - the gradient of auxiliary each field evaluated at the current point
2323 . t            - current time
2324 . x            - coordinates of the current point
2325 . n            - unit normal at the current point
2326 . numConstants - number of constant parameters
2327 . constants    - constant parameters
2328 - f0           - output values at the current point
2329 
2330   Level: intermediate
2331 
2332   Note:
2333   The calling sequence of `f1` is identical, and therefore omitted for brevity.
2334 
2335   We are using a first order FEM model for the weak form\:
2336   \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
2337 
2338 .seealso: `PetscDS`, `PetscDSGetBdResidual()`
2339 @*/
2340 PetscErrorCode PetscDSSetBdResidual(PetscDS ds, PetscInt f, void (*f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]), void (*f1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))
2341 {
2342   PetscFunctionBegin;
2343   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
2344   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2345   PetscCall(PetscWeakFormSetIndexBdResidual(ds->wf, NULL, 0, f, 0, 0, f0, 0, f1));
2346   PetscFunctionReturn(PETSC_SUCCESS);
2347 }
2348 
2349 /*@
2350   PetscDSHasBdJacobian - Indicates that boundary Jacobian functions have been set
2351 
2352   Not Collective
2353 
2354   Input Parameter:
2355 . ds - The `PetscDS`
2356 
2357   Output Parameter:
2358 . hasBdJac - flag that pointwise function for the boundary Jacobian has been set
2359 
2360   Level: intermediate
2361 
2362 .seealso: `PetscDS`, `PetscDSHasJacobian()`, `PetscDSSetBdJacobian()`, `PetscDSGetBdJacobian()`
2363 @*/
2364 PetscErrorCode PetscDSHasBdJacobian(PetscDS ds, PetscBool *hasBdJac)
2365 {
2366   PetscFunctionBegin;
2367   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
2368   PetscAssertPointer(hasBdJac, 2);
2369   PetscCall(PetscWeakFormHasBdJacobian(ds->wf, hasBdJac));
2370   PetscFunctionReturn(PETSC_SUCCESS);
2371 }
2372 
2373 /*@C
2374   PetscDSGetBdJacobian - Get the pointwise boundary Jacobian function for given test and basis field
2375 
2376   Not Collective
2377 
2378   Input Parameters:
2379 + ds - The `PetscDS`
2380 . f  - The test field number
2381 - g  - The field number
2382 
2383   Output Parameters:
2384 + g0 - integrand for the test and basis function term
2385 . g1 - integrand for the test function and basis function gradient term
2386 . g2 - integrand for the test function gradient and basis function term
2387 - g3 - integrand for the test function gradient and basis function gradient term
2388 
2389   Calling sequence of `g0`:
2390 + dim          - the spatial dimension
2391 . Nf           - the number of fields
2392 . NfAux        - the number of auxiliary fields
2393 . uOff         - the offset into u[] and u_t[] for each field
2394 . uOff_x       - the offset into u_x[] for each field
2395 . u            - each field evaluated at the current point
2396 . u_t          - the time derivative of each field evaluated at the current point
2397 . u_x          - the gradient of each field evaluated at the current point
2398 . aOff         - the offset into a[] and a_t[] for each auxiliary field
2399 . aOff_x       - the offset into a_x[] for each auxiliary field
2400 . a            - each auxiliary field evaluated at the current point
2401 . a_t          - the time derivative of each auxiliary field evaluated at the current point
2402 . a_x          - the gradient of auxiliary each field evaluated at the current point
2403 . t            - current time
2404 . u_tShift     - the multiplier a for dF/dU_t
2405 . x            - coordinates of the current point
2406 . n            - normal at the current point
2407 . numConstants - number of constant parameters
2408 . constants    - constant parameters
2409 - g0           - output values at the current point
2410 
2411   Level: intermediate
2412 
2413   Note:
2414   `g1`, `g2`, and `g3` have identical calling sequences to `g0` and are omitted for brevity.
2415 
2416   We are using a first order FEM model for the weak form\:
2417   \int_\Gamma \phi {\vec g}_0(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \cdot \hat n \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \hat n \cdot \nabla \psi
2418 
2419 .seealso: `PetscDS`, `PetscDSSetBdJacobian()`
2420 @*/
2421 PetscErrorCode PetscDSGetBdJacobian(PetscDS ds, PetscInt f, PetscInt g, void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]), void (**g1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (**g2)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (**g3)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))
2422 {
2423   PetscBdPointJac *tmp0, *tmp1, *tmp2, *tmp3;
2424   PetscInt         n0, n1, n2, n3;
2425 
2426   PetscFunctionBegin;
2427   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
2428   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);
2429   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);
2430   PetscCall(PetscWeakFormGetBdJacobian(ds->wf, NULL, 0, f, g, 0, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3));
2431   *g0 = tmp0 ? tmp0[0] : NULL;
2432   *g1 = tmp1 ? tmp1[0] : NULL;
2433   *g2 = tmp2 ? tmp2[0] : NULL;
2434   *g3 = tmp3 ? tmp3[0] : NULL;
2435   PetscFunctionReturn(PETSC_SUCCESS);
2436 }
2437 
2438 /*@C
2439   PetscDSSetBdJacobian - Set the pointwise boundary Jacobian function for given test and basis field
2440 
2441   Not Collective
2442 
2443   Input Parameters:
2444 + ds - The PetscDS
2445 . f  - The test field number
2446 . g  - The field number
2447 . g0 - integrand for the test and basis function term
2448 . g1 - integrand for the test function and basis function gradient term
2449 . g2 - integrand for the test function gradient and basis function term
2450 - g3 - integrand for the test function gradient and basis function gradient term
2451 
2452   Calling sequence of `g0`:
2453 + dim          - the spatial dimension
2454 . Nf           - the number of fields
2455 . NfAux        - the number of auxiliary fields
2456 . uOff         - the offset into u[] and u_t[] for each field
2457 . uOff_x       - the offset into u_x[] for each field
2458 . u            - each field evaluated at the current point
2459 . u_t          - the time derivative of each field evaluated at the current point
2460 . u_x          - the gradient of each field evaluated at the current point
2461 . aOff         - the offset into a[] and a_t[] for each auxiliary field
2462 . aOff_x       - the offset into a_x[] for each auxiliary field
2463 . a            - each auxiliary field evaluated at the current point
2464 . a_t          - the time derivative of each auxiliary field evaluated at the current point
2465 . a_x          - the gradient of auxiliary each field evaluated at the current point
2466 . t            - current time
2467 . u_tShift     - the multiplier a for dF/dU_t
2468 . x            - coordinates of the current point
2469 . n            - normal at the current point
2470 . numConstants - number of constant parameters
2471 . constants    - constant parameters
2472 - g0           - output values at the current point
2473 
2474   Level: intermediate
2475 
2476   Note:
2477   `g1`, `g2`, and `g3` have identical calling sequences to `g0` and are omitted for brevity.
2478 
2479   We are using a first order FEM model for the weak form\:
2480   \int_\Gamma \phi {\vec g}_0(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \cdot \hat n \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \hat n \cdot \nabla \psi
2481 
2482 .seealso: `PetscDS`, `PetscDSGetBdJacobian()`
2483 @*/
2484 PetscErrorCode PetscDSSetBdJacobian(PetscDS ds, PetscInt f, PetscInt g, void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]), void (*g1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (*g2)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (*g3)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))
2485 {
2486   PetscFunctionBegin;
2487   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
2488   if (g0) PetscValidFunction(g0, 4);
2489   if (g1) PetscValidFunction(g1, 5);
2490   if (g2) PetscValidFunction(g2, 6);
2491   if (g3) PetscValidFunction(g3, 7);
2492   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2493   PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", g);
2494   PetscCall(PetscWeakFormSetIndexBdJacobian(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3));
2495   PetscFunctionReturn(PETSC_SUCCESS);
2496 }
2497 
2498 /*@
2499   PetscDSHasBdJacobianPreconditioner - Signals that boundary Jacobian preconditioner functions have been set
2500 
2501   Not Collective
2502 
2503   Input Parameter:
2504 . ds - The `PetscDS`
2505 
2506   Output Parameter:
2507 . hasBdJacPre - flag that pointwise function for the boundary Jacobian preconditioner has been set
2508 
2509   Level: intermediate
2510 
2511 .seealso: `PetscDS`, `PetscDSHasJacobian()`, `PetscDSSetBdJacobian()`, `PetscDSGetBdJacobian()`
2512 @*/
2513 PetscErrorCode PetscDSHasBdJacobianPreconditioner(PetscDS ds, PetscBool *hasBdJacPre)
2514 {
2515   PetscFunctionBegin;
2516   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
2517   PetscAssertPointer(hasBdJacPre, 2);
2518   PetscCall(PetscWeakFormHasBdJacobianPreconditioner(ds->wf, hasBdJacPre));
2519   PetscFunctionReturn(PETSC_SUCCESS);
2520 }
2521 
2522 /*@C
2523   PetscDSGetBdJacobianPreconditioner - Get the pointwise boundary Jacobian preconditioner function for given test and basis field
2524 
2525   Not Collective; No Fortran Support
2526 
2527   Input Parameters:
2528 + ds - The `PetscDS`
2529 . f  - The test field number
2530 - g  - The field number
2531 
2532   Output Parameters:
2533 + g0 - integrand for the test and basis function term
2534 . g1 - integrand for the test function and basis function gradient term
2535 . g2 - integrand for the test function gradient and basis function term
2536 - g3 - integrand for the test function gradient and basis function gradient term
2537 
2538   Calling sequence of `g0`:
2539 + dim          - the spatial dimension
2540 . Nf           - the number of fields
2541 . NfAux        - the number of auxiliary fields
2542 . uOff         - the offset into u[] and u_t[] for each field
2543 . uOff_x       - the offset into u_x[] for each field
2544 . u            - each field evaluated at the current point
2545 . u_t          - the time derivative of each field evaluated at the current point
2546 . u_x          - the gradient of each field evaluated at the current point
2547 . aOff         - the offset into a[] and a_t[] for each auxiliary field
2548 . aOff_x       - the offset into a_x[] for each auxiliary field
2549 . a            - each auxiliary field evaluated at the current point
2550 . a_t          - the time derivative of each auxiliary field evaluated at the current point
2551 . a_x          - the gradient of auxiliary each field evaluated at the current point
2552 . t            - current time
2553 . u_tShift     - the multiplier a for dF/dU_t
2554 . x            - coordinates of the current point
2555 . n            - normal at the current point
2556 . numConstants - number of constant parameters
2557 . constants    - constant parameters
2558 - g0           - output values at the current point
2559 
2560   Level: intermediate
2561 
2562   Note:
2563   `g1`, `g2`, and `g3` have identical calling sequences to `g0` and are omitted for brevity.
2564 
2565   We are using a first order FEM model for the weak form\:
2566   \int_\Gamma \phi {\vec g}_0(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \cdot \hat n \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \hat n \cdot \nabla \psi
2567 
2568 .seealso: `PetscDS`, `PetscDSSetBdJacobianPreconditioner()`
2569 @*/
2570 PetscErrorCode PetscDSGetBdJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g, void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]), void (**g1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (**g2)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (**g3)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))
2571 {
2572   PetscBdPointJac *tmp0, *tmp1, *tmp2, *tmp3;
2573   PetscInt         n0, n1, n2, n3;
2574 
2575   PetscFunctionBegin;
2576   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
2577   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);
2578   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);
2579   PetscCall(PetscWeakFormGetBdJacobianPreconditioner(ds->wf, NULL, 0, f, g, 0, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3));
2580   *g0 = tmp0 ? tmp0[0] : NULL;
2581   *g1 = tmp1 ? tmp1[0] : NULL;
2582   *g2 = tmp2 ? tmp2[0] : NULL;
2583   *g3 = tmp3 ? tmp3[0] : NULL;
2584   PetscFunctionReturn(PETSC_SUCCESS);
2585 }
2586 
2587 /*@C
2588   PetscDSSetBdJacobianPreconditioner - Set the pointwise boundary Jacobian preconditioner function for given test and basis field
2589 
2590   Not Collective; No Fortran Support
2591 
2592   Input Parameters:
2593 + ds - The `PetscDS`
2594 . f  - The test field number
2595 . g  - The field number
2596 . g0 - integrand for the test and basis function term
2597 . g1 - integrand for the test function and basis function gradient term
2598 . g2 - integrand for the test function gradient and basis function term
2599 - g3 - integrand for the test function gradient and basis function gradient term
2600 
2601   Calling sequence of `g0':
2602 + dim          - the spatial dimension
2603 . Nf           - the number of fields
2604 . NfAux        - the number of auxiliary fields
2605 . uOff         - the offset into u[] and u_t[] for each field
2606 . uOff_x       - the offset into u_x[] for each field
2607 . u            - each field evaluated at the current point
2608 . u_t          - the time derivative of each field evaluated at the current point
2609 . u_x          - the gradient of each field evaluated at the current point
2610 . aOff         - the offset into a[] and a_t[] for each auxiliary field
2611 . aOff_x       - the offset into a_x[] for each auxiliary field
2612 . a            - each auxiliary field evaluated at the current point
2613 . a_t          - the time derivative of each auxiliary field evaluated at the current point
2614 . a_x          - the gradient of auxiliary each field evaluated at the current point
2615 . t            - current time
2616 . u_tShift     - the multiplier a for dF/dU_t
2617 . x            - coordinates of the current point
2618 . n            - normal at the current point
2619 . numConstants - number of constant parameters
2620 . constants    - constant parameters
2621 - g0           - output values at the current point
2622 
2623   Level: intermediate
2624 
2625   Note:
2626   `g1`, `g2`, and `g3` have identical calling sequences to `g0` and are omitted for brevity.
2627 
2628   We are using a first order FEM model for the weak form\:
2629   \int_\Gamma \phi {\vec g}_0(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \cdot \hat n \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \hat n \cdot \nabla \psi
2630 
2631 .seealso: `PetscDS`, `PetscDSGetBdJacobianPreconditioner()`
2632 @*/
2633 PetscErrorCode PetscDSSetBdJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g, void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]), void (*g1)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (*g2)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), void (*g3)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))
2634 {
2635   PetscFunctionBegin;
2636   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
2637   if (g0) PetscValidFunction(g0, 4);
2638   if (g1) PetscValidFunction(g1, 5);
2639   if (g2) PetscValidFunction(g2, 6);
2640   if (g3) PetscValidFunction(g3, 7);
2641   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2642   PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", g);
2643   PetscCall(PetscWeakFormSetIndexBdJacobianPreconditioner(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3));
2644   PetscFunctionReturn(PETSC_SUCCESS);
2645 }
2646 
2647 /*@C
2648   PetscDSGetExactSolution - Get the pointwise exact solution function for a given test field
2649 
2650   Not Collective
2651 
2652   Input Parameters:
2653 + prob - The PetscDS
2654 - f    - The test field number
2655 
2656   Output Parameters:
2657 + sol - exact solution for the test field
2658 - ctx - exact solution context
2659 
2660   Calling sequence of `exactSol`:
2661 + dim - the spatial dimension
2662 . t   - current time
2663 . x   - coordinates of the current point
2664 . Nc  - the number of field components
2665 . u   - the solution field evaluated at the current point
2666 - ctx - a user context
2667 
2668   Level: intermediate
2669 
2670 .seealso: `PetscDS`, `PetscDSSetExactSolution()`, `PetscDSGetExactSolutionTimeDerivative()`
2671 @*/
2672 PetscErrorCode PetscDSGetExactSolution(PetscDS prob, PetscInt f, PetscErrorCode (**sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void **ctx)
2673 {
2674   PetscFunctionBegin;
2675   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
2676   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);
2677   if (sol) {
2678     PetscAssertPointer(sol, 3);
2679     *sol = prob->exactSol[f];
2680   }
2681   if (ctx) {
2682     PetscAssertPointer(ctx, 4);
2683     *ctx = prob->exactCtx[f];
2684   }
2685   PetscFunctionReturn(PETSC_SUCCESS);
2686 }
2687 
2688 /*@C
2689   PetscDSSetExactSolution - Set the pointwise exact solution function for a given test field
2690 
2691   Not Collective
2692 
2693   Input Parameters:
2694 + prob - The `PetscDS`
2695 . f    - The test field number
2696 . sol  - solution function for the test fields
2697 - ctx  - solution context or `NULL`
2698 
2699   Calling sequence of `sol`:
2700 + dim - the spatial dimension
2701 . t   - current time
2702 . x   - coordinates of the current point
2703 . Nc  - the number of field components
2704 . u   - the solution field evaluated at the current point
2705 - ctx - a user context
2706 
2707   Level: intermediate
2708 
2709 .seealso: `PetscDS`, `PetscDSGetExactSolution()`
2710 @*/
2711 PetscErrorCode PetscDSSetExactSolution(PetscDS prob, PetscInt f, PetscErrorCode (*sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void *ctx)
2712 {
2713   PetscFunctionBegin;
2714   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
2715   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2716   PetscCall(PetscDSEnlarge_Static(prob, f + 1));
2717   if (sol) {
2718     PetscValidFunction(sol, 3);
2719     prob->exactSol[f] = sol;
2720   }
2721   if (ctx) {
2722     PetscValidFunction(ctx, 4);
2723     prob->exactCtx[f] = ctx;
2724   }
2725   PetscFunctionReturn(PETSC_SUCCESS);
2726 }
2727 
2728 /*@C
2729   PetscDSGetExactSolutionTimeDerivative - Get the pointwise time derivative of the exact solution function for a given test field
2730 
2731   Not Collective
2732 
2733   Input Parameters:
2734 + prob - The `PetscDS`
2735 - f    - The test field number
2736 
2737   Output Parameters:
2738 + sol - time derivative of the exact solution for the test field
2739 - ctx - time derivative of the exact solution context
2740 
2741   Calling sequence of `exactSol`:
2742 + dim - the spatial dimension
2743 . t   - current time
2744 . x   - coordinates of the current point
2745 . Nc  - the number of field components
2746 . u   - the solution field evaluated at the current point
2747 - ctx - a user context
2748 
2749   Level: intermediate
2750 
2751 .seealso: `PetscDS`, `PetscDSSetExactSolutionTimeDerivative()`, `PetscDSGetExactSolution()`
2752 @*/
2753 PetscErrorCode PetscDSGetExactSolutionTimeDerivative(PetscDS prob, PetscInt f, PetscErrorCode (**sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void **ctx)
2754 {
2755   PetscFunctionBegin;
2756   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
2757   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);
2758   if (sol) {
2759     PetscAssertPointer(sol, 3);
2760     *sol = prob->exactSol_t[f];
2761   }
2762   if (ctx) {
2763     PetscAssertPointer(ctx, 4);
2764     *ctx = prob->exactCtx_t[f];
2765   }
2766   PetscFunctionReturn(PETSC_SUCCESS);
2767 }
2768 
2769 /*@C
2770   PetscDSSetExactSolutionTimeDerivative - Set the pointwise time derivative of the exact solution function for a given test field
2771 
2772   Not Collective
2773 
2774   Input Parameters:
2775 + prob - The `PetscDS`
2776 . f    - The test field number
2777 . sol  - time derivative of the solution function for the test fields
2778 - ctx  - time derivative of the solution context or `NULL`
2779 
2780   Calling sequence of `sol`:
2781 + dim - the spatial dimension
2782 . t   - current time
2783 . x   - coordinates of the current point
2784 . Nc  - the number of field components
2785 . u   - the solution field evaluated at the current point
2786 - ctx - a user context
2787 
2788   Level: intermediate
2789 
2790 .seealso: `PetscDS`, `PetscDSGetExactSolutionTimeDerivative()`, `PetscDSSetExactSolution()`
2791 @*/
2792 PetscErrorCode PetscDSSetExactSolutionTimeDerivative(PetscDS prob, PetscInt f, PetscErrorCode (*sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void *ctx)
2793 {
2794   PetscFunctionBegin;
2795   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
2796   PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2797   PetscCall(PetscDSEnlarge_Static(prob, f + 1));
2798   if (sol) {
2799     PetscValidFunction(sol, 3);
2800     prob->exactSol_t[f] = sol;
2801   }
2802   if (ctx) {
2803     PetscValidFunction(ctx, 4);
2804     prob->exactCtx_t[f] = ctx;
2805   }
2806   PetscFunctionReturn(PETSC_SUCCESS);
2807 }
2808 
2809 /*@C
2810   PetscDSGetConstants - Returns the array of constants passed to point functions
2811 
2812   Not Collective
2813 
2814   Input Parameter:
2815 . prob - The `PetscDS` object
2816 
2817   Output Parameters:
2818 + numConstants - The number of constants
2819 - constants    - The array of constants, NULL if there are none
2820 
2821   Level: intermediate
2822 
2823 .seealso: `PetscDS`, `PetscDSSetConstants()`, `PetscDSCreate()`
2824 @*/
2825 PetscErrorCode PetscDSGetConstants(PetscDS prob, PetscInt *numConstants, const PetscScalar *constants[])
2826 {
2827   PetscFunctionBegin;
2828   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
2829   if (numConstants) {
2830     PetscAssertPointer(numConstants, 2);
2831     *numConstants = prob->numConstants;
2832   }
2833   if (constants) {
2834     PetscAssertPointer(constants, 3);
2835     *constants = prob->constants;
2836   }
2837   PetscFunctionReturn(PETSC_SUCCESS);
2838 }
2839 
2840 /*@C
2841   PetscDSSetConstants - Set the array of constants passed to point functions
2842 
2843   Not Collective
2844 
2845   Input Parameters:
2846 + prob         - The `PetscDS` object
2847 . numConstants - The number of constants
2848 - constants    - The array of constants, NULL if there are none
2849 
2850   Level: intermediate
2851 
2852 .seealso: `PetscDS`, `PetscDSGetConstants()`, `PetscDSCreate()`
2853 @*/
2854 PetscErrorCode PetscDSSetConstants(PetscDS prob, PetscInt numConstants, PetscScalar constants[])
2855 {
2856   PetscFunctionBegin;
2857   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
2858   if (numConstants != prob->numConstants) {
2859     PetscCall(PetscFree(prob->constants));
2860     prob->numConstants = numConstants;
2861     if (prob->numConstants) {
2862       PetscCall(PetscMalloc1(prob->numConstants, &prob->constants));
2863     } else {
2864       prob->constants = NULL;
2865     }
2866   }
2867   if (prob->numConstants) {
2868     PetscAssertPointer(constants, 3);
2869     PetscCall(PetscArraycpy(prob->constants, constants, prob->numConstants));
2870   }
2871   PetscFunctionReturn(PETSC_SUCCESS);
2872 }
2873 
2874 /*@
2875   PetscDSGetFieldIndex - Returns the index of the given field
2876 
2877   Not Collective
2878 
2879   Input Parameters:
2880 + prob - The `PetscDS` object
2881 - disc - The discretization object
2882 
2883   Output Parameter:
2884 . f - The field number
2885 
2886   Level: beginner
2887 
2888 .seealso: `PetscDS`, `PetscGetDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
2889 @*/
2890 PetscErrorCode PetscDSGetFieldIndex(PetscDS prob, PetscObject disc, PetscInt *f)
2891 {
2892   PetscInt g;
2893 
2894   PetscFunctionBegin;
2895   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
2896   PetscAssertPointer(f, 3);
2897   *f = -1;
2898   for (g = 0; g < prob->Nf; ++g) {
2899     if (disc == prob->disc[g]) break;
2900   }
2901   PetscCheck(g != prob->Nf, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONG, "Field not found in PetscDS.");
2902   *f = g;
2903   PetscFunctionReturn(PETSC_SUCCESS);
2904 }
2905 
2906 /*@
2907   PetscDSGetFieldSize - Returns the size of the given field in the full space basis
2908 
2909   Not Collective
2910 
2911   Input Parameters:
2912 + prob - The `PetscDS` object
2913 - f    - The field number
2914 
2915   Output Parameter:
2916 . size - The size
2917 
2918   Level: beginner
2919 
2920 .seealso: `PetscDS`, `PetscDSGetFieldOffset()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
2921 @*/
2922 PetscErrorCode PetscDSGetFieldSize(PetscDS prob, PetscInt f, PetscInt *size)
2923 {
2924   PetscFunctionBegin;
2925   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
2926   PetscAssertPointer(size, 3);
2927   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);
2928   PetscCall(PetscDSSetUp(prob));
2929   *size = prob->Nb[f];
2930   PetscFunctionReturn(PETSC_SUCCESS);
2931 }
2932 
2933 /*@
2934   PetscDSGetFieldOffset - Returns the offset of the given field in the full space basis
2935 
2936   Not Collective
2937 
2938   Input Parameters:
2939 + prob - The `PetscDS` object
2940 - f    - The field number
2941 
2942   Output Parameter:
2943 . off - The offset
2944 
2945   Level: beginner
2946 
2947 .seealso: `PetscDS`, `PetscDSGetFieldSize()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
2948 @*/
2949 PetscErrorCode PetscDSGetFieldOffset(PetscDS prob, PetscInt f, PetscInt *off)
2950 {
2951   PetscInt size, g;
2952 
2953   PetscFunctionBegin;
2954   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
2955   PetscAssertPointer(off, 3);
2956   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);
2957   *off = 0;
2958   for (g = 0; g < f; ++g) {
2959     PetscCall(PetscDSGetFieldSize(prob, g, &size));
2960     *off += size;
2961   }
2962   PetscFunctionReturn(PETSC_SUCCESS);
2963 }
2964 
2965 /*@
2966   PetscDSGetFieldOffsetCohesive - Returns the offset of the given field in the full space basis on a cohesive cell
2967 
2968   Not Collective
2969 
2970   Input Parameters:
2971 + ds - The `PetscDS` object
2972 - f  - The field number
2973 
2974   Output Parameter:
2975 . off - The offset
2976 
2977   Level: beginner
2978 
2979 .seealso: `PetscDS`, `PetscDSGetFieldSize()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
2980 @*/
2981 PetscErrorCode PetscDSGetFieldOffsetCohesive(PetscDS ds, PetscInt f, PetscInt *off)
2982 {
2983   PetscInt size, g;
2984 
2985   PetscFunctionBegin;
2986   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
2987   PetscAssertPointer(off, 3);
2988   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);
2989   *off = 0;
2990   for (g = 0; g < f; ++g) {
2991     PetscBool cohesive;
2992 
2993     PetscCall(PetscDSGetCohesive(ds, g, &cohesive));
2994     PetscCall(PetscDSGetFieldSize(ds, g, &size));
2995     *off += cohesive ? size : size * 2;
2996   }
2997   PetscFunctionReturn(PETSC_SUCCESS);
2998 }
2999 
3000 /*@
3001   PetscDSGetDimensions - Returns the size of the approximation space for each field on an evaluation point
3002 
3003   Not Collective
3004 
3005   Input Parameter:
3006 . prob - The `PetscDS` object
3007 
3008   Output Parameter:
3009 . dimensions - The number of dimensions
3010 
3011   Level: beginner
3012 
3013 .seealso: `PetscDS`, `PetscDSGetComponentOffsets()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3014 @*/
3015 PetscErrorCode PetscDSGetDimensions(PetscDS prob, PetscInt *dimensions[])
3016 {
3017   PetscFunctionBegin;
3018   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
3019   PetscCall(PetscDSSetUp(prob));
3020   PetscAssertPointer(dimensions, 2);
3021   *dimensions = prob->Nb;
3022   PetscFunctionReturn(PETSC_SUCCESS);
3023 }
3024 
3025 /*@
3026   PetscDSGetComponents - Returns the number of components for each field on an evaluation point
3027 
3028   Not Collective
3029 
3030   Input Parameter:
3031 . prob - The `PetscDS` object
3032 
3033   Output Parameter:
3034 . components - The number of components
3035 
3036   Level: beginner
3037 
3038 .seealso: `PetscDS`, `PetscDSGetComponentOffsets()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3039 @*/
3040 PetscErrorCode PetscDSGetComponents(PetscDS prob, PetscInt *components[])
3041 {
3042   PetscFunctionBegin;
3043   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
3044   PetscCall(PetscDSSetUp(prob));
3045   PetscAssertPointer(components, 2);
3046   *components = prob->Nc;
3047   PetscFunctionReturn(PETSC_SUCCESS);
3048 }
3049 
3050 /*@
3051   PetscDSGetComponentOffset - Returns the offset of the given field on an evaluation point
3052 
3053   Not Collective
3054 
3055   Input Parameters:
3056 + prob - The `PetscDS` object
3057 - f    - The field number
3058 
3059   Output Parameter:
3060 . off - The offset
3061 
3062   Level: beginner
3063 
3064 .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3065 @*/
3066 PetscErrorCode PetscDSGetComponentOffset(PetscDS prob, PetscInt f, PetscInt *off)
3067 {
3068   PetscFunctionBegin;
3069   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
3070   PetscAssertPointer(off, 3);
3071   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);
3072   PetscCall(PetscDSSetUp(prob));
3073   *off = prob->off[f];
3074   PetscFunctionReturn(PETSC_SUCCESS);
3075 }
3076 
3077 /*@
3078   PetscDSGetComponentOffsets - Returns the offset of each field on an evaluation point
3079 
3080   Not Collective
3081 
3082   Input Parameter:
3083 . prob - The `PetscDS` object
3084 
3085   Output Parameter:
3086 . offsets - The offsets
3087 
3088   Level: beginner
3089 
3090 .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3091 @*/
3092 PetscErrorCode PetscDSGetComponentOffsets(PetscDS prob, PetscInt *offsets[])
3093 {
3094   PetscFunctionBegin;
3095   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
3096   PetscAssertPointer(offsets, 2);
3097   PetscCall(PetscDSSetUp(prob));
3098   *offsets = prob->off;
3099   PetscFunctionReturn(PETSC_SUCCESS);
3100 }
3101 
3102 /*@
3103   PetscDSGetComponentDerivativeOffsets - Returns the offset of each field derivative on an evaluation point
3104 
3105   Not Collective
3106 
3107   Input Parameter:
3108 . prob - The `PetscDS` object
3109 
3110   Output Parameter:
3111 . offsets - The offsets
3112 
3113   Level: beginner
3114 
3115 .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3116 @*/
3117 PetscErrorCode PetscDSGetComponentDerivativeOffsets(PetscDS prob, PetscInt *offsets[])
3118 {
3119   PetscFunctionBegin;
3120   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
3121   PetscAssertPointer(offsets, 2);
3122   PetscCall(PetscDSSetUp(prob));
3123   *offsets = prob->offDer;
3124   PetscFunctionReturn(PETSC_SUCCESS);
3125 }
3126 
3127 /*@
3128   PetscDSGetComponentOffsetsCohesive - Returns the offset of each field on an evaluation point
3129 
3130   Not Collective
3131 
3132   Input Parameters:
3133 + ds - The `PetscDS` object
3134 - s  - The cohesive side, 0 for negative, 1 for positive, 2 for cohesive
3135 
3136   Output Parameter:
3137 . offsets - The offsets
3138 
3139   Level: beginner
3140 
3141 .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3142 @*/
3143 PetscErrorCode PetscDSGetComponentOffsetsCohesive(PetscDS ds, PetscInt s, PetscInt *offsets[])
3144 {
3145   PetscFunctionBegin;
3146   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
3147   PetscAssertPointer(offsets, 3);
3148   PetscCheck(ds->isCohesive, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cohesive offsets are only valid for a cohesive DS");
3149   PetscCheck(!(s < 0) && !(s > 2), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cohesive side %" PetscInt_FMT " is not in [0, 2]", s);
3150   PetscCall(PetscDSSetUp(ds));
3151   *offsets = ds->offCohesive[s];
3152   PetscFunctionReturn(PETSC_SUCCESS);
3153 }
3154 
3155 /*@
3156   PetscDSGetComponentDerivativeOffsetsCohesive - Returns the offset of each field derivative on an evaluation point
3157 
3158   Not Collective
3159 
3160   Input Parameters:
3161 + ds - The `PetscDS` object
3162 - s  - The cohesive side, 0 for negative, 1 for positive, 2 for cohesive
3163 
3164   Output Parameter:
3165 . offsets - The offsets
3166 
3167   Level: beginner
3168 
3169 .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3170 @*/
3171 PetscErrorCode PetscDSGetComponentDerivativeOffsetsCohesive(PetscDS ds, PetscInt s, PetscInt *offsets[])
3172 {
3173   PetscFunctionBegin;
3174   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
3175   PetscAssertPointer(offsets, 3);
3176   PetscCheck(ds->isCohesive, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cohesive offsets are only valid for a cohesive DS");
3177   PetscCheck(!(s < 0) && !(s > 2), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cohesive side %" PetscInt_FMT " is not in [0, 2]", s);
3178   PetscCall(PetscDSSetUp(ds));
3179   *offsets = ds->offDerCohesive[s];
3180   PetscFunctionReturn(PETSC_SUCCESS);
3181 }
3182 
3183 /*@C
3184   PetscDSGetTabulation - Return the basis tabulation at quadrature points for the volume discretization
3185 
3186   Not Collective
3187 
3188   Input Parameter:
3189 . prob - The `PetscDS` object
3190 
3191   Output Parameter:
3192 . T - The basis function and derivatives tabulation at quadrature points for each field
3193 
3194   Level: intermediate
3195 
3196 .seealso: `PetscDS`, `PetscTabulation`, `PetscDSCreate()`
3197 @*/
3198 PetscErrorCode PetscDSGetTabulation(PetscDS prob, PetscTabulation *T[])
3199 {
3200   PetscFunctionBegin;
3201   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
3202   PetscAssertPointer(T, 2);
3203   PetscCall(PetscDSSetUp(prob));
3204   *T = prob->T;
3205   PetscFunctionReturn(PETSC_SUCCESS);
3206 }
3207 
3208 /*@C
3209   PetscDSGetFaceTabulation - Return the basis tabulation at quadrature points on the faces
3210 
3211   Not Collective
3212 
3213   Input Parameter:
3214 . prob - The `PetscDS` object
3215 
3216   Output Parameter:
3217 . Tf - The basis function and derivative tabulation on each local face at quadrature points for each and field
3218 
3219   Level: intermediate
3220 
3221 .seealso: `PetscTabulation`, `PetscDS`, `PetscDSGetTabulation()`, `PetscDSCreate()`
3222 @*/
3223 PetscErrorCode PetscDSGetFaceTabulation(PetscDS prob, PetscTabulation *Tf[])
3224 {
3225   PetscFunctionBegin;
3226   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
3227   PetscAssertPointer(Tf, 2);
3228   PetscCall(PetscDSSetUp(prob));
3229   *Tf = prob->Tf;
3230   PetscFunctionReturn(PETSC_SUCCESS);
3231 }
3232 
3233 PetscErrorCode PetscDSGetEvaluationArrays(PetscDS prob, PetscScalar **u, PetscScalar **u_t, PetscScalar **u_x)
3234 {
3235   PetscFunctionBegin;
3236   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
3237   PetscCall(PetscDSSetUp(prob));
3238   if (u) {
3239     PetscAssertPointer(u, 2);
3240     *u = prob->u;
3241   }
3242   if (u_t) {
3243     PetscAssertPointer(u_t, 3);
3244     *u_t = prob->u_t;
3245   }
3246   if (u_x) {
3247     PetscAssertPointer(u_x, 4);
3248     *u_x = prob->u_x;
3249   }
3250   PetscFunctionReturn(PETSC_SUCCESS);
3251 }
3252 
3253 PetscErrorCode PetscDSGetWeakFormArrays(PetscDS prob, PetscScalar **f0, PetscScalar **f1, PetscScalar **g0, PetscScalar **g1, PetscScalar **g2, PetscScalar **g3)
3254 {
3255   PetscFunctionBegin;
3256   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
3257   PetscCall(PetscDSSetUp(prob));
3258   if (f0) {
3259     PetscAssertPointer(f0, 2);
3260     *f0 = prob->f0;
3261   }
3262   if (f1) {
3263     PetscAssertPointer(f1, 3);
3264     *f1 = prob->f1;
3265   }
3266   if (g0) {
3267     PetscAssertPointer(g0, 4);
3268     *g0 = prob->g0;
3269   }
3270   if (g1) {
3271     PetscAssertPointer(g1, 5);
3272     *g1 = prob->g1;
3273   }
3274   if (g2) {
3275     PetscAssertPointer(g2, 6);
3276     *g2 = prob->g2;
3277   }
3278   if (g3) {
3279     PetscAssertPointer(g3, 7);
3280     *g3 = prob->g3;
3281   }
3282   PetscFunctionReturn(PETSC_SUCCESS);
3283 }
3284 
3285 PetscErrorCode PetscDSGetWorkspace(PetscDS prob, PetscReal **x, PetscScalar **basisReal, PetscScalar **basisDerReal, PetscScalar **testReal, PetscScalar **testDerReal)
3286 {
3287   PetscFunctionBegin;
3288   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
3289   PetscCall(PetscDSSetUp(prob));
3290   if (x) {
3291     PetscAssertPointer(x, 2);
3292     *x = prob->x;
3293   }
3294   if (basisReal) {
3295     PetscAssertPointer(basisReal, 3);
3296     *basisReal = prob->basisReal;
3297   }
3298   if (basisDerReal) {
3299     PetscAssertPointer(basisDerReal, 4);
3300     *basisDerReal = prob->basisDerReal;
3301   }
3302   if (testReal) {
3303     PetscAssertPointer(testReal, 5);
3304     *testReal = prob->testReal;
3305   }
3306   if (testDerReal) {
3307     PetscAssertPointer(testDerReal, 6);
3308     *testDerReal = prob->testDerReal;
3309   }
3310   PetscFunctionReturn(PETSC_SUCCESS);
3311 }
3312 
3313 /*@C
3314   PetscDSAddBoundary - Add a boundary condition to the model.
3315 
3316   Collective
3317 
3318   Input Parameters:
3319 + ds       - The PetscDS object
3320 . type     - The type of condition, e.g. `DM_BC_ESSENTIAL`/`DM_BC_ESSENTIAL_FIELD` (Dirichlet), or `DM_BC_NATURAL` (Neumann)
3321 . name     - The BC name
3322 . label    - The label defining constrained points
3323 . Nv       - The number of `DMLabel` values for constrained points
3324 . values   - An array of label values for constrained points
3325 . field    - The field to constrain
3326 . Nc       - The number of constrained field components (0 will constrain all fields)
3327 . comps    - An array of constrained component numbers
3328 . bcFunc   - A pointwise function giving boundary values
3329 . bcFunc_t - A pointwise function giving the time derivative of the boundary values, or NULL
3330 - ctx      - An optional user context for bcFunc
3331 
3332   Output Parameter:
3333 . bd - The boundary number
3334 
3335   Options Database Keys:
3336 + -bc_<boundary name> <num>      - Overrides the boundary ids
3337 - -bc_<boundary name>_comp <num> - Overrides the boundary components
3338 
3339   Level: developer
3340 
3341   Note:
3342   Both `bcFunc` and `bcFunc_t` will depend on the boundary condition type. If the type if `DM_BC_ESSENTIAL`, then the calling sequence is\:
3343 
3344 $ void bcFunc(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar bcval[])
3345 
3346   If the type is `DM_BC_ESSENTIAL_FIELD` or other _FIELD value, then the calling sequence is\:
3347 .vb
3348   void bcFunc(PetscInt dim, PetscInt Nf, PetscInt NfAux,
3349               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
3350               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
3351               PetscReal time, const PetscReal x[], PetscScalar bcval[])
3352 .ve
3353 + dim - the spatial dimension
3354 . Nf - the number of fields
3355 . uOff - the offset into u[] and u_t[] for each field
3356 . uOff_x - the offset into u_x[] for each field
3357 . u - each field evaluated at the current point
3358 . u_t - the time derivative of each field evaluated at the current point
3359 . u_x - the gradient of each field evaluated at the current point
3360 . aOff - the offset into a[] and a_t[] for each auxiliary field
3361 . aOff_x - the offset into a_x[] for each auxiliary field
3362 . a - each auxiliary field evaluated at the current point
3363 . a_t - the time derivative of each auxiliary field evaluated at the current point
3364 . a_x - the gradient of auxiliary each field evaluated at the current point
3365 . t - current time
3366 . x - coordinates of the current point
3367 . numConstants - number of constant parameters
3368 . constants - constant parameters
3369 - bcval - output values at the current point
3370 
3371   Notes:
3372   The pointwise functions are used to provide boundary values for essential boundary
3373   conditions. In FEM, they are acting upon by dual basis functionals to generate FEM
3374   coefficients which are fixed. Natural boundary conditions signal to PETSc that boundary
3375   integrals should be performed, using the kernels from `PetscDSSetBdResidual()`.
3376 
3377 .seealso: `PetscDS`, `PetscWeakForm`, `DMLabel`, `DMBoundaryConditionType`, `PetscDSAddBoundaryByName()`, `PetscDSGetBoundary()`, `PetscDSSetResidual()`, `PetscDSSetBdResidual()`
3378 @*/
3379 PetscErrorCode PetscDSAddBoundary(PetscDS ds, DMBoundaryConditionType type, const char name[], DMLabel label, PetscInt Nv, const PetscInt values[], PetscInt field, PetscInt Nc, const PetscInt comps[], void (*bcFunc)(void), void (*bcFunc_t)(void), void *ctx, PetscInt *bd)
3380 {
3381   DSBoundary  head = ds->boundary, b;
3382   PetscInt    n    = 0;
3383   const char *lname;
3384 
3385   PetscFunctionBegin;
3386   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
3387   PetscValidLogicalCollectiveEnum(ds, type, 2);
3388   PetscAssertPointer(name, 3);
3389   PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 4);
3390   PetscValidLogicalCollectiveInt(ds, Nv, 5);
3391   PetscValidLogicalCollectiveInt(ds, field, 7);
3392   PetscValidLogicalCollectiveInt(ds, Nc, 8);
3393   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);
3394   if (Nc > 0) {
3395     PetscInt *fcomps;
3396     PetscInt  c;
3397 
3398     PetscCall(PetscDSGetComponents(ds, &fcomps));
3399     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);
3400     for (c = 0; c < Nc; ++c) {
3401       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);
3402     }
3403   }
3404   PetscCall(PetscNew(&b));
3405   PetscCall(PetscStrallocpy(name, (char **)&b->name));
3406   PetscCall(PetscWeakFormCreate(PETSC_COMM_SELF, &b->wf));
3407   PetscCall(PetscWeakFormSetNumFields(b->wf, ds->Nf));
3408   PetscCall(PetscMalloc1(Nv, &b->values));
3409   if (Nv) PetscCall(PetscArraycpy(b->values, values, Nv));
3410   PetscCall(PetscMalloc1(Nc, &b->comps));
3411   if (Nc) PetscCall(PetscArraycpy(b->comps, comps, Nc));
3412   PetscCall(PetscObjectGetName((PetscObject)label, &lname));
3413   PetscCall(PetscStrallocpy(lname, (char **)&b->lname));
3414   b->type   = type;
3415   b->label  = label;
3416   b->Nv     = Nv;
3417   b->field  = field;
3418   b->Nc     = Nc;
3419   b->func   = bcFunc;
3420   b->func_t = bcFunc_t;
3421   b->ctx    = ctx;
3422   b->next   = NULL;
3423   /* Append to linked list so that we can preserve the order */
3424   if (!head) ds->boundary = b;
3425   while (head) {
3426     if (!head->next) {
3427       head->next = b;
3428       head       = b;
3429     }
3430     head = head->next;
3431     ++n;
3432   }
3433   if (bd) {
3434     PetscAssertPointer(bd, 13);
3435     *bd = n;
3436   }
3437   PetscFunctionReturn(PETSC_SUCCESS);
3438 }
3439 
3440 // PetscClangLinter pragma ignore: -fdoc-section-header-unknown
3441 /*@C
3442   PetscDSAddBoundaryByName - Add a boundary condition to the model.
3443 
3444   Collective
3445 
3446   Input Parameters:
3447 + ds       - The `PetscDS` object
3448 . type     - The type of condition, e.g. `DM_BC_ESSENTIAL`/`DM_BC_ESSENTIAL_FIELD` (Dirichlet), or `DM_BC_NATURAL` (Neumann)
3449 . name     - The BC name
3450 . lname    - The naem of the label defining constrained points
3451 . Nv       - The number of `DMLabel` values for constrained points
3452 . values   - An array of label values for constrained points
3453 . field    - The field to constrain
3454 . Nc       - The number of constrained field components (0 will constrain all fields)
3455 . comps    - An array of constrained component numbers
3456 . bcFunc   - A pointwise function giving boundary values
3457 . bcFunc_t - A pointwise function giving the time derivative of the boundary values, or NULL
3458 - ctx      - An optional user context for bcFunc
3459 
3460   Output Parameter:
3461 . bd - The boundary number
3462 
3463   Options Database Keys:
3464 + -bc_<boundary name> <num>      - Overrides the boundary ids
3465 - -bc_<boundary name>_comp <num> - Overrides the boundary components
3466 
3467   Calling Sequence of `bcFunc` and `bcFunc_t`:
3468   If the type is `DM_BC_ESSENTIAL`
3469 .vb
3470   void bcFunc(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar bcval[])
3471 .ve
3472   If the type is `DM_BC_ESSENTIAL_FIELD` or other _FIELD value,
3473 .vb
3474   void bcFunc(PetscInt dim, PetscInt Nf, PetscInt NfAux,
3475               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
3476               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
3477               PetscReal time, const PetscReal x[], PetscScalar bcval[])
3478 .ve
3479 + dim - the spatial dimension
3480 . Nf - the number of fields
3481 . uOff - the offset into u[] and u_t[] for each field
3482 . uOff_x - the offset into u_x[] for each field
3483 . u - each field evaluated at the current point
3484 . u_t - the time derivative of each field evaluated at the current point
3485 . u_x - the gradient of each field evaluated at the current point
3486 . aOff - the offset into a[] and a_t[] for each auxiliary field
3487 . aOff_x - the offset into a_x[] for each auxiliary field
3488 . a - each auxiliary field evaluated at the current point
3489 . a_t - the time derivative of each auxiliary field evaluated at the current point
3490 . a_x - the gradient of auxiliary each field evaluated at the current point
3491 . t - current time
3492 . x - coordinates of the current point
3493 . numConstants - number of constant parameters
3494 . constants - constant parameters
3495 - bcval - output values at the current point
3496 
3497   Level: developer
3498 
3499   Notes:
3500   The pointwise functions are used to provide boundary values for essential boundary
3501   conditions. In FEM, they are acting upon by dual basis functionals to generate FEM
3502   coefficients which are fixed. Natural boundary conditions signal to PETSc that boundary
3503   integrals should be performed, using the kernels from `PetscDSSetBdResidual()`.
3504 
3505   This function should only be used with `DMFOREST` currently, since labels cannot be defined before the underlying `DMPLEX` is built.
3506 
3507 .seealso: `PetscDS`, `PetscWeakForm`, `DMLabel`, `DMBoundaryConditionType`, `PetscDSAddBoundary()`, `PetscDSGetBoundary()`, `PetscDSSetResidual()`, `PetscDSSetBdResidual()`
3508 @*/
3509 PetscErrorCode PetscDSAddBoundaryByName(PetscDS ds, DMBoundaryConditionType type, const char name[], const char lname[], PetscInt Nv, const PetscInt values[], PetscInt field, PetscInt Nc, const PetscInt comps[], void (*bcFunc)(void), void (*bcFunc_t)(void), void *ctx, PetscInt *bd)
3510 {
3511   DSBoundary head = ds->boundary, b;
3512   PetscInt   n    = 0;
3513 
3514   PetscFunctionBegin;
3515   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
3516   PetscValidLogicalCollectiveEnum(ds, type, 2);
3517   PetscAssertPointer(name, 3);
3518   PetscAssertPointer(lname, 4);
3519   PetscValidLogicalCollectiveInt(ds, Nv, 5);
3520   PetscValidLogicalCollectiveInt(ds, field, 7);
3521   PetscValidLogicalCollectiveInt(ds, Nc, 8);
3522   PetscCall(PetscNew(&b));
3523   PetscCall(PetscStrallocpy(name, (char **)&b->name));
3524   PetscCall(PetscWeakFormCreate(PETSC_COMM_SELF, &b->wf));
3525   PetscCall(PetscWeakFormSetNumFields(b->wf, ds->Nf));
3526   PetscCall(PetscMalloc1(Nv, &b->values));
3527   if (Nv) PetscCall(PetscArraycpy(b->values, values, Nv));
3528   PetscCall(PetscMalloc1(Nc, &b->comps));
3529   if (Nc) PetscCall(PetscArraycpy(b->comps, comps, Nc));
3530   PetscCall(PetscStrallocpy(lname, (char **)&b->lname));
3531   b->type   = type;
3532   b->label  = NULL;
3533   b->Nv     = Nv;
3534   b->field  = field;
3535   b->Nc     = Nc;
3536   b->func   = bcFunc;
3537   b->func_t = bcFunc_t;
3538   b->ctx    = ctx;
3539   b->next   = NULL;
3540   /* Append to linked list so that we can preserve the order */
3541   if (!head) ds->boundary = b;
3542   while (head) {
3543     if (!head->next) {
3544       head->next = b;
3545       head       = b;
3546     }
3547     head = head->next;
3548     ++n;
3549   }
3550   if (bd) {
3551     PetscAssertPointer(bd, 13);
3552     *bd = n;
3553   }
3554   PetscFunctionReturn(PETSC_SUCCESS);
3555 }
3556 
3557 /*@C
3558   PetscDSUpdateBoundary - Change a boundary condition for the model.
3559 
3560   Input Parameters:
3561 + ds       - The `PetscDS` object
3562 . bd       - The boundary condition number
3563 . type     - The type of condition, e.g. `DM_BC_ESSENTIAL`/`DM_BC_ESSENTIAL_FIELD` (Dirichlet), or `DM_BC_NATURAL` (Neumann)
3564 . name     - The BC name
3565 . label    - The label defining constrained points
3566 . Nv       - The number of `DMLabel` ids for constrained points
3567 . values   - An array of ids for constrained points
3568 . field    - The field to constrain
3569 . Nc       - The number of constrained field components
3570 . comps    - An array of constrained component numbers
3571 . bcFunc   - A pointwise function giving boundary values
3572 . bcFunc_t - A pointwise function giving the time derivative of the boundary values, or NULL
3573 - ctx      - An optional user context for bcFunc
3574 
3575   Level: developer
3576 
3577   Notes:
3578   The pointwise functions are used to provide boundary values for essential boundary
3579   conditions. In FEM, they are acting upon by dual basis functionals to generate FEM
3580   coefficients which are fixed. Natural boundary conditions signal to PETSc that boundary
3581   integrals should be performed, using the kernels from `PetscDSSetBdResidual()`.
3582 
3583   The boundary condition number is the order in which it was registered. The user can get the number of boundary conditions from `PetscDSGetNumBoundary()`.
3584   See `PetscDSAddBoundary()` for a description of the calling sequences for the callbacks.
3585 
3586 .seealso: `PetscDS`, `PetscWeakForm`, `DMBoundaryConditionType`, `PetscDSAddBoundary()`, `PetscDSGetBoundary()`, `PetscDSGetNumBoundary()`, `DMLabel`
3587 @*/
3588 PetscErrorCode PetscDSUpdateBoundary(PetscDS ds, PetscInt bd, DMBoundaryConditionType type, const char name[], DMLabel label, PetscInt Nv, const PetscInt values[], PetscInt field, PetscInt Nc, const PetscInt comps[], void (*bcFunc)(void), void (*bcFunc_t)(void), void *ctx)
3589 {
3590   DSBoundary b = ds->boundary;
3591   PetscInt   n = 0;
3592 
3593   PetscFunctionBegin;
3594   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
3595   while (b) {
3596     if (n == bd) break;
3597     b = b->next;
3598     ++n;
3599   }
3600   PetscCheck(b, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %" PetscInt_FMT " is not in [0, %" PetscInt_FMT ")", bd, n);
3601   if (name) {
3602     PetscCall(PetscFree(b->name));
3603     PetscCall(PetscStrallocpy(name, (char **)&b->name));
3604   }
3605   b->type = type;
3606   if (label) {
3607     const char *name;
3608 
3609     b->label = label;
3610     PetscCall(PetscFree(b->lname));
3611     PetscCall(PetscObjectGetName((PetscObject)label, &name));
3612     PetscCall(PetscStrallocpy(name, (char **)&b->lname));
3613   }
3614   if (Nv >= 0) {
3615     b->Nv = Nv;
3616     PetscCall(PetscFree(b->values));
3617     PetscCall(PetscMalloc1(Nv, &b->values));
3618     if (Nv) PetscCall(PetscArraycpy(b->values, values, Nv));
3619   }
3620   if (field >= 0) b->field = field;
3621   if (Nc >= 0) {
3622     b->Nc = Nc;
3623     PetscCall(PetscFree(b->comps));
3624     PetscCall(PetscMalloc1(Nc, &b->comps));
3625     if (Nc) PetscCall(PetscArraycpy(b->comps, comps, Nc));
3626   }
3627   if (bcFunc) b->func = bcFunc;
3628   if (bcFunc_t) b->func_t = bcFunc_t;
3629   if (ctx) b->ctx = ctx;
3630   PetscFunctionReturn(PETSC_SUCCESS);
3631 }
3632 
3633 /*@
3634   PetscDSGetNumBoundary - Get the number of registered BC
3635 
3636   Input Parameter:
3637 . ds - The `PetscDS` object
3638 
3639   Output Parameter:
3640 . numBd - The number of BC
3641 
3642   Level: intermediate
3643 
3644 .seealso: `PetscDS`, `PetscDSAddBoundary()`, `PetscDSGetBoundary()`
3645 @*/
3646 PetscErrorCode PetscDSGetNumBoundary(PetscDS ds, PetscInt *numBd)
3647 {
3648   DSBoundary b = ds->boundary;
3649 
3650   PetscFunctionBegin;
3651   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
3652   PetscAssertPointer(numBd, 2);
3653   *numBd = 0;
3654   while (b) {
3655     ++(*numBd);
3656     b = b->next;
3657   }
3658   PetscFunctionReturn(PETSC_SUCCESS);
3659 }
3660 
3661 /*@C
3662   PetscDSGetBoundary - Gets a boundary condition to the model
3663 
3664   Input Parameters:
3665 + ds - The `PetscDS` object
3666 - bd - The BC number
3667 
3668   Output Parameters:
3669 + wf     - The `PetscWeakForm` holding the pointwise functions
3670 . type   - The type of condition, e.g. `DM_BC_ESSENTIAL`/`DM_BC_ESSENTIAL_FIELD` (Dirichlet), or `DM_BC_NATURAL` (Neumann)
3671 . name   - The BC name
3672 . label  - The label defining constrained points
3673 . Nv     - The number of `DMLabel` ids for constrained points
3674 . values - An array of ids for constrained points
3675 . field  - The field to constrain
3676 . Nc     - The number of constrained field components
3677 . comps  - An array of constrained component numbers
3678 . func   - A pointwise function giving boundary values
3679 . func_t - A pointwise function giving the time derivative of the boundary values
3680 - ctx    - An optional user context for bcFunc
3681 
3682   Options Database Keys:
3683 + -bc_<boundary name> <num>      - Overrides the boundary ids
3684 - -bc_<boundary name>_comp <num> - Overrides the boundary components
3685 
3686   Level: developer
3687 
3688 .seealso: `PetscDS`, `PetscWeakForm`, `DMBoundaryConditionType`, `PetscDSAddBoundary()`, `DMLabel`
3689 @*/
3690 PetscErrorCode PetscDSGetBoundary(PetscDS ds, PetscInt bd, PetscWeakForm *wf, DMBoundaryConditionType *type, const char *name[], DMLabel *label, PetscInt *Nv, const PetscInt *values[], PetscInt *field, PetscInt *Nc, const PetscInt *comps[], void (**func)(void), void (**func_t)(void), void **ctx)
3691 {
3692   DSBoundary b = ds->boundary;
3693   PetscInt   n = 0;
3694 
3695   PetscFunctionBegin;
3696   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
3697   while (b) {
3698     if (n == bd) break;
3699     b = b->next;
3700     ++n;
3701   }
3702   PetscCheck(b, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %" PetscInt_FMT " is not in [0, %" PetscInt_FMT ")", bd, n);
3703   if (wf) {
3704     PetscAssertPointer(wf, 3);
3705     *wf = b->wf;
3706   }
3707   if (type) {
3708     PetscAssertPointer(type, 4);
3709     *type = b->type;
3710   }
3711   if (name) {
3712     PetscAssertPointer(name, 5);
3713     *name = b->name;
3714   }
3715   if (label) {
3716     PetscAssertPointer(label, 6);
3717     *label = b->label;
3718   }
3719   if (Nv) {
3720     PetscAssertPointer(Nv, 7);
3721     *Nv = b->Nv;
3722   }
3723   if (values) {
3724     PetscAssertPointer(values, 8);
3725     *values = b->values;
3726   }
3727   if (field) {
3728     PetscAssertPointer(field, 9);
3729     *field = b->field;
3730   }
3731   if (Nc) {
3732     PetscAssertPointer(Nc, 10);
3733     *Nc = b->Nc;
3734   }
3735   if (comps) {
3736     PetscAssertPointer(comps, 11);
3737     *comps = b->comps;
3738   }
3739   if (func) {
3740     PetscAssertPointer(func, 12);
3741     *func = b->func;
3742   }
3743   if (func_t) {
3744     PetscAssertPointer(func_t, 13);
3745     *func_t = b->func_t;
3746   }
3747   if (ctx) {
3748     PetscAssertPointer(ctx, 14);
3749     *ctx = b->ctx;
3750   }
3751   PetscFunctionReturn(PETSC_SUCCESS);
3752 }
3753 
3754 static PetscErrorCode DSBoundaryDuplicate_Internal(DSBoundary b, DSBoundary *bNew)
3755 {
3756   PetscFunctionBegin;
3757   PetscCall(PetscNew(bNew));
3758   PetscCall(PetscWeakFormCreate(PETSC_COMM_SELF, &(*bNew)->wf));
3759   PetscCall(PetscWeakFormCopy(b->wf, (*bNew)->wf));
3760   PetscCall(PetscStrallocpy(b->name, (char **)&((*bNew)->name)));
3761   PetscCall(PetscStrallocpy(b->lname, (char **)&((*bNew)->lname)));
3762   (*bNew)->type  = b->type;
3763   (*bNew)->label = b->label;
3764   (*bNew)->Nv    = b->Nv;
3765   PetscCall(PetscMalloc1(b->Nv, &(*bNew)->values));
3766   PetscCall(PetscArraycpy((*bNew)->values, b->values, b->Nv));
3767   (*bNew)->field = b->field;
3768   (*bNew)->Nc    = b->Nc;
3769   PetscCall(PetscMalloc1(b->Nc, &(*bNew)->comps));
3770   PetscCall(PetscArraycpy((*bNew)->comps, b->comps, b->Nc));
3771   (*bNew)->func   = b->func;
3772   (*bNew)->func_t = b->func_t;
3773   (*bNew)->ctx    = b->ctx;
3774   PetscFunctionReturn(PETSC_SUCCESS);
3775 }
3776 
3777 /*@
3778   PetscDSCopyBoundary - Copy all boundary condition objects to the new problem
3779 
3780   Not Collective
3781 
3782   Input Parameters:
3783 + ds        - The source `PetscDS` object
3784 . numFields - The number of selected fields, or `PETSC_DEFAULT` for all fields
3785 - fields    - The selected fields, or NULL for all fields
3786 
3787   Output Parameter:
3788 . newds - The target `PetscDS`, now with a copy of the boundary conditions
3789 
3790   Level: intermediate
3791 
3792 .seealso: `PetscDS`, `DMBoundary`, `PetscDSCopyEquations()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
3793 @*/
3794 PetscErrorCode PetscDSCopyBoundary(PetscDS ds, PetscInt numFields, const PetscInt fields[], PetscDS newds)
3795 {
3796   DSBoundary b, *lastnext;
3797 
3798   PetscFunctionBegin;
3799   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
3800   PetscValidHeaderSpecific(newds, PETSCDS_CLASSID, 4);
3801   if (ds == newds) PetscFunctionReturn(PETSC_SUCCESS);
3802   PetscCall(PetscDSDestroyBoundary(newds));
3803   lastnext = &(newds->boundary);
3804   for (b = ds->boundary; b; b = b->next) {
3805     DSBoundary bNew;
3806     PetscInt   fieldNew = -1;
3807 
3808     if (numFields > 0 && fields) {
3809       PetscInt f;
3810 
3811       for (f = 0; f < numFields; ++f)
3812         if (b->field == fields[f]) break;
3813       if (f == numFields) continue;
3814       fieldNew = f;
3815     }
3816     PetscCall(DSBoundaryDuplicate_Internal(b, &bNew));
3817     bNew->field = fieldNew < 0 ? b->field : fieldNew;
3818     *lastnext   = bNew;
3819     lastnext    = &(bNew->next);
3820   }
3821   PetscFunctionReturn(PETSC_SUCCESS);
3822 }
3823 
3824 /*@
3825   PetscDSDestroyBoundary - Remove all `DMBoundary` objects from the `PetscDS`
3826 
3827   Not Collective
3828 
3829   Input Parameter:
3830 . ds - The `PetscDS` object
3831 
3832   Level: intermediate
3833 
3834 .seealso: `PetscDS`, `DMBoundary`, `PetscDSCopyBoundary()`, `PetscDSCopyEquations()`
3835 @*/
3836 PetscErrorCode PetscDSDestroyBoundary(PetscDS ds)
3837 {
3838   DSBoundary next = ds->boundary;
3839 
3840   PetscFunctionBegin;
3841   while (next) {
3842     DSBoundary b = next;
3843 
3844     next = b->next;
3845     PetscCall(PetscWeakFormDestroy(&b->wf));
3846     PetscCall(PetscFree(b->name));
3847     PetscCall(PetscFree(b->lname));
3848     PetscCall(PetscFree(b->values));
3849     PetscCall(PetscFree(b->comps));
3850     PetscCall(PetscFree(b));
3851   }
3852   PetscFunctionReturn(PETSC_SUCCESS);
3853 }
3854 
3855 /*@
3856   PetscDSSelectDiscretizations - Copy discretizations to the new problem with different field layout
3857 
3858   Not Collective
3859 
3860   Input Parameters:
3861 + prob      - The `PetscDS` object
3862 . numFields - Number of new fields
3863 - fields    - Old field number for each new field
3864 
3865   Output Parameter:
3866 . newprob - The `PetscDS` copy
3867 
3868   Level: intermediate
3869 
3870 .seealso: `PetscDS`, `PetscDSSelectEquations()`, `PetscDSCopyBoundary()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
3871 @*/
3872 PetscErrorCode PetscDSSelectDiscretizations(PetscDS prob, PetscInt numFields, const PetscInt fields[], PetscDS newprob)
3873 {
3874   PetscInt Nf, Nfn, fn;
3875 
3876   PetscFunctionBegin;
3877   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
3878   if (fields) PetscAssertPointer(fields, 3);
3879   PetscValidHeaderSpecific(newprob, PETSCDS_CLASSID, 4);
3880   PetscCall(PetscDSGetNumFields(prob, &Nf));
3881   PetscCall(PetscDSGetNumFields(newprob, &Nfn));
3882   numFields = numFields < 0 ? Nf : numFields;
3883   for (fn = 0; fn < numFields; ++fn) {
3884     const PetscInt f = fields ? fields[fn] : fn;
3885     PetscObject    disc;
3886 
3887     if (f >= Nf) continue;
3888     PetscCall(PetscDSGetDiscretization(prob, f, &disc));
3889     PetscCall(PetscDSSetDiscretization(newprob, fn, disc));
3890   }
3891   PetscFunctionReturn(PETSC_SUCCESS);
3892 }
3893 
3894 /*@
3895   PetscDSSelectEquations - Copy pointwise function pointers to the new problem with different field layout
3896 
3897   Not Collective
3898 
3899   Input Parameters:
3900 + prob      - The `PetscDS` object
3901 . numFields - Number of new fields
3902 - fields    - Old field number for each new field
3903 
3904   Output Parameter:
3905 . newprob - The `PetscDS` copy
3906 
3907   Level: intermediate
3908 
3909 .seealso: `PetscDS`, `PetscDSSelectDiscretizations()`, `PetscDSCopyBoundary()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
3910 @*/
3911 PetscErrorCode PetscDSSelectEquations(PetscDS prob, PetscInt numFields, const PetscInt fields[], PetscDS newprob)
3912 {
3913   PetscInt Nf, Nfn, fn, gn;
3914 
3915   PetscFunctionBegin;
3916   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
3917   if (fields) PetscAssertPointer(fields, 3);
3918   PetscValidHeaderSpecific(newprob, PETSCDS_CLASSID, 4);
3919   PetscCall(PetscDSGetNumFields(prob, &Nf));
3920   PetscCall(PetscDSGetNumFields(newprob, &Nfn));
3921   PetscCheck(numFields <= Nfn, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_SIZ, "Number of fields %" PetscInt_FMT " to transfer must not be greater then the total number of fields %" PetscInt_FMT, numFields, Nfn);
3922   for (fn = 0; fn < numFields; ++fn) {
3923     const PetscInt   f = fields ? fields[fn] : fn;
3924     PetscPointFunc   obj;
3925     PetscPointFunc   f0, f1;
3926     PetscBdPointFunc f0Bd, f1Bd;
3927     PetscRiemannFunc r;
3928 
3929     if (f >= Nf) continue;
3930     PetscCall(PetscDSGetObjective(prob, f, &obj));
3931     PetscCall(PetscDSGetResidual(prob, f, &f0, &f1));
3932     PetscCall(PetscDSGetBdResidual(prob, f, &f0Bd, &f1Bd));
3933     PetscCall(PetscDSGetRiemannSolver(prob, f, &r));
3934     PetscCall(PetscDSSetObjective(newprob, fn, obj));
3935     PetscCall(PetscDSSetResidual(newprob, fn, f0, f1));
3936     PetscCall(PetscDSSetBdResidual(newprob, fn, f0Bd, f1Bd));
3937     PetscCall(PetscDSSetRiemannSolver(newprob, fn, r));
3938     for (gn = 0; gn < numFields; ++gn) {
3939       const PetscInt  g = fields ? fields[gn] : gn;
3940       PetscPointJac   g0, g1, g2, g3;
3941       PetscPointJac   g0p, g1p, g2p, g3p;
3942       PetscBdPointJac g0Bd, g1Bd, g2Bd, g3Bd;
3943 
3944       if (g >= Nf) continue;
3945       PetscCall(PetscDSGetJacobian(prob, f, g, &g0, &g1, &g2, &g3));
3946       PetscCall(PetscDSGetJacobianPreconditioner(prob, f, g, &g0p, &g1p, &g2p, &g3p));
3947       PetscCall(PetscDSGetBdJacobian(prob, f, g, &g0Bd, &g1Bd, &g2Bd, &g3Bd));
3948       PetscCall(PetscDSSetJacobian(newprob, fn, gn, g0, g1, g2, g3));
3949       PetscCall(PetscDSSetJacobianPreconditioner(newprob, fn, gn, g0p, g1p, g2p, g3p));
3950       PetscCall(PetscDSSetBdJacobian(newprob, fn, gn, g0Bd, g1Bd, g2Bd, g3Bd));
3951     }
3952   }
3953   PetscFunctionReturn(PETSC_SUCCESS);
3954 }
3955 
3956 /*@
3957   PetscDSCopyEquations - Copy all pointwise function pointers to another `PetscDS`
3958 
3959   Not Collective
3960 
3961   Input Parameter:
3962 . prob - The `PetscDS` object
3963 
3964   Output Parameter:
3965 . newprob - The `PetscDS` copy
3966 
3967   Level: intermediate
3968 
3969 .seealso: `PetscDS`, `PetscDSCopyBoundary()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
3970 @*/
3971 PetscErrorCode PetscDSCopyEquations(PetscDS prob, PetscDS newprob)
3972 {
3973   PetscWeakForm wf, newwf;
3974   PetscInt      Nf, Ng;
3975 
3976   PetscFunctionBegin;
3977   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
3978   PetscValidHeaderSpecific(newprob, PETSCDS_CLASSID, 2);
3979   PetscCall(PetscDSGetNumFields(prob, &Nf));
3980   PetscCall(PetscDSGetNumFields(newprob, &Ng));
3981   PetscCheck(Nf == Ng, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_SIZ, "Number of fields must match %" PetscInt_FMT " != %" PetscInt_FMT, Nf, Ng);
3982   PetscCall(PetscDSGetWeakForm(prob, &wf));
3983   PetscCall(PetscDSGetWeakForm(newprob, &newwf));
3984   PetscCall(PetscWeakFormCopy(wf, newwf));
3985   PetscFunctionReturn(PETSC_SUCCESS);
3986 }
3987 
3988 /*@
3989   PetscDSCopyConstants - Copy all constants to another `PetscDS`
3990 
3991   Not Collective
3992 
3993   Input Parameter:
3994 . prob - The `PetscDS` object
3995 
3996   Output Parameter:
3997 . newprob - The `PetscDS` copy
3998 
3999   Level: intermediate
4000 
4001 .seealso: `PetscDS`, `PetscDSCopyBoundary()`, `PetscDSCopyEquations()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
4002 @*/
4003 PetscErrorCode PetscDSCopyConstants(PetscDS prob, PetscDS newprob)
4004 {
4005   PetscInt           Nc;
4006   const PetscScalar *constants;
4007 
4008   PetscFunctionBegin;
4009   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
4010   PetscValidHeaderSpecific(newprob, PETSCDS_CLASSID, 2);
4011   PetscCall(PetscDSGetConstants(prob, &Nc, &constants));
4012   PetscCall(PetscDSSetConstants(newprob, Nc, (PetscScalar *)constants));
4013   PetscFunctionReturn(PETSC_SUCCESS);
4014 }
4015 
4016 /*@
4017   PetscDSCopyExactSolutions - Copy all exact solutions to another `PetscDS`
4018 
4019   Not Collective
4020 
4021   Input Parameter:
4022 . ds - The `PetscDS` object
4023 
4024   Output Parameter:
4025 . newds - The `PetscDS` copy
4026 
4027   Level: intermediate
4028 
4029 .seealso: `PetscDS`, `PetscDSCopyBoundary()`, `PetscDSCopyEquations()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
4030 @*/
4031 PetscErrorCode PetscDSCopyExactSolutions(PetscDS ds, PetscDS newds)
4032 {
4033   PetscSimplePointFunc sol;
4034   void                *ctx;
4035   PetscInt             Nf, f;
4036 
4037   PetscFunctionBegin;
4038   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
4039   PetscValidHeaderSpecific(newds, PETSCDS_CLASSID, 2);
4040   PetscCall(PetscDSGetNumFields(ds, &Nf));
4041   for (f = 0; f < Nf; ++f) {
4042     PetscCall(PetscDSGetExactSolution(ds, f, &sol, &ctx));
4043     PetscCall(PetscDSSetExactSolution(newds, f, sol, ctx));
4044     PetscCall(PetscDSGetExactSolutionTimeDerivative(ds, f, &sol, &ctx));
4045     PetscCall(PetscDSSetExactSolutionTimeDerivative(newds, f, sol, ctx));
4046   }
4047   PetscFunctionReturn(PETSC_SUCCESS);
4048 }
4049 
4050 PetscErrorCode PetscDSCopy(PetscDS ds, DM dmNew, PetscDS dsNew)
4051 {
4052   DSBoundary b;
4053   PetscInt   cdim, Nf, f, d;
4054   PetscBool  isCohesive;
4055   void      *ctx;
4056 
4057   PetscFunctionBegin;
4058   PetscCall(PetscDSCopyConstants(ds, dsNew));
4059   PetscCall(PetscDSCopyExactSolutions(ds, dsNew));
4060   PetscCall(PetscDSSelectDiscretizations(ds, PETSC_DETERMINE, NULL, dsNew));
4061   PetscCall(PetscDSCopyEquations(ds, dsNew));
4062   PetscCall(PetscDSGetNumFields(ds, &Nf));
4063   for (f = 0; f < Nf; ++f) {
4064     PetscCall(PetscDSGetContext(ds, f, &ctx));
4065     PetscCall(PetscDSSetContext(dsNew, f, ctx));
4066     PetscCall(PetscDSGetCohesive(ds, f, &isCohesive));
4067     PetscCall(PetscDSSetCohesive(dsNew, f, isCohesive));
4068     PetscCall(PetscDSGetJetDegree(ds, f, &d));
4069     PetscCall(PetscDSSetJetDegree(dsNew, f, d));
4070   }
4071   if (Nf) {
4072     PetscCall(PetscDSGetCoordinateDimension(ds, &cdim));
4073     PetscCall(PetscDSSetCoordinateDimension(dsNew, cdim));
4074   }
4075   PetscCall(PetscDSCopyBoundary(ds, PETSC_DETERMINE, NULL, dsNew));
4076   for (b = dsNew->boundary; b; b = b->next) {
4077     PetscCall(DMGetLabel(dmNew, b->lname, &b->label));
4078     /* Do not check if label exists here, since p4est calls this for the reference tree which does not have the labels */
4079     //PetscCheck(b->label,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Label %s missing in new DM", name);
4080   }
4081   PetscFunctionReturn(PETSC_SUCCESS);
4082 }
4083 
4084 PetscErrorCode PetscDSGetHeightSubspace(PetscDS prob, PetscInt height, PetscDS *subprob)
4085 {
4086   PetscInt dim, Nf, f;
4087 
4088   PetscFunctionBegin;
4089   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
4090   PetscAssertPointer(subprob, 3);
4091   if (height == 0) {
4092     *subprob = prob;
4093     PetscFunctionReturn(PETSC_SUCCESS);
4094   }
4095   PetscCall(PetscDSGetNumFields(prob, &Nf));
4096   PetscCall(PetscDSGetSpatialDimension(prob, &dim));
4097   PetscCheck(height <= dim, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_OUTOFRANGE, "DS can only handle height in [0, %" PetscInt_FMT "], not %" PetscInt_FMT, dim, height);
4098   if (!prob->subprobs) PetscCall(PetscCalloc1(dim, &prob->subprobs));
4099   if (!prob->subprobs[height - 1]) {
4100     PetscInt cdim;
4101 
4102     PetscCall(PetscDSCreate(PetscObjectComm((PetscObject)prob), &prob->subprobs[height - 1]));
4103     PetscCall(PetscDSGetCoordinateDimension(prob, &cdim));
4104     PetscCall(PetscDSSetCoordinateDimension(prob->subprobs[height - 1], cdim));
4105     for (f = 0; f < Nf; ++f) {
4106       PetscFE      subfe;
4107       PetscObject  obj;
4108       PetscClassId id;
4109 
4110       PetscCall(PetscDSGetDiscretization(prob, f, &obj));
4111       PetscCall(PetscObjectGetClassId(obj, &id));
4112       if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetHeightSubspace((PetscFE)obj, height, &subfe));
4113       else SETERRQ(PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONG, "Unsupported discretization type for field %" PetscInt_FMT, f);
4114       PetscCall(PetscDSSetDiscretization(prob->subprobs[height - 1], f, (PetscObject)subfe));
4115     }
4116   }
4117   *subprob = prob->subprobs[height - 1];
4118   PetscFunctionReturn(PETSC_SUCCESS);
4119 }
4120 
4121 PetscErrorCode PetscDSPermuteQuadPoint(PetscDS ds, PetscInt ornt, PetscInt field, PetscInt q, PetscInt *qperm)
4122 {
4123   IS              permIS;
4124   PetscQuadrature quad;
4125   DMPolytopeType  ct;
4126   const PetscInt *perm;
4127   PetscInt        Na, Nq;
4128 
4129   PetscFunctionBeginHot;
4130   PetscCall(PetscFEGetQuadrature((PetscFE)ds->disc[field], &quad));
4131   PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL));
4132   PetscCall(PetscQuadratureGetCellType(quad, &ct));
4133   PetscCheck(q >= 0 && q < Nq, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Quadrature point %" PetscInt_FMT " is not in [0, %" PetscInt_FMT ")", q, Nq);
4134   Na = DMPolytopeTypeGetNumArrangments(ct) / 2;
4135   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);
4136   if (!ds->quadPerm[(PetscInt)ct]) PetscCall(PetscQuadratureComputePermutations(quad, NULL, &ds->quadPerm[(PetscInt)ct]));
4137   permIS = ds->quadPerm[(PetscInt)ct][ornt + Na];
4138   PetscCall(ISGetIndices(permIS, &perm));
4139   *qperm = perm[q];
4140   PetscCall(ISRestoreIndices(permIS, &perm));
4141   PetscFunctionReturn(PETSC_SUCCESS);
4142 }
4143 
4144 PetscErrorCode PetscDSGetDiscType_Internal(PetscDS ds, PetscInt f, PetscDiscType *disctype)
4145 {
4146   PetscObject  obj;
4147   PetscClassId id;
4148   PetscInt     Nf;
4149 
4150   PetscFunctionBegin;
4151   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
4152   PetscAssertPointer(disctype, 3);
4153   *disctype = PETSC_DISC_NONE;
4154   PetscCall(PetscDSGetNumFields(ds, &Nf));
4155   PetscCheck(f < Nf, PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_SIZ, "Field %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, Nf);
4156   PetscCall(PetscDSGetDiscretization(ds, f, &obj));
4157   if (obj) {
4158     PetscCall(PetscObjectGetClassId(obj, &id));
4159     if (id == PETSCFE_CLASSID) *disctype = PETSC_DISC_FE;
4160     else *disctype = PETSC_DISC_FV;
4161   }
4162   PetscFunctionReturn(PETSC_SUCCESS);
4163 }
4164 
4165 static PetscErrorCode PetscDSDestroy_Basic(PetscDS ds)
4166 {
4167   PetscFunctionBegin;
4168   PetscCall(PetscFree(ds->data));
4169   PetscFunctionReturn(PETSC_SUCCESS);
4170 }
4171 
4172 static PetscErrorCode PetscDSInitialize_Basic(PetscDS ds)
4173 {
4174   PetscFunctionBegin;
4175   ds->ops->setfromoptions = NULL;
4176   ds->ops->setup          = NULL;
4177   ds->ops->view           = NULL;
4178   ds->ops->destroy        = PetscDSDestroy_Basic;
4179   PetscFunctionReturn(PETSC_SUCCESS);
4180 }
4181 
4182 /*MC
4183   PETSCDSBASIC = "basic" - A discrete system with pointwise residual and boundary residual functions
4184 
4185   Level: intermediate
4186 
4187 .seealso: `PetscDSType`, `PetscDSCreate()`, `PetscDSSetType()`
4188 M*/
4189 
4190 PETSC_EXTERN PetscErrorCode PetscDSCreate_Basic(PetscDS ds)
4191 {
4192   PetscDS_Basic *b;
4193 
4194   PetscFunctionBegin;
4195   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
4196   PetscCall(PetscNew(&b));
4197   ds->data = b;
4198 
4199   PetscCall(PetscDSInitialize_Basic(ds));
4200   PetscFunctionReturn(PETSC_SUCCESS);
4201 }
4202