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