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