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