Actual source code: dtds.c

petsc-main 2021-04-20
Report Typos and Errors
  1: #include <petsc/private/petscdsimpl.h>

  3: PetscClassId PETSCDS_CLASSID = 0;

  5: PetscFunctionList PetscDSList              = NULL;
  6: PetscBool         PetscDSRegisterAllCalled = PETSC_FALSE;

  8: /* A PetscDS (Discrete System) encodes a set of equations posed in a discrete space, which represents a set of
  9:    nonlinear continuum equations. The equations can have multiple fields, each field having a different
 10:    discretization. In addition, different pieces of the domain can have different field combinations and equations.

 12:    The DS provides the user a description of the approximation space on any given cell. It also gives pointwise
 13:    functions representing the equations.

 15:    Each field is associated with a label, marking the cells on which it is supported. Note that a field can be
 16:    supported on the closure of a cell not in the label due to overlap of the boundary of neighboring cells. The DM
 17:    then creates a DS for each set of cells with identical approximation spaces. When assembling, the user asks for
 18:    the space associated with a given cell. DMPlex uses the labels associated with each DS in the default integration loop.
 19: */

 21: /*@C
 22:   PetscDSRegister - Adds a new PetscDS implementation

 24:   Not Collective

 26:   Input Parameters:
 27: + name        - The name of a new user-defined creation routine
 28: - create_func - The creation routine itself

 30:   Notes:
 31:   PetscDSRegister() may be called multiple times to add several user-defined PetscDSs

 33:   Sample usage:
 34: .vb
 35:     PetscDSRegister("my_ds", MyPetscDSCreate);
 36: .ve

 38:   Then, your PetscDS type can be chosen with the procedural interface via
 39: .vb
 40:     PetscDSCreate(MPI_Comm, PetscDS *);
 41:     PetscDSSetType(PetscDS, "my_ds");
 42: .ve
 43:    or at runtime via the option
 44: .vb
 45:     -petscds_type my_ds
 46: .ve

 48:   Level: advanced

 50:    Not available from Fortran

 52: .seealso: PetscDSRegisterAll(), PetscDSRegisterDestroy()

 54: @*/
 55: PetscErrorCode PetscDSRegister(const char sname[], PetscErrorCode (*function)(PetscDS))
 56: {

 60:   PetscFunctionListAdd(&PetscDSList, sname, function);
 61:   return(0);
 62: }

 64: /*@C
 65:   PetscDSSetType - Builds a particular PetscDS

 67:   Collective on prob

 69:   Input Parameters:
 70: + prob - The PetscDS object
 71: - name - The kind of system

 73:   Options Database Key:
 74: . -petscds_type <type> - Sets the PetscDS type; use -help for a list of available types

 76:   Level: intermediate

 78:    Not available from Fortran

 80: .seealso: PetscDSGetType(), PetscDSCreate()
 81: @*/
 82: PetscErrorCode PetscDSSetType(PetscDS prob, PetscDSType name)
 83: {
 84:   PetscErrorCode (*r)(PetscDS);
 85:   PetscBool      match;

 90:   PetscObjectTypeCompare((PetscObject) prob, name, &match);
 91:   if (match) return(0);

 93:   PetscDSRegisterAll();
 94:   PetscFunctionListFind(PetscDSList, name, &r);
 95:   if (!r) SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscDS type: %s", name);

 97:   if (prob->ops->destroy) {
 98:     (*prob->ops->destroy)(prob);
 99:     prob->ops->destroy = NULL;
100:   }
101:   (*r)(prob);
102:   PetscObjectChangeTypeName((PetscObject) prob, name);
103:   return(0);
104: }

106: /*@C
107:   PetscDSGetType - Gets the PetscDS type name (as a string) from the object.

109:   Not Collective

111:   Input Parameter:
112: . prob  - The PetscDS

114:   Output Parameter:
115: . name - The PetscDS type name

117:   Level: intermediate

119:    Not available from Fortran

121: .seealso: PetscDSSetType(), PetscDSCreate()
122: @*/
123: PetscErrorCode PetscDSGetType(PetscDS prob, PetscDSType *name)
124: {

130:   PetscDSRegisterAll();
131:   *name = ((PetscObject) prob)->type_name;
132:   return(0);
133: }

135: static PetscErrorCode PetscDSView_Ascii(PetscDS prob, PetscViewer viewer)
136: {
137:   PetscViewerFormat  format;
138:   const PetscScalar *constants;
139:   PetscInt           numConstants, f;
140:   PetscErrorCode     ierr;

143:   PetscViewerGetFormat(viewer, &format);
144:   PetscViewerASCIIPrintf(viewer, "Discrete System with %d fields\n", prob->Nf);
145:   PetscViewerASCIIPushTab(viewer);
146:   PetscViewerASCIIPrintf(viewer, "  cell total dim %D total comp %D\n", prob->totDim, prob->totComp);
147:   if (prob->isHybrid) {PetscViewerASCIIPrintf(viewer, "  hybrid cell\n");}
148:   for (f = 0; f < prob->Nf; ++f) {
149:     DSBoundary      b;
150:     PetscObject     obj;
151:     PetscClassId    id;
152:     PetscQuadrature q;
153:     const char     *name;
154:     PetscInt        Nc, Nq, Nqc;

156:     PetscDSGetDiscretization(prob, f, &obj);
157:     PetscObjectGetClassId(obj, &id);
158:     PetscObjectGetName(obj, &name);
159:     PetscViewerASCIIPrintf(viewer, "Field %s", name ? name : "<unknown>");
160:     PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
161:     if (id == PETSCFE_CLASSID)      {
162:       PetscFEGetNumComponents((PetscFE) obj, &Nc);
163:       PetscFEGetQuadrature((PetscFE) obj, &q);
164:       PetscViewerASCIIPrintf(viewer, " FEM");
165:     } else if (id == PETSCFV_CLASSID) {
166:       PetscFVGetNumComponents((PetscFV) obj, &Nc);
167:       PetscFVGetQuadrature((PetscFV) obj, &q);
168:       PetscViewerASCIIPrintf(viewer, " FVM");
169:     }
170:     else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", f);
171:     if (Nc > 1) {PetscViewerASCIIPrintf(viewer, "%D components", Nc);}
172:     else        {PetscViewerASCIIPrintf(viewer, "%D component ", Nc);}
173:     if (prob->implicit[f]) {PetscViewerASCIIPrintf(viewer, " (implicit)");}
174:     else                   {PetscViewerASCIIPrintf(viewer, " (explicit)");}
175:     if (q) {
176:       PetscQuadratureGetData(q, NULL, &Nqc, &Nq, NULL, NULL);
177:       PetscViewerASCIIPrintf(viewer, " (Nq %D Nqc %D)", Nq, Nqc);
178:     }
179:     PetscViewerASCIIPrintf(viewer, " %D-jet", prob->jetDegree[f]);
180:     PetscViewerASCIIPrintf(viewer, "\n");
181:     PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
182:     PetscViewerASCIIPushTab(viewer);
183:     if (id == PETSCFE_CLASSID)      {PetscFEView((PetscFE) obj, viewer);}
184:     else if (id == PETSCFV_CLASSID) {PetscFVView((PetscFV) obj, viewer);}
185:     PetscViewerASCIIPopTab(viewer);

187:     for (b = prob->boundary; b; b = b->next) {
188:       const char *name;
189:       PetscInt    c, i;

191:       if (b->field != f) continue;
192:       PetscViewerASCIIPushTab(viewer);
193:       PetscViewerASCIIPrintf(viewer, "Boundary %s (%s) %s\n", b->name, b->lname, DMBoundaryConditionTypes[b->type]);
194:       if (!b->Nc) {
195:         PetscViewerASCIIPrintf(viewer, "  all components\n");
196:       } else {
197:         PetscViewerASCIIPrintf(viewer, "  components: ");
198:         PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
199:         for (c = 0; c < b->Nc; ++c) {
200:           if (c > 0) {PetscViewerASCIIPrintf(viewer, ", ");}
201:           PetscViewerASCIIPrintf(viewer, "%D", b->comps[c]);
202:         }
203:         PetscViewerASCIIPrintf(viewer, "\n");
204:         PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
205:       }
206:       PetscViewerASCIIPrintf(viewer, "  values: ");
207:       PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
208:       for (i = 0; i < b->Nv; ++i) {
209:         if (i > 0) {PetscViewerASCIIPrintf(viewer, ", ");}
210:         PetscViewerASCIIPrintf(viewer, "%D", b->values[i]);
211:       }
212:       PetscViewerASCIIPrintf(viewer, "\n");
213:       PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
214:       if (b->func) {
215:         PetscDLAddr(b->func, &name);
216:         if (name) {PetscViewerASCIIPrintf(viewer, "  func: %s\n", name);}
217:         else      {PetscViewerASCIIPrintf(viewer, "  func: %p\n", b->func);}
218:       }
219:       if (b->func_t) {
220:         PetscDLAddr(b->func_t, &name);
221:         if (name) {PetscViewerASCIIPrintf(viewer, "  func: %s\n", name);}
222:         else      {PetscViewerASCIIPrintf(viewer, "  func: %p\n", b->func_t);}
223:       }
224:       PetscWeakFormView(b->wf, viewer);
225:       PetscViewerASCIIPopTab(viewer);
226:     }
227:   }
228:   PetscDSGetConstants(prob, &numConstants, &constants);
229:   if (numConstants) {
230:     PetscViewerASCIIPrintf(viewer, "%D constants\n", numConstants);
231:     PetscViewerASCIIPushTab(viewer);
232:     for (f = 0; f < numConstants; ++f) {PetscViewerASCIIPrintf(viewer, "%g\n", (double) PetscRealPart(constants[f]));}
233:     PetscViewerASCIIPopTab(viewer);
234:   }
235:   PetscWeakFormView(prob->wf, viewer);
236:   PetscViewerASCIIPopTab(viewer);
237:   return(0);
238: }

240: /*@C
241:    PetscDSViewFromOptions - View from Options

243:    Collective on PetscDS

245:    Input Parameters:
246: +  A - the PetscDS object
247: .  obj - Optional object
248: -  name - command line option

250:    Level: intermediate
251: .seealso:  PetscDS, PetscDSView, PetscObjectViewFromOptions(), PetscDSCreate()
252: @*/
253: PetscErrorCode  PetscDSViewFromOptions(PetscDS A,PetscObject obj,const char name[])
254: {

259:   PetscObjectViewFromOptions((PetscObject)A,obj,name);
260:   return(0);
261: }

263: /*@C
264:   PetscDSView - Views a PetscDS

266:   Collective on prob

268:   Input Parameter:
269: + prob - the PetscDS object to view
270: - v  - the viewer

272:   Level: developer

274: .seealso PetscDSDestroy()
275: @*/
276: PetscErrorCode PetscDSView(PetscDS prob, PetscViewer v)
277: {
278:   PetscBool      iascii;

283:   if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) prob), &v);}
285:   PetscObjectTypeCompare((PetscObject) v, PETSCVIEWERASCII, &iascii);
286:   if (iascii) {PetscDSView_Ascii(prob, v);}
287:   if (prob->ops->view) {(*prob->ops->view)(prob, v);}
288:   return(0);
289: }

291: /*@
292:   PetscDSSetFromOptions - sets parameters in a PetscDS from the options database

294:   Collective on prob

296:   Input Parameter:
297: . prob - the PetscDS object to set options for

299:   Options Database:
300: + -petscds_type <type>     : Set the DS type
301: . -petscds_view <view opt> : View the DS
302: . -petscds_jac_pre         : Turn formation of a separate Jacobian preconditioner on and off
303: . -bc_<name> <ids>         : Specify a list of label ids for a boundary condition
304: - -bc_<name>_comp <comps>  : Specify a list of field components to constrain for a boundary condition

306:   Level: developer

308: .seealso PetscDSView()
309: @*/
310: PetscErrorCode PetscDSSetFromOptions(PetscDS prob)
311: {
312:   DSBoundary     b;
313:   const char    *defaultType;
314:   char           name[256];
315:   PetscBool      flg;

320:   if (!((PetscObject) prob)->type_name) {
321:     defaultType = PETSCDSBASIC;
322:   } else {
323:     defaultType = ((PetscObject) prob)->type_name;
324:   }
325:   PetscDSRegisterAll();

327:   PetscObjectOptionsBegin((PetscObject) prob);
328:   for (b = prob->boundary; b; b = b->next) {
329:     char       optname[1024];
330:     PetscInt   ids[1024], len = 1024;
331:     PetscBool  flg;

333:     PetscSNPrintf(optname, sizeof(optname), "-bc_%s", b->name);
334:     PetscMemzero(ids, sizeof(ids));
335:     PetscOptionsIntArray(optname, "List of boundary IDs", "", ids, &len, &flg);
336:     if (flg) {
337:       b->Nv = len;
338:       PetscFree(b->values);
339:       PetscMalloc1(len, &b->values);
340:       PetscArraycpy(b->values, ids, len);
341:       PetscWeakFormRewriteKeys(b->wf, b->label, len, b->values);
342:     }
343:     len = 1024;
344:     PetscSNPrintf(optname, sizeof(optname), "-bc_%s_comp", b->name);
345:     PetscMemzero(ids, sizeof(ids));
346:     PetscOptionsIntArray(optname, "List of boundary field components", "", ids, &len, &flg);
347:     if (flg) {
348:       b->Nc = len;
349:       PetscFree(b->comps);
350:       PetscMalloc1(len, &b->comps);
351:       PetscArraycpy(b->comps, ids, len);
352:     }
353:   }
354:   PetscOptionsFList("-petscds_type", "Discrete System", "PetscDSSetType", PetscDSList, defaultType, name, 256, &flg);
355:   if (flg) {
356:     PetscDSSetType(prob, name);
357:   } else if (!((PetscObject) prob)->type_name) {
358:     PetscDSSetType(prob, defaultType);
359:   }
360:   PetscOptionsBool("-petscds_jac_pre", "Discrete System", "PetscDSUseJacobianPreconditioner", prob->useJacPre, &prob->useJacPre, &flg);
361:   if (prob->ops->setfromoptions) {(*prob->ops->setfromoptions)(prob);}
362:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
363:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) prob);
364:   PetscOptionsEnd();
365:   if (prob->Nf) {PetscDSViewFromOptions(prob, NULL, "-petscds_view");}
366:   return(0);
367: }

369: /*@C
370:   PetscDSSetUp - Construct data structures for the PetscDS

372:   Collective on prob

374:   Input Parameter:
375: . prob - the PetscDS object to setup

377:   Level: developer

379: .seealso PetscDSView(), PetscDSDestroy()
380: @*/
381: PetscErrorCode PetscDSSetUp(PetscDS prob)
382: {
383:   const PetscInt Nf   = prob->Nf;
384:   PetscBool      hasH = PETSC_FALSE;
385:   PetscInt       dim, dimEmbed, NbMax = 0, NcMax = 0, NqMax = 0, NsMax = 1, f;

390:   if (prob->setup) return(0);
391:   /* Calculate sizes */
392:   PetscDSGetSpatialDimension(prob, &dim);
393:   PetscDSGetCoordinateDimension(prob, &dimEmbed);
394:   prob->totDim = prob->totComp = 0;
395:   PetscMalloc2(Nf,&prob->Nc,Nf,&prob->Nb);
396:   PetscCalloc2(Nf+1,&prob->off,Nf+1,&prob->offDer);
397:   PetscMalloc2(Nf,&prob->T,Nf,&prob->Tf);
398:   for (f = 0; f < Nf; ++f) {
399:     PetscObject     obj;
400:     PetscClassId    id;
401:     PetscQuadrature q = NULL;
402:     PetscInt        Nq = 0, Nb, Nc;

404:     PetscDSGetDiscretization(prob, f, &obj);
405:     if (prob->jetDegree[f] > 1) hasH = PETSC_TRUE;
406:     if (!obj) {
407:       /* Empty mesh */
408:       Nb = Nc = 0;
409:       prob->T[f] = prob->Tf[f] = NULL;
410:     } else {
411:       PetscObjectGetClassId(obj, &id);
412:       if (id == PETSCFE_CLASSID)      {
413:         PetscFE fe = (PetscFE) obj;

415:         PetscFEGetQuadrature(fe, &q);
416:         PetscFEGetDimension(fe, &Nb);
417:         PetscFEGetNumComponents(fe, &Nc);
418:         PetscFEGetCellTabulation(fe, prob->jetDegree[f], &prob->T[f]);
419:         PetscFEGetFaceTabulation(fe, prob->jetDegree[f], &prob->Tf[f]);
420:       } else if (id == PETSCFV_CLASSID) {
421:         PetscFV fv = (PetscFV) obj;

423:         PetscFVGetQuadrature(fv, &q);
424:         PetscFVGetNumComponents(fv, &Nc);
425:         Nb   = Nc;
426:         PetscFVGetCellTabulation(fv, &prob->T[f]);
427:         /* TODO: should PetscFV also have face tabulation? Otherwise there will be a null pointer in prob->basisFace */
428:       } else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
429:     }
430:     prob->Nc[f]       = Nc;
431:     prob->Nb[f]       = Nb;
432:     prob->off[f+1]    = Nc     + prob->off[f];
433:     prob->offDer[f+1] = Nc*dim + prob->offDer[f];
434:     if (q) {PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);}
435:     NqMax          = PetscMax(NqMax, Nq);
436:     NbMax          = PetscMax(NbMax, Nb);
437:     NcMax          = PetscMax(NcMax, Nc);
438:     prob->totDim  += Nb;
439:     prob->totComp += Nc;
440:     /* There are two faces for all fields but the cohesive field on a hybrid cell */
441:     if (prob->isHybrid && (f < Nf-1)) prob->totDim += Nb;
442:   }
443:   /* Allocate works space */
444:   if (prob->isHybrid) NsMax = 2;
445:   PetscMalloc3(NsMax*prob->totComp,&prob->u,NsMax*prob->totComp,&prob->u_t,NsMax*prob->totComp*dimEmbed + (hasH ? NsMax*prob->totComp*dimEmbed*dimEmbed : 0),&prob->u_x);
446:   PetscMalloc5(dimEmbed,&prob->x,NbMax*NcMax,&prob->basisReal,NbMax*NcMax*dimEmbed,&prob->basisDerReal,NbMax*NcMax,&prob->testReal,NbMax*NcMax*dimEmbed,&prob->testDerReal);
447:   PetscMalloc6(NsMax*NqMax*NcMax,&prob->f0,NsMax*NqMax*NcMax*dimEmbed,&prob->f1,
448:                       NsMax*NsMax*NqMax*NcMax*NcMax,&prob->g0,NsMax*NsMax*NqMax*NcMax*NcMax*dimEmbed,&prob->g1,
449:                       NsMax*NsMax*NqMax*NcMax*NcMax*dimEmbed,&prob->g2,NsMax*NsMax*NqMax*NcMax*NcMax*dimEmbed*dimEmbed,&prob->g3);
450:   if (prob->ops->setup) {(*prob->ops->setup)(prob);}
451:   prob->setup = PETSC_TRUE;
452:   return(0);
453: }

455: static PetscErrorCode PetscDSDestroyStructs_Static(PetscDS prob)
456: {

460:   PetscFree2(prob->Nc,prob->Nb);
461:   PetscFree2(prob->off,prob->offDer);
462:   PetscFree2(prob->T,prob->Tf);
463:   PetscFree3(prob->u,prob->u_t,prob->u_x);
464:   PetscFree5(prob->x,prob->basisReal, prob->basisDerReal,prob->testReal,prob->testDerReal);
465:   PetscFree6(prob->f0,prob->f1,prob->g0,prob->g1,prob->g2,prob->g3);
466:   return(0);
467: }

469: static PetscErrorCode PetscDSEnlarge_Static(PetscDS prob, PetscInt NfNew)
470: {
471:   PetscObject      *tmpd;
472:   PetscBool        *tmpi;
473:   PetscInt         *tmpk;
474:   PetscPointFunc   *tmpup;
475:   PetscSimplePointFunc *tmpexactSol,  *tmpexactSol_t;
476:   void                **tmpexactCtx, **tmpexactCtx_t;
477:   void            **tmpctx;
478:   PetscInt          Nf = prob->Nf, f;
479:   PetscErrorCode    ierr;

482:   if (Nf >= NfNew) return(0);
483:   prob->setup = PETSC_FALSE;
484:   PetscDSDestroyStructs_Static(prob);
485:   PetscMalloc3(NfNew, &tmpd, NfNew, &tmpi, NfNew, &tmpk);
486:   for (f = 0; f < Nf; ++f) {tmpd[f] = prob->disc[f]; tmpi[f] = prob->implicit[f]; tmpk[f] = prob->jetDegree[f];}
487:   for (f = Nf; f < NfNew; ++f) {tmpd[f] = NULL; tmpi[f] = PETSC_TRUE; tmpk[f] = 1;}
488:   PetscFree3(prob->disc, prob->implicit, prob->jetDegree);
489:   PetscWeakFormSetNumFields(prob->wf, NfNew);
490:   prob->Nf        = NfNew;
491:   prob->disc      = tmpd;
492:   prob->implicit  = tmpi;
493:   prob->jetDegree = tmpk;
494:   PetscCalloc2(NfNew, &tmpup, NfNew, &tmpctx);
495:   for (f = 0; f < Nf; ++f) tmpup[f] = prob->update[f];
496:   for (f = 0; f < Nf; ++f) tmpctx[f] = prob->ctx[f];
497:   for (f = Nf; f < NfNew; ++f) tmpup[f] = NULL;
498:   for (f = Nf; f < NfNew; ++f) tmpctx[f] = NULL;
499:   PetscFree2(prob->update, prob->ctx);
500:   prob->update = tmpup;
501:   prob->ctx = tmpctx;
502:   PetscCalloc4(NfNew, &tmpexactSol, NfNew, &tmpexactCtx, NfNew, &tmpexactSol_t, NfNew, &tmpexactCtx_t);
503:   for (f = 0; f < Nf; ++f) tmpexactSol[f] = prob->exactSol[f];
504:   for (f = 0; f < Nf; ++f) tmpexactCtx[f] = prob->exactCtx[f];
505:   for (f = 0; f < Nf; ++f) tmpexactSol_t[f] = prob->exactSol_t[f];
506:   for (f = 0; f < Nf; ++f) tmpexactCtx_t[f] = prob->exactCtx_t[f];
507:   for (f = Nf; f < NfNew; ++f) tmpexactSol[f] = NULL;
508:   for (f = Nf; f < NfNew; ++f) tmpexactCtx[f] = NULL;
509:   for (f = Nf; f < NfNew; ++f) tmpexactSol_t[f] = NULL;
510:   for (f = Nf; f < NfNew; ++f) tmpexactCtx_t[f] = NULL;
511:   PetscFree4(prob->exactSol, prob->exactCtx, prob->exactSol_t, prob->exactCtx_t);
512:   prob->exactSol = tmpexactSol;
513:   prob->exactCtx = tmpexactCtx;
514:   prob->exactSol_t = tmpexactSol_t;
515:   prob->exactCtx_t = tmpexactCtx_t;
516:   return(0);
517: }

519: /*@
520:   PetscDSDestroy - Destroys a PetscDS object

522:   Collective on prob

524:   Input Parameter:
525: . prob - the PetscDS object to destroy

527:   Level: developer

529: .seealso PetscDSView()
530: @*/
531: PetscErrorCode PetscDSDestroy(PetscDS *ds)
532: {
533:   PetscInt       f;

537:   if (!*ds) return(0);

540:   if (--((PetscObject)(*ds))->refct > 0) {*ds = NULL; return(0);}
541:   ((PetscObject) (*ds))->refct = 0;
542:   if ((*ds)->subprobs) {
543:     PetscInt dim, d;

545:     PetscDSGetSpatialDimension(*ds, &dim);
546:     for (d = 0; d < dim; ++d) {PetscDSDestroy(&(*ds)->subprobs[d]);}
547:   }
548:   PetscFree((*ds)->subprobs);
549:   PetscDSDestroyStructs_Static(*ds);
550:   for (f = 0; f < (*ds)->Nf; ++f) {
551:     PetscObjectDereference((*ds)->disc[f]);
552:   }
553:   PetscFree3((*ds)->disc, (*ds)->implicit, (*ds)->jetDegree);
554:   PetscWeakFormDestroy(&(*ds)->wf);
555:   PetscFree2((*ds)->update,(*ds)->ctx);
556:   PetscFree4((*ds)->exactSol,(*ds)->exactCtx,(*ds)->exactSol_t,(*ds)->exactCtx_t);
557:   if ((*ds)->ops->destroy) {(*(*ds)->ops->destroy)(*ds);}
558:   PetscDSDestroyBoundary(*ds);
559:   PetscFree((*ds)->constants);
560:   PetscHeaderDestroy(ds);
561:   return(0);
562: }

564: /*@
565:   PetscDSCreate - Creates an empty PetscDS object. The type can then be set with PetscDSSetType().

567:   Collective

569:   Input Parameter:
570: . comm - The communicator for the PetscDS object

572:   Output Parameter:
573: . ds   - The PetscDS object

575:   Level: beginner

577: .seealso: PetscDSSetType(), PETSCDSBASIC
578: @*/
579: PetscErrorCode PetscDSCreate(MPI_Comm comm, PetscDS *ds)
580: {
581:   PetscDS        p;

586:   *ds  = NULL;
587:   PetscDSInitializePackage();

589:   PetscHeaderCreate(p, PETSCDS_CLASSID, "PetscDS", "Discrete System", "PetscDS", comm, PetscDSDestroy, PetscDSView);

591:   p->Nf           = 0;
592:   p->setup        = PETSC_FALSE;
593:   p->numConstants = 0;
594:   p->constants    = NULL;
595:   p->dimEmbed     = -1;
596:   p->useJacPre    = PETSC_TRUE;
597:   PetscWeakFormCreate(comm, &p->wf);

599:   *ds = p;
600:   return(0);
601: }

603: /*@
604:   PetscDSGetNumFields - Returns the number of fields in the DS

606:   Not collective

608:   Input Parameter:
609: . prob - The PetscDS object

611:   Output Parameter:
612: . Nf - The number of fields

614:   Level: beginner

616: .seealso: PetscDSGetSpatialDimension(), PetscDSCreate()
617: @*/
618: PetscErrorCode PetscDSGetNumFields(PetscDS prob, PetscInt *Nf)
619: {
623:   *Nf = prob->Nf;
624:   return(0);
625: }

627: /*@
628:   PetscDSGetSpatialDimension - Returns the spatial dimension of the DS, meaning the topological dimension of the discretizations

630:   Not collective

632:   Input Parameter:
633: . prob - The PetscDS object

635:   Output Parameter:
636: . dim - The spatial dimension

638:   Level: beginner

640: .seealso: PetscDSGetCoordinateDimension(), PetscDSGetNumFields(), PetscDSCreate()
641: @*/
642: PetscErrorCode PetscDSGetSpatialDimension(PetscDS prob, PetscInt *dim)
643: {

649:   *dim = 0;
650:   if (prob->Nf) {
651:     PetscObject  obj;
652:     PetscClassId id;

654:     PetscDSGetDiscretization(prob, 0, &obj);
655:     if (obj) {
656:       PetscObjectGetClassId(obj, &id);
657:       if (id == PETSCFE_CLASSID)      {PetscFEGetSpatialDimension((PetscFE) obj, dim);}
658:       else if (id == PETSCFV_CLASSID) {PetscFVGetSpatialDimension((PetscFV) obj, dim);}
659:       else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", 0);
660:     }
661:   }
662:   return(0);
663: }

665: /*@
666:   PetscDSGetCoordinateDimension - Returns the coordinate dimension of the DS, meaning the dimension of the space into which the discretiaztions are embedded

668:   Not collective

670:   Input Parameter:
671: . prob - The PetscDS object

673:   Output Parameter:
674: . dimEmbed - The coordinate dimension

676:   Level: beginner

678: .seealso: PetscDSSetCoordinateDimension(), PetscDSGetSpatialDimension(), PetscDSGetNumFields(), PetscDSCreate()
679: @*/
680: PetscErrorCode PetscDSGetCoordinateDimension(PetscDS prob, PetscInt *dimEmbed)
681: {
685:   if (prob->dimEmbed < 0) SETERRQ(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONGSTATE, "No coordinate dimension set for this DS");
686:   *dimEmbed = prob->dimEmbed;
687:   return(0);
688: }

690: /*@
691:   PetscDSSetCoordinateDimension - Set the coordinate dimension of the DS, meaning the dimension of the space into which the discretiaztions are embedded

693:   Logically collective on prob

695:   Input Parameters:
696: + prob - The PetscDS object
697: - dimEmbed - The coordinate dimension

699:   Level: beginner

701: .seealso: PetscDSGetCoordinateDimension(), PetscDSGetSpatialDimension(), PetscDSGetNumFields(), PetscDSCreate()
702: @*/
703: PetscErrorCode PetscDSSetCoordinateDimension(PetscDS prob, PetscInt dimEmbed)
704: {
707:   if (dimEmbed < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Coordinate dimension must be non-negative, not %D", dimEmbed);
708:   prob->dimEmbed = dimEmbed;
709:   return(0);
710: }

712: /*@
713:   PetscDSGetHybrid - Returns the flag for a hybrid (cohesive) cell

715:   Not collective

717:   Input Parameter:
718: . prob - The PetscDS object

720:   Output Parameter:
721: . isHybrid - The flag

723:   Level: developer

725: .seealso: PetscDSSetHybrid(), PetscDSCreate()
726: @*/
727: PetscErrorCode PetscDSGetHybrid(PetscDS prob, PetscBool *isHybrid)
728: {
732:   *isHybrid = prob->isHybrid;
733:   return(0);
734: }

736: /*@
737:   PetscDSSetHybrid - Set the flag for a hybrid (cohesive) cell

739:   Not collective

741:   Input Parameters:
742: + prob - The PetscDS object
743: - isHybrid - The flag

745:   Level: developer

747: .seealso: PetscDSGetHybrid(), PetscDSCreate()
748: @*/
749: PetscErrorCode PetscDSSetHybrid(PetscDS prob, PetscBool isHybrid)
750: {
753:   prob->isHybrid = isHybrid;
754:   return(0);
755: }

757: /*@
758:   PetscDSGetTotalDimension - Returns the total size of the approximation space for this system

760:   Not collective

762:   Input Parameter:
763: . prob - The PetscDS object

765:   Output Parameter:
766: . dim - The total problem dimension

768:   Level: beginner

770: .seealso: PetscDSGetNumFields(), PetscDSCreate()
771: @*/
772: PetscErrorCode PetscDSGetTotalDimension(PetscDS prob, PetscInt *dim)
773: {

778:   PetscDSSetUp(prob);
780:   *dim = prob->totDim;
781:   return(0);
782: }

784: /*@
785:   PetscDSGetTotalComponents - Returns the total number of components in this system

787:   Not collective

789:   Input Parameter:
790: . prob - The PetscDS object

792:   Output Parameter:
793: . dim - The total number of components

795:   Level: beginner

797: .seealso: PetscDSGetNumFields(), PetscDSCreate()
798: @*/
799: PetscErrorCode PetscDSGetTotalComponents(PetscDS prob, PetscInt *Nc)
800: {

805:   PetscDSSetUp(prob);
807:   *Nc = prob->totComp;
808:   return(0);
809: }

811: /*@
812:   PetscDSGetDiscretization - Returns the discretization object for the given field

814:   Not collective

816:   Input Parameters:
817: + prob - The PetscDS object
818: - f - The field number

820:   Output Parameter:
821: . disc - The discretization object

823:   Level: beginner

825: .seealso: PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
826: @*/
827: PetscErrorCode PetscDSGetDiscretization(PetscDS prob, PetscInt f, PetscObject *disc)
828: {
832:   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);
833:   *disc = prob->disc[f];
834:   return(0);
835: }

837: /*@
838:   PetscDSSetDiscretization - Sets the discretization object for the given field

840:   Not collective

842:   Input Parameters:
843: + prob - The PetscDS object
844: . f - The field number
845: - disc - The discretization object

847:   Level: beginner

849: .seealso: PetscDSGetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
850: @*/
851: PetscErrorCode PetscDSSetDiscretization(PetscDS prob, PetscInt f, PetscObject disc)
852: {

858:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
859:   PetscDSEnlarge_Static(prob, f+1);
860:   PetscObjectDereference(prob->disc[f]);
861:   prob->disc[f] = disc;
862:   PetscObjectReference(disc);
863:   if (disc) {
864:     PetscClassId id;

866:     PetscObjectGetClassId(disc, &id);
867:     if (id == PETSCFE_CLASSID) {
868:       PetscDSSetImplicit(prob, f, PETSC_TRUE);
869:     } else if (id == PETSCFV_CLASSID) {
870:       PetscDSSetImplicit(prob, f, PETSC_FALSE);
871:     }
872:     PetscDSSetJetDegree(prob, f, 1);
873:   }
874:   return(0);
875: }

877: /*@
878:   PetscDSGetWeakForm - Returns the weak form object

880:   Not collective

882:   Input Parameter:
883: . ds - The PetscDS object

885:   Output Parameter:
886: . wf - The weak form object

888:   Level: beginner

890: .seealso: PetscDSSetWeakForm(), PetscDSGetNumFields(), PetscDSCreate()
891: @*/
892: PetscErrorCode PetscDSGetWeakForm(PetscDS ds, PetscWeakForm *wf)
893: {
897:   *wf = ds->wf;
898:   return(0);
899: }

901: /*@
902:   PetscDSSetWeakForm - Sets the weak form object

904:   Not collective

906:   Input Parameters:
907: + ds - The PetscDS object
908: - wf - The weak form object

910:   Level: beginner

912: .seealso: PetscDSGetWeakForm(), PetscDSGetNumFields(), PetscDSCreate()
913: @*/
914: PetscErrorCode PetscDSSetWeakForm(PetscDS ds, PetscWeakForm wf)
915: {

921:   PetscObjectDereference((PetscObject) ds->wf);
922:   ds->wf = wf;
923:   PetscObjectReference((PetscObject) wf);
924:   PetscWeakFormSetNumFields(wf, ds->Nf);
925:   return(0);
926: }

928: /*@
929:   PetscDSAddDiscretization - Adds a discretization object

931:   Not collective

933:   Input Parameters:
934: + prob - The PetscDS object
935: - disc - The boundary discretization object

937:   Level: beginner

939: .seealso: PetscDSGetDiscretization(), PetscDSSetDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
940: @*/
941: PetscErrorCode PetscDSAddDiscretization(PetscDS prob, PetscObject disc)
942: {

946:   PetscDSSetDiscretization(prob, prob->Nf, disc);
947:   return(0);
948: }

950: /*@
951:   PetscDSGetQuadrature - Returns the quadrature, which must agree for all fields in the DS

953:   Not collective

955:   Input Parameter:
956: . prob - The PetscDS object

958:   Output Parameter:
959: . q - The quadrature object

961: Level: intermediate

963: .seealso: PetscDSSetImplicit(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
964: @*/
965: PetscErrorCode PetscDSGetQuadrature(PetscDS prob, PetscQuadrature *q)
966: {
967:   PetscObject    obj;
968:   PetscClassId   id;

972:   *q = NULL;
973:   if (!prob->Nf) return(0);
974:   PetscDSGetDiscretization(prob, 0, &obj);
975:   PetscObjectGetClassId(obj, &id);
976:   if      (id == PETSCFE_CLASSID) {PetscFEGetQuadrature((PetscFE) obj, q);}
977:   else if (id == PETSCFV_CLASSID) {PetscFVGetQuadrature((PetscFV) obj, q);}
978:   else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", 0);
979:   return(0);
980: }

982: /*@
983:   PetscDSGetImplicit - Returns the flag for implicit solve for this field. This is just a guide for IMEX

985:   Not collective

987:   Input Parameters:
988: + prob - The PetscDS object
989: - f - The field number

991:   Output Parameter:
992: . implicit - The flag indicating what kind of solve to use for this field

994:   Level: developer

996: .seealso: PetscDSSetImplicit(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
997: @*/
998: PetscErrorCode PetscDSGetImplicit(PetscDS prob, PetscInt f, PetscBool *implicit)
999: {
1003:   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);
1004:   *implicit = prob->implicit[f];
1005:   return(0);
1006: }

1008: /*@
1009:   PetscDSSetImplicit - Set the flag for implicit solve for this field. This is just a guide for IMEX

1011:   Not collective

1013:   Input Parameters:
1014: + prob - The PetscDS object
1015: . f - The field number
1016: - implicit - The flag indicating what kind of solve to use for this field

1018:   Level: developer

1020: .seealso: PetscDSGetImplicit(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
1021: @*/
1022: PetscErrorCode PetscDSSetImplicit(PetscDS prob, PetscInt f, PetscBool implicit)
1023: {
1026:   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);
1027:   prob->implicit[f] = implicit;
1028:   return(0);
1029: }

1031: /*@
1032:   PetscDSGetJetDegree - Returns the highest derivative for this field equation, or the k-jet that the discretization needs to tabulate.

1034:   Not collective

1036:   Input Parameters:
1037: + ds - The PetscDS object
1038: - f  - The field number

1040:   Output Parameter:
1041: . k  - The highest derivative we need to tabulate

1043:   Level: developer

1045: .seealso: PetscDSSetJetDegree(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
1046: @*/
1047: PetscErrorCode PetscDSGetJetDegree(PetscDS ds, PetscInt f, PetscInt *k)
1048: {
1052:   if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
1053:   *k = ds->jetDegree[f];
1054:   return(0);
1055: }

1057: /*@
1058:   PetscDSSetJetDegree - Set the highest derivative for this field equation, or the k-jet that the discretization needs to tabulate.

1060:   Not collective

1062:   Input Parameters:
1063: + ds - The PetscDS object
1064: . f  - The field number
1065: - k  - The highest derivative we need to tabulate

1067:   Level: developer

1069: .seealso: PetscDSGetJetDegree(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
1070: @*/
1071: PetscErrorCode PetscDSSetJetDegree(PetscDS ds, PetscInt f, PetscInt k)
1072: {
1075:   if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
1076:   ds->jetDegree[f] = k;
1077:   return(0);
1078: }

1080: PetscErrorCode PetscDSGetObjective(PetscDS ds, PetscInt f,
1081:                                    void (**obj)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1082:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1083:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1084:                                                 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[]))
1085: {
1086:   PetscPointFunc *tmp;
1087:   PetscInt        n;
1088:   PetscErrorCode  ierr;

1093:   if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
1094:   PetscWeakFormGetObjective(ds->wf, NULL, 0, f, &n, &tmp);
1095:   *obj = tmp ? tmp[0] : NULL;
1096:   return(0);
1097: }

1099: PetscErrorCode PetscDSSetObjective(PetscDS ds, PetscInt f,
1100:                                    void (*obj)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1101:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1102:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1103:                                                PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[]))
1104: {

1110:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1111:   PetscWeakFormSetIndexObjective(ds->wf, NULL, 0, f, 0, obj);
1112:   return(0);
1113: }

1115: /*@C
1116:   PetscDSGetResidual - Get the pointwise residual function for a given test field

1118:   Not collective

1120:   Input Parameters:
1121: + ds - The PetscDS
1122: - f  - The test field number

1124:   Output Parameters:
1125: + f0 - integrand for the test function term
1126: - f1 - integrand for the test function gradient term

1128:   Note: We are using a first order FEM model for the weak form:

1130:   \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)

1132: The calling sequence for the callbacks f0 and f1 is given by:

1134: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1135: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1136: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1137: $    PetscReal t, const PetscReal x[], PetscScalar f0[])

1139: + dim - the spatial dimension
1140: . Nf - the number of fields
1141: . uOff - the offset into u[] and u_t[] for each field
1142: . uOff_x - the offset into u_x[] for each field
1143: . u - each field evaluated at the current point
1144: . u_t - the time derivative of each field evaluated at the current point
1145: . u_x - the gradient of each field evaluated at the current point
1146: . aOff - the offset into a[] and a_t[] for each auxiliary field
1147: . aOff_x - the offset into a_x[] for each auxiliary field
1148: . a - each auxiliary field evaluated at the current point
1149: . a_t - the time derivative of each auxiliary field evaluated at the current point
1150: . a_x - the gradient of auxiliary each field evaluated at the current point
1151: . t - current time
1152: . x - coordinates of the current point
1153: . numConstants - number of constant parameters
1154: . constants - constant parameters
1155: - f0 - output values at the current point

1157:   Level: intermediate

1159: .seealso: PetscDSSetResidual()
1160: @*/
1161: PetscErrorCode PetscDSGetResidual(PetscDS ds, PetscInt f,
1162:                                   void (**f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1163:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1164:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1165:                                               PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]),
1166:                                   void (**f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1167:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1168:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1169:                                               PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
1170: {
1171:   PetscPointFunc *tmp0, *tmp1;
1172:   PetscInt        n0, n1;
1173:   PetscErrorCode  ierr;

1177:   if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
1178:   PetscWeakFormGetResidual(ds->wf, NULL, 0, f, &n0, &tmp0, &n1, &tmp1);
1179:   *f0  = tmp0 ? tmp0[0] : NULL;
1180:   *f1  = tmp1 ? tmp1[0] : NULL;
1181:   return(0);
1182: }

1184: /*@C
1185:   PetscDSSetResidual - Set the pointwise residual function for a given test field

1187:   Not collective

1189:   Input Parameters:
1190: + ds - The PetscDS
1191: . f  - The test field number
1192: . f0 - integrand for the test function term
1193: - f1 - integrand for the test function gradient term

1195:   Note: We are using a first order FEM model for the weak form:

1197:   \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)

1199: The calling sequence for the callbacks f0 and f1 is given by:

1201: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1202: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1203: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1204: $    PetscReal t, const PetscReal x[], PetscScalar f0[])

1206: + dim - the spatial dimension
1207: . Nf - the number of fields
1208: . uOff - the offset into u[] and u_t[] for each field
1209: . uOff_x - the offset into u_x[] for each field
1210: . u - each field evaluated at the current point
1211: . u_t - the time derivative of each field evaluated at the current point
1212: . u_x - the gradient of each field evaluated at the current point
1213: . aOff - the offset into a[] and a_t[] for each auxiliary field
1214: . aOff_x - the offset into a_x[] for each auxiliary field
1215: . a - each auxiliary field evaluated at the current point
1216: . a_t - the time derivative of each auxiliary field evaluated at the current point
1217: . a_x - the gradient of auxiliary each field evaluated at the current point
1218: . t - current time
1219: . x - coordinates of the current point
1220: . numConstants - number of constant parameters
1221: . constants - constant parameters
1222: - f0 - output values at the current point

1224:   Level: intermediate

1226: .seealso: PetscDSGetResidual()
1227: @*/
1228: PetscErrorCode PetscDSSetResidual(PetscDS ds, PetscInt f,
1229:                                   void (*f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1230:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1231:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1232:                                              PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]),
1233:                                   void (*f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1234:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1235:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1236:                                              PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
1237: {

1244:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1245:   PetscWeakFormSetIndexResidual(ds->wf, NULL, 0, f, 0, f0, 0, f1);
1246:   return(0);
1247: }

1249: /*@C
1250:   PetscDSHasJacobian - Signals that Jacobian functions have been set

1252:   Not collective

1254:   Input Parameter:
1255: . prob - The PetscDS

1257:   Output Parameter:
1258: . hasJac - flag that pointwise function for the Jacobian has been set

1260:   Level: intermediate

1262: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1263: @*/
1264: PetscErrorCode PetscDSHasJacobian(PetscDS ds, PetscBool *hasJac)
1265: {

1270:   PetscWeakFormHasJacobian(ds->wf, hasJac);
1271:   return(0);
1272: }

1274: /*@C
1275:   PetscDSGetJacobian - Get the pointwise Jacobian function for given test and basis field

1277:   Not collective

1279:   Input Parameters:
1280: + ds - The PetscDS
1281: . f  - The test field number
1282: - g  - The field number

1284:   Output Parameters:
1285: + g0 - integrand for the test and basis function term
1286: . g1 - integrand for the test function and basis function gradient term
1287: . g2 - integrand for the test function gradient and basis function term
1288: - g3 - integrand for the test function gradient and basis function gradient term

1290:   Note: We are using a first order FEM model for the weak form:

1292:   \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

1294: The calling sequence for the callbacks g0, g1, g2 and g3 is given by:

1296: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1297: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1298: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1299: $    PetscReal t, const PetscReal u_tShift, const PetscReal x[], PetscScalar g0[])

1301: + dim - the spatial dimension
1302: . Nf - the number of fields
1303: . uOff - the offset into u[] and u_t[] for each field
1304: . uOff_x - the offset into u_x[] for each field
1305: . u - each field evaluated at the current point
1306: . u_t - the time derivative of each field evaluated at the current point
1307: . u_x - the gradient of each field evaluated at the current point
1308: . aOff - the offset into a[] and a_t[] for each auxiliary field
1309: . aOff_x - the offset into a_x[] for each auxiliary field
1310: . a - each auxiliary field evaluated at the current point
1311: . a_t - the time derivative of each auxiliary field evaluated at the current point
1312: . a_x - the gradient of auxiliary each field evaluated at the current point
1313: . t - current time
1314: . u_tShift - the multiplier a for dF/dU_t
1315: . x - coordinates of the current point
1316: . numConstants - number of constant parameters
1317: . constants - constant parameters
1318: - g0 - output values at the current point

1320:   Level: intermediate

1322: .seealso: PetscDSSetJacobian()
1323: @*/
1324: PetscErrorCode PetscDSGetJacobian(PetscDS ds, PetscInt f, PetscInt g,
1325:                                   void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1326:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1327:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1328:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1329:                                   void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1330:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1331:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1332:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1333:                                   void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1334:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1335:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1336:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1337:                                   void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1338:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1339:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1340:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1341: {
1342:   PetscPointJac *tmp0, *tmp1, *tmp2, *tmp3;
1343:   PetscInt       n0, n1, n2, n3;

1348:   if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
1349:   if ((g < 0) || (g >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, ds->Nf);
1350:   PetscWeakFormGetJacobian(ds->wf, NULL, 0, f, g, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3);
1351:   *g0  = tmp0 ? tmp0[0] : NULL;
1352:   *g1  = tmp1 ? tmp1[0] : NULL;
1353:   *g2  = tmp2 ? tmp2[0] : NULL;
1354:   *g3  = tmp3 ? tmp3[0] : NULL;
1355:   return(0);
1356: }

1358: /*@C
1359:   PetscDSSetJacobian - Set the pointwise Jacobian function for given test and basis fields

1361:   Not collective

1363:   Input Parameters:
1364: + ds - The PetscDS
1365: . f  - The test field number
1366: . g  - The field number
1367: . g0 - integrand for the test and basis function term
1368: . g1 - integrand for the test function and basis function gradient term
1369: . g2 - integrand for the test function gradient and basis function term
1370: - g3 - integrand for the test function gradient and basis function gradient term

1372:   Note: We are using a first order FEM model for the weak form:

1374:   \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

1376: The calling sequence for the callbacks g0, g1, g2 and g3 is given by:

1378: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1379: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1380: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1381: $    PetscReal t, const PetscReal x[], PetscScalar g0[])

1383: + dim - the spatial dimension
1384: . Nf - the number of fields
1385: . uOff - the offset into u[] and u_t[] for each field
1386: . uOff_x - the offset into u_x[] for each field
1387: . u - each field evaluated at the current point
1388: . u_t - the time derivative of each field evaluated at the current point
1389: . u_x - the gradient of each field evaluated at the current point
1390: . aOff - the offset into a[] and a_t[] for each auxiliary field
1391: . aOff_x - the offset into a_x[] for each auxiliary field
1392: . a - each auxiliary field evaluated at the current point
1393: . a_t - the time derivative of each auxiliary field evaluated at the current point
1394: . a_x - the gradient of auxiliary each field evaluated at the current point
1395: . t - current time
1396: . u_tShift - the multiplier a for dF/dU_t
1397: . x - coordinates of the current point
1398: . numConstants - number of constant parameters
1399: . constants - constant parameters
1400: - g0 - output values at the current point

1402:   Level: intermediate

1404: .seealso: PetscDSGetJacobian()
1405: @*/
1406: PetscErrorCode PetscDSSetJacobian(PetscDS ds, PetscInt f, PetscInt g,
1407:                                   void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1408:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1409:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1410:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1411:                                   void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1412:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1413:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1414:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1415:                                   void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1416:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1417:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1418:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1419:                                   void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1420:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1421:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1422:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1423: {

1432:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1433:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
1434:   PetscWeakFormSetIndexJacobian(ds->wf, NULL, 0, f, g, 0, g0, 0, g1, 0, g2, 0, g3);
1435:   return(0);
1436: }

1438: /*@C
1439:   PetscDSUseJacobianPreconditioner - Whether to construct a Jacobian preconditioner

1441:   Not collective

1443:   Input Parameters:
1444: + prob - The PetscDS
1445: - useJacPre - flag that enables construction of a Jacobian preconditioner

1447:   Level: intermediate

1449: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1450: @*/
1451: PetscErrorCode PetscDSUseJacobianPreconditioner(PetscDS prob, PetscBool useJacPre)
1452: {
1455:   prob->useJacPre = useJacPre;
1456:   return(0);
1457: }

1459: /*@C
1460:   PetscDSHasJacobianPreconditioner - Signals that a Jacobian preconditioner matrix has been set

1462:   Not collective

1464:   Input Parameter:
1465: . prob - The PetscDS

1467:   Output Parameter:
1468: . hasJacPre - flag that pointwise function for Jacobian preconditioner matrix has been set

1470:   Level: intermediate

1472: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1473: @*/
1474: PetscErrorCode PetscDSHasJacobianPreconditioner(PetscDS ds, PetscBool *hasJacPre)
1475: {

1480:   *hasJacPre = PETSC_FALSE;
1481:   if (!ds->useJacPre) return(0);
1482:   PetscWeakFormHasJacobianPreconditioner(ds->wf, hasJacPre);
1483:   return(0);
1484: }

1486: /*@C
1487:   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.

1489:   Not collective

1491:   Input Parameters:
1492: + ds - The PetscDS
1493: . f  - The test field number
1494: - g  - The field number

1496:   Output Parameters:
1497: + g0 - integrand for the test and basis function term
1498: . g1 - integrand for the test function and basis function gradient term
1499: . g2 - integrand for the test function gradient and basis function term
1500: - g3 - integrand for the test function gradient and basis function gradient term

1502:   Note: We are using a first order FEM model for the weak form:

1504:   \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

1506: The calling sequence for the callbacks g0, g1, g2 and g3 is given by:

1508: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1509: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1510: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1511: $    PetscReal t, const PetscReal u_tShift, const PetscReal x[], PetscScalar g0[])

1513: + dim - the spatial dimension
1514: . Nf - the number of fields
1515: . uOff - the offset into u[] and u_t[] for each field
1516: . uOff_x - the offset into u_x[] for each field
1517: . u - each field evaluated at the current point
1518: . u_t - the time derivative of each field evaluated at the current point
1519: . u_x - the gradient of each field evaluated at the current point
1520: . aOff - the offset into a[] and a_t[] for each auxiliary field
1521: . aOff_x - the offset into a_x[] for each auxiliary field
1522: . a - each auxiliary field evaluated at the current point
1523: . a_t - the time derivative of each auxiliary field evaluated at the current point
1524: . a_x - the gradient of auxiliary each field evaluated at the current point
1525: . t - current time
1526: . u_tShift - the multiplier a for dF/dU_t
1527: . x - coordinates of the current point
1528: . numConstants - number of constant parameters
1529: . constants - constant parameters
1530: - g0 - output values at the current point

1532:   Level: intermediate

1534: .seealso: PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1535: @*/
1536: PetscErrorCode PetscDSGetJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g,
1537:                                   void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1538:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1539:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1540:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1541:                                   void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1542:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1543:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1544:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1545:                                   void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1546:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1547:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1548:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1549:                                   void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1550:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1551:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1552:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1553: {
1554:   PetscPointJac *tmp0, *tmp1, *tmp2, *tmp3;
1555:   PetscInt       n0, n1, n2, n3;

1560:   if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
1561:   if ((g < 0) || (g >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, ds->Nf);
1562:   PetscWeakFormGetJacobianPreconditioner(ds->wf, NULL, 0, f, g, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3);
1563:   *g0  = tmp0 ? tmp0[0] : NULL;
1564:   *g1  = tmp1 ? tmp1[0] : NULL;
1565:   *g2  = tmp2 ? tmp2[0] : NULL;
1566:   *g3  = tmp3 ? tmp3[0] : NULL;
1567:   return(0);
1568: }

1570: /*@C
1571:   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.

1573:   Not collective

1575:   Input Parameters:
1576: + ds - The PetscDS
1577: . f  - The test field number
1578: . g  - The field number
1579: . g0 - integrand for the test and basis function term
1580: . g1 - integrand for the test function and basis function gradient term
1581: . g2 - integrand for the test function gradient and basis function term
1582: - g3 - integrand for the test function gradient and basis function gradient term

1584:   Note: We are using a first order FEM model for the weak form:

1586:   \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

1588: The calling sequence for the callbacks g0, g1, g2 and g3 is given by:

1590: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1591: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1592: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1593: $    PetscReal t, const PetscReal x[], PetscScalar g0[])

1595: + dim - the spatial dimension
1596: . Nf - the number of fields
1597: . uOff - the offset into u[] and u_t[] for each field
1598: . uOff_x - the offset into u_x[] for each field
1599: . u - each field evaluated at the current point
1600: . u_t - the time derivative of each field evaluated at the current point
1601: . u_x - the gradient of each field evaluated at the current point
1602: . aOff - the offset into a[] and a_t[] for each auxiliary field
1603: . aOff_x - the offset into a_x[] for each auxiliary field
1604: . a - each auxiliary field evaluated at the current point
1605: . a_t - the time derivative of each auxiliary field evaluated at the current point
1606: . a_x - the gradient of auxiliary each field evaluated at the current point
1607: . t - current time
1608: . u_tShift - the multiplier a for dF/dU_t
1609: . x - coordinates of the current point
1610: . numConstants - number of constant parameters
1611: . constants - constant parameters
1612: - g0 - output values at the current point

1614:   Level: intermediate

1616: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobian()
1617: @*/
1618: PetscErrorCode PetscDSSetJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g,
1619:                                   void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1620:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1621:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1622:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1623:                                   void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1624:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1625:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1626:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1627:                                   void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1628:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1629:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1630:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1631:                                   void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1632:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1633:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1634:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1635: {

1644:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1645:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
1646:   PetscWeakFormSetIndexJacobianPreconditioner(ds->wf, NULL, 0, f, g, 0, g0, 0, g1, 0, g2, 0, g3);
1647:   return(0);
1648: }

1650: /*@C
1651:   PetscDSHasDynamicJacobian - Signals that a dynamic Jacobian, dF/du_t, has been set

1653:   Not collective

1655:   Input Parameter:
1656: . ds - The PetscDS

1658:   Output Parameter:
1659: . hasDynJac - flag that pointwise function for dynamic Jacobian has been set

1661:   Level: intermediate

1663: .seealso: PetscDSGetDynamicJacobian(), PetscDSSetDynamicJacobian(), PetscDSGetJacobian()
1664: @*/
1665: PetscErrorCode PetscDSHasDynamicJacobian(PetscDS ds, PetscBool *hasDynJac)
1666: {

1671:   PetscWeakFormHasDynamicJacobian(ds->wf, hasDynJac);
1672:   return(0);
1673: }

1675: /*@C
1676:   PetscDSGetDynamicJacobian - Get the pointwise dynamic Jacobian, dF/du_t, function for given test and basis field

1678:   Not collective

1680:   Input Parameters:
1681: + ds - The PetscDS
1682: . f  - The test field number
1683: - g  - The field number

1685:   Output Parameters:
1686: + g0 - integrand for the test and basis function term
1687: . g1 - integrand for the test function and basis function gradient term
1688: . g2 - integrand for the test function gradient and basis function term
1689: - g3 - integrand for the test function gradient and basis function gradient term

1691:   Note: We are using a first order FEM model for the weak form:

1693:   \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

1695: The calling sequence for the callbacks g0, g1, g2 and g3 is given by:

1697: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1698: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1699: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1700: $    PetscReal t, const PetscReal u_tShift, const PetscReal x[], PetscScalar g0[])

1702: + dim - the spatial dimension
1703: . Nf - the number of fields
1704: . uOff - the offset into u[] and u_t[] for each field
1705: . uOff_x - the offset into u_x[] for each field
1706: . u - each field evaluated at the current point
1707: . u_t - the time derivative of each field evaluated at the current point
1708: . u_x - the gradient of each field evaluated at the current point
1709: . aOff - the offset into a[] and a_t[] for each auxiliary field
1710: . aOff_x - the offset into a_x[] for each auxiliary field
1711: . a - each auxiliary field evaluated at the current point
1712: . a_t - the time derivative of each auxiliary field evaluated at the current point
1713: . a_x - the gradient of auxiliary each field evaluated at the current point
1714: . t - current time
1715: . u_tShift - the multiplier a for dF/dU_t
1716: . x - coordinates of the current point
1717: . numConstants - number of constant parameters
1718: . constants - constant parameters
1719: - g0 - output values at the current point

1721:   Level: intermediate

1723: .seealso: PetscDSSetJacobian()
1724: @*/
1725: PetscErrorCode PetscDSGetDynamicJacobian(PetscDS ds, PetscInt f, PetscInt g,
1726:                                          void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1727:                                                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1728:                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1729:                                                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1730:                                          void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1731:                                                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1732:                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1733:                                                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1734:                                          void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1735:                                                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1736:                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1737:                                                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1738:                                          void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1739:                                                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1740:                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1741:                                                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1742: {
1743:   PetscPointJac *tmp0, *tmp1, *tmp2, *tmp3;
1744:   PetscInt       n0, n1, n2, n3;

1749:   if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
1750:   if ((g < 0) || (g >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, ds->Nf);
1751:   PetscWeakFormGetDynamicJacobian(ds->wf, NULL, 0, f, g, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3);
1752:   *g0  = tmp0 ? tmp0[0] : NULL;
1753:   *g1  = tmp1 ? tmp1[0] : NULL;
1754:   *g2  = tmp2 ? tmp2[0] : NULL;
1755:   *g3  = tmp3 ? tmp3[0] : NULL;
1756:   return(0);
1757: }

1759: /*@C
1760:   PetscDSSetDynamicJacobian - Set the pointwise dynamic Jacobian, dF/du_t, function for given test and basis fields

1762:   Not collective

1764:   Input Parameters:
1765: + ds - The PetscDS
1766: . f  - The test field number
1767: . g  - The field number
1768: . g0 - integrand for the test and basis function term
1769: . g1 - integrand for the test function and basis function gradient term
1770: . g2 - integrand for the test function gradient and basis function term
1771: - g3 - integrand for the test function gradient and basis function gradient term

1773:   Note: We are using a first order FEM model for the weak form:

1775:   \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

1777: The calling sequence for the callbacks g0, g1, g2 and g3 is given by:

1779: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1780: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1781: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1782: $    PetscReal t, const PetscReal x[], PetscScalar g0[])

1784: + dim - the spatial dimension
1785: . Nf - the number of fields
1786: . uOff - the offset into u[] and u_t[] for each field
1787: . uOff_x - the offset into u_x[] for each field
1788: . u - each field evaluated at the current point
1789: . u_t - the time derivative of each field evaluated at the current point
1790: . u_x - the gradient of each field evaluated at the current point
1791: . aOff - the offset into a[] and a_t[] for each auxiliary field
1792: . aOff_x - the offset into a_x[] for each auxiliary field
1793: . a - each auxiliary field evaluated at the current point
1794: . a_t - the time derivative of each auxiliary field evaluated at the current point
1795: . a_x - the gradient of auxiliary each field evaluated at the current point
1796: . t - current time
1797: . u_tShift - the multiplier a for dF/dU_t
1798: . x - coordinates of the current point
1799: . numConstants - number of constant parameters
1800: . constants - constant parameters
1801: - g0 - output values at the current point

1803:   Level: intermediate

1805: .seealso: PetscDSGetJacobian()
1806: @*/
1807: PetscErrorCode PetscDSSetDynamicJacobian(PetscDS ds, PetscInt f, PetscInt g,
1808:                                          void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1809:                                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1810:                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1811:                                                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1812:                                          void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1813:                                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1814:                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1815:                                                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1816:                                          void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1817:                                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1818:                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1819:                                                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1820:                                          void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1821:                                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1822:                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1823:                                                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1824: {

1833:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1834:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
1835:   PetscWeakFormSetIndexDynamicJacobian(ds->wf, NULL, 0, f, g, 0, g0, 0, g1, 0, g2, 0, g3);
1836:   return(0);
1837: }

1839: /*@C
1840:   PetscDSGetRiemannSolver - Returns the Riemann solver for the given field

1842:   Not collective

1844:   Input Arguments:
1845: + ds - The PetscDS object
1846: - f  - The field number

1848:   Output Argument:
1849: . r    - Riemann solver

1851:   Calling sequence for r:

1853: $ r(PetscInt dim, PetscInt Nf, const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscScalar flux[], void *ctx)

1855: + dim  - The spatial dimension
1856: . Nf   - The number of fields
1857: . x    - The coordinates at a point on the interface
1858: . n    - The normal vector to the interface
1859: . uL   - The state vector to the left of the interface
1860: . uR   - The state vector to the right of the interface
1861: . flux - output array of flux through the interface
1862: . numConstants - number of constant parameters
1863: . constants - constant parameters
1864: - ctx  - optional user context

1866:   Level: intermediate

1868: .seealso: PetscDSSetRiemannSolver()
1869: @*/
1870: PetscErrorCode PetscDSGetRiemannSolver(PetscDS ds, PetscInt f,
1871:                                        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))
1872: {
1873:   PetscRiemannFunc *tmp;
1874:   PetscInt          n;
1875:   PetscErrorCode    ierr;

1880:   if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
1881:   PetscWeakFormGetRiemannSolver(ds->wf, NULL, 0, f, &n, &tmp);
1882:   *r   = tmp ? tmp[0] : NULL;
1883:   return(0);
1884: }

1886: /*@C
1887:   PetscDSSetRiemannSolver - Sets the Riemann solver for the given field

1889:   Not collective

1891:   Input Arguments:
1892: + ds - The PetscDS object
1893: . f  - The field number
1894: - r  - Riemann solver

1896:   Calling sequence for r:

1898: $ r(PetscInt dim, PetscInt Nf, const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscScalar flux[], void *ctx)

1900: + dim  - The spatial dimension
1901: . Nf   - The number of fields
1902: . x    - The coordinates at a point on the interface
1903: . n    - The normal vector to the interface
1904: . uL   - The state vector to the left of the interface
1905: . uR   - The state vector to the right of the interface
1906: . flux - output array of flux through the interface
1907: . numConstants - number of constant parameters
1908: . constants - constant parameters
1909: - ctx  - optional user context

1911:   Level: intermediate

1913: .seealso: PetscDSGetRiemannSolver()
1914: @*/
1915: PetscErrorCode PetscDSSetRiemannSolver(PetscDS ds, PetscInt f,
1916:                                        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))
1917: {

1923:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1924:   PetscWeakFormSetIndexRiemannSolver(ds->wf, NULL, 0, f, 0, r);
1925:   return(0);
1926: }

1928: /*@C
1929:   PetscDSGetUpdate - Get the pointwise update function for a given field

1931:   Not collective

1933:   Input Parameters:
1934: + ds - The PetscDS
1935: - f  - The field number

1937:   Output Parameters:
1938: . update - update function

1940:   Note: The calling sequence for the callback update is given by:

1942: $ update(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1943: $        const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1944: $        const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1945: $        PetscReal t, const PetscReal x[], PetscScalar uNew[])

1947: + dim - the spatial dimension
1948: . Nf - the number of fields
1949: . uOff - the offset into u[] and u_t[] for each field
1950: . uOff_x - the offset into u_x[] for each field
1951: . u - each field evaluated at the current point
1952: . u_t - the time derivative of each field evaluated at the current point
1953: . u_x - the gradient of each field evaluated at the current point
1954: . aOff - the offset into a[] and a_t[] for each auxiliary field
1955: . aOff_x - the offset into a_x[] for each auxiliary field
1956: . a - each auxiliary field evaluated at the current point
1957: . a_t - the time derivative of each auxiliary field evaluated at the current point
1958: . a_x - the gradient of auxiliary each field evaluated at the current point
1959: . t - current time
1960: . x - coordinates of the current point
1961: - uNew - new value for field at the current point

1963:   Level: intermediate

1965: .seealso: PetscDSSetUpdate(), PetscDSSetResidual()
1966: @*/
1967: PetscErrorCode PetscDSGetUpdate(PetscDS ds, PetscInt f,
1968:                                   void (**update)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1969:                                                   const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1970:                                                   const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1971:                                                   PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar uNew[]))
1972: {
1975:   if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
1977:   return(0);
1978: }

1980: /*@C
1981:   PetscDSSetUpdate - Set the pointwise update function for a given field

1983:   Not collective

1985:   Input Parameters:
1986: + ds     - The PetscDS
1987: . f      - The field number
1988: - update - update function

1990:   Note: The calling sequence for the callback update is given by:

1992: $ update(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1993: $        const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1994: $        const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1995: $        PetscReal t, const PetscReal x[], PetscScalar uNew[])

1997: + dim - the spatial dimension
1998: . Nf - the number of fields
1999: . uOff - the offset into u[] and u_t[] for each field
2000: . uOff_x - the offset into u_x[] for each field
2001: . u - each field evaluated at the current point
2002: . u_t - the time derivative of each field evaluated at the current point
2003: . u_x - the gradient of each field evaluated at the current point
2004: . aOff - the offset into a[] and a_t[] for each auxiliary field
2005: . aOff_x - the offset into a_x[] for each auxiliary field
2006: . a - each auxiliary field evaluated at the current point
2007: . a_t - the time derivative of each auxiliary field evaluated at the current point
2008: . a_x - the gradient of auxiliary each field evaluated at the current point
2009: . t - current time
2010: . x - coordinates of the current point
2011: - uNew - new field values at the current point

2013:   Level: intermediate

2015: .seealso: PetscDSGetResidual()
2016: @*/
2017: PetscErrorCode PetscDSSetUpdate(PetscDS ds, PetscInt f,
2018:                                 void (*update)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2019:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2020:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2021:                                                PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar uNew[]))
2022: {

2028:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2029:   PetscDSEnlarge_Static(ds, f+1);
2030:   ds->update[f] = update;
2031:   return(0);
2032: }

2034: PetscErrorCode PetscDSGetContext(PetscDS ds, PetscInt f, void **ctx)
2035: {
2038:   if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
2040:   *ctx = ds->ctx[f];
2041:   return(0);
2042: }

2044: PetscErrorCode PetscDSSetContext(PetscDS ds, PetscInt f, void *ctx)
2045: {

2050:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2051:   PetscDSEnlarge_Static(ds, f+1);
2052:   ds->ctx[f] = ctx;
2053:   return(0);
2054: }

2056: /*@C
2057:   PetscDSGetBdResidual - Get the pointwise boundary residual function for a given test field

2059:   Not collective

2061:   Input Parameters:
2062: + ds - The PetscDS
2063: - f  - The test field number

2065:   Output Parameters:
2066: + f0 - boundary integrand for the test function term
2067: - f1 - boundary integrand for the test function gradient term

2069:   Note: We are using a first order FEM model for the weak form:

2071:   \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

2073: The calling sequence for the callbacks f0 and f1 is given by:

2075: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2076: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2077: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2078: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[])

2080: + dim - the spatial dimension
2081: . Nf - the number of fields
2082: . uOff - the offset into u[] and u_t[] for each field
2083: . uOff_x - the offset into u_x[] for each field
2084: . u - each field evaluated at the current point
2085: . u_t - the time derivative of each field evaluated at the current point
2086: . u_x - the gradient of each field evaluated at the current point
2087: . aOff - the offset into a[] and a_t[] for each auxiliary field
2088: . aOff_x - the offset into a_x[] for each auxiliary field
2089: . a - each auxiliary field evaluated at the current point
2090: . a_t - the time derivative of each auxiliary field evaluated at the current point
2091: . a_x - the gradient of auxiliary each field evaluated at the current point
2092: . t - current time
2093: . x - coordinates of the current point
2094: . n - unit normal at the current point
2095: . numConstants - number of constant parameters
2096: . constants - constant parameters
2097: - f0 - output values at the current point

2099:   Level: intermediate

2101: .seealso: PetscDSSetBdResidual()
2102: @*/
2103: PetscErrorCode PetscDSGetBdResidual(PetscDS ds, PetscInt f,
2104:                                     void (**f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2105:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2106:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2107:                                                 PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]),
2108:                                     void (**f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2109:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2110:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2111:                                                 PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
2112: {
2113:   PetscBdPointFunc *tmp0, *tmp1;
2114:   PetscInt          n0, n1;
2115:   PetscErrorCode    ierr;

2119:   if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
2120:   PetscWeakFormGetBdResidual(ds->wf, NULL, 0, f, &n0, &tmp0, &n1, &tmp1);
2121:   *f0  = tmp0 ? tmp0[0] : NULL;
2122:   *f1  = tmp1 ? tmp1[0] : NULL;
2123:   return(0);
2124: }

2126: /*@C
2127:   PetscDSSetBdResidual - Get the pointwise boundary residual function for a given test field

2129:   Not collective

2131:   Input Parameters:
2132: + ds - The PetscDS
2133: . f  - The test field number
2134: . f0 - boundary integrand for the test function term
2135: - f1 - boundary integrand for the test function gradient term

2137:   Note: We are using a first order FEM model for the weak form:

2139:   \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

2141: The calling sequence for the callbacks f0 and f1 is given by:

2143: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2144: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2145: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2146: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[])

2148: + dim - the spatial dimension
2149: . Nf - the number of fields
2150: . uOff - the offset into u[] and u_t[] for each field
2151: . uOff_x - the offset into u_x[] for each field
2152: . u - each field evaluated at the current point
2153: . u_t - the time derivative of each field evaluated at the current point
2154: . u_x - the gradient of each field evaluated at the current point
2155: . aOff - the offset into a[] and a_t[] for each auxiliary field
2156: . aOff_x - the offset into a_x[] for each auxiliary field
2157: . a - each auxiliary field evaluated at the current point
2158: . a_t - the time derivative of each auxiliary field evaluated at the current point
2159: . a_x - the gradient of auxiliary each field evaluated at the current point
2160: . t - current time
2161: . x - coordinates of the current point
2162: . n - unit normal at the current point
2163: . numConstants - number of constant parameters
2164: . constants - constant parameters
2165: - f0 - output values at the current point

2167:   Level: intermediate

2169: .seealso: PetscDSGetBdResidual()
2170: @*/
2171: PetscErrorCode PetscDSSetBdResidual(PetscDS ds, PetscInt f,
2172:                                     void (*f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2173:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2174:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2175:                                                PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]),
2176:                                     void (*f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2177:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2178:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2179:                                                PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
2180: {

2185:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2186:   PetscWeakFormSetIndexBdResidual(ds->wf, NULL, 0, f, 0, f0, 0, f1);
2187:   return(0);
2188: }

2190: /*@
2191:   PetscDSHasBdJacobian - Signals that boundary Jacobian functions have been set

2193:   Not collective

2195:   Input Parameter:
2196: . ds - The PetscDS

2198:   Output Parameter:
2199: . hasBdJac - flag that pointwise function for the boundary Jacobian has been set

2201:   Level: intermediate

2203: .seealso: PetscDSHasJacobian(), PetscDSSetBdJacobian(), PetscDSGetBdJacobian()
2204: @*/
2205: PetscErrorCode PetscDSHasBdJacobian(PetscDS ds, PetscBool *hasBdJac)
2206: {

2212:   PetscWeakFormHasBdJacobian(ds->wf, hasBdJac);
2213:   return(0);
2214: }

2216: /*@C
2217:   PetscDSGetBdJacobian - Get the pointwise boundary Jacobian function for given test and basis field

2219:   Not collective

2221:   Input Parameters:
2222: + ds - The PetscDS
2223: . f  - The test field number
2224: - g  - The field number

2226:   Output Parameters:
2227: + g0 - integrand for the test and basis function term
2228: . g1 - integrand for the test function and basis function gradient term
2229: . g2 - integrand for the test function gradient and basis function term
2230: - g3 - integrand for the test function gradient and basis function gradient term

2232:   Note: We are using a first order FEM model for the weak form:

2234:   \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

2236: The calling sequence for the callbacks g0, g1, g2 and g3 is given by:

2238: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2239: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2240: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2241: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar g0[])

2243: + dim - the spatial dimension
2244: . Nf - the number of fields
2245: . uOff - the offset into u[] and u_t[] for each field
2246: . uOff_x - the offset into u_x[] for each field
2247: . u - each field evaluated at the current point
2248: . u_t - the time derivative of each field evaluated at the current point
2249: . u_x - the gradient of each field evaluated at the current point
2250: . aOff - the offset into a[] and a_t[] for each auxiliary field
2251: . aOff_x - the offset into a_x[] for each auxiliary field
2252: . a - each auxiliary field evaluated at the current point
2253: . a_t - the time derivative of each auxiliary field evaluated at the current point
2254: . a_x - the gradient of auxiliary each field evaluated at the current point
2255: . t - current time
2256: . u_tShift - the multiplier a for dF/dU_t
2257: . x - coordinates of the current point
2258: . n - normal at the current point
2259: . numConstants - number of constant parameters
2260: . constants - constant parameters
2261: - g0 - output values at the current point

2263:   Level: intermediate

2265: .seealso: PetscDSSetBdJacobian()
2266: @*/
2267: PetscErrorCode PetscDSGetBdJacobian(PetscDS ds, PetscInt f, PetscInt g,
2268:                                     void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2269:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2270:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2271:                                                 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
2272:                                     void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2273:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2274:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2275:                                                 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
2276:                                     void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2277:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2278:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2279:                                                 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
2280:                                     void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2281:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2282:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2283:                                                 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
2284: {
2285:   PetscBdPointJac *tmp0, *tmp1, *tmp2, *tmp3;
2286:   PetscInt         n0, n1, n2, n3;
2287:   PetscErrorCode   ierr;

2291:   if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
2292:   if ((g < 0) || (g >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, ds->Nf);
2293:   PetscWeakFormGetBdJacobian(ds->wf, NULL, 0, f, g, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3);
2294:   *g0  = tmp0 ? tmp0[0] : NULL;
2295:   *g1  = tmp1 ? tmp1[0] : NULL;
2296:   *g2  = tmp2 ? tmp2[0] : NULL;
2297:   *g3  = tmp3 ? tmp3[0] : NULL;
2298:   return(0);
2299: }

2301: /*@C
2302:   PetscDSSetBdJacobian - Set the pointwise boundary Jacobian function for given test and basis field

2304:   Not collective

2306:   Input Parameters:
2307: + ds - The PetscDS
2308: . f  - The test field number
2309: . g  - The field number
2310: . g0 - integrand for the test and basis function term
2311: . g1 - integrand for the test function and basis function gradient term
2312: . g2 - integrand for the test function gradient and basis function term
2313: - g3 - integrand for the test function gradient and basis function gradient term

2315:   Note: We are using a first order FEM model for the weak form:

2317:   \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

2319: The calling sequence for the callbacks g0, g1, g2 and g3 is given by:

2321: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2322: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2323: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2324: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar g0[])

2326: + dim - the spatial dimension
2327: . Nf - the number of fields
2328: . uOff - the offset into u[] and u_t[] for each field
2329: . uOff_x - the offset into u_x[] for each field
2330: . u - each field evaluated at the current point
2331: . u_t - the time derivative of each field evaluated at the current point
2332: . u_x - the gradient of each field evaluated at the current point
2333: . aOff - the offset into a[] and a_t[] for each auxiliary field
2334: . aOff_x - the offset into a_x[] for each auxiliary field
2335: . a - each auxiliary field evaluated at the current point
2336: . a_t - the time derivative of each auxiliary field evaluated at the current point
2337: . a_x - the gradient of auxiliary each field evaluated at the current point
2338: . t - current time
2339: . u_tShift - the multiplier a for dF/dU_t
2340: . x - coordinates of the current point
2341: . n - normal at the current point
2342: . numConstants - number of constant parameters
2343: . constants - constant parameters
2344: - g0 - output values at the current point

2346:   Level: intermediate

2348: .seealso: PetscDSGetBdJacobian()
2349: @*/
2350: PetscErrorCode PetscDSSetBdJacobian(PetscDS ds, PetscInt f, PetscInt g,
2351:                                     void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2352:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2353:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2354:                                                PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
2355:                                     void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2356:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2357:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2358:                                                PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
2359:                                     void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2360:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2361:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2362:                                                PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
2363:                                     void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2364:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2365:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2366:                                                PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
2367: {

2376:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2377:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
2378:   PetscWeakFormSetIndexBdJacobian(ds->wf, NULL, 0, f, g, 0, g0, 0, g1, 0, g2, 0, g3);
2379:   return(0);
2380: }

2382: /*@
2383:   PetscDSHasBdJacobianPreconditioner - Signals that boundary Jacobian preconditioner functions have been set

2385:   Not collective

2387:   Input Parameter:
2388: . ds - The PetscDS

2390:   Output Parameter:
2391: . hasBdJac - flag that pointwise function for the boundary Jacobian preconditioner has been set

2393:   Level: intermediate

2395: .seealso: PetscDSHasJacobian(), PetscDSSetBdJacobian(), PetscDSGetBdJacobian()
2396: @*/
2397: PetscErrorCode PetscDSHasBdJacobianPreconditioner(PetscDS ds, PetscBool *hasBdJacPre)
2398: {

2404:   PetscWeakFormHasBdJacobianPreconditioner(ds->wf, hasBdJacPre);
2405:   return(0);
2406: }

2408: /*@C
2409:   PetscDSGetBdJacobianPreconditioner - Get the pointwise boundary Jacobian preconditioner function for given test and basis field

2411:   Not collective

2413:   Input Parameters:
2414: + ds - The PetscDS
2415: . f  - The test field number
2416: - g  - The field number

2418:   Output Parameters:
2419: + g0 - integrand for the test and basis function term
2420: . g1 - integrand for the test function and basis function gradient term
2421: . g2 - integrand for the test function gradient and basis function term
2422: - g3 - integrand for the test function gradient and basis function gradient term

2424:   Note: We are using a first order FEM model for the weak form:

2426:   \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

2428: The calling sequence for the callbacks g0, g1, g2 and g3 is given by:

2430: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2431: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2432: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2433: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])

2435: + dim - the spatial dimension
2436: . Nf - the number of fields
2437: . NfAux - the number of auxiliary fields
2438: . uOff - the offset into u[] and u_t[] for each field
2439: . uOff_x - the offset into u_x[] for each field
2440: . u - each field evaluated at the current point
2441: . u_t - the time derivative of each field evaluated at the current point
2442: . u_x - the gradient of each field evaluated at the current point
2443: . aOff - the offset into a[] and a_t[] for each auxiliary field
2444: . aOff_x - the offset into a_x[] for each auxiliary field
2445: . a - each auxiliary field evaluated at the current point
2446: . a_t - the time derivative of each auxiliary field evaluated at the current point
2447: . a_x - the gradient of auxiliary each field evaluated at the current point
2448: . t - current time
2449: . u_tShift - the multiplier a for dF/dU_t
2450: . x - coordinates of the current point
2451: . n - normal at the current point
2452: . numConstants - number of constant parameters
2453: . constants - constant parameters
2454: - g0 - output values at the current point

2456:   This is not yet available in Fortran.

2458:   Level: intermediate

2460: .seealso: PetscDSSetBdJacobianPreconditioner()
2461: @*/
2462: PetscErrorCode PetscDSGetBdJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g,
2463:                                                   void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2464:                                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2465:                                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2466:                                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
2467:                                                   void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2468:                                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2469:                                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2470:                                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
2471:                                                   void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2472:                                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2473:                                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2474:                                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
2475:                                                   void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2476:                                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2477:                                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2478:                                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
2479: {
2480:   PetscBdPointJac *tmp0, *tmp1, *tmp2, *tmp3;
2481:   PetscInt         n0, n1, n2, n3;
2482:   PetscErrorCode   ierr;

2486:   if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
2487:   if ((g < 0) || (g >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, ds->Nf);
2488:   PetscWeakFormGetBdJacobianPreconditioner(ds->wf, NULL, 0, f, g, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3);
2489:   *g0  = tmp0 ? tmp0[0] : NULL;
2490:   *g1  = tmp1 ? tmp1[0] : NULL;
2491:   *g2  = tmp2 ? tmp2[0] : NULL;
2492:   *g3  = tmp3 ? tmp3[0] : NULL;
2493:   return(0);
2494: }

2496: /*@C
2497:   PetscDSSetBdJacobianPreconditioner - Set the pointwise boundary Jacobian preconditioner function for given test and basis field

2499:   Not collective

2501:   Input Parameters:
2502: + ds - The PetscDS
2503: . f  - The test field number
2504: . g  - The field number
2505: . g0 - integrand for the test and basis function term
2506: . g1 - integrand for the test function and basis function gradient term
2507: . g2 - integrand for the test function gradient and basis function term
2508: - g3 - integrand for the test function gradient and basis function gradient term

2510:   Note: We are using a first order FEM model for the weak form:

2512:   \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

2514: The calling sequence for the callbacks g0, g1, g2 and g3 is given by:

2516: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2517: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2518: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2519: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])

2521: + dim - the spatial dimension
2522: . Nf - the number of fields
2523: . NfAux - the number of auxiliary fields
2524: . uOff - the offset into u[] and u_t[] for each field
2525: . uOff_x - the offset into u_x[] for each field
2526: . u - each field evaluated at the current point
2527: . u_t - the time derivative of each field evaluated at the current point
2528: . u_x - the gradient of each field evaluated at the current point
2529: . aOff - the offset into a[] and a_t[] for each auxiliary field
2530: . aOff_x - the offset into a_x[] for each auxiliary field
2531: . a - each auxiliary field evaluated at the current point
2532: . a_t - the time derivative of each auxiliary field evaluated at the current point
2533: . a_x - the gradient of auxiliary each field evaluated at the current point
2534: . t - current time
2535: . u_tShift - the multiplier a for dF/dU_t
2536: . x - coordinates of the current point
2537: . n - normal at the current point
2538: . numConstants - number of constant parameters
2539: . constants - constant parameters
2540: - g0 - output values at the current point

2542:   This is not yet available in Fortran.

2544:   Level: intermediate

2546: .seealso: PetscDSGetBdJacobianPreconditioner()
2547: @*/
2548: PetscErrorCode PetscDSSetBdJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g,
2549:                                                   void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2550:                                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2551:                                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2552:                                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
2553:                                                   void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2554:                                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2555:                                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2556:                                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
2557:                                                   void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2558:                                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2559:                                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2560:                                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
2561:                                                   void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2562:                                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2563:                                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2564:                                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
2565: {

2574:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2575:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
2576:   PetscWeakFormSetIndexBdJacobianPreconditioner(ds->wf, NULL, 0, f, g, 0, g0, 0, g1, 0, g2, 0, g3);
2577:   return(0);
2578: }

2580: /*@C
2581:   PetscDSGetExactSolution - Get the pointwise exact solution function for a given test field

2583:   Not collective

2585:   Input Parameters:
2586: + prob - The PetscDS
2587: - f    - The test field number

2589:   Output Parameter:
2590: + exactSol - exact solution for the test field
2591: - exactCtx - exact solution context

2593:   Note: The calling sequence for the solution functions is given by:

2595: $ sol(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx)

2597: + dim - the spatial dimension
2598: . t - current time
2599: . x - coordinates of the current point
2600: . Nc - the number of field components
2601: . u - the solution field evaluated at the current point
2602: - ctx - a user context

2604:   Level: intermediate

2606: .seealso: PetscDSSetExactSolution(), PetscDSGetExactSolutionTimeDerivative()
2607: @*/
2608: PetscErrorCode PetscDSGetExactSolution(PetscDS prob, PetscInt f, PetscErrorCode (**sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void **ctx)
2609: {
2612:   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);
2615:   return(0);
2616: }

2618: /*@C
2619:   PetscDSSetExactSolution - Set the pointwise exact solution function for a given test field

2621:   Not collective

2623:   Input Parameters:
2624: + prob - The PetscDS
2625: . f    - The test field number
2626: . sol  - solution function for the test fields
2627: - ctx  - solution context or NULL

2629:   Note: The calling sequence for solution functions is given by:

2631: $ sol(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx)

2633: + dim - the spatial dimension
2634: . t - current time
2635: . x - coordinates of the current point
2636: . Nc - the number of field components
2637: . u - the solution field evaluated at the current point
2638: - ctx - a user context

2640:   Level: intermediate

2642: .seealso: PetscDSGetExactSolution()
2643: @*/
2644: PetscErrorCode PetscDSSetExactSolution(PetscDS prob, PetscInt f, PetscErrorCode (*sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void *ctx)
2645: {

2650:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2651:   PetscDSEnlarge_Static(prob, f+1);
2654:   return(0);
2655: }

2657: /*@C
2658:   PetscDSGetExactSolutionTimeDerivative - Get the pointwise time derivative of the exact solution function for a given test field

2660:   Not collective

2662:   Input Parameters:
2663: + prob - The PetscDS
2664: - f    - The test field number

2666:   Output Parameter:
2667: + exactSol - time derivative of the exact solution for the test field
2668: - exactCtx - time derivative of the exact solution context

2670:   Note: The calling sequence for the solution functions is given by:

2672: $ sol(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx)

2674: + dim - the spatial dimension
2675: . t - current time
2676: . x - coordinates of the current point
2677: . Nc - the number of field components
2678: . u - the solution field evaluated at the current point
2679: - ctx - a user context

2681:   Level: intermediate

2683: .seealso: PetscDSSetExactSolutionTimeDerivative(), PetscDSGetExactSolution()
2684: @*/
2685: PetscErrorCode PetscDSGetExactSolutionTimeDerivative(PetscDS prob, PetscInt f, PetscErrorCode (**sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void **ctx)
2686: {
2689:   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);
2692:   return(0);
2693: }

2695: /*@C
2696:   PetscDSSetExactSolutionTimeDerivative - Set the pointwise time derivative of the exact solution function for a given test field

2698:   Not collective

2700:   Input Parameters:
2701: + prob - The PetscDS
2702: . f    - The test field number
2703: . sol  - time derivative of the solution function for the test fields
2704: - ctx  - time derivative of the solution context or NULL

2706:   Note: The calling sequence for solution functions is given by:

2708: $ sol(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx)

2710: + dim - the spatial dimension
2711: . t - current time
2712: . x - coordinates of the current point
2713: . Nc - the number of field components
2714: . u - the solution field evaluated at the current point
2715: - ctx - a user context

2717:   Level: intermediate

2719: .seealso: PetscDSGetExactSolutionTimeDerivative(), PetscDSSetExactSolution()
2720: @*/
2721: PetscErrorCode PetscDSSetExactSolutionTimeDerivative(PetscDS prob, PetscInt f, PetscErrorCode (*sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void *ctx)
2722: {

2727:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2728:   PetscDSEnlarge_Static(prob, f+1);
2731:   return(0);
2732: }

2734: /*@C
2735:   PetscDSGetConstants - Returns the array of constants passed to point functions

2737:   Not collective

2739:   Input Parameter:
2740: . prob - The PetscDS object

2742:   Output Parameters:
2743: + numConstants - The number of constants
2744: - constants    - The array of constants, NULL if there are none

2746:   Level: intermediate

2748: .seealso: PetscDSSetConstants(), PetscDSCreate()
2749: @*/
2750: PetscErrorCode PetscDSGetConstants(PetscDS prob, PetscInt *numConstants, const PetscScalar *constants[])
2751: {
2756:   return(0);
2757: }

2759: /*@C
2760:   PetscDSSetConstants - Set the array of constants passed to point functions

2762:   Not collective

2764:   Input Parameters:
2765: + prob         - The PetscDS object
2766: . numConstants - The number of constants
2767: - constants    - The array of constants, NULL if there are none

2769:   Level: intermediate

2771: .seealso: PetscDSGetConstants(), PetscDSCreate()
2772: @*/
2773: PetscErrorCode PetscDSSetConstants(PetscDS prob, PetscInt numConstants, PetscScalar constants[])
2774: {

2779:   if (numConstants != prob->numConstants) {
2780:     PetscFree(prob->constants);
2781:     prob->numConstants = numConstants;
2782:     if (prob->numConstants) {
2783:       PetscMalloc1(prob->numConstants, &prob->constants);
2784:     } else {
2785:       prob->constants = NULL;
2786:     }
2787:   }
2788:   if (prob->numConstants) {
2790:     PetscArraycpy(prob->constants, constants, prob->numConstants);
2791:   }
2792:   return(0);
2793: }

2795: /*@
2796:   PetscDSGetFieldIndex - Returns the index of the given field

2798:   Not collective

2800:   Input Parameters:
2801: + prob - The PetscDS object
2802: - disc - The discretization object

2804:   Output Parameter:
2805: . f - The field number

2807:   Level: beginner

2809: .seealso: PetscGetDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
2810: @*/
2811: PetscErrorCode PetscDSGetFieldIndex(PetscDS prob, PetscObject disc, PetscInt *f)
2812: {
2813:   PetscInt g;

2818:   *f = -1;
2819:   for (g = 0; g < prob->Nf; ++g) {if (disc == prob->disc[g]) break;}
2820:   if (g == prob->Nf) SETERRQ(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Field not found in PetscDS.");
2821:   *f = g;
2822:   return(0);
2823: }

2825: /*@
2826:   PetscDSGetFieldSize - Returns the size of the given field in the full space basis

2828:   Not collective

2830:   Input Parameters:
2831: + prob - The PetscDS object
2832: - f - The field number

2834:   Output Parameter:
2835: . size - The size

2837:   Level: beginner

2839: .seealso: PetscDSGetFieldOffset(), PetscDSGetNumFields(), PetscDSCreate()
2840: @*/
2841: PetscErrorCode PetscDSGetFieldSize(PetscDS prob, PetscInt f, PetscInt *size)
2842: {

2848:   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);
2849:   PetscDSSetUp(prob);
2850:   *size = prob->Nb[f];
2851:   return(0);
2852: }

2854: /*@
2855:   PetscDSGetFieldOffset - Returns the offset of the given field in the full space basis

2857:   Not collective

2859:   Input Parameters:
2860: + prob - The PetscDS object
2861: - f - The field number

2863:   Output Parameter:
2864: . off - The offset

2866:   Level: beginner

2868: .seealso: PetscDSGetFieldSize(), PetscDSGetNumFields(), PetscDSCreate()
2869: @*/
2870: PetscErrorCode PetscDSGetFieldOffset(PetscDS prob, PetscInt f, PetscInt *off)
2871: {
2872:   PetscInt       size, g;

2878:   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);
2879:   *off = 0;
2880:   for (g = 0; g < f; ++g) {
2881:     PetscDSGetFieldSize(prob, g, &size);
2882:     *off += size;
2883:   }
2884:   return(0);
2885: }

2887: /*@
2888:   PetscDSGetDimensions - Returns the size of the approximation space for each field on an evaluation point

2890:   Not collective

2892:   Input Parameter:
2893: . prob - The PetscDS object

2895:   Output Parameter:
2896: . dimensions - The number of dimensions

2898:   Level: beginner

2900: .seealso: PetscDSGetComponentOffsets(), PetscDSGetNumFields(), PetscDSCreate()
2901: @*/
2902: PetscErrorCode PetscDSGetDimensions(PetscDS prob, PetscInt *dimensions[])
2903: {

2908:   PetscDSSetUp(prob);
2910:   *dimensions = prob->Nb;
2911:   return(0);
2912: }

2914: /*@
2915:   PetscDSGetComponents - Returns the number of components for each field on an evaluation point

2917:   Not collective

2919:   Input Parameter:
2920: . prob - The PetscDS object

2922:   Output Parameter:
2923: . components - The number of components

2925:   Level: beginner

2927: .seealso: PetscDSGetComponentOffsets(), PetscDSGetNumFields(), PetscDSCreate()
2928: @*/
2929: PetscErrorCode PetscDSGetComponents(PetscDS prob, PetscInt *components[])
2930: {

2935:   PetscDSSetUp(prob);
2937:   *components = prob->Nc;
2938:   return(0);
2939: }

2941: /*@
2942:   PetscDSGetComponentOffset - Returns the offset of the given field on an evaluation point

2944:   Not collective

2946:   Input Parameters:
2947: + prob - The PetscDS object
2948: - f - The field number

2950:   Output Parameter:
2951: . off - The offset

2953:   Level: beginner

2955: .seealso: PetscDSGetNumFields(), PetscDSCreate()
2956: @*/
2957: PetscErrorCode PetscDSGetComponentOffset(PetscDS prob, PetscInt f, PetscInt *off)
2958: {
2962:   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);
2963:   *off = prob->off[f];
2964:   return(0);
2965: }

2967: /*@
2968:   PetscDSGetComponentOffsets - Returns the offset of each field on an evaluation point

2970:   Not collective

2972:   Input Parameter:
2973: . prob - The PetscDS object

2975:   Output Parameter:
2976: . offsets - The offsets

2978:   Level: beginner

2980: .seealso: PetscDSGetNumFields(), PetscDSCreate()
2981: @*/
2982: PetscErrorCode PetscDSGetComponentOffsets(PetscDS prob, PetscInt *offsets[])
2983: {
2987:   *offsets = prob->off;
2988:   return(0);
2989: }

2991: /*@
2992:   PetscDSGetComponentDerivativeOffsets - Returns the offset of each field derivative on an evaluation point

2994:   Not collective

2996:   Input Parameter:
2997: . prob - The PetscDS object

2999:   Output Parameter:
3000: . offsets - The offsets

3002:   Level: beginner

3004: .seealso: PetscDSGetNumFields(), PetscDSCreate()
3005: @*/
3006: PetscErrorCode PetscDSGetComponentDerivativeOffsets(PetscDS prob, PetscInt *offsets[])
3007: {
3011:   *offsets = prob->offDer;
3012:   return(0);
3013: }

3015: /*@C
3016:   PetscDSGetTabulation - Return the basis tabulation at quadrature points for the volume discretization

3018:   Not collective

3020:   Input Parameter:
3021: . prob - The PetscDS object

3023:   Output Parameter:
3024: . T - The basis function and derivatives tabulation at quadrature points for each field

3026:   Level: intermediate

3028: .seealso: PetscDSCreate()
3029: @*/
3030: PetscErrorCode PetscDSGetTabulation(PetscDS prob, PetscTabulation *T[])
3031: {

3037:   PetscDSSetUp(prob);
3038:   *T = prob->T;
3039:   return(0);
3040: }

3042: /*@C
3043:   PetscDSGetFaceTabulation - Return the basis tabulation at quadrature points on the faces

3045:   Not collective

3047:   Input Parameter:
3048: . prob - The PetscDS object

3050:   Output Parameter:
3051: . Tf - The basis function and derviative tabulation on each local face at quadrature points for each and field

3053:   Level: intermediate

3055: .seealso: PetscDSGetTabulation(), PetscDSCreate()
3056: @*/
3057: PetscErrorCode PetscDSGetFaceTabulation(PetscDS prob, PetscTabulation *Tf[])
3058: {

3064:   PetscDSSetUp(prob);
3065:   *Tf = prob->Tf;
3066:   return(0);
3067: }

3069: PetscErrorCode PetscDSGetEvaluationArrays(PetscDS prob, PetscScalar **u, PetscScalar **u_t, PetscScalar **u_x)
3070: {

3075:   PetscDSSetUp(prob);
3079:   return(0);
3080: }

3082: PetscErrorCode PetscDSGetWeakFormArrays(PetscDS prob, PetscScalar **f0, PetscScalar **f1, PetscScalar **g0, PetscScalar **g1, PetscScalar **g2, PetscScalar **g3)
3083: {

3088:   PetscDSSetUp(prob);
3095:   return(0);
3096: }

3098: PetscErrorCode PetscDSGetWorkspace(PetscDS prob, PetscReal **x, PetscScalar **basisReal, PetscScalar **basisDerReal, PetscScalar **testReal, PetscScalar **testDerReal)
3099: {

3104:   PetscDSSetUp(prob);
3110:   return(0);
3111: }

3113: /*@C
3114:   PetscDSAddBoundary - Add a boundary condition to the model. The pointwise functions are used to provide boundary values for essential boundary conditions. In FEM, they are acting upon by dual basis functionals to generate FEM coefficients which are fixed. Natural boundary conditions signal to PETSc that boundary integrals should be performaed, using the kernels from PetscDSSetBdResidual().

3116:   Collective on ds

3118:   Input Parameters:
3119: + ds       - The PetscDS object
3120: . type     - The type of condition, e.g. DM_BC_ESSENTIAL/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
3121: . name     - The BC name
3122: . label    - The label defining constrained points
3123: . Nv       - The number of DMLabel values for constrained points
3124: . values   - An array of label values for constrained points
3125: . field    - The field to constrain
3126: . Nc       - The number of constrained field components (0 will constrain all fields)
3127: . comps    - An array of constrained component numbers
3128: . bcFunc   - A pointwise function giving boundary values
3129: . bcFunc_t - A pointwise function giving the time derviative of the boundary values, or NULL
3130: - ctx      - An optional user context for bcFunc

3132:   Output Parameters:
3133: - bd       - The boundary number

3135:   Options Database Keys:
3136: + -bc_<boundary name> <num> - Overrides the boundary ids
3137: - -bc_<boundary name>_comp <num> - Overrides the boundary components

3139:   Note:
3140:   Both bcFunc abd bcFunc_t will depend on the boundary condition type. If the type if DM_BC_ESSENTIAL, Then the calling sequence is:

3142: $ bcFunc(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar bcval[])

3144:   If the type is DM_BC_ESSENTIAL_FIELD or other _FIELD value, then the calling sequence is:

3146: $ bcFunc(PetscInt dim, PetscInt Nf, PetscInt NfAux,
3147: $        const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
3148: $        const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
3149: $        PetscReal time, const PetscReal x[], PetscScalar bcval[])

3151: + dim - the spatial dimension
3152: . Nf - the number of fields
3153: . uOff - the offset into u[] and u_t[] for each field
3154: . uOff_x - the offset into u_x[] for each field
3155: . u - each field evaluated at the current point
3156: . u_t - the time derivative of each field evaluated at the current point
3157: . u_x - the gradient of each field evaluated at the current point
3158: . aOff - the offset into a[] and a_t[] for each auxiliary field
3159: . aOff_x - the offset into a_x[] for each auxiliary field
3160: . a - each auxiliary field evaluated at the current point
3161: . a_t - the time derivative of each auxiliary field evaluated at the current point
3162: . a_x - the gradient of auxiliary each field evaluated at the current point
3163: . t - current time
3164: . x - coordinates of the current point
3165: . numConstants - number of constant parameters
3166: . constants - constant parameters
3167: - bcval - output values at the current point

3169:   Level: developer

3171: .seealso: PetscDSAddBoundaryByName(), PetscDSGetBoundary(), PetscDSSetResidual(), PetscDSSetBdResidual()
3172: @*/
3173: PetscErrorCode PetscDSAddBoundary(PetscDS ds, DMBoundaryConditionType type, const char name[], DMLabel label, PetscInt Nv, const PetscInt values[], PetscInt field, PetscInt Nc, const PetscInt comps[], void (*bcFunc)(void), void (*bcFunc_t)(void), void *ctx, PetscInt *bd)
3174: {
3175:   DSBoundary     head = ds->boundary, b;
3176:   PetscInt       n    = 0;
3177:   const char    *lname;

3188:   PetscNew(&b);
3189:   PetscStrallocpy(name, (char **) &b->name);
3190:   PetscWeakFormCreate(PETSC_COMM_SELF, &b->wf);
3191:   PetscWeakFormSetNumFields(b->wf, ds->Nf);
3192:   PetscMalloc1(Nv, &b->values);
3193:   if (Nv) {PetscArraycpy(b->values, values, Nv);}
3194:   PetscMalloc1(Nc, &b->comps);
3195:   if (Nc) {PetscArraycpy(b->comps, comps, Nc);}
3196:   PetscObjectGetName((PetscObject) label, &lname);
3197:   PetscStrallocpy(lname, (char **) &b->lname);
3198:   b->type   = type;
3199:   b->label  = label;
3200:   b->Nv     = Nv;
3201:   b->field  = field;
3202:   b->Nc     = Nc;
3203:   b->func   = bcFunc;
3204:   b->func_t = bcFunc_t;
3205:   b->ctx    = ctx;
3206:   b->next   = NULL;
3207:   /* Append to linked list so that we can preserve the order */
3208:   if (!head) ds->boundary = b;
3209:   while (head) {
3210:     if (!head->next) {
3211:       head->next = b;
3212:       head       = b;
3213:     }
3214:     head = head->next;
3215:     ++n;
3216:   }
3218:   return(0);
3219: }

3221: /*@C
3222:   PetscDSAddBoundaryByName - Add a boundary condition to the model. The pointwise functions are used to provide boundary values for essential boundary conditions. In FEM, they are acting upon by dual basis functionals to generate FEM coefficients which are fixed. Natural boundary conditions signal to PETSc that boundary integrals should be performaed, using the kernels from PetscDSSetBdResidual().

3224:   Collective on ds

3226:   Input Parameters:
3227: + ds       - The PetscDS object
3228: . type     - The type of condition, e.g. DM_BC_ESSENTIAL/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
3229: . name     - The BC name
3230: . lname    - The naem of the label defining constrained points
3231: . Nv       - The number of DMLabel values for constrained points
3232: . values   - An array of label values for constrained points
3233: . field    - The field to constrain
3234: . Nc       - The number of constrained field components (0 will constrain all fields)
3235: . comps    - An array of constrained component numbers
3236: . bcFunc   - A pointwise function giving boundary values
3237: . bcFunc_t - A pointwise function giving the time derviative of the boundary values, or NULL
3238: - ctx      - An optional user context for bcFunc

3240:   Output Parameters:
3241: - bd       - The boundary number

3243:   Options Database Keys:
3244: + -bc_<boundary name> <num> - Overrides the boundary ids
3245: - -bc_<boundary name>_comp <num> - Overrides the boundary components

3247:   Note:
3248:   This function should only be used with DMForest currently, since labels cannot be defined before the underlygin Plex is built.

3250:   Both bcFunc abd bcFunc_t will depend on the boundary condition type. If the type if DM_BC_ESSENTIAL, Then the calling sequence is:

3252: $ bcFunc(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar bcval[])

3254:   If the type is DM_BC_ESSENTIAL_FIELD or other _FIELD value, then the calling sequence is:

3256: $ bcFunc(PetscInt dim, PetscInt Nf, PetscInt NfAux,
3257: $        const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
3258: $        const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
3259: $        PetscReal time, const PetscReal x[], PetscScalar bcval[])

3261: + dim - the spatial dimension
3262: . Nf - the number of fields
3263: . uOff - the offset into u[] and u_t[] for each field
3264: . uOff_x - the offset into u_x[] for each field
3265: . u - each field evaluated at the current point
3266: . u_t - the time derivative of each field evaluated at the current point
3267: . u_x - the gradient of each field evaluated at the current point
3268: . aOff - the offset into a[] and a_t[] for each auxiliary field
3269: . aOff_x - the offset into a_x[] for each auxiliary field
3270: . a - each auxiliary field evaluated at the current point
3271: . a_t - the time derivative of each auxiliary field evaluated at the current point
3272: . a_x - the gradient of auxiliary each field evaluated at the current point
3273: . t - current time
3274: . x - coordinates of the current point
3275: . numConstants - number of constant parameters
3276: . constants - constant parameters
3277: - bcval - output values at the current point

3279:   Level: developer

3281: .seealso: PetscDSAddBoundary(), PetscDSGetBoundary(), PetscDSSetResidual(), PetscDSSetBdResidual()
3282: @*/
3283: PetscErrorCode PetscDSAddBoundaryByName(PetscDS ds, DMBoundaryConditionType type, const char name[], const char lname[], PetscInt Nv, const PetscInt values[], PetscInt field, PetscInt Nc, const PetscInt comps[], void (*bcFunc)(void), void (*bcFunc_t)(void), void *ctx, PetscInt *bd)
3284: {
3285:   DSBoundary     head = ds->boundary, b;
3286:   PetscInt       n    = 0;

3297:   PetscNew(&b);
3298:   PetscStrallocpy(name, (char **) &b->name);
3299:   PetscWeakFormCreate(PETSC_COMM_SELF, &b->wf);
3300:   PetscWeakFormSetNumFields(b->wf, ds->Nf);
3301:   PetscMalloc1(Nv, &b->values);
3302:   if (Nv) {PetscArraycpy(b->values, values, Nv);}
3303:   PetscMalloc1(Nc, &b->comps);
3304:   if (Nc) {PetscArraycpy(b->comps, comps, Nc);}
3305:   PetscStrallocpy(lname, (char **) &b->lname);
3306:   b->type   = type;
3307:   b->label  = NULL;
3308:   b->Nv     = Nv;
3309:   b->field  = field;
3310:   b->Nc     = Nc;
3311:   b->func   = bcFunc;
3312:   b->func_t = bcFunc_t;
3313:   b->ctx    = ctx;
3314:   b->next   = NULL;
3315:   /* Append to linked list so that we can preserve the order */
3316:   if (!head) ds->boundary = b;
3317:   while (head) {
3318:     if (!head->next) {
3319:       head->next = b;
3320:       head       = b;
3321:     }
3322:     head = head->next;
3323:     ++n;
3324:   }
3326:   return(0);
3327: }

3329: /*@C
3330:   PetscDSUpdateBoundary - Change a boundary condition for the model. The pointwise functions are used to provide boundary values for essential boundary conditions. In FEM, they are acting upon by dual basis functionals to generate FEM coefficients which are fixed. Natural boundary conditions signal to PETSc that boundary integrals should be performaed, using the kernels from PetscDSSetBdResidual().

3332:   Input Parameters:
3333: + ds       - The PetscDS object
3334: . bd       - The boundary condition number
3335: . type     - The type of condition, e.g. DM_BC_ESSENTIAL/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
3336: . name     - The BC name
3337: . label    - The label defining constrained points
3338: . Nv       - The number of DMLabel ids for constrained points
3339: . values   - An array of ids for constrained points
3340: . field    - The field to constrain
3341: . Nc       - The number of constrained field components
3342: . comps    - An array of constrained component numbers
3343: . bcFunc   - A pointwise function giving boundary values
3344: . bcFunc_t - A pointwise function giving the time derviative of the boundary values, or NULL
3345: - ctx      - An optional user context for bcFunc

3347:   Note:
3348:   The boundary condition number is the order in which it was registered. The user can get the number of boundary conditions from PetscDSGetNumBoundary(). See PetscDSAddBoundary() for a description of the calling sequences for the callbacks.

3350:   Level: developer

3352: .seealso: PetscDSAddBoundary(), PetscDSGetBoundary(), PetscDSGetNumBoundary()
3353: @*/
3354: PetscErrorCode PetscDSUpdateBoundary(PetscDS ds, PetscInt bd, DMBoundaryConditionType type, const char name[], DMLabel label, PetscInt Nv, const PetscInt values[], PetscInt field, PetscInt Nc, const PetscInt comps[], void (*bcFunc)(void), void (*bcFunc_t)(void), void *ctx)
3355: {
3356:   DSBoundary     b = ds->boundary;
3357:   PetscInt       n = 0;

3362:   while (b) {
3363:     if (n == bd) break;
3364:     b = b->next;
3365:     ++n;
3366:   }
3367:   if (!b) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %d is not in [0, %d)", bd, n);
3368:   if (name) {
3369:     PetscFree(b->name);
3370:     PetscStrallocpy(name, (char **) &b->name);
3371:   }
3372:   b->type = type;
3373:   if (label) {
3374:     const char *name;

3376:     b->label = label;
3377:     PetscFree(b->lname);
3378:     PetscObjectGetName((PetscObject) label, &name);
3379:     PetscStrallocpy(name, (char **) &b->lname);
3380:   }
3381:   if (Nv >= 0) {
3382:     b->Nv = Nv;
3383:     PetscFree(b->values);
3384:     PetscMalloc1(Nv, &b->values);
3385:     if (Nv) {PetscArraycpy(b->values, values, Nv);}
3386:   }
3387:   if (field >= 0) b->field = field;
3388:   if (Nc >= 0) {
3389:     b->Nc = Nc;
3390:     PetscFree(b->comps);
3391:     PetscMalloc1(Nc, &b->comps);
3392:     if (Nc) {PetscArraycpy(b->comps, comps, Nc);}
3393:   }
3394:   if (bcFunc)   b->func   = bcFunc;
3395:   if (bcFunc_t) b->func_t = bcFunc_t;
3396:   if (ctx)      b->ctx    = ctx;
3397:   return(0);
3398: }

3400: /*@
3401:   PetscDSGetNumBoundary - Get the number of registered BC

3403:   Input Parameters:
3404: . ds - The PetscDS object

3406:   Output Parameters:
3407: . numBd - The number of BC

3409:   Level: intermediate

3411: .seealso: PetscDSAddBoundary(), PetscDSGetBoundary()
3412: @*/
3413: PetscErrorCode PetscDSGetNumBoundary(PetscDS ds, PetscInt *numBd)
3414: {
3415:   DSBoundary b = ds->boundary;

3420:   *numBd = 0;
3421:   while (b) {++(*numBd); b = b->next;}
3422:   return(0);
3423: }

3425: /*@C
3426:   PetscDSGetBoundary - Gets a boundary condition to the model

3428:   Input Parameters:
3429: + ds          - The PetscDS object
3430: - bd          - The BC number

3432:   Output Parameters:
3433: + wf       - The PetscWeakForm holding the pointwise functions
3434: . type     - The type of condition, e.g. DM_BC_ESSENTIAL/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
3435: . name     - The BC name
3436: . label    - The label defining constrained points
3437: . Nv       - The number of DMLabel ids for constrained points
3438: . values   - An array of ids for constrained points
3439: . field    - The field to constrain
3440: . Nc       - The number of constrained field components
3441: . comps    - An array of constrained component numbers
3442: . bcFunc   - A pointwise function giving boundary values
3443: . bcFunc_t - A pointwise function giving the time derviative of the boundary values
3444: - ctx      - An optional user context for bcFunc

3446:   Options Database Keys:
3447: + -bc_<boundary name> <num> - Overrides the boundary ids
3448: - -bc_<boundary name>_comp <num> - Overrides the boundary components

3450:   Level: developer

3452: .seealso: PetscDSAddBoundary()
3453: @*/
3454: PetscErrorCode PetscDSGetBoundary(PetscDS ds, PetscInt bd, PetscWeakForm *wf, DMBoundaryConditionType *type, const char *name[], DMLabel *label, PetscInt *Nv, const PetscInt *values[], PetscInt *field, PetscInt *Nc, const PetscInt *comps[], void (**func)(void), void (**func_t)(void), void **ctx)
3455: {
3456:   DSBoundary b = ds->boundary;
3457:   PetscInt   n = 0;

3461:   while (b) {
3462:     if (n == bd) break;
3463:     b = b->next;
3464:     ++n;
3465:   }
3466:   if (!b) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %d is not in [0, %d)", bd, n);
3467:   if (wf) {
3469:     *wf = b->wf;
3470:   }
3471:   if (type) {
3473:     *type = b->type;
3474:   }
3475:   if (name) {
3477:     *name = b->name;
3478:   }
3479:   if (label) {
3481:     *label = b->label;
3482:   }
3483:   if (Nv) {
3485:     *Nv = b->Nv;
3486:   }
3487:   if (values) {
3489:     *values = b->values;
3490:   }
3491:   if (field) {
3493:     *field = b->field;
3494:   }
3495:   if (Nc) {
3497:     *Nc = b->Nc;
3498:   }
3499:   if (comps) {
3501:     *comps = b->comps;
3502:   }
3503:   if (func) {
3505:     *func = b->func;
3506:   }
3507:   if (func_t) {
3509:     *func_t = b->func_t;
3510:   }
3511:   if (ctx) {
3513:     *ctx = b->ctx;
3514:   }
3515:   return(0);
3516: }

3518: static PetscErrorCode DSBoundaryDuplicate_Internal(DSBoundary b, DSBoundary *bNew)
3519: {

3523:   PetscNew(bNew);
3524:   PetscWeakFormCreate(PETSC_COMM_SELF, &(*bNew)->wf);
3525:   PetscWeakFormCopy(b->wf, (*bNew)->wf);
3526:   PetscStrallocpy(b->name,(char **) &((*bNew)->name));
3527:   PetscStrallocpy(b->lname,(char **) &((*bNew)->lname));
3528:   (*bNew)->type   = b->type;
3529:   (*bNew)->label  = b->label;
3530:   (*bNew)->Nv     = b->Nv;
3531:   PetscMalloc1(b->Nv, &(*bNew)->values);
3532:   PetscArraycpy((*bNew)->values, b->values, b->Nv);
3533:   (*bNew)->field  = b->field;
3534:   (*bNew)->Nc     = b->Nc;
3535:   PetscMalloc1(b->Nc, &(*bNew)->comps);
3536:   PetscArraycpy((*bNew)->comps, b->comps, b->Nc);
3537:   (*bNew)->func   = b->func;
3538:   (*bNew)->func_t = b->func_t;
3539:   (*bNew)->ctx    = b->ctx;
3540:   return(0);
3541: }

3543: /*@
3544:   PetscDSCopyBoundary - Copy all boundary condition objects to the new problem

3546:   Not collective

3548:   Input Parameters:
3549: + ds        - The source PetscDS object
3550: . numFields - The number of selected fields, or PETSC_DEFAULT for all fields
3551: - fields    - The selected fields, or NULL for all fields

3553:   Output Parameter:
3554: . newds     - The target PetscDS, now with a copy of the boundary conditions

3556:   Level: intermediate

3558: .seealso: PetscDSCopyEquations(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3559: @*/
3560: PetscErrorCode PetscDSCopyBoundary(PetscDS ds, PetscInt numFields, const PetscInt fields[], PetscDS newds)
3561: {
3562:   DSBoundary     b, *lastnext;

3568:   if (ds == newds) return(0);
3569:   PetscDSDestroyBoundary(newds);
3570:   lastnext = &(newds->boundary);
3571:   for (b = ds->boundary; b; b = b->next) {
3572:     DSBoundary bNew;
3573:     PetscInt   fieldNew = -1;

3575:     if (numFields > 0 && fields) {
3576:       PetscInt f;

3578:       for (f = 0; f < numFields; ++f) if (b->field == fields[f]) break;
3579:       if (f == numFields) continue;
3580:       fieldNew = f;
3581:     }
3582:     DSBoundaryDuplicate_Internal(b, &bNew);
3583:     bNew->field = fieldNew < 0 ? b->field : fieldNew;
3584:     *lastnext = bNew;
3585:     lastnext  = &(bNew->next);
3586:   }
3587:   return(0);
3588: }

3590: /*@
3591:   PetscDSDestroyBoundary - Remove all DMBoundary objects from the PetscDS

3593:   Not collective

3595:   Input Parameter:
3596: . ds - The PetscDS object

3598:   Level: intermediate

3600: .seealso: PetscDSCopyBoundary(), PetscDSCopyEquations()
3601: @*/
3602: PetscErrorCode PetscDSDestroyBoundary(PetscDS ds)
3603: {
3604:   DSBoundary     next = ds->boundary;

3608:   while (next) {
3609:     DSBoundary b = next;

3611:     next = b->next;
3612:     PetscWeakFormDestroy(&b->wf);
3613:     PetscFree(b->name);
3614:     PetscFree(b->lname);
3615:     PetscFree(b->values);
3616:     PetscFree(b->comps);
3617:     PetscFree(b);
3618:   }
3619:   return(0);
3620: }

3622: /*@
3623:   PetscDSSelectDiscretizations - Copy discretizations to the new problem with different field layout

3625:   Not collective

3627:   Input Parameter:
3628: + prob - The PetscDS object
3629: . numFields - Number of new fields
3630: - fields - Old field number for each new field

3632:   Output Parameter:
3633: . newprob - The PetscDS copy

3635:   Level: intermediate

3637: .seealso: PetscDSSelectEquations(), PetscDSCopyBoundary(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3638: @*/
3639: PetscErrorCode PetscDSSelectDiscretizations(PetscDS prob, PetscInt numFields, const PetscInt fields[], PetscDS newprob)
3640: {
3641:   PetscInt       Nf, Nfn, fn;

3648:   PetscDSGetNumFields(prob, &Nf);
3649:   PetscDSGetNumFields(newprob, &Nfn);
3650:   numFields = numFields < 0 ? Nf : numFields;
3651:   for (fn = 0; fn < numFields; ++fn) {
3652:     const PetscInt f = fields ? fields[fn] : fn;
3653:     PetscObject    disc;

3655:     if (f >= Nf) continue;
3656:     PetscDSGetDiscretization(prob, f, &disc);
3657:     PetscDSSetDiscretization(newprob, fn, disc);
3658:   }
3659:   return(0);
3660: }

3662: /*@
3663:   PetscDSSelectEquations - Copy pointwise function pointers to the new problem with different field layout

3665:   Not collective

3667:   Input Parameter:
3668: + prob - The PetscDS object
3669: . numFields - Number of new fields
3670: - fields - Old field number for each new field

3672:   Output Parameter:
3673: . newprob - The PetscDS copy

3675:   Level: intermediate

3677: .seealso: PetscDSSelectDiscretizations(), PetscDSCopyBoundary(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3678: @*/
3679: PetscErrorCode PetscDSSelectEquations(PetscDS prob, PetscInt numFields, const PetscInt fields[], PetscDS newprob)
3680: {
3681:   PetscInt       Nf, Nfn, fn, gn;

3688:   PetscDSGetNumFields(prob, &Nf);
3689:   PetscDSGetNumFields(newprob, &Nfn);
3690:   if (numFields > Nfn) SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_SIZ, "Number of fields %D to transfer must not be greater then the total number of fields %D", numFields, Nfn);
3691:   for (fn = 0; fn < numFields; ++fn) {
3692:     const PetscInt   f = fields ? fields[fn] : fn;
3693:     PetscPointFunc   obj;
3694:     PetscPointFunc   f0, f1;
3695:     PetscBdPointFunc f0Bd, f1Bd;
3696:     PetscRiemannFunc r;

3698:     if (f >= Nf) continue;
3699:     PetscDSGetObjective(prob, f, &obj);
3700:     PetscDSGetResidual(prob, f, &f0, &f1);
3701:     PetscDSGetBdResidual(prob, f, &f0Bd, &f1Bd);
3702:     PetscDSGetRiemannSolver(prob, f, &r);
3703:     PetscDSSetObjective(newprob, fn, obj);
3704:     PetscDSSetResidual(newprob, fn, f0, f1);
3705:     PetscDSSetBdResidual(newprob, fn, f0Bd, f1Bd);
3706:     PetscDSSetRiemannSolver(newprob, fn, r);
3707:     for (gn = 0; gn < numFields; ++gn) {
3708:       const PetscInt  g = fields ? fields[gn] : gn;
3709:       PetscPointJac   g0, g1, g2, g3;
3710:       PetscPointJac   g0p, g1p, g2p, g3p;
3711:       PetscBdPointJac g0Bd, g1Bd, g2Bd, g3Bd;

3713:       if (g >= Nf) continue;
3714:       PetscDSGetJacobian(prob, f, g, &g0, &g1, &g2, &g3);
3715:       PetscDSGetJacobianPreconditioner(prob, f, g, &g0p, &g1p, &g2p, &g3p);
3716:       PetscDSGetBdJacobian(prob, f, g, &g0Bd, &g1Bd, &g2Bd, &g3Bd);
3717:       PetscDSSetJacobian(newprob, fn, gn, g0, g1, g2, g3);
3718:       PetscDSSetJacobianPreconditioner(newprob, fn, gn, g0p, g1p, g2p, g3p);
3719:       PetscDSSetBdJacobian(newprob, fn, gn, g0Bd, g1Bd, g2Bd, g3Bd);
3720:     }
3721:   }
3722:   return(0);
3723: }

3725: /*@
3726:   PetscDSCopyEquations - Copy all pointwise function pointers to the new problem

3728:   Not collective

3730:   Input Parameter:
3731: . prob - The PetscDS object

3733:   Output Parameter:
3734: . newprob - The PetscDS copy

3736:   Level: intermediate

3738: .seealso: PetscDSCopyBoundary(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3739: @*/
3740: PetscErrorCode PetscDSCopyEquations(PetscDS prob, PetscDS newprob)
3741: {
3742:   PetscInt       Nf, Ng;

3748:   PetscDSGetNumFields(prob, &Nf);
3749:   PetscDSGetNumFields(newprob, &Ng);
3750:   if (Nf != Ng) SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_SIZ, "Number of fields must match %D != %D", Nf, Ng);
3751:   PetscDSSelectEquations(prob, Nf, NULL, newprob);
3752:   return(0);
3753: }

3755: /*@
3756:   PetscDSCopyConstants - Copy all constants to the new problem

3758:   Not collective

3760:   Input Parameter:
3761: . prob - The PetscDS object

3763:   Output Parameter:
3764: . newprob - The PetscDS copy

3766:   Level: intermediate

3768: .seealso: PetscDSCopyBoundary(), PetscDSCopyEquations(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3769: @*/
3770: PetscErrorCode PetscDSCopyConstants(PetscDS prob, PetscDS newprob)
3771: {
3772:   PetscInt           Nc;
3773:   const PetscScalar *constants;
3774:   PetscErrorCode     ierr;

3779:   PetscDSGetConstants(prob, &Nc, &constants);
3780:   PetscDSSetConstants(newprob, Nc, (PetscScalar *) constants);
3781:   return(0);
3782: }

3784: /*@
3785:   PetscDSCopyExactSolutions - Copy all exact solutions to the new problem

3787:   Not collective

3789:   Input Parameter:
3790: . ds - The PetscDS object

3792:   Output Parameter:
3793: . newds - The PetscDS copy

3795:   Level: intermediate

3797: .seealso: PetscDSCopyBoundary(), PetscDSCopyEquations(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3798: @*/
3799: PetscErrorCode PetscDSCopyExactSolutions(PetscDS ds, PetscDS newds)
3800: {
3801:   PetscSimplePointFunc sol;
3802:   void                *ctx;
3803:   PetscInt             Nf, f;
3804:   PetscErrorCode       ierr;

3809:   PetscDSGetNumFields(ds, &Nf);
3810:   for (f = 0; f < Nf; ++f) {
3811:     PetscDSGetExactSolution(ds,    f, &sol, &ctx);
3812:     PetscDSSetExactSolution(newds, f,  sol,  ctx);
3813:     PetscDSGetExactSolutionTimeDerivative(ds,    f, &sol, &ctx);
3814:     PetscDSSetExactSolutionTimeDerivative(newds, f,  sol,  ctx);
3815:   }
3816:   return(0);
3817: }

3819: PetscErrorCode PetscDSGetHeightSubspace(PetscDS prob, PetscInt height, PetscDS *subprob)
3820: {
3821:   PetscInt       dim, Nf, f;

3827:   if (height == 0) {*subprob = prob; return(0);}
3828:   PetscDSGetNumFields(prob, &Nf);
3829:   PetscDSGetSpatialDimension(prob, &dim);
3830:   if (height > dim) SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_OUTOFRANGE, "DS can only handle height in [0, %D], not %D", dim, height);
3831:   if (!prob->subprobs) {PetscCalloc1(dim, &prob->subprobs);}
3832:   if (!prob->subprobs[height-1]) {
3833:     PetscInt cdim;

3835:     PetscDSCreate(PetscObjectComm((PetscObject) prob), &prob->subprobs[height-1]);
3836:     PetscDSGetCoordinateDimension(prob, &cdim);
3837:     PetscDSSetCoordinateDimension(prob->subprobs[height-1], cdim);
3838:     for (f = 0; f < Nf; ++f) {
3839:       PetscFE      subfe;
3840:       PetscObject  obj;
3841:       PetscClassId id;

3843:       PetscDSGetDiscretization(prob, f, &obj);
3844:       PetscObjectGetClassId(obj, &id);
3845:       if (id == PETSCFE_CLASSID) {PetscFEGetHeightSubspace((PetscFE) obj, height, &subfe);}
3846:       else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unsupported discretization type for field %d", f);
3847:       PetscDSSetDiscretization(prob->subprobs[height-1], f, (PetscObject) subfe);
3848:     }
3849:   }
3850:   *subprob = prob->subprobs[height-1];
3851:   return(0);
3852: }

3854: PetscErrorCode PetscDSGetDiscType_Internal(PetscDS ds, PetscInt f, PetscDiscType *disctype)
3855: {
3856:   PetscObject    obj;
3857:   PetscClassId   id;
3858:   PetscInt       Nf;

3864:   *disctype = PETSC_DISC_NONE;
3865:   PetscDSGetNumFields(ds, &Nf);
3866:   if (f >= Nf) SETERRQ2(PetscObjectComm((PetscObject) ds), PETSC_ERR_ARG_SIZ, "Field %D must be in [0, %D)", f, Nf);
3867:   PetscDSGetDiscretization(ds, f, &obj);
3868:   if (obj) {
3869:     PetscObjectGetClassId(obj, &id);
3870:     if (id == PETSCFE_CLASSID) *disctype = PETSC_DISC_FE;
3871:     else                       *disctype = PETSC_DISC_FV;
3872:   }
3873:   return(0);
3874: }

3876: static PetscErrorCode PetscDSDestroy_Basic(PetscDS ds)
3877: {

3881:   PetscFree(ds->data);
3882:   return(0);
3883: }

3885: static PetscErrorCode PetscDSInitialize_Basic(PetscDS ds)
3886: {
3888:   ds->ops->setfromoptions = NULL;
3889:   ds->ops->setup          = NULL;
3890:   ds->ops->view           = NULL;
3891:   ds->ops->destroy        = PetscDSDestroy_Basic;
3892:   return(0);
3893: }

3895: /*MC
3896:   PETSCDSBASIC = "basic" - A discrete system with pointwise residual and boundary residual functions

3898:   Level: intermediate

3900: .seealso: PetscDSType, PetscDSCreate(), PetscDSSetType()
3901: M*/

3903: PETSC_EXTERN PetscErrorCode PetscDSCreate_Basic(PetscDS ds)
3904: {
3905:   PetscDS_Basic *b;

3910:   PetscNewLog(ds, &b);
3911:   ds->data = b;

3913:   PetscDSInitialize_Basic(ds);
3914:   return(0);
3915: }