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