xref: /petsc/src/dm/dt/interface/dtds.c (revision f17e87946373dc9caaa8e0926271963e0f9af11e)
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 /*@C
9   PetscDSRegister - Adds a new PetscDS implementation
10 
11   Not Collective
12 
13   Input Parameters:
14 + name        - The name of a new user-defined creation routine
15 - create_func - The creation routine itself
16 
17   Notes:
18   PetscDSRegister() may be called multiple times to add several user-defined PetscDSs
19 
20   Sample usage:
21 .vb
22     PetscDSRegister("my_ds", MyPetscDSCreate);
23 .ve
24 
25   Then, your PetscDS type can be chosen with the procedural interface via
26 .vb
27     PetscDSCreate(MPI_Comm, PetscDS *);
28     PetscDSSetType(PetscDS, "my_ds");
29 .ve
30    or at runtime via the option
31 .vb
32     -petscds_type my_ds
33 .ve
34 
35   Level: advanced
36 
37 .keywords: PetscDS, register
38 .seealso: PetscDSRegisterAll(), PetscDSRegisterDestroy()
39 
40 @*/
41 PetscErrorCode PetscDSRegister(const char sname[], PetscErrorCode (*function)(PetscDS))
42 {
43   PetscErrorCode ierr;
44 
45   PetscFunctionBegin;
46   ierr = PetscFunctionListAdd(&PetscDSList, sname, function);CHKERRQ(ierr);
47   PetscFunctionReturn(0);
48 }
49 
50 /*@C
51   PetscDSSetType - Builds a particular PetscDS
52 
53   Collective on PetscDS
54 
55   Input Parameters:
56 + prob - The PetscDS object
57 - name - The kind of system
58 
59   Options Database Key:
60 . -petscds_type <type> - Sets the PetscDS type; use -help for a list of available types
61 
62   Level: intermediate
63 
64 .keywords: PetscDS, set, type
65 .seealso: PetscDSGetType(), PetscDSCreate()
66 @*/
67 PetscErrorCode PetscDSSetType(PetscDS prob, PetscDSType name)
68 {
69   PetscErrorCode (*r)(PetscDS);
70   PetscBool      match;
71   PetscErrorCode ierr;
72 
73   PetscFunctionBegin;
74   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
75   ierr = PetscObjectTypeCompare((PetscObject) prob, name, &match);CHKERRQ(ierr);
76   if (match) PetscFunctionReturn(0);
77 
78   ierr = PetscDSRegisterAll();CHKERRQ(ierr);
79   ierr = PetscFunctionListFind(PetscDSList, name, &r);CHKERRQ(ierr);
80   if (!r) SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscDS type: %s", name);
81 
82   if (prob->ops->destroy) {
83     ierr             = (*prob->ops->destroy)(prob);CHKERRQ(ierr);
84     prob->ops->destroy = NULL;
85   }
86   ierr = (*r)(prob);CHKERRQ(ierr);
87   ierr = PetscObjectChangeTypeName((PetscObject) prob, name);CHKERRQ(ierr);
88   PetscFunctionReturn(0);
89 }
90 
91 /*@C
92   PetscDSGetType - Gets the PetscDS type name (as a string) from the object.
93 
94   Not Collective
95 
96   Input Parameter:
97 . prob  - The PetscDS
98 
99   Output Parameter:
100 . name - The PetscDS type name
101 
102   Level: intermediate
103 
104 .keywords: PetscDS, get, type, name
105 .seealso: PetscDSSetType(), PetscDSCreate()
106 @*/
107 PetscErrorCode PetscDSGetType(PetscDS prob, PetscDSType *name)
108 {
109   PetscErrorCode ierr;
110 
111   PetscFunctionBegin;
112   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
113   PetscValidPointer(name, 2);
114   ierr = PetscDSRegisterAll();CHKERRQ(ierr);
115   *name = ((PetscObject) prob)->type_name;
116   PetscFunctionReturn(0);
117 }
118 
119 static PetscErrorCode PetscDSView_Ascii(PetscDS prob, PetscViewer viewer)
120 {
121   PetscViewerFormat  format;
122   const PetscScalar *constants;
123   PetscInt           numConstants, f;
124   PetscErrorCode     ierr;
125 
126   PetscFunctionBegin;
127   ierr = PetscViewerGetFormat(viewer, &format);CHKERRQ(ierr);
128   ierr = PetscViewerASCIIPrintf(viewer, "Discrete System with %d fields\n", prob->Nf);CHKERRQ(ierr);
129   ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
130   for (f = 0; f < prob->Nf; ++f) {
131     PetscObject  obj;
132     PetscClassId id;
133     const char  *name;
134     PetscInt     Nc;
135 
136     ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr);
137     ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr);
138     ierr = PetscObjectGetName(obj, &name);CHKERRQ(ierr);
139     ierr = PetscViewerASCIIPrintf(viewer, "Field %s", name ? name : "<unknown>");CHKERRQ(ierr);
140     if (id == PETSCFE_CLASSID)      {
141       ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);
142       ierr = PetscViewerASCIIPrintf(viewer, " FEM");CHKERRQ(ierr);
143     } else if (id == PETSCFV_CLASSID) {
144       ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);
145       ierr = PetscViewerASCIIPrintf(viewer, " FVM");CHKERRQ(ierr);
146     }
147     else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", f);
148     if (Nc > 1) {ierr = PetscViewerASCIIPrintf(viewer, "%D components", Nc);CHKERRQ(ierr);}
149     else        {ierr = PetscViewerASCIIPrintf(viewer, "%D component ", Nc);CHKERRQ(ierr);}
150     if (prob->implicit[f]) {ierr = PetscViewerASCIIPrintf(viewer, " (implicit)");CHKERRQ(ierr);}
151     else                   {ierr = PetscViewerASCIIPrintf(viewer, " (explicit)");CHKERRQ(ierr);}
152     if (prob->adjacency[f*2+0]) {
153       if (prob->adjacency[f*2+1]) {ierr = PetscViewerASCIIPrintf(viewer, " (adj FVM++)");CHKERRQ(ierr);}
154       else                        {ierr = PetscViewerASCIIPrintf(viewer, " (adj FVM)");CHKERRQ(ierr);}
155     } else {
156       if (prob->adjacency[f*2+1]) {ierr = PetscViewerASCIIPrintf(viewer, " (adj FEM)");CHKERRQ(ierr);}
157       else                        {ierr = PetscViewerASCIIPrintf(viewer, " (adj FUNKY)");CHKERRQ(ierr);}
158     }
159     ierr = PetscViewerASCIIPrintf(viewer, "\n");CHKERRQ(ierr);
160     if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
161       if (id == PETSCFE_CLASSID)      {ierr = PetscFEView((PetscFE) obj, viewer);CHKERRQ(ierr);}
162       else if (id == PETSCFV_CLASSID) {ierr = PetscFVView((PetscFV) obj, viewer);CHKERRQ(ierr);}
163     }
164   }
165   ierr = PetscDSGetConstants(prob, &numConstants, &constants);CHKERRQ(ierr);
166   if (numConstants) {
167     ierr = PetscViewerASCIIPrintf(viewer, "%D constants\n", numConstants);CHKERRQ(ierr);
168     ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
169     for (f = 0; f < numConstants; ++f) {ierr = PetscViewerASCIIPrintf(viewer, "%g\n", (double) PetscRealPart(constants[f]));CHKERRQ(ierr);}
170     ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
171   }
172   ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
173   PetscFunctionReturn(0);
174 }
175 
176 /*@C
177   PetscDSView - Views a PetscDS
178 
179   Collective on PetscDS
180 
181   Input Parameter:
182 + prob - the PetscDS object to view
183 - v  - the viewer
184 
185   Level: developer
186 
187 .seealso PetscDSDestroy()
188 @*/
189 PetscErrorCode PetscDSView(PetscDS prob, PetscViewer v)
190 {
191   PetscBool      iascii;
192   PetscErrorCode ierr;
193 
194   PetscFunctionBegin;
195   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
196   if (!v) {ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) prob), &v);CHKERRQ(ierr);}
197   else    {PetscValidHeaderSpecific(v, PETSC_VIEWER_CLASSID, 2);}
198   ierr = PetscObjectTypeCompare((PetscObject) v, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr);
199   if (iascii) {ierr = PetscDSView_Ascii(prob, v);CHKERRQ(ierr);}
200   if (prob->ops->view) {ierr = (*prob->ops->view)(prob, v);CHKERRQ(ierr);}
201   PetscFunctionReturn(0);
202 }
203 
204 /*@
205   PetscDSSetFromOptions - sets parameters in a PetscDS from the options database
206 
207   Collective on PetscDS
208 
209   Input Parameter:
210 . prob - the PetscDS object to set options for
211 
212   Options Database:
213 
214   Level: developer
215 
216 .seealso PetscDSView()
217 @*/
218 PetscErrorCode PetscDSSetFromOptions(PetscDS prob)
219 {
220   DSBoundary     b;
221   const char    *defaultType;
222   char           name[256];
223   PetscBool      flg;
224   PetscErrorCode ierr;
225 
226   PetscFunctionBegin;
227   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
228   if (!((PetscObject) prob)->type_name) {
229     defaultType = PETSCDSBASIC;
230   } else {
231     defaultType = ((PetscObject) prob)->type_name;
232   }
233   ierr = PetscDSRegisterAll();CHKERRQ(ierr);
234 
235   ierr = PetscObjectOptionsBegin((PetscObject) prob);CHKERRQ(ierr);
236   for (b = prob->boundary; b; b = b->next) {
237     char       optname[1024];
238     PetscInt   ids[1024], len = 1024;
239     PetscBool  flg;
240 
241     ierr = PetscSNPrintf(optname, sizeof(optname), "-bc_%s", b->name);CHKERRQ(ierr);
242     ierr = PetscMemzero(ids, sizeof(ids));CHKERRQ(ierr);
243     ierr = PetscOptionsIntArray(optname, "List of boundary IDs", "", ids, &len, &flg);CHKERRQ(ierr);
244     if (flg) {
245       b->numids = len;
246       ierr = PetscFree(b->ids);CHKERRQ(ierr);
247       ierr = PetscMalloc1(len, &b->ids);CHKERRQ(ierr);
248       ierr = PetscMemcpy(b->ids, ids, len*sizeof(PetscInt));CHKERRQ(ierr);
249     }
250     len = 1024;
251     ierr = PetscSNPrintf(optname, sizeof(optname), "-bc_%s_comp", b->name);CHKERRQ(ierr);
252     ierr = PetscMemzero(ids, sizeof(ids));CHKERRQ(ierr);
253     ierr = PetscOptionsIntArray(optname, "List of boundary field components", "", ids, &len, &flg);CHKERRQ(ierr);
254     if (flg) {
255       b->numcomps = len;
256       ierr = PetscFree(b->comps);CHKERRQ(ierr);
257       ierr = PetscMalloc1(len, &b->comps);CHKERRQ(ierr);
258       ierr = PetscMemcpy(b->comps, ids, len*sizeof(PetscInt));CHKERRQ(ierr);
259     }
260   }
261   ierr = PetscOptionsFList("-petscds_type", "Discrete System", "PetscDSSetType", PetscDSList, defaultType, name, 256, &flg);CHKERRQ(ierr);
262   if (flg) {
263     ierr = PetscDSSetType(prob, name);CHKERRQ(ierr);
264   } else if (!((PetscObject) prob)->type_name) {
265     ierr = PetscDSSetType(prob, defaultType);CHKERRQ(ierr);
266   }
267   if (prob->ops->setfromoptions) {ierr = (*prob->ops->setfromoptions)(prob);CHKERRQ(ierr);}
268   /* process any options handlers added with PetscObjectAddOptionsHandler() */
269   ierr = PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) prob);CHKERRQ(ierr);
270   ierr = PetscOptionsEnd();CHKERRQ(ierr);
271   ierr = PetscDSViewFromOptions(prob, NULL, "-petscds_view");CHKERRQ(ierr);
272   PetscFunctionReturn(0);
273 }
274 
275 /*@C
276   PetscDSSetUp - Construct data structures for the PetscDS
277 
278   Collective on PetscDS
279 
280   Input Parameter:
281 . prob - the PetscDS object to setup
282 
283   Level: developer
284 
285 .seealso PetscDSView(), PetscDSDestroy()
286 @*/
287 PetscErrorCode PetscDSSetUp(PetscDS prob)
288 {
289   const PetscInt Nf = prob->Nf;
290   PetscInt       dim, work, NcMax = 0, NqMax = 0, f;
291   PetscErrorCode ierr;
292 
293   PetscFunctionBegin;
294   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
295   if (prob->setup) PetscFunctionReturn(0);
296   /* Calculate sizes */
297   ierr = PetscDSGetSpatialDimension(prob, &dim);CHKERRQ(ierr);
298   prob->totDim = prob->totComp = 0;
299   ierr = PetscMalloc2(Nf,&prob->Nc,Nf,&prob->Nb);CHKERRQ(ierr);
300   ierr = PetscCalloc2(Nf+1,&prob->off,Nf+1,&prob->offDer);CHKERRQ(ierr);
301   ierr = PetscMalloc4(Nf,&prob->basis,Nf,&prob->basisDer,Nf,&prob->basisFace,Nf,&prob->basisDerFace);CHKERRQ(ierr);
302   for (f = 0; f < Nf; ++f) {
303     PetscObject     obj;
304     PetscClassId    id;
305     PetscQuadrature q;
306     PetscInt        Nq = 0, Nb, Nc;
307 
308     ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr);
309     ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr);
310     if (id == PETSCFE_CLASSID)      {
311       PetscFE fe = (PetscFE) obj;
312 
313       ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr);
314       ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr);
315       ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr);
316       ierr = PetscFEGetDefaultTabulation(fe, &prob->basis[f], &prob->basisDer[f], NULL);CHKERRQ(ierr);
317       ierr = PetscFEGetFaceTabulation(fe, &prob->basisFace[f], &prob->basisDerFace[f], NULL);CHKERRQ(ierr);
318     } else if (id == PETSCFV_CLASSID) {
319       PetscFV fv = (PetscFV) obj;
320 
321       ierr = PetscFVGetQuadrature(fv, &q);CHKERRQ(ierr);
322       ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr);
323       Nb   = Nc;
324       ierr = PetscFVGetDefaultTabulation(fv, &prob->basis[f], &prob->basisDer[f], NULL);CHKERRQ(ierr);
325       /* TODO: should PetscFV also have face tabulation? Otherwise there will be a null pointer in prob->basisFace */
326     } else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
327     prob->Nc[f]       = Nc;
328     prob->Nb[f]       = Nb;
329     prob->off[f+1]    = Nc     + prob->off[f];
330     prob->offDer[f+1] = Nc*dim + prob->offDer[f];
331     if (q) {ierr = PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);CHKERRQ(ierr);}
332     NqMax          = PetscMax(NqMax, Nq);
333     NcMax          = PetscMax(NcMax, Nc);
334     prob->totDim  += Nb;
335     prob->totComp += Nc;
336   }
337   work = PetscMax(prob->totComp*dim, PetscSqr(NcMax*dim));
338   /* Allocate works space */
339   ierr = PetscMalloc5(prob->totComp,&prob->u,prob->totComp,&prob->u_t,prob->totComp*dim,&prob->u_x,dim,&prob->x,work,&prob->refSpaceDer);CHKERRQ(ierr);
340   ierr = PetscMalloc6(NqMax*NcMax,&prob->f0,NqMax*NcMax*dim,&prob->f1,NqMax*NcMax*NcMax,&prob->g0,NqMax*NcMax*NcMax*dim,&prob->g1,NqMax*NcMax*NcMax*dim,&prob->g2,NqMax*NcMax*NcMax*dim*dim,&prob->g3);CHKERRQ(ierr);
341   if (prob->ops->setup) {ierr = (*prob->ops->setup)(prob);CHKERRQ(ierr);}
342   prob->setup = PETSC_TRUE;
343   PetscFunctionReturn(0);
344 }
345 
346 static PetscErrorCode PetscDSDestroyStructs_Static(PetscDS prob)
347 {
348   PetscErrorCode ierr;
349 
350   PetscFunctionBegin;
351   ierr = PetscFree2(prob->Nc,prob->Nb);CHKERRQ(ierr);
352   ierr = PetscFree2(prob->off,prob->offDer);CHKERRQ(ierr);
353   ierr = PetscFree4(prob->basis,prob->basisDer,prob->basisFace,prob->basisDerFace);CHKERRQ(ierr);
354   ierr = PetscFree5(prob->u,prob->u_t,prob->u_x,prob->x,prob->refSpaceDer);CHKERRQ(ierr);
355   ierr = PetscFree6(prob->f0,prob->f1,prob->g0,prob->g1,prob->g2,prob->g3);CHKERRQ(ierr);
356   PetscFunctionReturn(0);
357 }
358 
359 static PetscErrorCode PetscDSEnlarge_Static(PetscDS prob, PetscInt NfNew)
360 {
361   PetscObject      *tmpd;
362   PetscBool        *tmpi, *tmpa;
363   PetscPointFunc   *tmpobj, *tmpf, *tmpup;
364   PetscPointJac    *tmpg, *tmpgp, *tmpgt;
365   PetscBdPointFunc *tmpfbd;
366   PetscBdPointJac  *tmpgbd;
367   PetscRiemannFunc *tmpr;
368   PetscSimplePointFunc *tmpexactSol;
369   void            **tmpctx;
370   PetscInt          Nf = prob->Nf, f, i;
371   PetscErrorCode    ierr;
372 
373   PetscFunctionBegin;
374   if (Nf >= NfNew) PetscFunctionReturn(0);
375   prob->setup = PETSC_FALSE;
376   ierr = PetscDSDestroyStructs_Static(prob);CHKERRQ(ierr);
377   ierr = PetscMalloc3(NfNew, &tmpd, NfNew, &tmpi, NfNew*2, &tmpa);CHKERRQ(ierr);
378   for (f = 0; f < Nf; ++f) {tmpd[f] = prob->disc[f]; tmpi[f] = prob->implicit[f]; for (i = 0; i < 2; ++i) tmpa[f*2+i] = prob->adjacency[f*2+i];}
379   for (f = Nf; f < NfNew; ++f) {tmpd[f] = NULL; tmpi[f] = PETSC_TRUE; tmpa[f*2+0] = PETSC_FALSE; tmpa[f*2+1] = PETSC_TRUE;}
380   ierr = PetscFree3(prob->disc, prob->implicit, prob->adjacency);CHKERRQ(ierr);
381   prob->Nf        = NfNew;
382   prob->disc      = tmpd;
383   prob->implicit  = tmpi;
384   prob->adjacency = tmpa;
385   ierr = PetscCalloc7(NfNew, &tmpobj, NfNew*2, &tmpf, NfNew*NfNew*4, &tmpg, NfNew*NfNew*4, &tmpgp, NfNew*NfNew*4, &tmpgt, NfNew, &tmpr, NfNew, &tmpctx);CHKERRQ(ierr);
386   ierr = PetscCalloc1(NfNew, &tmpup);CHKERRQ(ierr);
387   for (f = 0; f < Nf; ++f) tmpobj[f] = prob->obj[f];
388   for (f = 0; f < Nf*2; ++f) tmpf[f] = prob->f[f];
389   for (f = 0; f < Nf*Nf*4; ++f) tmpg[f] = prob->g[f];
390   for (f = 0; f < Nf*Nf*4; ++f) tmpgp[f] = prob->gp[f];
391   for (f = 0; f < Nf; ++f) tmpr[f] = prob->r[f];
392   for (f = 0; f < Nf; ++f) tmpup[f] = prob->update[f];
393   for (f = 0; f < Nf; ++f) tmpctx[f] = prob->ctx[f];
394   for (f = Nf; f < NfNew; ++f) tmpobj[f] = NULL;
395   for (f = Nf*2; f < NfNew*2; ++f) tmpf[f] = NULL;
396   for (f = Nf*Nf*4; f < NfNew*NfNew*4; ++f) tmpg[f] = NULL;
397   for (f = Nf*Nf*4; f < NfNew*NfNew*4; ++f) tmpgp[f] = NULL;
398   for (f = Nf*Nf*4; f < NfNew*NfNew*4; ++f) tmpgt[f] = NULL;
399   for (f = Nf; f < NfNew; ++f) tmpr[f] = NULL;
400   for (f = Nf; f < NfNew; ++f) tmpup[f] = NULL;
401   for (f = Nf; f < NfNew; ++f) tmpctx[f] = NULL;
402   ierr = PetscFree7(prob->obj, prob->f, prob->g, prob->gp, prob->gt, prob->r, prob->ctx);CHKERRQ(ierr);
403   ierr = PetscFree(prob->update);CHKERRQ(ierr);
404   prob->obj = tmpobj;
405   prob->f   = tmpf;
406   prob->g   = tmpg;
407   prob->gp  = tmpgp;
408   prob->gt  = tmpgt;
409   prob->r   = tmpr;
410   prob->update = tmpup;
411   prob->ctx = tmpctx;
412   ierr = PetscCalloc3(NfNew*2, &tmpfbd, NfNew*NfNew*4, &tmpgbd, NfNew, &tmpexactSol);CHKERRQ(ierr);
413   for (f = 0; f < Nf*2; ++f) tmpfbd[f] = prob->fBd[f];
414   for (f = 0; f < Nf*Nf*4; ++f) tmpgbd[f] = prob->gBd[f];
415   for (f = 0; f < Nf; ++f) tmpexactSol[f] = prob->exactSol[f];
416   for (f = Nf*2; f < NfNew*2; ++f) tmpfbd[f] = NULL;
417   for (f = Nf*Nf*4; f < NfNew*NfNew*4; ++f) tmpgbd[f] = NULL;
418   for (f = Nf; f < NfNew; ++f) tmpexactSol[f] = NULL;
419   ierr = PetscFree3(prob->fBd, prob->gBd, prob->exactSol);CHKERRQ(ierr);
420   prob->fBd = tmpfbd;
421   prob->gBd = tmpgbd;
422   prob->exactSol = tmpexactSol;
423   PetscFunctionReturn(0);
424 }
425 
426 /*@
427   PetscDSDestroy - Destroys a PetscDS object
428 
429   Collective on PetscDS
430 
431   Input Parameter:
432 . prob - the PetscDS object to destroy
433 
434   Level: developer
435 
436 .seealso PetscDSView()
437 @*/
438 PetscErrorCode PetscDSDestroy(PetscDS *prob)
439 {
440   PetscInt       f;
441   DSBoundary     next;
442   PetscErrorCode ierr;
443 
444   PetscFunctionBegin;
445   if (!*prob) PetscFunctionReturn(0);
446   PetscValidHeaderSpecific((*prob), PETSCDS_CLASSID, 1);
447 
448   if (--((PetscObject)(*prob))->refct > 0) {*prob = 0; PetscFunctionReturn(0);}
449   ((PetscObject) (*prob))->refct = 0;
450   ierr = PetscDSDestroyStructs_Static(*prob);CHKERRQ(ierr);
451   for (f = 0; f < (*prob)->Nf; ++f) {
452     ierr = PetscObjectDereference((*prob)->disc[f]);CHKERRQ(ierr);
453   }
454   ierr = PetscFree3((*prob)->disc, (*prob)->implicit, (*prob)->adjacency);CHKERRQ(ierr);
455   ierr = PetscFree7((*prob)->obj,(*prob)->f,(*prob)->g,(*prob)->gp,(*prob)->gt,(*prob)->r,(*prob)->ctx);CHKERRQ(ierr);
456   ierr = PetscFree((*prob)->update);CHKERRQ(ierr);
457   ierr = PetscFree3((*prob)->fBd,(*prob)->gBd,(*prob)->exactSol);CHKERRQ(ierr);
458   if ((*prob)->ops->destroy) {ierr = (*(*prob)->ops->destroy)(*prob);CHKERRQ(ierr);}
459   next = (*prob)->boundary;
460   while (next) {
461     DSBoundary b = next;
462 
463     next = b->next;
464     ierr = PetscFree(b->comps);CHKERRQ(ierr);
465     ierr = PetscFree(b->ids);CHKERRQ(ierr);
466     ierr = PetscFree(b->name);CHKERRQ(ierr);
467     ierr = PetscFree(b->labelname);CHKERRQ(ierr);
468     ierr = PetscFree(b);CHKERRQ(ierr);
469   }
470   ierr = PetscFree((*prob)->constants);CHKERRQ(ierr);
471   ierr = PetscHeaderDestroy(prob);CHKERRQ(ierr);
472   PetscFunctionReturn(0);
473 }
474 
475 /*@
476   PetscDSCreate - Creates an empty PetscDS object. The type can then be set with PetscDSSetType().
477 
478   Collective on MPI_Comm
479 
480   Input Parameter:
481 . comm - The communicator for the PetscDS object
482 
483   Output Parameter:
484 . prob - The PetscDS object
485 
486   Level: beginner
487 
488 .seealso: PetscDSSetType(), PETSCDSBASIC
489 @*/
490 PetscErrorCode PetscDSCreate(MPI_Comm comm, PetscDS *prob)
491 {
492   PetscDS   p;
493   PetscErrorCode ierr;
494 
495   PetscFunctionBegin;
496   PetscValidPointer(prob, 2);
497   *prob  = NULL;
498   ierr = PetscDSInitializePackage();CHKERRQ(ierr);
499 
500   ierr = PetscHeaderCreate(p, PETSCDS_CLASSID, "PetscDS", "Discrete System", "PetscDS", comm, PetscDSDestroy, PetscDSView);CHKERRQ(ierr);
501 
502   p->Nf    = 0;
503   p->setup = PETSC_FALSE;
504   p->numConstants = 0;
505   p->constants    = NULL;
506   p->dimEmbed     = -1;
507 
508   *prob = p;
509   PetscFunctionReturn(0);
510 }
511 
512 /*@
513   PetscDSGetNumFields - Returns the number of fields in the DS
514 
515   Not collective
516 
517   Input Parameter:
518 . prob - The PetscDS object
519 
520   Output Parameter:
521 . Nf - The number of fields
522 
523   Level: beginner
524 
525 .seealso: PetscDSGetSpatialDimension(), PetscDSCreate()
526 @*/
527 PetscErrorCode PetscDSGetNumFields(PetscDS prob, PetscInt *Nf)
528 {
529   PetscFunctionBegin;
530   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
531   PetscValidPointer(Nf, 2);
532   *Nf = prob->Nf;
533   PetscFunctionReturn(0);
534 }
535 
536 /*@
537   PetscDSGetSpatialDimension - Returns the spatial dimension of the DS, meaning the topological dimension of the discretizations
538 
539   Not collective
540 
541   Input Parameter:
542 . prob - The PetscDS object
543 
544   Output Parameter:
545 . dim - The spatial dimension
546 
547   Level: beginner
548 
549 .seealso: PetscDSGetCoordinateDimension(), PetscDSGetNumFields(), PetscDSCreate()
550 @*/
551 PetscErrorCode PetscDSGetSpatialDimension(PetscDS prob, PetscInt *dim)
552 {
553   PetscErrorCode ierr;
554 
555   PetscFunctionBegin;
556   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
557   PetscValidPointer(dim, 2);
558   *dim = 0;
559   if (prob->Nf) {
560     PetscObject  obj;
561     PetscClassId id;
562 
563     ierr = PetscDSGetDiscretization(prob, 0, &obj);CHKERRQ(ierr);
564     ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr);
565     if (id == PETSCFE_CLASSID)      {ierr = PetscFEGetSpatialDimension((PetscFE) obj, dim);CHKERRQ(ierr);}
566     else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetSpatialDimension((PetscFV) obj, dim);CHKERRQ(ierr);}
567     else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", 0);
568   }
569   PetscFunctionReturn(0);
570 }
571 
572 /*@
573   PetscDSGetCoordinateDimension - Returns the coordinate dimension of the DS, meaning the dimension of the space into which the discretiaztions are embedded
574 
575   Not collective
576 
577   Input Parameter:
578 . prob - The PetscDS object
579 
580   Output Parameter:
581 . dimEmbed - The coordinate dimension
582 
583   Level: beginner
584 
585 .seealso: PetscDSSetCoordinateDimension(), PetscDSGetSpatialDimension(), PetscDSGetNumFields(), PetscDSCreate()
586 @*/
587 PetscErrorCode PetscDSGetCoordinateDimension(PetscDS prob, PetscInt *dimEmbed)
588 {
589   PetscFunctionBegin;
590   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
591   PetscValidPointer(dimEmbed, 2);
592   if (prob->dimEmbed < 0) SETERRQ(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONGSTATE, "No coordinate dimension set for this DS");
593   *dimEmbed = prob->dimEmbed;
594   PetscFunctionReturn(0);
595 }
596 
597 /*@
598   PetscDSSetCoordinateDimension - Set the coordinate dimension of the DS, meaning the dimension of the space into which the discretiaztions are embedded
599 
600   Not collective
601 
602   Input Parameters:
603 + prob - The PetscDS object
604 - dimEmbed - The coordinate dimension
605 
606   Level: beginner
607 
608 .seealso: PetscDSGetCoordinateDimension(), PetscDSGetSpatialDimension(), PetscDSGetNumFields(), PetscDSCreate()
609 @*/
610 PetscErrorCode PetscDSSetCoordinateDimension(PetscDS prob, PetscInt dimEmbed)
611 {
612   PetscFunctionBegin;
613   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
614   prob->dimEmbed = dimEmbed;
615   PetscFunctionReturn(0);
616 }
617 
618 /*@
619   PetscDSGetTotalDimension - Returns the total size of the approximation space for this system
620 
621   Not collective
622 
623   Input Parameter:
624 . prob - The PetscDS object
625 
626   Output Parameter:
627 . dim - The total problem dimension
628 
629   Level: beginner
630 
631 .seealso: PetscDSGetNumFields(), PetscDSCreate()
632 @*/
633 PetscErrorCode PetscDSGetTotalDimension(PetscDS prob, PetscInt *dim)
634 {
635   PetscErrorCode ierr;
636 
637   PetscFunctionBegin;
638   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
639   ierr = PetscDSSetUp(prob);CHKERRQ(ierr);
640   PetscValidPointer(dim, 2);
641   *dim = prob->totDim;
642   PetscFunctionReturn(0);
643 }
644 
645 /*@
646   PetscDSGetTotalComponents - Returns the total number of components in this system
647 
648   Not collective
649 
650   Input Parameter:
651 . prob - The PetscDS object
652 
653   Output Parameter:
654 . dim - The total number of components
655 
656   Level: beginner
657 
658 .seealso: PetscDSGetNumFields(), PetscDSCreate()
659 @*/
660 PetscErrorCode PetscDSGetTotalComponents(PetscDS prob, PetscInt *Nc)
661 {
662   PetscErrorCode ierr;
663 
664   PetscFunctionBegin;
665   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
666   ierr = PetscDSSetUp(prob);CHKERRQ(ierr);
667   PetscValidPointer(Nc, 2);
668   *Nc = prob->totComp;
669   PetscFunctionReturn(0);
670 }
671 
672 /*@
673   PetscDSGetDiscretization - Returns the discretization object for the given field
674 
675   Not collective
676 
677   Input Parameters:
678 + prob - The PetscDS object
679 - f - The field number
680 
681   Output Parameter:
682 . disc - The discretization object
683 
684   Level: beginner
685 
686 .seealso: PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
687 @*/
688 PetscErrorCode PetscDSGetDiscretization(PetscDS prob, PetscInt f, PetscObject *disc)
689 {
690   PetscFunctionBegin;
691   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
692   PetscValidPointer(disc, 3);
693   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
694   *disc = prob->disc[f];
695   PetscFunctionReturn(0);
696 }
697 
698 /*@
699   PetscDSSetDiscretization - Sets the discretization object for the given field
700 
701   Not collective
702 
703   Input Parameters:
704 + prob - The PetscDS object
705 . f - The field number
706 - disc - The discretization object
707 
708   Level: beginner
709 
710 .seealso: PetscDSGetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
711 @*/
712 PetscErrorCode PetscDSSetDiscretization(PetscDS prob, PetscInt f, PetscObject disc)
713 {
714   PetscErrorCode ierr;
715 
716   PetscFunctionBegin;
717   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
718   PetscValidPointer(disc, 3);
719   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
720   ierr = PetscDSEnlarge_Static(prob, f+1);CHKERRQ(ierr);
721   if (prob->disc[f]) {ierr = PetscObjectDereference(prob->disc[f]);CHKERRQ(ierr);}
722   prob->disc[f] = disc;
723   ierr = PetscObjectReference(disc);CHKERRQ(ierr);
724   {
725     PetscClassId id;
726 
727     ierr = PetscObjectGetClassId(disc, &id);CHKERRQ(ierr);
728     if (id == PETSCFV_CLASSID) {
729       prob->implicit[f]      = PETSC_FALSE;
730       prob->adjacency[f*2+0] = PETSC_TRUE;
731       prob->adjacency[f*2+1] = PETSC_FALSE;
732     }
733   }
734   PetscFunctionReturn(0);
735 }
736 
737 /*@
738   PetscDSAddDiscretization - Adds a discretization object
739 
740   Not collective
741 
742   Input Parameters:
743 + prob - The PetscDS object
744 - disc - The boundary discretization object
745 
746   Level: beginner
747 
748 .seealso: PetscDSGetDiscretization(), PetscDSSetDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
749 @*/
750 PetscErrorCode PetscDSAddDiscretization(PetscDS prob, PetscObject disc)
751 {
752   PetscErrorCode ierr;
753 
754   PetscFunctionBegin;
755   ierr = PetscDSSetDiscretization(prob, prob->Nf, disc);CHKERRQ(ierr);
756   PetscFunctionReturn(0);
757 }
758 
759 /*@
760   PetscDSGetImplicit - Returns the flag for implicit solve for this field. This is just a guide for IMEX
761 
762   Not collective
763 
764   Input Parameters:
765 + prob - The PetscDS object
766 - f - The field number
767 
768   Output Parameter:
769 . implicit - The flag indicating what kind of solve to use for this field
770 
771   Level: developer
772 
773 .seealso: PetscDSSetImplicit(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
774 @*/
775 PetscErrorCode PetscDSGetImplicit(PetscDS prob, PetscInt f, PetscBool *implicit)
776 {
777   PetscFunctionBegin;
778   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
779   PetscValidPointer(implicit, 3);
780   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
781   *implicit = prob->implicit[f];
782   PetscFunctionReturn(0);
783 }
784 
785 /*@
786   PetscDSSetImplicit - Set the flag for implicit solve for this field. This is just a guide for IMEX
787 
788   Not collective
789 
790   Input Parameters:
791 + prob - The PetscDS object
792 . f - The field number
793 - implicit - The flag indicating what kind of solve to use for this field
794 
795   Level: developer
796 
797 .seealso: PetscDSGetImplicit(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
798 @*/
799 PetscErrorCode PetscDSSetImplicit(PetscDS prob, PetscInt f, PetscBool implicit)
800 {
801   PetscFunctionBegin;
802   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
803   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
804   prob->implicit[f] = implicit;
805   PetscFunctionReturn(0);
806 }
807 
808 /*@
809   PetscDSGetAdjacency - Returns the flags for determining variable influence
810 
811   Not collective
812 
813   Input Parameters:
814 + prob - The PetscDS object
815 - f - The field number
816 
817   Output Parameter:
818 + useCone    - Flag for variable influence starting with the cone operation
819 - useClosure - Flag for variable influence using transitive closure
820 
821   Note: See the discussion in DMPlexGetAdjacencyUseCone() and DMPlexGetAdjacencyUseClosure()
822 
823   Level: developer
824 
825 .seealso: PetscDSSetAdjacency(), DMPlexGetAdjacencyUseCone(), DMPlexGetAdjacencyUseClosure(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
826 @*/
827 PetscErrorCode PetscDSGetAdjacency(PetscDS prob, PetscInt f, PetscBool *useCone, PetscBool *useClosure)
828 {
829   PetscFunctionBegin;
830   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
831   PetscValidPointer(useCone, 3);
832   PetscValidPointer(useClosure, 4);
833   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
834   *useCone    = prob->adjacency[f*2+0];
835   *useClosure = prob->adjacency[f*2+1];
836   PetscFunctionReturn(0);
837 }
838 
839 /*@
840   PetscDSSetAdjacency - Set the flags for determining variable influence
841 
842   Not collective
843 
844   Input Parameters:
845 + prob - The PetscDS object
846 . f - The field number
847 . useCone    - Flag for variable influence starting with the cone operation
848 - useClosure - Flag for variable influence using transitive closure
849 
850   Note: See the discussion in DMPlexGetAdjacencyUseCone() and DMPlexGetAdjacencyUseClosure()
851 
852   Level: developer
853 
854 .seealso: PetscDSGetAdjacency(), DMPlexGetAdjacencyUseCone(), DMPlexGetAdjacencyUseClosure(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
855 @*/
856 PetscErrorCode PetscDSSetAdjacency(PetscDS prob, PetscInt f, PetscBool useCone, PetscBool useClosure)
857 {
858   PetscFunctionBegin;
859   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
860   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
861   prob->adjacency[f*2+0] = useCone;
862   prob->adjacency[f*2+1] = useClosure;
863   PetscFunctionReturn(0);
864 }
865 
866 PetscErrorCode PetscDSGetObjective(PetscDS prob, PetscInt f,
867                                    void (**obj)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
868                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
869                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
870                                                 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[]))
871 {
872   PetscFunctionBegin;
873   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
874   PetscValidPointer(obj, 2);
875   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
876   *obj = prob->obj[f];
877   PetscFunctionReturn(0);
878 }
879 
880 PetscErrorCode PetscDSSetObjective(PetscDS prob, PetscInt f,
881                                    void (*obj)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
882                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
883                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
884                                                PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[]))
885 {
886   PetscErrorCode ierr;
887 
888   PetscFunctionBegin;
889   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
890   if (obj) PetscValidFunction(obj, 2);
891   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
892   ierr = PetscDSEnlarge_Static(prob, f+1);CHKERRQ(ierr);
893   prob->obj[f] = obj;
894   PetscFunctionReturn(0);
895 }
896 
897 /*@C
898   PetscDSGetResidual - Get the pointwise residual function for a given test field
899 
900   Not collective
901 
902   Input Parameters:
903 + prob - The PetscDS
904 - f    - The test field number
905 
906   Output Parameters:
907 + f0 - integrand for the test function term
908 - f1 - integrand for the test function gradient term
909 
910   Note: We are using a first order FEM model for the weak form:
911 
912   \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)
913 
914 The calling sequence for the callbacks f0 and f1 is given by:
915 
916 $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
917 $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
918 $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
919 $    PetscReal t, const PetscReal x[], PetscScalar f0[])
920 
921 + dim - the spatial dimension
922 . Nf - the number of fields
923 . uOff - the offset into u[] and u_t[] for each field
924 . uOff_x - the offset into u_x[] for each field
925 . u - each field evaluated at the current point
926 . u_t - the time derivative of each field evaluated at the current point
927 . u_x - the gradient of each field evaluated at the current point
928 . aOff - the offset into a[] and a_t[] for each auxiliary field
929 . aOff_x - the offset into a_x[] for each auxiliary field
930 . a - each auxiliary field evaluated at the current point
931 . a_t - the time derivative of each auxiliary field evaluated at the current point
932 . a_x - the gradient of auxiliary each field evaluated at the current point
933 . t - current time
934 . x - coordinates of the current point
935 . numConstants - number of constant parameters
936 . constants - constant parameters
937 - f0 - output values at the current point
938 
939   Level: intermediate
940 
941 .seealso: PetscDSSetResidual()
942 @*/
943 PetscErrorCode PetscDSGetResidual(PetscDS prob, PetscInt f,
944                                   void (**f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
945                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
946                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
947                                               PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]),
948                                   void (**f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
949                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
950                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
951                                               PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
952 {
953   PetscFunctionBegin;
954   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
955   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
956   if (f0) {PetscValidPointer(f0, 3); *f0 = prob->f[f*2+0];}
957   if (f1) {PetscValidPointer(f1, 4); *f1 = prob->f[f*2+1];}
958   PetscFunctionReturn(0);
959 }
960 
961 /*@C
962   PetscDSSetResidual - Set the pointwise residual function for a given test field
963 
964   Not collective
965 
966   Input Parameters:
967 + prob - The PetscDS
968 . f    - The test field number
969 . f0 - integrand for the test function term
970 - f1 - integrand for the test function gradient term
971 
972   Note: We are using a first order FEM model for the weak form:
973 
974   \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)
975 
976 The calling sequence for the callbacks f0 and f1 is given by:
977 
978 $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
979 $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
980 $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
981 $    PetscReal t, const PetscReal x[], PetscScalar f0[])
982 
983 + dim - the spatial dimension
984 . Nf - the number of fields
985 . uOff - the offset into u[] and u_t[] for each field
986 . uOff_x - the offset into u_x[] for each field
987 . u - each field evaluated at the current point
988 . u_t - the time derivative of each field evaluated at the current point
989 . u_x - the gradient of each field evaluated at the current point
990 . aOff - the offset into a[] and a_t[] for each auxiliary field
991 . aOff_x - the offset into a_x[] for each auxiliary field
992 . a - each auxiliary field evaluated at the current point
993 . a_t - the time derivative of each auxiliary field evaluated at the current point
994 . a_x - the gradient of auxiliary each field evaluated at the current point
995 . t - current time
996 . x - coordinates of the current point
997 . numConstants - number of constant parameters
998 . constants - constant parameters
999 - f0 - output values at the current point
1000 
1001   Level: intermediate
1002 
1003 .seealso: PetscDSGetResidual()
1004 @*/
1005 PetscErrorCode PetscDSSetResidual(PetscDS prob, PetscInt f,
1006                                   void (*f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1007                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1008                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1009                                              PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]),
1010                                   void (*f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1011                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1012                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1013                                              PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
1014 {
1015   PetscErrorCode ierr;
1016 
1017   PetscFunctionBegin;
1018   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
1019   if (f0) PetscValidFunction(f0, 3);
1020   if (f1) PetscValidFunction(f1, 4);
1021   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1022   ierr = PetscDSEnlarge_Static(prob, f+1);CHKERRQ(ierr);
1023   prob->f[f*2+0] = f0;
1024   prob->f[f*2+1] = f1;
1025   PetscFunctionReturn(0);
1026 }
1027 
1028 /*@C
1029   PetscDSHasJacobian - Signals that Jacobian functions have been set
1030 
1031   Not collective
1032 
1033   Input Parameter:
1034 . prob - The PetscDS
1035 
1036   Output Parameter:
1037 . hasJac - flag that pointwise function for the Jacobian has been set
1038 
1039   Level: intermediate
1040 
1041 .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1042 @*/
1043 PetscErrorCode PetscDSHasJacobian(PetscDS prob, PetscBool *hasJac)
1044 {
1045   PetscInt f, g, h;
1046 
1047   PetscFunctionBegin;
1048   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
1049   *hasJac = PETSC_FALSE;
1050   for (f = 0; f < prob->Nf; ++f) {
1051     for (g = 0; g < prob->Nf; ++g) {
1052       for (h = 0; h < 4; ++h) {
1053         if (prob->g[(f*prob->Nf + g)*4+h]) *hasJac = PETSC_TRUE;
1054       }
1055     }
1056   }
1057   PetscFunctionReturn(0);
1058 }
1059 
1060 /*@C
1061   PetscDSGetJacobian - Get the pointwise Jacobian function for given test and basis field
1062 
1063   Not collective
1064 
1065   Input Parameters:
1066 + prob - The PetscDS
1067 . f    - The test field number
1068 - g    - The field number
1069 
1070   Output Parameters:
1071 + g0 - integrand for the test and basis function term
1072 . g1 - integrand for the test function and basis function gradient term
1073 . g2 - integrand for the test function gradient and basis function term
1074 - g3 - integrand for the test function gradient and basis function gradient term
1075 
1076   Note: We are using a first order FEM model for the weak form:
1077 
1078   \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
1079 
1080 The calling sequence for the callbacks g0, g1, g2 and g3 is given by:
1081 
1082 $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1083 $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1084 $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1085 $    PetscReal t, const PetscReal u_tShift, const PetscReal x[], PetscScalar g0[])
1086 
1087 + dim - the spatial dimension
1088 . Nf - the number of fields
1089 . uOff - the offset into u[] and u_t[] for each field
1090 . uOff_x - the offset into u_x[] for each field
1091 . u - each field evaluated at the current point
1092 . u_t - the time derivative of each field evaluated at the current point
1093 . u_x - the gradient of each field evaluated at the current point
1094 . aOff - the offset into a[] and a_t[] for each auxiliary field
1095 . aOff_x - the offset into a_x[] for each auxiliary field
1096 . a - each auxiliary field evaluated at the current point
1097 . a_t - the time derivative of each auxiliary field evaluated at the current point
1098 . a_x - the gradient of auxiliary each field evaluated at the current point
1099 . t - current time
1100 . u_tShift - the multiplier a for dF/dU_t
1101 . x - coordinates of the current point
1102 . numConstants - number of constant parameters
1103 . constants - constant parameters
1104 - g0 - output values at the current point
1105 
1106   Level: intermediate
1107 
1108 .seealso: PetscDSSetJacobian()
1109 @*/
1110 PetscErrorCode PetscDSGetJacobian(PetscDS prob, PetscInt f, PetscInt g,
1111                                   void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1112                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1113                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1114                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1115                                   void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1116                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1117                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1118                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1119                                   void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1120                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1121                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1122                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1123                                   void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1124                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1125                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1126                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1127 {
1128   PetscFunctionBegin;
1129   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
1130   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
1131   if ((g < 0) || (g >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, prob->Nf);
1132   if (g0) {PetscValidPointer(g0, 4); *g0 = prob->g[(f*prob->Nf + g)*4+0];}
1133   if (g1) {PetscValidPointer(g1, 5); *g1 = prob->g[(f*prob->Nf + g)*4+1];}
1134   if (g2) {PetscValidPointer(g2, 6); *g2 = prob->g[(f*prob->Nf + g)*4+2];}
1135   if (g3) {PetscValidPointer(g3, 7); *g3 = prob->g[(f*prob->Nf + g)*4+3];}
1136   PetscFunctionReturn(0);
1137 }
1138 
1139 /*@C
1140   PetscDSSetJacobian - Set the pointwise Jacobian function for given test and basis fields
1141 
1142   Not collective
1143 
1144   Input Parameters:
1145 + prob - The PetscDS
1146 . f    - The test field number
1147 . g    - The field number
1148 . g0 - integrand for the test and basis function term
1149 . g1 - integrand for the test function and basis function gradient term
1150 . g2 - integrand for the test function gradient and basis function term
1151 - g3 - integrand for the test function gradient and basis function gradient term
1152 
1153   Note: We are using a first order FEM model for the weak form:
1154 
1155   \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
1156 
1157 The calling sequence for the callbacks g0, g1, g2 and g3 is given by:
1158 
1159 $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1160 $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1161 $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1162 $    PetscReal t, const PetscReal x[], PetscScalar g0[])
1163 
1164 + dim - the spatial dimension
1165 . Nf - the number of fields
1166 . uOff - the offset into u[] and u_t[] for each field
1167 . uOff_x - the offset into u_x[] for each field
1168 . u - each field evaluated at the current point
1169 . u_t - the time derivative of each field evaluated at the current point
1170 . u_x - the gradient of each field evaluated at the current point
1171 . aOff - the offset into a[] and a_t[] for each auxiliary field
1172 . aOff_x - the offset into a_x[] for each auxiliary field
1173 . a - each auxiliary field evaluated at the current point
1174 . a_t - the time derivative of each auxiliary field evaluated at the current point
1175 . a_x - the gradient of auxiliary each field evaluated at the current point
1176 . t - current time
1177 . u_tShift - the multiplier a for dF/dU_t
1178 . x - coordinates of the current point
1179 . numConstants - number of constant parameters
1180 . constants - constant parameters
1181 - g0 - output values at the current point
1182 
1183   Level: intermediate
1184 
1185 .seealso: PetscDSGetJacobian()
1186 @*/
1187 PetscErrorCode PetscDSSetJacobian(PetscDS prob, PetscInt f, PetscInt g,
1188                                   void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1189                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1190                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1191                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1192                                   void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1193                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1194                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1195                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1196                                   void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1197                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1198                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1199                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1200                                   void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1201                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1202                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1203                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1204 {
1205   PetscErrorCode ierr;
1206 
1207   PetscFunctionBegin;
1208   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
1209   if (g0) PetscValidFunction(g0, 4);
1210   if (g1) PetscValidFunction(g1, 5);
1211   if (g2) PetscValidFunction(g2, 6);
1212   if (g3) PetscValidFunction(g3, 7);
1213   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1214   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
1215   ierr = PetscDSEnlarge_Static(prob, PetscMax(f, g)+1);CHKERRQ(ierr);
1216   prob->g[(f*prob->Nf + g)*4+0] = g0;
1217   prob->g[(f*prob->Nf + g)*4+1] = g1;
1218   prob->g[(f*prob->Nf + g)*4+2] = g2;
1219   prob->g[(f*prob->Nf + g)*4+3] = g3;
1220   PetscFunctionReturn(0);
1221 }
1222 
1223 /*@C
1224   PetscDSHasJacobianPreconditioner - Signals that a Jacobian preconditioner matrix has been set
1225 
1226   Not collective
1227 
1228   Input Parameter:
1229 . prob - The PetscDS
1230 
1231   Output Parameter:
1232 . hasJacPre - flag that pointwise function for Jacobian preconditioner matrix has been set
1233 
1234   Level: intermediate
1235 
1236 .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1237 @*/
1238 PetscErrorCode PetscDSHasJacobianPreconditioner(PetscDS prob, PetscBool *hasJacPre)
1239 {
1240   PetscInt f, g, h;
1241 
1242   PetscFunctionBegin;
1243   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
1244   *hasJacPre = PETSC_FALSE;
1245   for (f = 0; f < prob->Nf; ++f) {
1246     for (g = 0; g < prob->Nf; ++g) {
1247       for (h = 0; h < 4; ++h) {
1248         if (prob->gp[(f*prob->Nf + g)*4+h]) *hasJacPre = PETSC_TRUE;
1249       }
1250     }
1251   }
1252   PetscFunctionReturn(0);
1253 }
1254 
1255 /*@C
1256   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.
1257 
1258   Not collective
1259 
1260   Input Parameters:
1261 + prob - The PetscDS
1262 . f    - The test field number
1263 - g    - The field number
1264 
1265   Output Parameters:
1266 + g0 - integrand for the test and basis function term
1267 . g1 - integrand for the test function and basis function gradient term
1268 . g2 - integrand for the test function gradient and basis function term
1269 - g3 - integrand for the test function gradient and basis function gradient term
1270 
1271   Note: We are using a first order FEM model for the weak form:
1272 
1273   \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
1274 
1275 The calling sequence for the callbacks g0, g1, g2 and g3 is given by:
1276 
1277 $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1278 $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1279 $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1280 $    PetscReal t, const PetscReal u_tShift, const PetscReal x[], PetscScalar g0[])
1281 
1282 + dim - the spatial dimension
1283 . Nf - the number of fields
1284 . uOff - the offset into u[] and u_t[] for each field
1285 . uOff_x - the offset into u_x[] for each field
1286 . u - each field evaluated at the current point
1287 . u_t - the time derivative of each field evaluated at the current point
1288 . u_x - the gradient of each field evaluated at the current point
1289 . aOff - the offset into a[] and a_t[] for each auxiliary field
1290 . aOff_x - the offset into a_x[] for each auxiliary field
1291 . a - each auxiliary field evaluated at the current point
1292 . a_t - the time derivative of each auxiliary field evaluated at the current point
1293 . a_x - the gradient of auxiliary each field evaluated at the current point
1294 . t - current time
1295 . u_tShift - the multiplier a for dF/dU_t
1296 . x - coordinates of the current point
1297 . numConstants - number of constant parameters
1298 . constants - constant parameters
1299 - g0 - output values at the current point
1300 
1301   Level: intermediate
1302 
1303 .seealso: PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1304 @*/
1305 PetscErrorCode PetscDSGetJacobianPreconditioner(PetscDS prob, PetscInt f, PetscInt g,
1306                                   void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1307                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1308                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1309                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1310                                   void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1311                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1312                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1313                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1314                                   void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1315                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1316                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1317                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1318                                   void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1319                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1320                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1321                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1322 {
1323   PetscFunctionBegin;
1324   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
1325   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
1326   if ((g < 0) || (g >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, prob->Nf);
1327   if (g0) {PetscValidPointer(g0, 4); *g0 = prob->gp[(f*prob->Nf + g)*4+0];}
1328   if (g1) {PetscValidPointer(g1, 5); *g1 = prob->gp[(f*prob->Nf + g)*4+1];}
1329   if (g2) {PetscValidPointer(g2, 6); *g2 = prob->gp[(f*prob->Nf + g)*4+2];}
1330   if (g3) {PetscValidPointer(g3, 7); *g3 = prob->gp[(f*prob->Nf + g)*4+3];}
1331   PetscFunctionReturn(0);
1332 }
1333 
1334 /*@C
1335   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.
1336 
1337   Not collective
1338 
1339   Input Parameters:
1340 + prob - The PetscDS
1341 . f    - The test field number
1342 . g    - The field number
1343 . g0 - integrand for the test and basis function term
1344 . g1 - integrand for the test function and basis function gradient term
1345 . g2 - integrand for the test function gradient and basis function term
1346 - g3 - integrand for the test function gradient and basis function gradient term
1347 
1348   Note: We are using a first order FEM model for the weak form:
1349 
1350   \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
1351 
1352 The calling sequence for the callbacks g0, g1, g2 and g3 is given by:
1353 
1354 $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1355 $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1356 $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1357 $    PetscReal t, const PetscReal x[], PetscScalar g0[])
1358 
1359 + dim - the spatial dimension
1360 . Nf - the number of fields
1361 . uOff - the offset into u[] and u_t[] for each field
1362 . uOff_x - the offset into u_x[] for each field
1363 . u - each field evaluated at the current point
1364 . u_t - the time derivative of each field evaluated at the current point
1365 . u_x - the gradient of each field evaluated at the current point
1366 . aOff - the offset into a[] and a_t[] for each auxiliary field
1367 . aOff_x - the offset into a_x[] for each auxiliary field
1368 . a - each auxiliary field evaluated at the current point
1369 . a_t - the time derivative of each auxiliary field evaluated at the current point
1370 . a_x - the gradient of auxiliary each field evaluated at the current point
1371 . t - current time
1372 . u_tShift - the multiplier a for dF/dU_t
1373 . x - coordinates of the current point
1374 . numConstants - number of constant parameters
1375 . constants - constant parameters
1376 - g0 - output values at the current point
1377 
1378   Level: intermediate
1379 
1380 .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobian()
1381 @*/
1382 PetscErrorCode PetscDSSetJacobianPreconditioner(PetscDS prob, PetscInt f, PetscInt g,
1383                                   void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1384                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1385                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1386                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1387                                   void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1388                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1389                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1390                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1391                                   void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1392                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1393                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1394                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1395                                   void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1396                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1397                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1398                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1399 {
1400   PetscErrorCode ierr;
1401 
1402   PetscFunctionBegin;
1403   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
1404   if (g0) PetscValidFunction(g0, 4);
1405   if (g1) PetscValidFunction(g1, 5);
1406   if (g2) PetscValidFunction(g2, 6);
1407   if (g3) PetscValidFunction(g3, 7);
1408   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1409   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
1410   ierr = PetscDSEnlarge_Static(prob, PetscMax(f, g)+1);CHKERRQ(ierr);
1411   prob->gp[(f*prob->Nf + g)*4+0] = g0;
1412   prob->gp[(f*prob->Nf + g)*4+1] = g1;
1413   prob->gp[(f*prob->Nf + g)*4+2] = g2;
1414   prob->gp[(f*prob->Nf + g)*4+3] = g3;
1415   PetscFunctionReturn(0);
1416 }
1417 
1418 /*@C
1419   PetscDSHasDynamicJacobian - Signals that a dynamic Jacobian, dF/du_t, has been set
1420 
1421   Not collective
1422 
1423   Input Parameter:
1424 . prob - The PetscDS
1425 
1426   Output Parameter:
1427 . hasDynJac - flag that pointwise function for dynamic Jacobian has been set
1428 
1429   Level: intermediate
1430 
1431 .seealso: PetscDSGetDynamicJacobian(), PetscDSSetDynamicJacobian(), PetscDSGetJacobian()
1432 @*/
1433 PetscErrorCode PetscDSHasDynamicJacobian(PetscDS prob, PetscBool *hasDynJac)
1434 {
1435   PetscInt f, g, h;
1436 
1437   PetscFunctionBegin;
1438   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
1439   *hasDynJac = PETSC_FALSE;
1440   for (f = 0; f < prob->Nf; ++f) {
1441     for (g = 0; g < prob->Nf; ++g) {
1442       for (h = 0; h < 4; ++h) {
1443         if (prob->gt[(f*prob->Nf + g)*4+h]) *hasDynJac = PETSC_TRUE;
1444       }
1445     }
1446   }
1447   PetscFunctionReturn(0);
1448 }
1449 
1450 /*@C
1451   PetscDSGetDynamicJacobian - Get the pointwise dynamic Jacobian, dF/du_t, function for given test and basis field
1452 
1453   Not collective
1454 
1455   Input Parameters:
1456 + prob - The PetscDS
1457 . f    - The test field number
1458 - g    - The field number
1459 
1460   Output Parameters:
1461 + g0 - integrand for the test and basis function term
1462 . g1 - integrand for the test function and basis function gradient term
1463 . g2 - integrand for the test function gradient and basis function term
1464 - g3 - integrand for the test function gradient and basis function gradient term
1465 
1466   Note: We are using a first order FEM model for the weak form:
1467 
1468   \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
1469 
1470 The calling sequence for the callbacks g0, g1, g2 and g3 is given by:
1471 
1472 $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1473 $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1474 $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1475 $    PetscReal t, const PetscReal u_tShift, const PetscReal x[], PetscScalar g0[])
1476 
1477 + dim - the spatial dimension
1478 . Nf - the number of fields
1479 . uOff - the offset into u[] and u_t[] for each field
1480 . uOff_x - the offset into u_x[] for each field
1481 . u - each field evaluated at the current point
1482 . u_t - the time derivative of each field evaluated at the current point
1483 . u_x - the gradient of each field evaluated at the current point
1484 . aOff - the offset into a[] and a_t[] for each auxiliary field
1485 . aOff_x - the offset into a_x[] for each auxiliary field
1486 . a - each auxiliary field evaluated at the current point
1487 . a_t - the time derivative of each auxiliary field evaluated at the current point
1488 . a_x - the gradient of auxiliary each field evaluated at the current point
1489 . t - current time
1490 . u_tShift - the multiplier a for dF/dU_t
1491 . x - coordinates of the current point
1492 . numConstants - number of constant parameters
1493 . constants - constant parameters
1494 - g0 - output values at the current point
1495 
1496   Level: intermediate
1497 
1498 .seealso: PetscDSSetJacobian()
1499 @*/
1500 PetscErrorCode PetscDSGetDynamicJacobian(PetscDS prob, PetscInt f, PetscInt g,
1501                                          void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1502                                                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1503                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1504                                                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1505                                          void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1506                                                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1507                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1508                                                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1509                                          void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1510                                                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1511                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1512                                                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1513                                          void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1514                                                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1515                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1516                                                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1517 {
1518   PetscFunctionBegin;
1519   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
1520   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
1521   if ((g < 0) || (g >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, prob->Nf);
1522   if (g0) {PetscValidPointer(g0, 4); *g0 = prob->gt[(f*prob->Nf + g)*4+0];}
1523   if (g1) {PetscValidPointer(g1, 5); *g1 = prob->gt[(f*prob->Nf + g)*4+1];}
1524   if (g2) {PetscValidPointer(g2, 6); *g2 = prob->gt[(f*prob->Nf + g)*4+2];}
1525   if (g3) {PetscValidPointer(g3, 7); *g3 = prob->gt[(f*prob->Nf + g)*4+3];}
1526   PetscFunctionReturn(0);
1527 }
1528 
1529 /*@C
1530   PetscDSSetDynamicJacobian - Set the pointwise dynamic Jacobian, dF/du_t, function for given test and basis fields
1531 
1532   Not collective
1533 
1534   Input Parameters:
1535 + prob - The PetscDS
1536 . f    - The test field number
1537 . g    - The field number
1538 . g0 - integrand for the test and basis function term
1539 . g1 - integrand for the test function and basis function gradient term
1540 . g2 - integrand for the test function gradient and basis function term
1541 - g3 - integrand for the test function gradient and basis function gradient term
1542 
1543   Note: We are using a first order FEM model for the weak form:
1544 
1545   \int_\Omega \phi g_0(u, u_t, \nabla u, x, t) \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi + \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
1546 
1547 The calling sequence for the callbacks g0, g1, g2 and g3 is given by:
1548 
1549 $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1550 $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1551 $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1552 $    PetscReal t, const PetscReal x[], PetscScalar g0[])
1553 
1554 + dim - the spatial dimension
1555 . Nf - the number of fields
1556 . uOff - the offset into u[] and u_t[] for each field
1557 . uOff_x - the offset into u_x[] for each field
1558 . u - each field evaluated at the current point
1559 . u_t - the time derivative of each field evaluated at the current point
1560 . u_x - the gradient of each field evaluated at the current point
1561 . aOff - the offset into a[] and a_t[] for each auxiliary field
1562 . aOff_x - the offset into a_x[] for each auxiliary field
1563 . a - each auxiliary field evaluated at the current point
1564 . a_t - the time derivative of each auxiliary field evaluated at the current point
1565 . a_x - the gradient of auxiliary each field evaluated at the current point
1566 . t - current time
1567 . u_tShift - the multiplier a for dF/dU_t
1568 . x - coordinates of the current point
1569 . numConstants - number of constant parameters
1570 . constants - constant parameters
1571 - g0 - output values at the current point
1572 
1573   Level: intermediate
1574 
1575 .seealso: PetscDSGetJacobian()
1576 @*/
1577 PetscErrorCode PetscDSSetDynamicJacobian(PetscDS prob, PetscInt f, PetscInt g,
1578                                          void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1579                                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1580                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1581                                                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1582                                          void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1583                                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1584                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1585                                                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1586                                          void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1587                                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1588                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1589                                                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1590                                          void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1591                                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1592                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1593                                                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1594 {
1595   PetscErrorCode ierr;
1596 
1597   PetscFunctionBegin;
1598   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
1599   if (g0) PetscValidFunction(g0, 4);
1600   if (g1) PetscValidFunction(g1, 5);
1601   if (g2) PetscValidFunction(g2, 6);
1602   if (g3) PetscValidFunction(g3, 7);
1603   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1604   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
1605   ierr = PetscDSEnlarge_Static(prob, PetscMax(f, g)+1);CHKERRQ(ierr);
1606   prob->gt[(f*prob->Nf + g)*4+0] = g0;
1607   prob->gt[(f*prob->Nf + g)*4+1] = g1;
1608   prob->gt[(f*prob->Nf + g)*4+2] = g2;
1609   prob->gt[(f*prob->Nf + g)*4+3] = g3;
1610   PetscFunctionReturn(0);
1611 }
1612 
1613 /*@C
1614   PetscDSGetRiemannSolver - Returns the Riemann solver for the given field
1615 
1616   Not collective
1617 
1618   Input Arguments:
1619 + prob - The PetscDS object
1620 - f    - The field number
1621 
1622   Output Argument:
1623 . r    - Riemann solver
1624 
1625   Calling sequence for r:
1626 
1627 $ r(PetscInt dim, PetscInt Nf, const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscScalar flux[], void *ctx)
1628 
1629 + dim  - The spatial dimension
1630 . Nf   - The number of fields
1631 . x    - The coordinates at a point on the interface
1632 . n    - The normal vector to the interface
1633 . uL   - The state vector to the left of the interface
1634 . uR   - The state vector to the right of the interface
1635 . flux - output array of flux through the interface
1636 . numConstants - number of constant parameters
1637 . constants - constant parameters
1638 - ctx  - optional user context
1639 
1640   Level: intermediate
1641 
1642 .seealso: PetscDSSetRiemannSolver()
1643 @*/
1644 PetscErrorCode PetscDSGetRiemannSolver(PetscDS prob, PetscInt f,
1645                                        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))
1646 {
1647   PetscFunctionBegin;
1648   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
1649   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
1650   PetscValidPointer(r, 3);
1651   *r = prob->r[f];
1652   PetscFunctionReturn(0);
1653 }
1654 
1655 /*@C
1656   PetscDSSetRiemannSolver - Sets the Riemann solver for the given field
1657 
1658   Not collective
1659 
1660   Input Arguments:
1661 + prob - The PetscDS object
1662 . f    - The field number
1663 - r    - Riemann solver
1664 
1665   Calling sequence for r:
1666 
1667 $ r(PetscInt dim, PetscInt Nf, const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscScalar flux[], void *ctx)
1668 
1669 + dim  - The spatial dimension
1670 . Nf   - The number of fields
1671 . x    - The coordinates at a point on the interface
1672 . n    - The normal vector to the interface
1673 . uL   - The state vector to the left of the interface
1674 . uR   - The state vector to the right of the interface
1675 . flux - output array of flux through the interface
1676 . numConstants - number of constant parameters
1677 . constants - constant parameters
1678 - ctx  - optional user context
1679 
1680   Level: intermediate
1681 
1682 .seealso: PetscDSGetRiemannSolver()
1683 @*/
1684 PetscErrorCode PetscDSSetRiemannSolver(PetscDS prob, PetscInt f,
1685                                        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))
1686 {
1687   PetscErrorCode ierr;
1688 
1689   PetscFunctionBegin;
1690   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
1691   if (r) PetscValidFunction(r, 3);
1692   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1693   ierr = PetscDSEnlarge_Static(prob, f+1);CHKERRQ(ierr);
1694   prob->r[f] = r;
1695   PetscFunctionReturn(0);
1696 }
1697 
1698 /*@C
1699   PetscDSGetUpdate - Get the pointwise update function for a given field
1700 
1701   Not collective
1702 
1703   Input Parameters:
1704 + prob - The PetscDS
1705 - f    - The field number
1706 
1707   Output Parameters:
1708 + update - update function
1709 
1710   Note: The calling sequence for the callback update is given by:
1711 
1712 $ update(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1713 $        const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1714 $        const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1715 $        PetscReal t, const PetscReal x[], PetscScalar uNew[])
1716 
1717 + dim - the spatial dimension
1718 . Nf - the number of fields
1719 . uOff - the offset into u[] and u_t[] for each field
1720 . uOff_x - the offset into u_x[] for each field
1721 . u - each field evaluated at the current point
1722 . u_t - the time derivative of each field evaluated at the current point
1723 . u_x - the gradient of each field evaluated at the current point
1724 . aOff - the offset into a[] and a_t[] for each auxiliary field
1725 . aOff_x - the offset into a_x[] for each auxiliary field
1726 . a - each auxiliary field evaluated at the current point
1727 . a_t - the time derivative of each auxiliary field evaluated at the current point
1728 . a_x - the gradient of auxiliary each field evaluated at the current point
1729 . t - current time
1730 . x - coordinates of the current point
1731 - uNew - new value for field at the current point
1732 
1733   Level: intermediate
1734 
1735 .seealso: PetscDSSetUpdate(), PetscDSSetResidual()
1736 @*/
1737 PetscErrorCode PetscDSGetUpdate(PetscDS prob, PetscInt f,
1738                                   void (**update)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1739                                                   const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1740                                                   const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1741                                                   PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar uNew[]))
1742 {
1743   PetscFunctionBegin;
1744   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
1745   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
1746   if (update) {PetscValidPointer(update, 3); *update = prob->update[f];}
1747   PetscFunctionReturn(0);
1748 }
1749 
1750 /*@C
1751   PetscDSSetUpdate - Set the pointwise update function for a given field
1752 
1753   Not collective
1754 
1755   Input Parameters:
1756 + prob   - The PetscDS
1757 . f      - The field number
1758 - update - update function
1759 
1760   Note: The calling sequence for the callback update is given by:
1761 
1762 $ update(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1763 $        const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1764 $        const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1765 $        PetscReal t, const PetscReal x[], PetscScalar uNew[])
1766 
1767 + dim - the spatial dimension
1768 . Nf - the number of fields
1769 . uOff - the offset into u[] and u_t[] for each field
1770 . uOff_x - the offset into u_x[] for each field
1771 . u - each field evaluated at the current point
1772 . u_t - the time derivative of each field evaluated at the current point
1773 . u_x - the gradient of each field evaluated at the current point
1774 . aOff - the offset into a[] and a_t[] for each auxiliary field
1775 . aOff_x - the offset into a_x[] for each auxiliary field
1776 . a - each auxiliary field evaluated at the current point
1777 . a_t - the time derivative of each auxiliary field evaluated at the current point
1778 . a_x - the gradient of auxiliary each field evaluated at the current point
1779 . t - current time
1780 . x - coordinates of the current point
1781 - uNew - new field values at the current point
1782 
1783   Level: intermediate
1784 
1785 .seealso: PetscDSGetResidual()
1786 @*/
1787 PetscErrorCode PetscDSSetUpdate(PetscDS prob, PetscInt f,
1788                                 void (*update)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1789                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1790                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1791                                                PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar uNew[]))
1792 {
1793   PetscErrorCode ierr;
1794 
1795   PetscFunctionBegin;
1796   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
1797   if (update) PetscValidFunction(update, 3);
1798   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1799   ierr = PetscDSEnlarge_Static(prob, f+1);CHKERRQ(ierr);
1800   prob->update[f] = update;
1801   PetscFunctionReturn(0);
1802 }
1803 
1804 PetscErrorCode PetscDSGetContext(PetscDS prob, PetscInt f, void **ctx)
1805 {
1806   PetscFunctionBegin;
1807   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
1808   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
1809   PetscValidPointer(ctx, 3);
1810   *ctx = prob->ctx[f];
1811   PetscFunctionReturn(0);
1812 }
1813 
1814 PetscErrorCode PetscDSSetContext(PetscDS prob, PetscInt f, void *ctx)
1815 {
1816   PetscErrorCode ierr;
1817 
1818   PetscFunctionBegin;
1819   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
1820   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1821   ierr = PetscDSEnlarge_Static(prob, f+1);CHKERRQ(ierr);
1822   prob->ctx[f] = ctx;
1823   PetscFunctionReturn(0);
1824 }
1825 
1826 /*@C
1827   PetscDSGetBdResidual - Get the pointwise boundary residual function for a given test field
1828 
1829   Not collective
1830 
1831   Input Parameters:
1832 + prob - The PetscDS
1833 - f    - The test field number
1834 
1835   Output Parameters:
1836 + f0 - boundary integrand for the test function term
1837 - f1 - boundary integrand for the test function gradient term
1838 
1839   Note: We are using a first order FEM model for the weak form:
1840 
1841   \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
1842 
1843 The calling sequence for the callbacks f0 and f1 is given by:
1844 
1845 $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1846 $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1847 $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1848 $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[])
1849 
1850 + dim - the spatial dimension
1851 . Nf - the number of fields
1852 . uOff - the offset into u[] and u_t[] for each field
1853 . uOff_x - the offset into u_x[] for each field
1854 . u - each field evaluated at the current point
1855 . u_t - the time derivative of each field evaluated at the current point
1856 . u_x - the gradient of each field evaluated at the current point
1857 . aOff - the offset into a[] and a_t[] for each auxiliary field
1858 . aOff_x - the offset into a_x[] for each auxiliary field
1859 . a - each auxiliary field evaluated at the current point
1860 . a_t - the time derivative of each auxiliary field evaluated at the current point
1861 . a_x - the gradient of auxiliary each field evaluated at the current point
1862 . t - current time
1863 . x - coordinates of the current point
1864 . n - unit normal at the current point
1865 . numConstants - number of constant parameters
1866 . constants - constant parameters
1867 - f0 - output values at the current point
1868 
1869   Level: intermediate
1870 
1871 .seealso: PetscDSSetBdResidual()
1872 @*/
1873 PetscErrorCode PetscDSGetBdResidual(PetscDS prob, PetscInt f,
1874                                     void (**f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1875                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1876                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1877                                                 PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]),
1878                                     void (**f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1879                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1880                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1881                                                 PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
1882 {
1883   PetscFunctionBegin;
1884   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
1885   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
1886   if (f0) {PetscValidPointer(f0, 3); *f0 = prob->fBd[f*2+0];}
1887   if (f1) {PetscValidPointer(f1, 4); *f1 = prob->fBd[f*2+1];}
1888   PetscFunctionReturn(0);
1889 }
1890 
1891 /*@C
1892   PetscDSSetBdResidual - Get the pointwise boundary residual function for a given test field
1893 
1894   Not collective
1895 
1896   Input Parameters:
1897 + prob - The PetscDS
1898 . f    - The test field number
1899 . f0 - boundary integrand for the test function term
1900 - f1 - boundary integrand for the test function gradient term
1901 
1902   Note: We are using a first order FEM model for the weak form:
1903 
1904   \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
1905 
1906 The calling sequence for the callbacks f0 and f1 is given by:
1907 
1908 $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1909 $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1910 $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1911 $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[])
1912 
1913 + dim - the spatial dimension
1914 . Nf - the number of fields
1915 . uOff - the offset into u[] and u_t[] for each field
1916 . uOff_x - the offset into u_x[] for each field
1917 . u - each field evaluated at the current point
1918 . u_t - the time derivative of each field evaluated at the current point
1919 . u_x - the gradient of each field evaluated at the current point
1920 . aOff - the offset into a[] and a_t[] for each auxiliary field
1921 . aOff_x - the offset into a_x[] for each auxiliary field
1922 . a - each auxiliary field evaluated at the current point
1923 . a_t - the time derivative of each auxiliary field evaluated at the current point
1924 . a_x - the gradient of auxiliary each field evaluated at the current point
1925 . t - current time
1926 . x - coordinates of the current point
1927 . n - unit normal at the current point
1928 . numConstants - number of constant parameters
1929 . constants - constant parameters
1930 - f0 - output values at the current point
1931 
1932   Level: intermediate
1933 
1934 .seealso: PetscDSGetBdResidual()
1935 @*/
1936 PetscErrorCode PetscDSSetBdResidual(PetscDS prob, PetscInt f,
1937                                     void (*f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1938                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1939                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1940                                                PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]),
1941                                     void (*f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1942                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1943                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1944                                                PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
1945 {
1946   PetscErrorCode ierr;
1947 
1948   PetscFunctionBegin;
1949   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
1950   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1951   ierr = PetscDSEnlarge_Static(prob, f+1);CHKERRQ(ierr);
1952   if (f0) {PetscValidFunction(f0, 3); prob->fBd[f*2+0] = f0;}
1953   if (f1) {PetscValidFunction(f1, 4); prob->fBd[f*2+1] = f1;}
1954   PetscFunctionReturn(0);
1955 }
1956 
1957 /*@C
1958   PetscDSGetBdJacobian - Get the pointwise boundary Jacobian function for given test and basis field
1959 
1960   Not collective
1961 
1962   Input Parameters:
1963 + prob - The PetscDS
1964 . f    - The test field number
1965 - g    - The field number
1966 
1967   Output Parameters:
1968 + g0 - integrand for the test and basis function term
1969 . g1 - integrand for the test function and basis function gradient term
1970 . g2 - integrand for the test function gradient and basis function term
1971 - g3 - integrand for the test function gradient and basis function gradient term
1972 
1973   Note: We are using a first order FEM model for the weak form:
1974 
1975   \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
1976 
1977 The calling sequence for the callbacks g0, g1, g2 and g3 is given by:
1978 
1979 $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1980 $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1981 $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1982 $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar g0[])
1983 
1984 + dim - the spatial dimension
1985 . Nf - the number of fields
1986 . uOff - the offset into u[] and u_t[] for each field
1987 . uOff_x - the offset into u_x[] for each field
1988 . u - each field evaluated at the current point
1989 . u_t - the time derivative of each field evaluated at the current point
1990 . u_x - the gradient of each field evaluated at the current point
1991 . aOff - the offset into a[] and a_t[] for each auxiliary field
1992 . aOff_x - the offset into a_x[] for each auxiliary field
1993 . a - each auxiliary field evaluated at the current point
1994 . a_t - the time derivative of each auxiliary field evaluated at the current point
1995 . a_x - the gradient of auxiliary each field evaluated at the current point
1996 . t - current time
1997 . u_tShift - the multiplier a for dF/dU_t
1998 . x - coordinates of the current point
1999 . n - normal at the current point
2000 . numConstants - number of constant parameters
2001 . constants - constant parameters
2002 - g0 - output values at the current point
2003 
2004   Level: intermediate
2005 
2006 .seealso: PetscDSSetBdJacobian()
2007 @*/
2008 PetscErrorCode PetscDSGetBdJacobian(PetscDS prob, PetscInt f, PetscInt g,
2009                                     void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2010                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2011                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2012                                                 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
2013                                     void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2014                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2015                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2016                                                 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
2017                                     void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2018                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2019                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2020                                                 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
2021                                     void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2022                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2023                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2024                                                 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
2025 {
2026   PetscFunctionBegin;
2027   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
2028   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
2029   if ((g < 0) || (g >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, prob->Nf);
2030   if (g0) {PetscValidPointer(g0, 4); *g0 = prob->gBd[(f*prob->Nf + g)*4+0];}
2031   if (g1) {PetscValidPointer(g1, 5); *g1 = prob->gBd[(f*prob->Nf + g)*4+1];}
2032   if (g2) {PetscValidPointer(g2, 6); *g2 = prob->gBd[(f*prob->Nf + g)*4+2];}
2033   if (g3) {PetscValidPointer(g3, 7); *g3 = prob->gBd[(f*prob->Nf + g)*4+3];}
2034   PetscFunctionReturn(0);
2035 }
2036 
2037 /*@C
2038   PetscDSSetBdJacobian - Set the pointwise boundary Jacobian function for given test and basis field
2039 
2040   Not collective
2041 
2042   Input Parameters:
2043 + prob - The PetscDS
2044 . f    - The test field number
2045 . g    - The field number
2046 . g0 - integrand for the test and basis function term
2047 . g1 - integrand for the test function and basis function gradient term
2048 . g2 - integrand for the test function gradient and basis function term
2049 - g3 - integrand for the test function gradient and basis function gradient term
2050 
2051   Note: We are using a first order FEM model for the weak form:
2052 
2053   \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
2054 
2055 The calling sequence for the callbacks g0, g1, g2 and g3 is given by:
2056 
2057 $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2058 $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2059 $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2060 $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar g0[])
2061 
2062 + dim - the spatial dimension
2063 . Nf - the number of fields
2064 . uOff - the offset into u[] and u_t[] for each field
2065 . uOff_x - the offset into u_x[] for each field
2066 . u - each field evaluated at the current point
2067 . u_t - the time derivative of each field evaluated at the current point
2068 . u_x - the gradient of each field evaluated at the current point
2069 . aOff - the offset into a[] and a_t[] for each auxiliary field
2070 . aOff_x - the offset into a_x[] for each auxiliary field
2071 . a - each auxiliary field evaluated at the current point
2072 . a_t - the time derivative of each auxiliary field evaluated at the current point
2073 . a_x - the gradient of auxiliary each field evaluated at the current point
2074 . t - current time
2075 . u_tShift - the multiplier a for dF/dU_t
2076 . x - coordinates of the current point
2077 . n - normal at the current point
2078 . numConstants - number of constant parameters
2079 . constants - constant parameters
2080 - g0 - output values at the current point
2081 
2082   Level: intermediate
2083 
2084 .seealso: PetscDSGetBdJacobian()
2085 @*/
2086 PetscErrorCode PetscDSSetBdJacobian(PetscDS prob, PetscInt f, PetscInt g,
2087                                     void (*g0)(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, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
2091                                     void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2092                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2093                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2094                                                PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
2095                                     void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2096                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2097                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2098                                                PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
2099                                     void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2100                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2101                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2102                                                PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
2103 {
2104   PetscErrorCode ierr;
2105 
2106   PetscFunctionBegin;
2107   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
2108   if (g0) PetscValidFunction(g0, 4);
2109   if (g1) PetscValidFunction(g1, 5);
2110   if (g2) PetscValidFunction(g2, 6);
2111   if (g3) PetscValidFunction(g3, 7);
2112   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2113   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
2114   ierr = PetscDSEnlarge_Static(prob, PetscMax(f, g)+1);CHKERRQ(ierr);
2115   prob->gBd[(f*prob->Nf + g)*4+0] = g0;
2116   prob->gBd[(f*prob->Nf + g)*4+1] = g1;
2117   prob->gBd[(f*prob->Nf + g)*4+2] = g2;
2118   prob->gBd[(f*prob->Nf + g)*4+3] = g3;
2119   PetscFunctionReturn(0);
2120 }
2121 
2122 /*@C
2123   PetscDSGetExactSolution - Get the pointwise exact solution function for a given test field
2124 
2125   Not collective
2126 
2127   Input Parameters:
2128 + prob - The PetscDS
2129 - f    - The test field number
2130 
2131   Output Parameter:
2132 . exactSol - exact solution for the test field
2133 
2134   Note: The calling sequence for the solution functions is given by:
2135 
2136 $ sol(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx)
2137 
2138 + dim - the spatial dimension
2139 . t - current time
2140 . x - coordinates of the current point
2141 . Nc - the number of field components
2142 . u - the solution field evaluated at the current point
2143 - ctx - a user context
2144 
2145   Level: intermediate
2146 
2147 .seealso: PetscDSSetExactSolution()
2148 @*/
2149 PetscErrorCode PetscDSGetExactSolution(PetscDS prob, PetscInt f, PetscErrorCode (**sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx))
2150 {
2151   PetscFunctionBegin;
2152   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
2153   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
2154   if (sol) {PetscValidPointer(sol, 3); *sol = prob->exactSol[f];}
2155   PetscFunctionReturn(0);
2156 }
2157 
2158 /*@C
2159   PetscDSSetExactSolution - Get the pointwise exact solution function for a given test field
2160 
2161   Not collective
2162 
2163   Input Parameters:
2164 + prob - The PetscDS
2165 . f    - The test field number
2166 - sol  - solution function for the test fields
2167 
2168   Note: The calling sequence for solution functions is given by:
2169 
2170 $ sol(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx)
2171 
2172 + dim - the spatial dimension
2173 . t - current time
2174 . x - coordinates of the current point
2175 . Nc - the number of field components
2176 . u - the solution field evaluated at the current point
2177 - ctx - a user context
2178 
2179   Level: intermediate
2180 
2181 .seealso: PetscDSGetExactSolution()
2182 @*/
2183 PetscErrorCode PetscDSSetExactSolution(PetscDS prob, PetscInt f, PetscErrorCode (*sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx))
2184 {
2185   PetscErrorCode ierr;
2186 
2187   PetscFunctionBegin;
2188   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
2189   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2190   ierr = PetscDSEnlarge_Static(prob, f+1);CHKERRQ(ierr);
2191   if (sol) {PetscValidFunction(sol, 3); prob->exactSol[f] = sol;}
2192   PetscFunctionReturn(0);
2193 }
2194 
2195 /*@C
2196   PetscDSGetConstants - Returns the array of constants passed to point functions
2197 
2198   Not collective
2199 
2200   Input Parameter:
2201 . prob - The PetscDS object
2202 
2203   Output Parameters:
2204 + numConstants - The number of constants
2205 - constants    - The array of constants, NULL if there are none
2206 
2207   Level: intermediate
2208 
2209 .seealso: PetscDSSetConstants(), PetscDSCreate()
2210 @*/
2211 PetscErrorCode PetscDSGetConstants(PetscDS prob, PetscInt *numConstants, const PetscScalar *constants[])
2212 {
2213   PetscFunctionBegin;
2214   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
2215   if (numConstants) {PetscValidPointer(numConstants, 2); *numConstants = prob->numConstants;}
2216   if (constants)    {PetscValidPointer(constants, 3);    *constants    = prob->constants;}
2217   PetscFunctionReturn(0);
2218 }
2219 
2220 /*@C
2221   PetscDSSetConstants - Set the array of constants passed to point functions
2222 
2223   Not collective
2224 
2225   Input Parameters:
2226 + prob         - The PetscDS object
2227 . numConstants - The number of constants
2228 - constants    - The array of constants, NULL if there are none
2229 
2230   Level: intermediate
2231 
2232 .seealso: PetscDSGetConstants(), PetscDSCreate()
2233 @*/
2234 PetscErrorCode PetscDSSetConstants(PetscDS prob, PetscInt numConstants, PetscScalar constants[])
2235 {
2236   PetscErrorCode ierr;
2237 
2238   PetscFunctionBegin;
2239   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
2240   if (numConstants != prob->numConstants) {
2241     ierr = PetscFree(prob->constants);CHKERRQ(ierr);
2242     prob->numConstants = numConstants;
2243     if (prob->numConstants) {
2244       ierr = PetscMalloc1(prob->numConstants, &prob->constants);CHKERRQ(ierr);
2245     } else {
2246       prob->constants = NULL;
2247     }
2248   }
2249   if (prob->numConstants) {
2250     PetscValidPointer(constants, 3);
2251     ierr = PetscMemcpy(prob->constants, constants, prob->numConstants * sizeof(PetscScalar));CHKERRQ(ierr);
2252   }
2253   PetscFunctionReturn(0);
2254 }
2255 
2256 /*@
2257   PetscDSGetFieldIndex - Returns the index of the given field
2258 
2259   Not collective
2260 
2261   Input Parameters:
2262 + prob - The PetscDS object
2263 - disc - The discretization object
2264 
2265   Output Parameter:
2266 . f - The field number
2267 
2268   Level: beginner
2269 
2270 .seealso: PetscGetDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
2271 @*/
2272 PetscErrorCode PetscDSGetFieldIndex(PetscDS prob, PetscObject disc, PetscInt *f)
2273 {
2274   PetscInt g;
2275 
2276   PetscFunctionBegin;
2277   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
2278   PetscValidPointer(f, 3);
2279   *f = -1;
2280   for (g = 0; g < prob->Nf; ++g) {if (disc == prob->disc[g]) break;}
2281   if (g == prob->Nf) SETERRQ(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Field not found in PetscDS.");
2282   *f = g;
2283   PetscFunctionReturn(0);
2284 }
2285 
2286 /*@
2287   PetscDSGetFieldSize - Returns the size of the given field in the full space basis
2288 
2289   Not collective
2290 
2291   Input Parameters:
2292 + prob - The PetscDS object
2293 - f - The field number
2294 
2295   Output Parameter:
2296 . size - The size
2297 
2298   Level: beginner
2299 
2300 .seealso: PetscDSGetFieldOffset(), PetscDSGetNumFields(), PetscDSCreate()
2301 @*/
2302 PetscErrorCode PetscDSGetFieldSize(PetscDS prob, PetscInt f, PetscInt *size)
2303 {
2304   PetscErrorCode ierr;
2305 
2306   PetscFunctionBegin;
2307   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
2308   PetscValidPointer(size, 3);
2309   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
2310   ierr = PetscDSSetUp(prob);CHKERRQ(ierr);
2311   *size = prob->Nb[f];
2312   PetscFunctionReturn(0);
2313 }
2314 
2315 /*@
2316   PetscDSGetFieldOffset - Returns the offset of the given field in the full space basis
2317 
2318   Not collective
2319 
2320   Input Parameters:
2321 + prob - The PetscDS object
2322 - f - The field number
2323 
2324   Output Parameter:
2325 . off - The offset
2326 
2327   Level: beginner
2328 
2329 .seealso: PetscDSGetFieldSize(), PetscDSGetNumFields(), PetscDSCreate()
2330 @*/
2331 PetscErrorCode PetscDSGetFieldOffset(PetscDS prob, PetscInt f, PetscInt *off)
2332 {
2333   PetscInt       size, g;
2334   PetscErrorCode ierr;
2335 
2336   PetscFunctionBegin;
2337   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
2338   PetscValidPointer(off, 3);
2339   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
2340   *off = 0;
2341   for (g = 0; g < f; ++g) {
2342     ierr = PetscDSGetFieldSize(prob, g, &size);CHKERRQ(ierr);
2343     *off += size;
2344   }
2345   PetscFunctionReturn(0);
2346 }
2347 
2348 /*@
2349   PetscDSGetDimensions - Returns the size of the approximation space for each field on an evaluation point
2350 
2351   Not collective
2352 
2353   Input Parameter:
2354 . prob - The PetscDS object
2355 
2356   Output Parameter:
2357 . dimensions - The number of dimensions
2358 
2359   Level: beginner
2360 
2361 .seealso: PetscDSGetComponentOffsets(), PetscDSGetNumFields(), PetscDSCreate()
2362 @*/
2363 PetscErrorCode PetscDSGetDimensions(PetscDS prob, PetscInt *dimensions[])
2364 {
2365   PetscErrorCode ierr;
2366 
2367   PetscFunctionBegin;
2368   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
2369   ierr = PetscDSSetUp(prob);CHKERRQ(ierr);
2370   PetscValidPointer(dimensions, 2);
2371   *dimensions = prob->Nb;
2372   PetscFunctionReturn(0);
2373 }
2374 
2375 /*@
2376   PetscDSGetComponents - Returns the number of components for each field on an evaluation point
2377 
2378   Not collective
2379 
2380   Input Parameter:
2381 . prob - The PetscDS object
2382 
2383   Output Parameter:
2384 . components - The number of components
2385 
2386   Level: beginner
2387 
2388 .seealso: PetscDSGetComponentOffsets(), PetscDSGetNumFields(), PetscDSCreate()
2389 @*/
2390 PetscErrorCode PetscDSGetComponents(PetscDS prob, PetscInt *components[])
2391 {
2392   PetscErrorCode ierr;
2393 
2394   PetscFunctionBegin;
2395   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
2396   ierr = PetscDSSetUp(prob);CHKERRQ(ierr);
2397   PetscValidPointer(components, 2);
2398   *components = prob->Nc;
2399   PetscFunctionReturn(0);
2400 }
2401 
2402 /*@
2403   PetscDSGetComponentOffset - Returns the offset of the given field on an evaluation point
2404 
2405   Not collective
2406 
2407   Input Parameters:
2408 + prob - The PetscDS object
2409 - f - The field number
2410 
2411   Output Parameter:
2412 . off - The offset
2413 
2414   Level: beginner
2415 
2416 .seealso: PetscDSGetNumFields(), PetscDSCreate()
2417 @*/
2418 PetscErrorCode PetscDSGetComponentOffset(PetscDS prob, PetscInt f, PetscInt *off)
2419 {
2420   PetscFunctionBegin;
2421   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
2422   PetscValidPointer(off, 3);
2423   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
2424   *off = prob->off[f];
2425   PetscFunctionReturn(0);
2426 }
2427 
2428 /*@
2429   PetscDSGetComponentOffsets - Returns the offset of each field on an evaluation point
2430 
2431   Not collective
2432 
2433   Input Parameter:
2434 . prob - The PetscDS object
2435 
2436   Output Parameter:
2437 . offsets - The offsets
2438 
2439   Level: beginner
2440 
2441 .seealso: PetscDSGetNumFields(), PetscDSCreate()
2442 @*/
2443 PetscErrorCode PetscDSGetComponentOffsets(PetscDS prob, PetscInt *offsets[])
2444 {
2445   PetscFunctionBegin;
2446   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
2447   PetscValidPointer(offsets, 2);
2448   *offsets = prob->off;
2449   PetscFunctionReturn(0);
2450 }
2451 
2452 /*@
2453   PetscDSGetComponentDerivativeOffsets - Returns the offset of each field derivative on an evaluation point
2454 
2455   Not collective
2456 
2457   Input Parameter:
2458 . prob - The PetscDS object
2459 
2460   Output Parameter:
2461 . offsets - The offsets
2462 
2463   Level: beginner
2464 
2465 .seealso: PetscDSGetNumFields(), PetscDSCreate()
2466 @*/
2467 PetscErrorCode PetscDSGetComponentDerivativeOffsets(PetscDS prob, PetscInt *offsets[])
2468 {
2469   PetscFunctionBegin;
2470   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
2471   PetscValidPointer(offsets, 2);
2472   *offsets = prob->offDer;
2473   PetscFunctionReturn(0);
2474 }
2475 
2476 /*@C
2477   PetscDSGetTabulation - Return the basis tabulation at quadrature points for the volume discretization
2478 
2479   Not collective
2480 
2481   Input Parameter:
2482 . prob - The PetscDS object
2483 
2484   Output Parameters:
2485 + basis - The basis function tabulation at quadrature points
2486 - basisDer - The basis function derivative tabulation at quadrature points
2487 
2488   Level: intermediate
2489 
2490 .seealso: PetscDSCreate()
2491 @*/
2492 PetscErrorCode PetscDSGetTabulation(PetscDS prob, PetscReal ***basis, PetscReal ***basisDer)
2493 {
2494   PetscErrorCode ierr;
2495 
2496   PetscFunctionBegin;
2497   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
2498   ierr = PetscDSSetUp(prob);CHKERRQ(ierr);
2499   if (basis)    {PetscValidPointer(basis, 2);    *basis    = prob->basis;}
2500   if (basisDer) {PetscValidPointer(basisDer, 3); *basisDer = prob->basisDer;}
2501   PetscFunctionReturn(0);
2502 }
2503 
2504 /*@C
2505   PetscDSGetFaceTabulation - Return the basis tabulation at quadrature points on the faces
2506 
2507   Not collective
2508 
2509   Input Parameter:
2510 . prob - The PetscDS object
2511 
2512   Output Parameters:
2513 + basisFace - The basis function tabulation at quadrature points
2514 - basisDerFace - The basis function derivative tabulation at quadrature points
2515 
2516   Level: intermediate
2517 
2518 .seealso: PetscDSGetTabulation(), PetscDSCreate()
2519 @*/
2520 PetscErrorCode PetscDSGetFaceTabulation(PetscDS prob, PetscReal ***basis, PetscReal ***basisDer)
2521 {
2522   PetscErrorCode ierr;
2523 
2524   PetscFunctionBegin;
2525   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
2526   ierr = PetscDSSetUp(prob);CHKERRQ(ierr);
2527   if (basis)    {PetscValidPointer(basis, 2);    *basis    = prob->basisFace;}
2528   if (basisDer) {PetscValidPointer(basisDer, 3); *basisDer = prob->basisDerFace;}
2529   PetscFunctionReturn(0);
2530 }
2531 
2532 PetscErrorCode PetscDSGetEvaluationArrays(PetscDS prob, PetscScalar **u, PetscScalar **u_t, PetscScalar **u_x)
2533 {
2534   PetscErrorCode ierr;
2535 
2536   PetscFunctionBegin;
2537   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
2538   ierr = PetscDSSetUp(prob);CHKERRQ(ierr);
2539   if (u)   {PetscValidPointer(u, 2);   *u   = prob->u;}
2540   if (u_t) {PetscValidPointer(u_t, 3); *u_t = prob->u_t;}
2541   if (u_x) {PetscValidPointer(u_x, 4); *u_x = prob->u_x;}
2542   PetscFunctionReturn(0);
2543 }
2544 
2545 PetscErrorCode PetscDSGetWeakFormArrays(PetscDS prob, PetscScalar **f0, PetscScalar **f1, PetscScalar **g0, PetscScalar **g1, PetscScalar **g2, PetscScalar **g3)
2546 {
2547   PetscErrorCode ierr;
2548 
2549   PetscFunctionBegin;
2550   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
2551   ierr = PetscDSSetUp(prob);CHKERRQ(ierr);
2552   if (f0) {PetscValidPointer(f0, 2); *f0 = prob->f0;}
2553   if (f1) {PetscValidPointer(f1, 3); *f1 = prob->f1;}
2554   if (g0) {PetscValidPointer(g0, 4); *g0 = prob->g0;}
2555   if (g1) {PetscValidPointer(g1, 5); *g1 = prob->g1;}
2556   if (g2) {PetscValidPointer(g2, 6); *g2 = prob->g2;}
2557   if (g3) {PetscValidPointer(g3, 7); *g3 = prob->g3;}
2558   PetscFunctionReturn(0);
2559 }
2560 
2561 PetscErrorCode PetscDSGetRefCoordArrays(PetscDS prob, PetscReal **x, PetscScalar **refSpaceDer)
2562 {
2563   PetscErrorCode ierr;
2564 
2565   PetscFunctionBegin;
2566   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
2567   ierr = PetscDSSetUp(prob);CHKERRQ(ierr);
2568   if (x)           {PetscValidPointer(x, 2);           *x           = prob->x;}
2569   if (refSpaceDer) {PetscValidPointer(refSpaceDer, 3); *refSpaceDer = prob->refSpaceDer;}
2570   PetscFunctionReturn(0);
2571 }
2572 
2573 /*@C
2574   PetscDSAddBoundary - Add a boundary condition to the model
2575 
2576   Input Parameters:
2577 + ds          - The PetscDS object
2578 . type        - The type of condition, e.g. DM_BC_ESSENTIAL/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
2579 . name        - The BC name
2580 . labelname   - The label defining constrained points
2581 . field       - The field to constrain
2582 . numcomps    - The number of constrained field components
2583 . comps       - An array of constrained component numbers
2584 . bcFunc      - A pointwise function giving boundary values
2585 . numids      - The number of DMLabel ids for constrained points
2586 . ids         - An array of ids for constrained points
2587 - ctx         - An optional user context for bcFunc
2588 
2589   Options Database Keys:
2590 + -bc_<boundary name> <num> - Overrides the boundary ids
2591 - -bc_<boundary name>_comp <num> - Overrides the boundary components
2592 
2593   Level: developer
2594 
2595 .seealso: PetscDSGetBoundary()
2596 @*/
2597 PetscErrorCode PetscDSAddBoundary(PetscDS ds, DMBoundaryConditionType type, const char name[], const char labelname[], PetscInt field, PetscInt numcomps, const PetscInt *comps, void (*bcFunc)(void), PetscInt numids, const PetscInt *ids, void *ctx)
2598 {
2599   DSBoundary     b;
2600   PetscErrorCode ierr;
2601 
2602   PetscFunctionBegin;
2603   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
2604   ierr = PetscNew(&b);CHKERRQ(ierr);
2605   ierr = PetscStrallocpy(name, (char **) &b->name);CHKERRQ(ierr);
2606   ierr = PetscStrallocpy(labelname, (char **) &b->labelname);CHKERRQ(ierr);
2607   ierr = PetscMalloc1(numcomps, &b->comps);CHKERRQ(ierr);
2608   if (numcomps) {ierr = PetscMemcpy(b->comps, comps, numcomps*sizeof(PetscInt));CHKERRQ(ierr);}
2609   ierr = PetscMalloc1(numids, &b->ids);CHKERRQ(ierr);
2610   if (numids) {ierr = PetscMemcpy(b->ids, ids, numids*sizeof(PetscInt));CHKERRQ(ierr);}
2611   b->type            = type;
2612   b->field           = field;
2613   b->numcomps        = numcomps;
2614   b->func            = bcFunc;
2615   b->numids          = numids;
2616   b->ctx             = ctx;
2617   b->next            = ds->boundary;
2618   ds->boundary       = b;
2619   PetscFunctionReturn(0);
2620 }
2621 
2622 /*@
2623   PetscDSGetNumBoundary - Get the number of registered BC
2624 
2625   Input Parameters:
2626 . ds - The PetscDS object
2627 
2628   Output Parameters:
2629 . numBd - The number of BC
2630 
2631   Level: intermediate
2632 
2633 .seealso: PetscDSAddBoundary(), PetscDSGetBoundary()
2634 @*/
2635 PetscErrorCode PetscDSGetNumBoundary(PetscDS ds, PetscInt *numBd)
2636 {
2637   DSBoundary b = ds->boundary;
2638 
2639   PetscFunctionBegin;
2640   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
2641   PetscValidPointer(numBd, 2);
2642   *numBd = 0;
2643   while (b) {++(*numBd); b = b->next;}
2644   PetscFunctionReturn(0);
2645 }
2646 
2647 /*@C
2648   PetscDSGetBoundary - Add a boundary condition to the model
2649 
2650   Input Parameters:
2651 + ds          - The PetscDS object
2652 - bd          - The BC number
2653 
2654   Output Parameters:
2655 + type        - The type of condition, e.g. DM_BC_ESSENTIAL/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
2656 . name        - The BC name
2657 . labelname   - The label defining constrained points
2658 . field       - The field to constrain
2659 . numcomps    - The number of constrained field components
2660 . comps       - An array of constrained component numbers
2661 . bcFunc      - A pointwise function giving boundary values
2662 . numids      - The number of DMLabel ids for constrained points
2663 . ids         - An array of ids for constrained points
2664 - ctx         - An optional user context for bcFunc
2665 
2666   Options Database Keys:
2667 + -bc_<boundary name> <num> - Overrides the boundary ids
2668 - -bc_<boundary name>_comp <num> - Overrides the boundary components
2669 
2670   Level: developer
2671 
2672 .seealso: PetscDSAddBoundary()
2673 @*/
2674 PetscErrorCode PetscDSGetBoundary(PetscDS ds, PetscInt bd, DMBoundaryConditionType *type, const char **name, const char **labelname, PetscInt *field, PetscInt *numcomps, const PetscInt **comps, void (**func)(void), PetscInt *numids, const PetscInt **ids, void **ctx)
2675 {
2676   DSBoundary b    = ds->boundary;
2677   PetscInt   n    = 0;
2678 
2679   PetscFunctionBegin;
2680   PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
2681   while (b) {
2682     if (n == bd) break;
2683     b = b->next;
2684     ++n;
2685   }
2686   if (!b) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %d is not in [0, %d)", bd, n);
2687   if (type) {
2688     PetscValidPointer(type, 3);
2689     *type = b->type;
2690   }
2691   if (name) {
2692     PetscValidPointer(name, 4);
2693     *name = b->name;
2694   }
2695   if (labelname) {
2696     PetscValidPointer(labelname, 5);
2697     *labelname = b->labelname;
2698   }
2699   if (field) {
2700     PetscValidPointer(field, 6);
2701     *field = b->field;
2702   }
2703   if (numcomps) {
2704     PetscValidPointer(numcomps, 7);
2705     *numcomps = b->numcomps;
2706   }
2707   if (comps) {
2708     PetscValidPointer(comps, 8);
2709     *comps = b->comps;
2710   }
2711   if (func) {
2712     PetscValidPointer(func, 9);
2713     *func = b->func;
2714   }
2715   if (numids) {
2716     PetscValidPointer(numids, 10);
2717     *numids = b->numids;
2718   }
2719   if (ids) {
2720     PetscValidPointer(ids, 11);
2721     *ids = b->ids;
2722   }
2723   if (ctx) {
2724     PetscValidPointer(ctx, 12);
2725     *ctx = b->ctx;
2726   }
2727   PetscFunctionReturn(0);
2728 }
2729 
2730 PetscErrorCode PetscDSCopyBoundary(PetscDS probA, PetscDS probB)
2731 {
2732   DSBoundary     b, next, *lastnext;
2733   PetscErrorCode ierr;
2734 
2735   PetscFunctionBegin;
2736   PetscValidHeaderSpecific(probA, PETSCDS_CLASSID, 1);
2737   PetscValidHeaderSpecific(probB, PETSCDS_CLASSID, 2);
2738   if (probA == probB) PetscFunctionReturn(0);
2739   next = probB->boundary;
2740   while (next) {
2741     DSBoundary b = next;
2742 
2743     next = b->next;
2744     ierr = PetscFree(b->comps);CHKERRQ(ierr);
2745     ierr = PetscFree(b->ids);CHKERRQ(ierr);
2746     ierr = PetscFree(b->name);CHKERRQ(ierr);
2747     ierr = PetscFree(b->labelname);CHKERRQ(ierr);
2748     ierr = PetscFree(b);CHKERRQ(ierr);
2749   }
2750   lastnext = &(probB->boundary);
2751   for (b = probA->boundary; b; b = b->next) {
2752     DSBoundary bNew;
2753 
2754     ierr = PetscNew(&bNew);CHKERRQ(ierr);
2755     bNew->numcomps = b->numcomps;
2756     ierr = PetscMalloc1(bNew->numcomps, &bNew->comps);CHKERRQ(ierr);
2757     ierr = PetscMemcpy(bNew->comps, b->comps, bNew->numcomps*sizeof(PetscInt));CHKERRQ(ierr);
2758     bNew->numids = b->numids;
2759     ierr = PetscMalloc1(bNew->numids, &bNew->ids);CHKERRQ(ierr);
2760     ierr = PetscMemcpy(bNew->ids, b->ids, bNew->numids*sizeof(PetscInt));CHKERRQ(ierr);
2761     ierr = PetscStrallocpy(b->labelname,(char **) &(bNew->labelname));CHKERRQ(ierr);
2762     ierr = PetscStrallocpy(b->name,(char **) &(bNew->name));CHKERRQ(ierr);
2763     bNew->ctx   = b->ctx;
2764     bNew->type  = b->type;
2765     bNew->field = b->field;
2766     bNew->func  = b->func;
2767 
2768     *lastnext = bNew;
2769     lastnext = &(bNew->next);
2770   }
2771   PetscFunctionReturn(0);
2772 }
2773 
2774 /*@
2775   PetscDSCopyEquations - Copy all pointwise function pointers to the new problem
2776 
2777   Not collective
2778 
2779   Input Parameter:
2780 . prob - The PetscDS object
2781 
2782   Output Parameter:
2783 . newprob - The PetscDS copy
2784 
2785   Level: intermediate
2786 
2787 .seealso: PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
2788 @*/
2789 PetscErrorCode PetscDSCopyEquations(PetscDS prob, PetscDS newprob)
2790 {
2791   PetscInt       Nf, Ng, f, g;
2792   PetscErrorCode ierr;
2793 
2794   PetscFunctionBegin;
2795   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
2796   PetscValidHeaderSpecific(newprob, PETSCDS_CLASSID, 2);
2797   ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr);
2798   ierr = PetscDSGetNumFields(newprob, &Ng);CHKERRQ(ierr);
2799   if (Nf != Ng) SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_SIZ, "Number of fields must match %D != %D", Nf, Ng);
2800   for (f = 0; f < Nf; ++f) {
2801     PetscPointFunc   obj;
2802     PetscPointFunc   f0, f1;
2803     PetscPointJac    g0, g1, g2, g3;
2804     PetscBdPointFunc f0Bd, f1Bd;
2805     PetscBdPointJac  g0Bd, g1Bd, g2Bd, g3Bd;
2806     PetscRiemannFunc r;
2807 
2808     ierr = PetscDSGetObjective(prob, f, &obj);CHKERRQ(ierr);
2809     ierr = PetscDSGetResidual(prob, f, &f0, &f1);CHKERRQ(ierr);
2810     ierr = PetscDSGetBdResidual(prob, f, &f0Bd, &f1Bd);CHKERRQ(ierr);
2811     ierr = PetscDSGetRiemannSolver(prob, f, &r);CHKERRQ(ierr);
2812     ierr = PetscDSSetObjective(newprob, f, obj);CHKERRQ(ierr);
2813     ierr = PetscDSSetResidual(newprob, f, f0, f1);CHKERRQ(ierr);
2814     ierr = PetscDSSetBdResidual(newprob, f, f0Bd, f1Bd);CHKERRQ(ierr);
2815     ierr = PetscDSSetRiemannSolver(newprob, f, r);CHKERRQ(ierr);
2816     for (g = 0; g < Nf; ++g) {
2817       ierr = PetscDSGetJacobian(prob, f, g, &g0, &g1, &g2, &g3);CHKERRQ(ierr);
2818       ierr = PetscDSGetBdJacobian(prob, f, g, &g0Bd, &g1Bd, &g2Bd, &g3Bd);CHKERRQ(ierr);
2819       ierr = PetscDSSetJacobian(newprob, f, g, g0, g1, g2, g3);CHKERRQ(ierr);
2820       ierr = PetscDSSetBdJacobian(newprob, f, g, g0Bd, g1Bd, g2Bd, g3Bd);CHKERRQ(ierr);
2821     }
2822   }
2823   PetscFunctionReturn(0);
2824 }
2825 
2826 static PetscErrorCode PetscDSDestroy_Basic(PetscDS prob)
2827 {
2828   PetscErrorCode      ierr;
2829 
2830   PetscFunctionBegin;
2831   ierr = PetscFree(prob->data);CHKERRQ(ierr);
2832   PetscFunctionReturn(0);
2833 }
2834 
2835 static PetscErrorCode PetscDSInitialize_Basic(PetscDS prob)
2836 {
2837   PetscFunctionBegin;
2838   prob->ops->setfromoptions = NULL;
2839   prob->ops->setup          = NULL;
2840   prob->ops->view           = NULL;
2841   prob->ops->destroy        = PetscDSDestroy_Basic;
2842   PetscFunctionReturn(0);
2843 }
2844 
2845 /*MC
2846   PETSCDSBASIC = "basic" - A discrete system with pointwise residual and boundary residual functions
2847 
2848   Level: intermediate
2849 
2850 .seealso: PetscDSType, PetscDSCreate(), PetscDSSetType()
2851 M*/
2852 
2853 PETSC_EXTERN PetscErrorCode PetscDSCreate_Basic(PetscDS prob)
2854 {
2855   PetscDS_Basic *b;
2856   PetscErrorCode      ierr;
2857 
2858   PetscFunctionBegin;
2859   PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
2860   ierr       = PetscNewLog(prob, &b);CHKERRQ(ierr);
2861   prob->data = b;
2862 
2863   ierr = PetscDSInitialize_Basic(prob);CHKERRQ(ierr);
2864   PetscFunctionReturn(0);
2865 }
2866