Actual source code: dtds.c

petsc-3.11.2 2019-05-18
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: .keywords: PetscDS, register
 53: .seealso: PetscDSRegisterAll(), PetscDSRegisterDestroy()

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

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

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

 68:   Collective on PetscDS

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

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

 77:   Level: intermediate

 79:    Not available from Fortran

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

 92:   PetscObjectTypeCompare((PetscObject) prob, name, &match);
 93:   if (match) return(0);

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

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

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

111:   Not Collective

113:   Input Parameter:
114: . prob  - The PetscDS

116:   Output Parameter:
117: . name - The PetscDS type name

119:   Level: intermediate

121:    Not available from Fortran

123: .keywords: PetscDS, get, type, name
124: .seealso: PetscDSSetType(), PetscDSCreate()
125: @*/
126: PetscErrorCode PetscDSGetType(PetscDS prob, PetscDSType *name)
127: {

133:   PetscDSRegisterAll();
134:   *name = ((PetscObject) prob)->type_name;
135:   return(0);
136: }

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

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

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

189:     for (b = prob->boundary; b; b = b->next) {
190:       PetscInt c, i;

192:       if (b->field != f) continue;
193:       PetscViewerASCIIPushTab(viewer);
194:       PetscViewerASCIIPrintf(viewer, "Boundary %s (%s) %s\n", b->name, b->labelname, DMBoundaryConditionTypes[b->type]);
195:       if (!b->numcomps) {
196:         PetscViewerASCIIPrintf(viewer, "  all components\n");
197:       } else {
198:         PetscViewerASCIIPrintf(viewer, "  components: ");
199:         PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
200:         for (c = 0; c < b->numcomps; ++c) {
201:           if (c > 0) {PetscViewerASCIIPrintf(viewer, ", ");}
202:           PetscViewerASCIIPrintf(viewer, "%D", b->comps[c]);
203:         }
204:         PetscViewerASCIIPrintf(viewer, "\n");
205:         PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
206:       }
207:       PetscViewerASCIIPrintf(viewer, "  ids: ");
208:       PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
209:       for (i = 0; i < b->numids; ++i) {
210:         if (i > 0) {PetscViewerASCIIPrintf(viewer, ", ");}
211:         PetscViewerASCIIPrintf(viewer, "%D", b->ids[i]);
212:       }
213:       PetscViewerASCIIPrintf(viewer, "\n");
214:       PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
215:       PetscViewerASCIIPopTab(viewer);
216:     }
217:   }
218:   PetscDSGetConstants(prob, &numConstants, &constants);
219:   if (numConstants) {
220:     PetscViewerASCIIPrintf(viewer, "%D constants\n", numConstants);
221:     PetscViewerASCIIPushTab(viewer);
222:     for (f = 0; f < numConstants; ++f) {PetscViewerASCIIPrintf(viewer, "%g\n", (double) PetscRealPart(constants[f]));}
223:     PetscViewerASCIIPopTab(viewer);
224:   }
225:   PetscViewerASCIIPopTab(viewer);
226:   return(0);
227: }

229: /*@C
230:   PetscDSView - Views a PetscDS

232:   Collective on PetscDS

234:   Input Parameter:
235: + prob - the PetscDS object to view
236: - v  - the viewer

238:   Level: developer

240: .seealso PetscDSDestroy()
241: @*/
242: PetscErrorCode PetscDSView(PetscDS prob, PetscViewer v)
243: {
244:   PetscBool      iascii;

249:   if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) prob), &v);}
251:   PetscObjectTypeCompare((PetscObject) v, PETSCVIEWERASCII, &iascii);
252:   if (iascii) {PetscDSView_Ascii(prob, v);}
253:   if (prob->ops->view) {(*prob->ops->view)(prob, v);}
254:   return(0);
255: }

257: /*@
258:   PetscDSSetFromOptions - sets parameters in a PetscDS from the options database

260:   Collective on PetscDS

262:   Input Parameter:
263: . prob - the PetscDS object to set options for

265:   Options Database:
266: + -petscds_type <type>     : Set the DS type
267: . -petscds_view <view opt> : View the DS
268: . -petscds_jac_pre         : Turn formation of a separate Jacobian preconditioner on and off
269: . -bc_<name> <ids>         : Specify a list of label ids for a boundary condition
270: - -bc_<name>_comp <comps>  : Specify a list of field components to constrain for a boundary condition

272:   Level: developer

274: .seealso PetscDSView()
275: @*/
276: PetscErrorCode PetscDSSetFromOptions(PetscDS prob)
277: {
278:   DSBoundary     b;
279:   const char    *defaultType;
280:   char           name[256];
281:   PetscBool      flg;

286:   if (!((PetscObject) prob)->type_name) {
287:     defaultType = PETSCDSBASIC;
288:   } else {
289:     defaultType = ((PetscObject) prob)->type_name;
290:   }
291:   PetscDSRegisterAll();

293:   PetscObjectOptionsBegin((PetscObject) prob);
294:   for (b = prob->boundary; b; b = b->next) {
295:     char       optname[1024];
296:     PetscInt   ids[1024], len = 1024;
297:     PetscBool  flg;

299:     PetscSNPrintf(optname, sizeof(optname), "-bc_%s", b->name);
300:     PetscMemzero(ids, sizeof(ids));
301:     PetscOptionsIntArray(optname, "List of boundary IDs", "", ids, &len, &flg);
302:     if (flg) {
303:       b->numids = len;
304:       PetscFree(b->ids);
305:       PetscMalloc1(len, &b->ids);
306:       PetscMemcpy(b->ids, ids, len*sizeof(PetscInt));
307:     }
308:     len = 1024;
309:     PetscSNPrintf(optname, sizeof(optname), "-bc_%s_comp", b->name);
310:     PetscMemzero(ids, sizeof(ids));
311:     PetscOptionsIntArray(optname, "List of boundary field components", "", ids, &len, &flg);
312:     if (flg) {
313:       b->numcomps = len;
314:       PetscFree(b->comps);
315:       PetscMalloc1(len, &b->comps);
316:       PetscMemcpy(b->comps, ids, len*sizeof(PetscInt));
317:     }
318:   }
319:   PetscOptionsFList("-petscds_type", "Discrete System", "PetscDSSetType", PetscDSList, defaultType, name, 256, &flg);
320:   if (flg) {
321:     PetscDSSetType(prob, name);
322:   } else if (!((PetscObject) prob)->type_name) {
323:     PetscDSSetType(prob, defaultType);
324:   }
325:   PetscOptionsBool("-petscds_jac_pre", "Discrete System", "PetscDSUseJacobianPreconditioner", prob->useJacPre, &prob->useJacPre, &flg);
326:   if (prob->ops->setfromoptions) {(*prob->ops->setfromoptions)(prob);}
327:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
328:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) prob);
329:   PetscOptionsEnd();
330:   if (prob->Nf) {PetscDSViewFromOptions(prob, NULL, "-petscds_view");}
331:   return(0);
332: }

334: /*@C
335:   PetscDSSetUp - Construct data structures for the PetscDS

337:   Collective on PetscDS

339:   Input Parameter:
340: . prob - the PetscDS object to setup

342:   Level: developer

344: .seealso PetscDSView(), PetscDSDestroy()
345: @*/
346: PetscErrorCode PetscDSSetUp(PetscDS prob)
347: {
348:   const PetscInt Nf = prob->Nf;
349:   PetscInt       dim, dimEmbed, work, NcMax = 0, NqMax = 0, NsMax = 1, f;

354:   if (prob->setup) return(0);
355:   /* Calculate sizes */
356:   PetscDSGetSpatialDimension(prob, &dim);
357:   PetscDSGetCoordinateDimension(prob, &dimEmbed);
358:   prob->totDim = prob->totComp = 0;
359:   PetscMalloc2(Nf,&prob->Nc,Nf,&prob->Nb);
360:   PetscCalloc2(Nf+1,&prob->off,Nf+1,&prob->offDer);
361:   PetscMalloc4(Nf,&prob->basis,Nf,&prob->basisDer,Nf,&prob->basisFace,Nf,&prob->basisDerFace);
362:   for (f = 0; f < Nf; ++f) {
363:     PetscObject     obj;
364:     PetscClassId    id;
365:     PetscQuadrature q;
366:     PetscInt        Nq = 0, Nb, Nc;

368:     PetscDSGetDiscretization(prob, f, &obj);
369:     PetscObjectGetClassId(obj, &id);
370:     if (id == PETSCFE_CLASSID)      {
371:       PetscFE fe = (PetscFE) obj;

373:       PetscFEGetQuadrature(fe, &q);
374:       PetscFEGetDimension(fe, &Nb);
375:       PetscFEGetNumComponents(fe, &Nc);
376:       PetscFEGetDefaultTabulation(fe, &prob->basis[f], &prob->basisDer[f], NULL);
377:       PetscFEGetFaceTabulation(fe, &prob->basisFace[f], &prob->basisDerFace[f], NULL);
378:     } else if (id == PETSCFV_CLASSID) {
379:       PetscFV fv = (PetscFV) obj;

381:       PetscFVGetQuadrature(fv, &q);
382:       PetscFVGetNumComponents(fv, &Nc);
383:       Nb   = Nc;
384:       PetscFVGetDefaultTabulation(fv, &prob->basis[f], &prob->basisDer[f], NULL);
385:       /* TODO: should PetscFV also have face tabulation? Otherwise there will be a null pointer in prob->basisFace */
386:     } else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
387:     prob->Nc[f]       = Nc;
388:     prob->Nb[f]       = Nb;
389:     prob->off[f+1]    = Nc     + prob->off[f];
390:     prob->offDer[f+1] = Nc*dim + prob->offDer[f];
391:     if (q) {PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);}
392:     NqMax          = PetscMax(NqMax, Nq);
393:     NcMax          = PetscMax(NcMax, Nc);
394:     prob->totDim  += Nb;
395:     prob->totComp += Nc;
396:     /* There are two faces for all fields but the cohesive field on a hybrid cell */
397:     if (prob->isHybrid && (f < Nf-1)) prob->totDim += Nb;
398:   }
399:   work = PetscMax(prob->totComp*dim, PetscSqr(NcMax*dim));
400:   /* Allocate works space */
401:   if (prob->isHybrid) NsMax = 2;
402:   PetscMalloc5(NsMax*prob->totComp,&prob->u,NsMax*prob->totComp,&prob->u_t,NsMax*prob->totComp*dimEmbed,&prob->u_x,dimEmbed,&prob->x,work,&prob->refSpaceDer);
403:   PetscMalloc6(NsMax*NqMax*NcMax,&prob->f0,NsMax*NqMax*NcMax*dim,&prob->f1,
404:                       NsMax*NsMax*NqMax*NcMax*NcMax,&prob->g0,NsMax*NsMax*NqMax*NcMax*NcMax*dim,&prob->g1,
405:                       NsMax*NsMax*NqMax*NcMax*NcMax*dim,&prob->g2,NsMax*NsMax*NqMax*NcMax*NcMax*dim*dim,&prob->g3);
406:   if (prob->ops->setup) {(*prob->ops->setup)(prob);}
407:   prob->setup = PETSC_TRUE;
408:   return(0);
409: }

411: static PetscErrorCode PetscDSDestroyStructs_Static(PetscDS prob)
412: {

416:   PetscFree2(prob->Nc,prob->Nb);
417:   PetscFree2(prob->off,prob->offDer);
418:   PetscFree4(prob->basis,prob->basisDer,prob->basisFace,prob->basisDerFace);
419:   PetscFree5(prob->u,prob->u_t,prob->u_x,prob->x,prob->refSpaceDer);
420:   PetscFree6(prob->f0,prob->f1,prob->g0,prob->g1,prob->g2,prob->g3);
421:   return(0);
422: }

424: static PetscErrorCode PetscDSEnlarge_Static(PetscDS prob, PetscInt NfNew)
425: {
426:   PetscObject      *tmpd;
427:   PetscBool        *tmpi;
428:   PetscPointFunc   *tmpobj, *tmpf, *tmpup;
429:   PetscPointJac    *tmpg, *tmpgp, *tmpgt;
430:   PetscBdPointFunc *tmpfbd;
431:   PetscBdPointJac  *tmpgbd;
432:   PetscRiemannFunc *tmpr;
433:   PetscSimplePointFunc *tmpexactSol;
434:   void            **tmpctx;
435:   PetscInt          Nf = prob->Nf, f;
436:   PetscErrorCode    ierr;

439:   if (Nf >= NfNew) return(0);
440:   prob->setup = PETSC_FALSE;
441:   PetscDSDestroyStructs_Static(prob);
442:   PetscMalloc2(NfNew, &tmpd, NfNew, &tmpi);
443:   for (f = 0; f < Nf; ++f) {tmpd[f] = prob->disc[f]; tmpi[f] = prob->implicit[f];}
444:   for (f = Nf; f < NfNew; ++f) {tmpd[f] = NULL; tmpi[f] = PETSC_TRUE;}
445:   PetscFree2(prob->disc, prob->implicit);
446:   prob->Nf        = NfNew;
447:   prob->disc      = tmpd;
448:   prob->implicit  = tmpi;
449:   PetscCalloc7(NfNew, &tmpobj, NfNew*2, &tmpf, NfNew*NfNew*4, &tmpg, NfNew*NfNew*4, &tmpgp, NfNew*NfNew*4, &tmpgt, NfNew, &tmpr, NfNew, &tmpctx);
450:   PetscCalloc1(NfNew, &tmpup);
451:   for (f = 0; f < Nf; ++f) tmpobj[f] = prob->obj[f];
452:   for (f = 0; f < Nf*2; ++f) tmpf[f] = prob->f[f];
453:   for (f = 0; f < Nf*Nf*4; ++f) tmpg[f] = prob->g[f];
454:   for (f = 0; f < Nf*Nf*4; ++f) tmpgp[f] = prob->gp[f];
455:   for (f = 0; f < Nf; ++f) tmpr[f] = prob->r[f];
456:   for (f = 0; f < Nf; ++f) tmpup[f] = prob->update[f];
457:   for (f = 0; f < Nf; ++f) tmpctx[f] = prob->ctx[f];
458:   for (f = Nf; f < NfNew; ++f) tmpobj[f] = NULL;
459:   for (f = Nf*2; f < NfNew*2; ++f) tmpf[f] = NULL;
460:   for (f = Nf*Nf*4; f < NfNew*NfNew*4; ++f) tmpg[f] = NULL;
461:   for (f = Nf*Nf*4; f < NfNew*NfNew*4; ++f) tmpgp[f] = NULL;
462:   for (f = Nf*Nf*4; f < NfNew*NfNew*4; ++f) tmpgt[f] = NULL;
463:   for (f = Nf; f < NfNew; ++f) tmpr[f] = NULL;
464:   for (f = Nf; f < NfNew; ++f) tmpup[f] = NULL;
465:   for (f = Nf; f < NfNew; ++f) tmpctx[f] = NULL;
466:   PetscFree7(prob->obj, prob->f, prob->g, prob->gp, prob->gt, prob->r, prob->ctx);
467:   PetscFree(prob->update);
468:   prob->obj = tmpobj;
469:   prob->f   = tmpf;
470:   prob->g   = tmpg;
471:   prob->gp  = tmpgp;
472:   prob->gt  = tmpgt;
473:   prob->r   = tmpr;
474:   prob->update = tmpup;
475:   prob->ctx = tmpctx;
476:   PetscCalloc3(NfNew*2, &tmpfbd, NfNew*NfNew*4, &tmpgbd, NfNew, &tmpexactSol);
477:   for (f = 0; f < Nf*2; ++f) tmpfbd[f] = prob->fBd[f];
478:   for (f = 0; f < Nf*Nf*4; ++f) tmpgbd[f] = prob->gBd[f];
479:   for (f = 0; f < Nf; ++f) tmpexactSol[f] = prob->exactSol[f];
480:   for (f = Nf*2; f < NfNew*2; ++f) tmpfbd[f] = NULL;
481:   for (f = Nf*Nf*4; f < NfNew*NfNew*4; ++f) tmpgbd[f] = NULL;
482:   for (f = Nf; f < NfNew; ++f) tmpexactSol[f] = NULL;
483:   PetscFree3(prob->fBd, prob->gBd, prob->exactSol);
484:   prob->fBd = tmpfbd;
485:   prob->gBd = tmpgbd;
486:   prob->exactSol = tmpexactSol;
487:   return(0);
488: }

490: /*@
491:   PetscDSDestroy - Destroys a PetscDS object

493:   Collective on PetscDS

495:   Input Parameter:
496: . prob - the PetscDS object to destroy

498:   Level: developer

500: .seealso PetscDSView()
501: @*/
502: PetscErrorCode PetscDSDestroy(PetscDS *prob)
503: {
504:   PetscInt       f;
505:   DSBoundary     next;

509:   if (!*prob) return(0);

512:   if (--((PetscObject)(*prob))->refct > 0) {*prob = 0; return(0);}
513:   ((PetscObject) (*prob))->refct = 0;
514:   if ((*prob)->subprobs) {
515:     PetscInt dim, d;

517:     PetscDSGetSpatialDimension(*prob, &dim);
518:     for (d = 0; d < dim; ++d) {PetscDSDestroy(&(*prob)->subprobs[d]);}
519:   }
520:   PetscFree((*prob)->subprobs);
521:   PetscDSDestroyStructs_Static(*prob);
522:   for (f = 0; f < (*prob)->Nf; ++f) {
523:     PetscObjectDereference((*prob)->disc[f]);
524:   }
525:   PetscFree2((*prob)->disc, (*prob)->implicit);
526:   PetscFree7((*prob)->obj,(*prob)->f,(*prob)->g,(*prob)->gp,(*prob)->gt,(*prob)->r,(*prob)->ctx);
527:   PetscFree((*prob)->update);
528:   PetscFree3((*prob)->fBd,(*prob)->gBd,(*prob)->exactSol);
529:   if ((*prob)->ops->destroy) {(*(*prob)->ops->destroy)(*prob);}
530:   next = (*prob)->boundary;
531:   while (next) {
532:     DSBoundary b = next;

534:     next = b->next;
535:     PetscFree(b->comps);
536:     PetscFree(b->ids);
537:     PetscFree(b->name);
538:     PetscFree(b->labelname);
539:     PetscFree(b);
540:   }
541:   PetscFree((*prob)->constants);
542:   PetscHeaderDestroy(prob);
543:   return(0);
544: }

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

549:   Collective on MPI_Comm

551:   Input Parameter:
552: . comm - The communicator for the PetscDS object

554:   Output Parameter:
555: . prob - The PetscDS object

557:   Level: beginner

559: .seealso: PetscDSSetType(), PETSCDSBASIC
560: @*/
561: PetscErrorCode PetscDSCreate(MPI_Comm comm, PetscDS *prob)
562: {
563:   PetscDS   p;

568:   *prob  = NULL;
569:   PetscDSInitializePackage();

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

573:   p->Nf           = 0;
574:   p->setup        = PETSC_FALSE;
575:   p->numConstants = 0;
576:   p->constants    = NULL;
577:   p->dimEmbed     = -1;
578:   p->useJacPre    = PETSC_TRUE;

580:   *prob = p;
581:   return(0);
582: }

584: /*@
585:   PetscDSGetNumFields - Returns the number of fields in the DS

587:   Not collective

589:   Input Parameter:
590: . prob - The PetscDS object

592:   Output Parameter:
593: . Nf - The number of fields

595:   Level: beginner

597: .seealso: PetscDSGetSpatialDimension(), PetscDSCreate()
598: @*/
599: PetscErrorCode PetscDSGetNumFields(PetscDS prob, PetscInt *Nf)
600: {
604:   *Nf = prob->Nf;
605:   return(0);
606: }

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

611:   Not collective

613:   Input Parameter:
614: . prob - The PetscDS object

616:   Output Parameter:
617: . dim - The spatial dimension

619:   Level: beginner

621: .seealso: PetscDSGetCoordinateDimension(), PetscDSGetNumFields(), PetscDSCreate()
622: @*/
623: PetscErrorCode PetscDSGetSpatialDimension(PetscDS prob, PetscInt *dim)
624: {

630:   *dim = 0;
631:   if (prob->Nf) {
632:     PetscObject  obj;
633:     PetscClassId id;

635:     PetscDSGetDiscretization(prob, 0, &obj);
636:     PetscObjectGetClassId(obj, &id);
637:     if (id == PETSCFE_CLASSID)      {PetscFEGetSpatialDimension((PetscFE) obj, dim);}
638:     else if (id == PETSCFV_CLASSID) {PetscFVGetSpatialDimension((PetscFV) obj, dim);}
639:     else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", 0);
640:   }
641:   return(0);
642: }

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

647:   Not collective

649:   Input Parameter:
650: . prob - The PetscDS object

652:   Output Parameter:
653: . dimEmbed - The coordinate dimension

655:   Level: beginner

657: .seealso: PetscDSSetCoordinateDimension(), PetscDSGetSpatialDimension(), PetscDSGetNumFields(), PetscDSCreate()
658: @*/
659: PetscErrorCode PetscDSGetCoordinateDimension(PetscDS prob, PetscInt *dimEmbed)
660: {
664:   if (prob->dimEmbed < 0) SETERRQ(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONGSTATE, "No coordinate dimension set for this DS");
665:   *dimEmbed = prob->dimEmbed;
666:   return(0);
667: }

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

672:   Logically collective on DS

674:   Input Parameters:
675: + prob - The PetscDS object
676: - dimEmbed - The coordinate dimension

678:   Level: beginner

680: .seealso: PetscDSGetCoordinateDimension(), PetscDSGetSpatialDimension(), PetscDSGetNumFields(), PetscDSCreate()
681: @*/
682: PetscErrorCode PetscDSSetCoordinateDimension(PetscDS prob, PetscInt dimEmbed)
683: {
686:   if (dimEmbed < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Coordinate dimension must be non-negative, not %D", dimEmbed);
687:   prob->dimEmbed = dimEmbed;
688:   return(0);
689: }

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

694:   Not collective

696:   Input Parameter:
697: . prob - The PetscDS object

699:   Output Parameter:
700: . isHybrid - The flag

702:   Level: developer

704: .seealso: PetscDSSetHybrid(), PetscDSCreate()
705: @*/
706: PetscErrorCode PetscDSGetHybrid(PetscDS prob, PetscBool *isHybrid)
707: {
711:   *isHybrid = prob->isHybrid;
712:   return(0);
713: }

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

718:   Not collective

720:   Input Parameters:
721: + prob - The PetscDS object
722: - isHybrid - The flag

724:   Level: developer

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

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

739:   Not collective

741:   Input Parameter:
742: . prob - The PetscDS object

744:   Output Parameter:
745: . dim - The total problem dimension

747:   Level: beginner

749: .seealso: PetscDSGetNumFields(), PetscDSCreate()
750: @*/
751: PetscErrorCode PetscDSGetTotalDimension(PetscDS prob, PetscInt *dim)
752: {

757:   PetscDSSetUp(prob);
759:   *dim = prob->totDim;
760:   return(0);
761: }

763: /*@
764:   PetscDSGetTotalComponents - Returns the total number of components in this system

766:   Not collective

768:   Input Parameter:
769: . prob - The PetscDS object

771:   Output Parameter:
772: . dim - The total number of components

774:   Level: beginner

776: .seealso: PetscDSGetNumFields(), PetscDSCreate()
777: @*/
778: PetscErrorCode PetscDSGetTotalComponents(PetscDS prob, PetscInt *Nc)
779: {

784:   PetscDSSetUp(prob);
786:   *Nc = prob->totComp;
787:   return(0);
788: }

790: /*@
791:   PetscDSGetDiscretization - Returns the discretization object for the given field

793:   Not collective

795:   Input Parameters:
796: + prob - The PetscDS object
797: - f - The field number

799:   Output Parameter:
800: . disc - The discretization object

802:   Level: beginner

804: .seealso: PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
805: @*/
806: PetscErrorCode PetscDSGetDiscretization(PetscDS prob, PetscInt f, PetscObject *disc)
807: {
811:   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);
812:   *disc = prob->disc[f];
813:   return(0);
814: }

816: /*@
817:   PetscDSSetDiscretization - Sets the discretization object for the given field

819:   Not collective

821:   Input Parameters:
822: + prob - The PetscDS object
823: . f - The field number
824: - disc - The discretization object

826:   Level: beginner

828: .seealso: PetscDSGetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
829: @*/
830: PetscErrorCode PetscDSSetDiscretization(PetscDS prob, PetscInt f, PetscObject disc)
831: {

837:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
838:   PetscDSEnlarge_Static(prob, f+1);
839:   PetscObjectDereference(prob->disc[f]);
840:   prob->disc[f] = disc;
841:   PetscObjectReference(disc);
842:   {
843:     PetscClassId id;

845:     PetscObjectGetClassId(disc, &id);
846:     if (id == PETSCFE_CLASSID) {
847:       PetscDSSetImplicit(prob, f, PETSC_TRUE);
848:     } else if (id == PETSCFV_CLASSID) {
849:       PetscDSSetImplicit(prob, f, PETSC_FALSE);
850:     }
851:   }
852:   return(0);
853: }

855: /*@
856:   PetscDSAddDiscretization - Adds a discretization object

858:   Not collective

860:   Input Parameters:
861: + prob - The PetscDS object
862: - disc - The boundary discretization object

864:   Level: beginner

866: .seealso: PetscDSGetDiscretization(), PetscDSSetDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
867: @*/
868: PetscErrorCode PetscDSAddDiscretization(PetscDS prob, PetscObject disc)
869: {

873:   PetscDSSetDiscretization(prob, prob->Nf, disc);
874:   return(0);
875: }

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

880:   Not collective

882:   Input Parameters:
883: + prob - The PetscDS object
884: - f - The field number

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

889:   Level: developer

891: .seealso: PetscDSSetImplicit(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
892: @*/
893: PetscErrorCode PetscDSGetImplicit(PetscDS prob, PetscInt f, PetscBool *implicit)
894: {
898:   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);
899:   *implicit = prob->implicit[f];
900:   return(0);
901: }

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

906:   Not collective

908:   Input Parameters:
909: + prob - The PetscDS object
910: . f - The field number
911: - implicit - The flag indicating what kind of solve to use for this field

913:   Level: developer

915: .seealso: PetscDSGetImplicit(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
916: @*/
917: PetscErrorCode PetscDSSetImplicit(PetscDS prob, PetscInt f, PetscBool implicit)
918: {
921:   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);
922:   prob->implicit[f] = implicit;
923:   return(0);
924: }

926: PetscErrorCode PetscDSGetObjective(PetscDS prob, PetscInt f,
927:                                    void (**obj)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
928:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
929:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
930:                                                 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[]))
931: {
935:   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
936:   *obj = prob->obj[f];
937:   return(0);
938: }

940: PetscErrorCode PetscDSSetObjective(PetscDS prob, PetscInt f,
941:                                    void (*obj)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
942:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
943:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
944:                                                PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[]))
945: {

951:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
952:   PetscDSEnlarge_Static(prob, f+1);
953:   prob->obj[f] = obj;
954:   return(0);
955: }

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

960:   Not collective

962:   Input Parameters:
963: + prob - The PetscDS
964: - f    - The test field number

966:   Output Parameters:
967: + f0 - integrand for the test function term
968: - f1 - integrand for the test function gradient term

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

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

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

976: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
977: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
978: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
979: $    PetscReal t, const PetscReal x[], PetscScalar f0[])

981: + dim - the spatial dimension
982: . Nf - the number of fields
983: . uOff - the offset into u[] and u_t[] for each field
984: . uOff_x - the offset into u_x[] for each field
985: . u - each field evaluated at the current point
986: . u_t - the time derivative of each field evaluated at the current point
987: . u_x - the gradient of each field evaluated at the current point
988: . aOff - the offset into a[] and a_t[] for each auxiliary field
989: . aOff_x - the offset into a_x[] for each auxiliary field
990: . a - each auxiliary field evaluated at the current point
991: . a_t - the time derivative of each auxiliary field evaluated at the current point
992: . a_x - the gradient of auxiliary each field evaluated at the current point
993: . t - current time
994: . x - coordinates of the current point
995: . numConstants - number of constant parameters
996: . constants - constant parameters
997: - f0 - output values at the current point

999:   Level: intermediate

1001: .seealso: PetscDSSetResidual()
1002: @*/
1003: PetscErrorCode PetscDSGetResidual(PetscDS prob, PetscInt f,
1004:                                   void (**f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1005:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1006:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1007:                                               PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]),
1008:                                   void (**f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1009:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1010:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1011:                                               PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
1012: {
1015:   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
1018:   return(0);
1019: }

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

1024:   Not collective

1026:   Input Parameters:
1027: + prob - The PetscDS
1028: . f    - The test field number
1029: . f0 - integrand for the test function term
1030: - f1 - integrand for the test function gradient term

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

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

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

1038: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1039: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1040: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1041: $    PetscReal t, const PetscReal x[], PetscScalar f0[])

1043: + dim - the spatial dimension
1044: . Nf - the number of fields
1045: . uOff - the offset into u[] and u_t[] for each field
1046: . uOff_x - the offset into u_x[] for each field
1047: . u - each field evaluated at the current point
1048: . u_t - the time derivative of each field evaluated at the current point
1049: . u_x - the gradient of each field evaluated at the current point
1050: . aOff - the offset into a[] and a_t[] for each auxiliary field
1051: . aOff_x - the offset into a_x[] for each auxiliary field
1052: . a - each auxiliary field evaluated at the current point
1053: . a_t - the time derivative of each auxiliary field evaluated at the current point
1054: . a_x - the gradient of auxiliary each field evaluated at the current point
1055: . t - current time
1056: . x - coordinates of the current point
1057: . numConstants - number of constant parameters
1058: . constants - constant parameters
1059: - f0 - output values at the current point

1061:   Level: intermediate

1063: .seealso: PetscDSGetResidual()
1064: @*/
1065: PetscErrorCode PetscDSSetResidual(PetscDS prob, PetscInt f,
1066:                                   void (*f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1067:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1068:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1069:                                              PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]),
1070:                                   void (*f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1071:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1072:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1073:                                              PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
1074: {

1081:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1082:   PetscDSEnlarge_Static(prob, f+1);
1083:   prob->f[f*2+0] = f0;
1084:   prob->f[f*2+1] = f1;
1085:   return(0);
1086: }

1088: /*@C
1089:   PetscDSHasJacobian - Signals that Jacobian functions have been set

1091:   Not collective

1093:   Input Parameter:
1094: . prob - The PetscDS

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

1099:   Level: intermediate

1101: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1102: @*/
1103: PetscErrorCode PetscDSHasJacobian(PetscDS prob, PetscBool *hasJac)
1104: {
1105:   PetscInt f, g, h;

1109:   *hasJac = PETSC_FALSE;
1110:   for (f = 0; f < prob->Nf; ++f) {
1111:     for (g = 0; g < prob->Nf; ++g) {
1112:       for (h = 0; h < 4; ++h) {
1113:         if (prob->g[(f*prob->Nf + g)*4+h]) *hasJac = PETSC_TRUE;
1114:       }
1115:     }
1116:   }
1117:   return(0);
1118: }

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

1123:   Not collective

1125:   Input Parameters:
1126: + prob - The PetscDS
1127: . f    - The test field number
1128: - g    - The field number

1130:   Output Parameters:
1131: + g0 - integrand for the test and basis function term
1132: . g1 - integrand for the test function and basis function gradient term
1133: . g2 - integrand for the test function gradient and basis function term
1134: - g3 - integrand for the test function gradient and basis function gradient term

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

1138:   \int_\Omega \phi g_0(u, u_t, \nabla u, x, t) \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi

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

1142: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1143: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1144: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1145: $    PetscReal t, const PetscReal u_tShift, const PetscReal x[], PetscScalar g0[])

1147: + dim - the spatial dimension
1148: . Nf - the number of fields
1149: . uOff - the offset into u[] and u_t[] for each field
1150: . uOff_x - the offset into u_x[] for each field
1151: . u - each field evaluated at the current point
1152: . u_t - the time derivative of each field evaluated at the current point
1153: . u_x - the gradient of each field evaluated at the current point
1154: . aOff - the offset into a[] and a_t[] for each auxiliary field
1155: . aOff_x - the offset into a_x[] for each auxiliary field
1156: . a - each auxiliary field evaluated at the current point
1157: . a_t - the time derivative of each auxiliary field evaluated at the current point
1158: . a_x - the gradient of auxiliary each field evaluated at the current point
1159: . t - current time
1160: . u_tShift - the multiplier a for dF/dU_t
1161: . x - coordinates of the current point
1162: . numConstants - number of constant parameters
1163: . constants - constant parameters
1164: - g0 - output values at the current point

1166:   Level: intermediate

1168: .seealso: PetscDSSetJacobian()
1169: @*/
1170: PetscErrorCode PetscDSGetJacobian(PetscDS prob, PetscInt f, PetscInt g,
1171:                                   void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1172:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1173:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1174:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1175:                                   void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1176:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1177:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1178:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1179:                                   void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1180:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1181:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1182:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1183:                                   void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1184:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1185:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1186:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1187: {
1190:   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);
1191:   if ((g < 0) || (g >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, prob->Nf);
1196:   return(0);
1197: }

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

1202:   Not collective

1204:   Input Parameters:
1205: + prob - The PetscDS
1206: . f    - The test field number
1207: . g    - The field number
1208: . g0 - integrand for the test and basis function term
1209: . g1 - integrand for the test function and basis function gradient term
1210: . g2 - integrand for the test function gradient and basis function term
1211: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

1219: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1220: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1221: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1222: $    PetscReal t, const PetscReal x[], PetscScalar g0[])

1224: + dim - the spatial dimension
1225: . Nf - the number of fields
1226: . uOff - the offset into u[] and u_t[] for each field
1227: . uOff_x - the offset into u_x[] for each field
1228: . u - each field evaluated at the current point
1229: . u_t - the time derivative of each field evaluated at the current point
1230: . u_x - the gradient of each field evaluated at the current point
1231: . aOff - the offset into a[] and a_t[] for each auxiliary field
1232: . aOff_x - the offset into a_x[] for each auxiliary field
1233: . a - each auxiliary field evaluated at the current point
1234: . a_t - the time derivative of each auxiliary field evaluated at the current point
1235: . a_x - the gradient of auxiliary each field evaluated at the current point
1236: . t - current time
1237: . u_tShift - the multiplier a for dF/dU_t
1238: . x - coordinates of the current point
1239: . numConstants - number of constant parameters
1240: . constants - constant parameters
1241: - g0 - output values at the current point

1243:   Level: intermediate

1245: .seealso: PetscDSGetJacobian()
1246: @*/
1247: PetscErrorCode PetscDSSetJacobian(PetscDS prob, PetscInt f, PetscInt g,
1248:                                   void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1249:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1250:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1251:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1252:                                   void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1253:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1254:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1255:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1256:                                   void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1257:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1258:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1259:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1260:                                   void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1261:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1262:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1263:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1264: {

1273:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1274:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
1275:   PetscDSEnlarge_Static(prob, PetscMax(f, g)+1);
1276:   prob->g[(f*prob->Nf + g)*4+0] = g0;
1277:   prob->g[(f*prob->Nf + g)*4+1] = g1;
1278:   prob->g[(f*prob->Nf + g)*4+2] = g2;
1279:   prob->g[(f*prob->Nf + g)*4+3] = g3;
1280:   return(0);
1281: }

1283: /*@C
1284:   PetscDSUseJacobianPreconditioner - Whether to construct a Jacobian preconditioner

1286:   Not collective

1288:   Input Parameters:
1289: + prob - The PetscDS
1290: - useJacPre - flag that enables construction of a Jacobian preconditioner

1292:   Level: intermediate

1294: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1295: @*/
1296: PetscErrorCode PetscDSUseJacobianPreconditioner(PetscDS prob, PetscBool useJacPre)
1297: {
1300:   prob->useJacPre = useJacPre;
1301:   return(0);
1302: }

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

1307:   Not collective

1309:   Input Parameter:
1310: . prob - The PetscDS

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

1315:   Level: intermediate

1317: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1318: @*/
1319: PetscErrorCode PetscDSHasJacobianPreconditioner(PetscDS prob, PetscBool *hasJacPre)
1320: {
1321:   PetscInt f, g, h;

1325:   *hasJacPre = PETSC_FALSE;
1326:   if (!prob->useJacPre) return(0);
1327:   for (f = 0; f < prob->Nf; ++f) {
1328:     for (g = 0; g < prob->Nf; ++g) {
1329:       for (h = 0; h < 4; ++h) {
1330:         if (prob->gp[(f*prob->Nf + g)*4+h]) *hasJacPre = PETSC_TRUE;
1331:       }
1332:     }
1333:   }
1334:   return(0);
1335: }

1337: /*@C
1338:   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.

1340:   Not collective

1342:   Input Parameters:
1343: + prob - The PetscDS
1344: . f    - The test field number
1345: - g    - The field number

1347:   Output Parameters:
1348: + g0 - integrand for the test and basis function term
1349: . g1 - integrand for the test function and basis function gradient term
1350: . g2 - integrand for the test function gradient and basis function term
1351: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

1359: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1360: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1361: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1362: $    PetscReal t, const PetscReal u_tShift, const PetscReal x[], PetscScalar g0[])

1364: + dim - the spatial dimension
1365: . Nf - the number of fields
1366: . uOff - the offset into u[] and u_t[] for each field
1367: . uOff_x - the offset into u_x[] for each field
1368: . u - each field evaluated at the current point
1369: . u_t - the time derivative of each field evaluated at the current point
1370: . u_x - the gradient of each field evaluated at the current point
1371: . aOff - the offset into a[] and a_t[] for each auxiliary field
1372: . aOff_x - the offset into a_x[] for each auxiliary field
1373: . a - each auxiliary field evaluated at the current point
1374: . a_t - the time derivative of each auxiliary field evaluated at the current point
1375: . a_x - the gradient of auxiliary each field evaluated at the current point
1376: . t - current time
1377: . u_tShift - the multiplier a for dF/dU_t
1378: . x - coordinates of the current point
1379: . numConstants - number of constant parameters
1380: . constants - constant parameters
1381: - g0 - output values at the current point

1383:   Level: intermediate

1385: .seealso: PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1386: @*/
1387: PetscErrorCode PetscDSGetJacobianPreconditioner(PetscDS prob, PetscInt f, PetscInt g,
1388:                                   void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1389:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1390:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1391:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1392:                                   void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1393:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1394:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1395:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1396:                                   void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1397:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1398:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1399:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1400:                                   void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1401:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1402:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1403:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1404: {
1407:   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);
1408:   if ((g < 0) || (g >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, prob->Nf);
1413:   return(0);
1414: }

1416: /*@C
1417:   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.

1419:   Not collective

1421:   Input Parameters:
1422: + prob - The PetscDS
1423: . f    - The test field number
1424: . g    - The field number
1425: . g0 - integrand for the test and basis function term
1426: . g1 - integrand for the test function and basis function gradient term
1427: . g2 - integrand for the test function gradient and basis function term
1428: - g3 - integrand for the test function gradient and basis function gradient term

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

1432:   \int_\Omega \phi g_0(u, u_t, \nabla u, x, t) \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi

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

1436: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1437: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1438: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1439: $    PetscReal t, const PetscReal x[], PetscScalar g0[])

1441: + dim - the spatial dimension
1442: . Nf - the number of fields
1443: . uOff - the offset into u[] and u_t[] for each field
1444: . uOff_x - the offset into u_x[] for each field
1445: . u - each field evaluated at the current point
1446: . u_t - the time derivative of each field evaluated at the current point
1447: . u_x - the gradient of each field evaluated at the current point
1448: . aOff - the offset into a[] and a_t[] for each auxiliary field
1449: . aOff_x - the offset into a_x[] for each auxiliary field
1450: . a - each auxiliary field evaluated at the current point
1451: . a_t - the time derivative of each auxiliary field evaluated at the current point
1452: . a_x - the gradient of auxiliary each field evaluated at the current point
1453: . t - current time
1454: . u_tShift - the multiplier a for dF/dU_t
1455: . x - coordinates of the current point
1456: . numConstants - number of constant parameters
1457: . constants - constant parameters
1458: - g0 - output values at the current point

1460:   Level: intermediate

1462: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobian()
1463: @*/
1464: PetscErrorCode PetscDSSetJacobianPreconditioner(PetscDS prob, PetscInt f, PetscInt g,
1465:                                   void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1466:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1467:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1468:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1469:                                   void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1470:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1471:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1472:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1473:                                   void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1474:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1475:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1476:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1477:                                   void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1478:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1479:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1480:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1481: {

1490:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1491:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
1492:   PetscDSEnlarge_Static(prob, PetscMax(f, g)+1);
1493:   prob->gp[(f*prob->Nf + g)*4+0] = g0;
1494:   prob->gp[(f*prob->Nf + g)*4+1] = g1;
1495:   prob->gp[(f*prob->Nf + g)*4+2] = g2;
1496:   prob->gp[(f*prob->Nf + g)*4+3] = g3;
1497:   return(0);
1498: }

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

1503:   Not collective

1505:   Input Parameter:
1506: . prob - The PetscDS

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

1511:   Level: intermediate

1513: .seealso: PetscDSGetDynamicJacobian(), PetscDSSetDynamicJacobian(), PetscDSGetJacobian()
1514: @*/
1515: PetscErrorCode PetscDSHasDynamicJacobian(PetscDS prob, PetscBool *hasDynJac)
1516: {
1517:   PetscInt f, g, h;

1521:   *hasDynJac = PETSC_FALSE;
1522:   for (f = 0; f < prob->Nf; ++f) {
1523:     for (g = 0; g < prob->Nf; ++g) {
1524:       for (h = 0; h < 4; ++h) {
1525:         if (prob->gt[(f*prob->Nf + g)*4+h]) *hasDynJac = PETSC_TRUE;
1526:       }
1527:     }
1528:   }
1529:   return(0);
1530: }

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

1535:   Not collective

1537:   Input Parameters:
1538: + prob - The PetscDS
1539: . f    - The test field number
1540: - g    - The field number

1542:   Output Parameters:
1543: + g0 - integrand for the test and basis function term
1544: . g1 - integrand for the test function and basis function gradient term
1545: . g2 - integrand for the test function gradient and basis function term
1546: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

1554: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1555: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1556: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1557: $    PetscReal t, const PetscReal u_tShift, const PetscReal x[], PetscScalar g0[])

1559: + dim - the spatial dimension
1560: . Nf - the number of fields
1561: . uOff - the offset into u[] and u_t[] for each field
1562: . uOff_x - the offset into u_x[] for each field
1563: . u - each field evaluated at the current point
1564: . u_t - the time derivative of each field evaluated at the current point
1565: . u_x - the gradient of each field evaluated at the current point
1566: . aOff - the offset into a[] and a_t[] for each auxiliary field
1567: . aOff_x - the offset into a_x[] for each auxiliary field
1568: . a - each auxiliary field evaluated at the current point
1569: . a_t - the time derivative of each auxiliary field evaluated at the current point
1570: . a_x - the gradient of auxiliary each field evaluated at the current point
1571: . t - current time
1572: . u_tShift - the multiplier a for dF/dU_t
1573: . x - coordinates of the current point
1574: . numConstants - number of constant parameters
1575: . constants - constant parameters
1576: - g0 - output values at the current point

1578:   Level: intermediate

1580: .seealso: PetscDSSetJacobian()
1581: @*/
1582: PetscErrorCode PetscDSGetDynamicJacobian(PetscDS prob, PetscInt f, PetscInt g,
1583:                                          void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1584:                                                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1585:                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1586:                                                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1587:                                          void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1588:                                                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1589:                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1590:                                                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1591:                                          void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1592:                                                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1593:                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1594:                                                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1595:                                          void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1596:                                                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1597:                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1598:                                                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1599: {
1602:   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);
1603:   if ((g < 0) || (g >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, prob->Nf);
1608:   return(0);
1609: }

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

1614:   Not collective

1616:   Input Parameters:
1617: + prob - The PetscDS
1618: . f    - The test field number
1619: . g    - The field number
1620: . g0 - integrand for the test and basis function term
1621: . g1 - integrand for the test function and basis function gradient term
1622: . g2 - integrand for the test function gradient and basis function term
1623: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

1631: $ g0(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, const PetscReal x[], PetscScalar g0[])

1636: + dim - the spatial dimension
1637: . Nf - the number of fields
1638: . uOff - the offset into u[] and u_t[] for each field
1639: . uOff_x - the offset into u_x[] for each field
1640: . u - each field evaluated at the current point
1641: . u_t - the time derivative of each field evaluated at the current point
1642: . u_x - the gradient of each field evaluated at the current point
1643: . aOff - the offset into a[] and a_t[] for each auxiliary field
1644: . aOff_x - the offset into a_x[] for each auxiliary field
1645: . a - each auxiliary field evaluated at the current point
1646: . a_t - the time derivative of each auxiliary field evaluated at the current point
1647: . a_x - the gradient of auxiliary each field evaluated at the current point
1648: . t - current time
1649: . u_tShift - the multiplier a for dF/dU_t
1650: . x - coordinates of the current point
1651: . numConstants - number of constant parameters
1652: . constants - constant parameters
1653: - g0 - output values at the current point

1655:   Level: intermediate

1657: .seealso: PetscDSGetJacobian()
1658: @*/
1659: PetscErrorCode PetscDSSetDynamicJacobian(PetscDS prob, PetscInt f, PetscInt g,
1660:                                          void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1661:                                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1662:                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1663:                                                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1664:                                          void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1665:                                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1666:                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1667:                                                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1668:                                          void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1669:                                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1670:                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1671:                                                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1672:                                          void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1673:                                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1674:                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1675:                                                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1676: {

1685:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1686:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
1687:   PetscDSEnlarge_Static(prob, PetscMax(f, g)+1);
1688:   prob->gt[(f*prob->Nf + g)*4+0] = g0;
1689:   prob->gt[(f*prob->Nf + g)*4+1] = g1;
1690:   prob->gt[(f*prob->Nf + g)*4+2] = g2;
1691:   prob->gt[(f*prob->Nf + g)*4+3] = g3;
1692:   return(0);
1693: }

1695: /*@C
1696:   PetscDSGetRiemannSolver - Returns the Riemann solver for the given field

1698:   Not collective

1700:   Input Arguments:
1701: + prob - The PetscDS object
1702: - f    - The field number

1704:   Output Argument:
1705: . r    - Riemann solver

1707:   Calling sequence for r:

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

1711: + dim  - The spatial dimension
1712: . Nf   - The number of fields
1713: . x    - The coordinates at a point on the interface
1714: . n    - The normal vector to the interface
1715: . uL   - The state vector to the left of the interface
1716: . uR   - The state vector to the right of the interface
1717: . flux - output array of flux through the interface
1718: . numConstants - number of constant parameters
1719: . constants - constant parameters
1720: - ctx  - optional user context

1722:   Level: intermediate

1724: .seealso: PetscDSSetRiemannSolver()
1725: @*/
1726: PetscErrorCode PetscDSGetRiemannSolver(PetscDS prob, PetscInt f,
1727:                                        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))
1728: {
1731:   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);
1733:   *r = prob->r[f];
1734:   return(0);
1735: }

1737: /*@C
1738:   PetscDSSetRiemannSolver - Sets the Riemann solver for the given field

1740:   Not collective

1742:   Input Arguments:
1743: + prob - The PetscDS object
1744: . f    - The field number
1745: - r    - Riemann solver

1747:   Calling sequence for r:

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

1751: + dim  - The spatial dimension
1752: . Nf   - The number of fields
1753: . x    - The coordinates at a point on the interface
1754: . n    - The normal vector to the interface
1755: . uL   - The state vector to the left of the interface
1756: . uR   - The state vector to the right of the interface
1757: . flux - output array of flux through the interface
1758: . numConstants - number of constant parameters
1759: . constants - constant parameters
1760: - ctx  - optional user context

1762:   Level: intermediate

1764: .seealso: PetscDSGetRiemannSolver()
1765: @*/
1766: PetscErrorCode PetscDSSetRiemannSolver(PetscDS prob, PetscInt f,
1767:                                        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))
1768: {

1774:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1775:   PetscDSEnlarge_Static(prob, f+1);
1776:   prob->r[f] = r;
1777:   return(0);
1778: }

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

1783:   Not collective

1785:   Input Parameters:
1786: + prob - The PetscDS
1787: - f    - The field number

1789:   Output Parameters:
1790: + update - update function

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

1794: $ update(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1795: $        const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1796: $        const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1797: $        PetscReal t, const PetscReal x[], PetscScalar uNew[])

1799: + dim - the spatial dimension
1800: . Nf - the number of fields
1801: . uOff - the offset into u[] and u_t[] for each field
1802: . uOff_x - the offset into u_x[] for each field
1803: . u - each field evaluated at the current point
1804: . u_t - the time derivative of each field evaluated at the current point
1805: . u_x - the gradient of each field evaluated at the current point
1806: . aOff - the offset into a[] and a_t[] for each auxiliary field
1807: . aOff_x - the offset into a_x[] for each auxiliary field
1808: . a - each auxiliary field evaluated at the current point
1809: . a_t - the time derivative of each auxiliary field evaluated at the current point
1810: . a_x - the gradient of auxiliary each field evaluated at the current point
1811: . t - current time
1812: . x - coordinates of the current point
1813: - uNew - new value for field at the current point

1815:   Level: intermediate

1817: .seealso: PetscDSSetUpdate(), PetscDSSetResidual()
1818: @*/
1819: PetscErrorCode PetscDSGetUpdate(PetscDS prob, PetscInt f,
1820:                                   void (**update)(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, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar uNew[]))
1824: {
1827:   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);
1829:   return(0);
1830: }

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

1835:   Not collective

1837:   Input Parameters:
1838: + prob   - The PetscDS
1839: . f      - The field number
1840: - update - update function

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

1844: $ update(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1845: $        const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1846: $        const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1847: $        PetscReal t, const PetscReal x[], PetscScalar uNew[])

1849: + dim - the spatial dimension
1850: . Nf - the number of fields
1851: . uOff - the offset into u[] and u_t[] for each field
1852: . uOff_x - the offset into u_x[] for each field
1853: . u - each field evaluated at the current point
1854: . u_t - the time derivative of each field evaluated at the current point
1855: . u_x - the gradient of each field evaluated at the current point
1856: . aOff - the offset into a[] and a_t[] for each auxiliary field
1857: . aOff_x - the offset into a_x[] for each auxiliary field
1858: . a - each auxiliary field evaluated at the current point
1859: . a_t - the time derivative of each auxiliary field evaluated at the current point
1860: . a_x - the gradient of auxiliary each field evaluated at the current point
1861: . t - current time
1862: . x - coordinates of the current point
1863: - uNew - new field values at the current point

1865:   Level: intermediate

1867: .seealso: PetscDSGetResidual()
1868: @*/
1869: PetscErrorCode PetscDSSetUpdate(PetscDS prob, PetscInt f,
1870:                                 void (*update)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1871:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1872:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1873:                                                PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar uNew[]))
1874: {

1880:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1881:   PetscDSEnlarge_Static(prob, f+1);
1882:   prob->update[f] = update;
1883:   return(0);
1884: }

1886: PetscErrorCode PetscDSGetContext(PetscDS prob, PetscInt f, void **ctx)
1887: {
1890:   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);
1892:   *ctx = prob->ctx[f];
1893:   return(0);
1894: }

1896: PetscErrorCode PetscDSSetContext(PetscDS prob, PetscInt f, void *ctx)
1897: {

1902:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1903:   PetscDSEnlarge_Static(prob, f+1);
1904:   prob->ctx[f] = ctx;
1905:   return(0);
1906: }

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

1911:   Not collective

1913:   Input Parameters:
1914: + prob - The PetscDS
1915: - f    - The test field number

1917:   Output Parameters:
1918: + f0 - boundary integrand for the test function term
1919: - f1 - boundary integrand for the test function gradient term

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

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

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

1927: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1928: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1929: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1930: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[])

1932: + dim - the spatial dimension
1933: . Nf - the number of fields
1934: . uOff - the offset into u[] and u_t[] for each field
1935: . uOff_x - the offset into u_x[] for each field
1936: . u - each field evaluated at the current point
1937: . u_t - the time derivative of each field evaluated at the current point
1938: . u_x - the gradient of each field evaluated at the current point
1939: . aOff - the offset into a[] and a_t[] for each auxiliary field
1940: . aOff_x - the offset into a_x[] for each auxiliary field
1941: . a - each auxiliary field evaluated at the current point
1942: . a_t - the time derivative of each auxiliary field evaluated at the current point
1943: . a_x - the gradient of auxiliary each field evaluated at the current point
1944: . t - current time
1945: . x - coordinates of the current point
1946: . n - unit normal at the current point
1947: . numConstants - number of constant parameters
1948: . constants - constant parameters
1949: - f0 - output values at the current point

1951:   Level: intermediate

1953: .seealso: PetscDSSetBdResidual()
1954: @*/
1955: PetscErrorCode PetscDSGetBdResidual(PetscDS prob, PetscInt f,
1956:                                     void (**f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1957:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1958:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1959:                                                 PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]),
1960:                                     void (**f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1961:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1962:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1963:                                                 PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
1964: {
1967:   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);
1970:   return(0);
1971: }

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

1976:   Not collective

1978:   Input Parameters:
1979: + prob - The PetscDS
1980: . f    - The test field number
1981: . f0 - boundary integrand for the test function term
1982: - f1 - boundary integrand for the test function gradient term

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

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

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

1990: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1991: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1992: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1993: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[])

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

2014:   Level: intermediate

2016: .seealso: PetscDSGetBdResidual()
2017: @*/
2018: PetscErrorCode PetscDSSetBdResidual(PetscDS prob, PetscInt f,
2019:                                     void (*f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2020:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2021:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2022:                                                PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]),
2023:                                     void (*f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2024:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2025:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2026:                                                PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
2027: {

2032:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2033:   PetscDSEnlarge_Static(prob, f+1);
2036:   return(0);
2037: }

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

2042:   Not collective

2044:   Input Parameters:
2045: + prob - The PetscDS
2046: . f    - The test field number
2047: - g    - The field number

2049:   Output Parameters:
2050: + g0 - integrand for the test and basis function term
2051: . g1 - integrand for the test function and basis function gradient term
2052: . g2 - integrand for the test function gradient and basis function term
2053: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

2061: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2062: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2063: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2064: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar g0[])

2066: + dim - the spatial dimension
2067: . Nf - the number of fields
2068: . uOff - the offset into u[] and u_t[] for each field
2069: . uOff_x - the offset into u_x[] for each field
2070: . u - each field evaluated at the current point
2071: . u_t - the time derivative of each field evaluated at the current point
2072: . u_x - the gradient of each field evaluated at the current point
2073: . aOff - the offset into a[] and a_t[] for each auxiliary field
2074: . aOff_x - the offset into a_x[] for each auxiliary field
2075: . a - each auxiliary field evaluated at the current point
2076: . a_t - the time derivative of each auxiliary field evaluated at the current point
2077: . a_x - the gradient of auxiliary each field evaluated at the current point
2078: . t - current time
2079: . u_tShift - the multiplier a for dF/dU_t
2080: . x - coordinates of the current point
2081: . n - normal at the current point
2082: . numConstants - number of constant parameters
2083: . constants - constant parameters
2084: - g0 - output values at the current point

2086:   Level: intermediate

2088: .seealso: PetscDSSetBdJacobian()
2089: @*/
2090: PetscErrorCode PetscDSGetBdJacobian(PetscDS prob, PetscInt f, PetscInt g,
2091:                                     void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2092:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2093:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2094:                                                 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
2095:                                     void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2096:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2097:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2098:                                                 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
2099:                                     void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2100:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2101:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2102:                                                 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
2103:                                     void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2104:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2105:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2106:                                                 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
2107: {
2110:   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);
2111:   if ((g < 0) || (g >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, prob->Nf);
2116:   return(0);
2117: }

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

2122:   Not collective

2124:   Input Parameters:
2125: + prob - The PetscDS
2126: . f    - The test field number
2127: . g    - The field number
2128: . g0 - integrand for the test and basis function term
2129: . g1 - integrand for the test function and basis function gradient term
2130: . g2 - integrand for the test function gradient and basis function term
2131: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

2139: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2140: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2141: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2142: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar g0[])

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

2164:   Level: intermediate

2166: .seealso: PetscDSGetBdJacobian()
2167: @*/
2168: PetscErrorCode PetscDSSetBdJacobian(PetscDS prob, PetscInt f, PetscInt g,
2169:                                     void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2170:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2171:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2172:                                                PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
2173:                                     void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2174:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2175:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2176:                                                PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
2177:                                     void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2178:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2179:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2180:                                                PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
2181:                                     void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2182:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2183:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2184:                                                PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
2185: {

2194:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2195:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
2196:   PetscDSEnlarge_Static(prob, PetscMax(f, g)+1);
2197:   prob->gBd[(f*prob->Nf + g)*4+0] = g0;
2198:   prob->gBd[(f*prob->Nf + g)*4+1] = g1;
2199:   prob->gBd[(f*prob->Nf + g)*4+2] = g2;
2200:   prob->gBd[(f*prob->Nf + g)*4+3] = g3;
2201:   return(0);
2202: }

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

2207:   Not collective

2209:   Input Parameters:
2210: + prob - The PetscDS
2211: - f    - The test field number

2213:   Output Parameter:
2214: . exactSol - exact solution for the test field

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

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

2220: + dim - the spatial dimension
2221: . t - current time
2222: . x - coordinates of the current point
2223: . Nc - the number of field components
2224: . u - the solution field evaluated at the current point
2225: - ctx - a user context

2227:   Level: intermediate

2229: .seealso: PetscDSSetExactSolution()
2230: @*/
2231: PetscErrorCode PetscDSGetExactSolution(PetscDS prob, PetscInt f, PetscErrorCode (**sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx))
2232: {
2235:   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);
2237:   return(0);
2238: }

2240: /*@C
2241:   PetscDSSetExactSolution - Get the pointwise exact solution function for a given test field

2243:   Not collective

2245:   Input Parameters:
2246: + prob - The PetscDS
2247: . f    - The test field number
2248: - sol  - solution function for the test fields

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

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

2254: + dim - the spatial dimension
2255: . t - current time
2256: . x - coordinates of the current point
2257: . Nc - the number of field components
2258: . u - the solution field evaluated at the current point
2259: - ctx - a user context

2261:   Level: intermediate

2263: .seealso: PetscDSGetExactSolution()
2264: @*/
2265: PetscErrorCode PetscDSSetExactSolution(PetscDS prob, PetscInt f, PetscErrorCode (*sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx))
2266: {

2271:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2272:   PetscDSEnlarge_Static(prob, f+1);
2274:   return(0);
2275: }

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

2280:   Not collective

2282:   Input Parameter:
2283: . prob - The PetscDS object

2285:   Output Parameters:
2286: + numConstants - The number of constants
2287: - constants    - The array of constants, NULL if there are none

2289:   Level: intermediate

2291: .seealso: PetscDSSetConstants(), PetscDSCreate()
2292: @*/
2293: PetscErrorCode PetscDSGetConstants(PetscDS prob, PetscInt *numConstants, const PetscScalar *constants[])
2294: {
2299:   return(0);
2300: }

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

2305:   Not collective

2307:   Input Parameters:
2308: + prob         - The PetscDS object
2309: . numConstants - The number of constants
2310: - constants    - The array of constants, NULL if there are none

2312:   Level: intermediate

2314: .seealso: PetscDSGetConstants(), PetscDSCreate()
2315: @*/
2316: PetscErrorCode PetscDSSetConstants(PetscDS prob, PetscInt numConstants, PetscScalar constants[])
2317: {

2322:   if (numConstants != prob->numConstants) {
2323:     PetscFree(prob->constants);
2324:     prob->numConstants = numConstants;
2325:     if (prob->numConstants) {
2326:       PetscMalloc1(prob->numConstants, &prob->constants);
2327:     } else {
2328:       prob->constants = NULL;
2329:     }
2330:   }
2331:   if (prob->numConstants) {
2333:     PetscMemcpy(prob->constants, constants, prob->numConstants * sizeof(PetscScalar));
2334:   }
2335:   return(0);
2336: }

2338: /*@
2339:   PetscDSGetFieldIndex - Returns the index of the given field

2341:   Not collective

2343:   Input Parameters:
2344: + prob - The PetscDS object
2345: - disc - The discretization object

2347:   Output Parameter:
2348: . f - The field number

2350:   Level: beginner

2352: .seealso: PetscGetDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
2353: @*/
2354: PetscErrorCode PetscDSGetFieldIndex(PetscDS prob, PetscObject disc, PetscInt *f)
2355: {
2356:   PetscInt g;

2361:   *f = -1;
2362:   for (g = 0; g < prob->Nf; ++g) {if (disc == prob->disc[g]) break;}
2363:   if (g == prob->Nf) SETERRQ(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Field not found in PetscDS.");
2364:   *f = g;
2365:   return(0);
2366: }

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

2371:   Not collective

2373:   Input Parameters:
2374: + prob - The PetscDS object
2375: - f - The field number

2377:   Output Parameter:
2378: . size - The size

2380:   Level: beginner

2382: .seealso: PetscDSGetFieldOffset(), PetscDSGetNumFields(), PetscDSCreate()
2383: @*/
2384: PetscErrorCode PetscDSGetFieldSize(PetscDS prob, PetscInt f, PetscInt *size)
2385: {

2391:   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);
2392:   PetscDSSetUp(prob);
2393:   *size = prob->Nb[f];
2394:   return(0);
2395: }

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

2400:   Not collective

2402:   Input Parameters:
2403: + prob - The PetscDS object
2404: - f - The field number

2406:   Output Parameter:
2407: . off - The offset

2409:   Level: beginner

2411: .seealso: PetscDSGetFieldSize(), PetscDSGetNumFields(), PetscDSCreate()
2412: @*/
2413: PetscErrorCode PetscDSGetFieldOffset(PetscDS prob, PetscInt f, PetscInt *off)
2414: {
2415:   PetscInt       size, g;

2421:   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);
2422:   *off = 0;
2423:   for (g = 0; g < f; ++g) {
2424:     PetscDSGetFieldSize(prob, g, &size);
2425:     *off += size;
2426:   }
2427:   return(0);
2428: }

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

2433:   Not collective

2435:   Input Parameter:
2436: . prob - The PetscDS object

2438:   Output Parameter:
2439: . dimensions - The number of dimensions

2441:   Level: beginner

2443: .seealso: PetscDSGetComponentOffsets(), PetscDSGetNumFields(), PetscDSCreate()
2444: @*/
2445: PetscErrorCode PetscDSGetDimensions(PetscDS prob, PetscInt *dimensions[])
2446: {

2451:   PetscDSSetUp(prob);
2453:   *dimensions = prob->Nb;
2454:   return(0);
2455: }

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

2460:   Not collective

2462:   Input Parameter:
2463: . prob - The PetscDS object

2465:   Output Parameter:
2466: . components - The number of components

2468:   Level: beginner

2470: .seealso: PetscDSGetComponentOffsets(), PetscDSGetNumFields(), PetscDSCreate()
2471: @*/
2472: PetscErrorCode PetscDSGetComponents(PetscDS prob, PetscInt *components[])
2473: {

2478:   PetscDSSetUp(prob);
2480:   *components = prob->Nc;
2481:   return(0);
2482: }

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

2487:   Not collective

2489:   Input Parameters:
2490: + prob - The PetscDS object
2491: - f - The field number

2493:   Output Parameter:
2494: . off - The offset

2496:   Level: beginner

2498: .seealso: PetscDSGetNumFields(), PetscDSCreate()
2499: @*/
2500: PetscErrorCode PetscDSGetComponentOffset(PetscDS prob, PetscInt f, PetscInt *off)
2501: {
2505:   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);
2506:   *off = prob->off[f];
2507:   return(0);
2508: }

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

2513:   Not collective

2515:   Input Parameter:
2516: . prob - The PetscDS object

2518:   Output Parameter:
2519: . offsets - The offsets

2521:   Level: beginner

2523: .seealso: PetscDSGetNumFields(), PetscDSCreate()
2524: @*/
2525: PetscErrorCode PetscDSGetComponentOffsets(PetscDS prob, PetscInt *offsets[])
2526: {
2530:   *offsets = prob->off;
2531:   return(0);
2532: }

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

2537:   Not collective

2539:   Input Parameter:
2540: . prob - The PetscDS object

2542:   Output Parameter:
2543: . offsets - The offsets

2545:   Level: beginner

2547: .seealso: PetscDSGetNumFields(), PetscDSCreate()
2548: @*/
2549: PetscErrorCode PetscDSGetComponentDerivativeOffsets(PetscDS prob, PetscInt *offsets[])
2550: {
2554:   *offsets = prob->offDer;
2555:   return(0);
2556: }

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

2561:   Not collective

2563:   Input Parameter:
2564: . prob - The PetscDS object

2566:   Output Parameters:
2567: + basis - The basis function tabulation at quadrature points
2568: - basisDer - The basis function derivative tabulation at quadrature points

2570:   Level: intermediate

2572: .seealso: PetscDSCreate()
2573: @*/
2574: PetscErrorCode PetscDSGetTabulation(PetscDS prob, PetscReal ***basis, PetscReal ***basisDer)
2575: {

2580:   PetscDSSetUp(prob);
2583:   return(0);
2584: }

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

2589:   Not collective

2591:   Input Parameter:
2592: . prob - The PetscDS object

2594:   Output Parameters:
2595: + basisFace - The basis function tabulation at quadrature points
2596: - basisDerFace - The basis function derivative tabulation at quadrature points

2598:   Level: intermediate

2600: .seealso: PetscDSGetTabulation(), PetscDSCreate()
2601: @*/
2602: PetscErrorCode PetscDSGetFaceTabulation(PetscDS prob, PetscReal ***basis, PetscReal ***basisDer)
2603: {

2608:   PetscDSSetUp(prob);
2611:   return(0);
2612: }

2614: PetscErrorCode PetscDSGetEvaluationArrays(PetscDS prob, PetscScalar **u, PetscScalar **u_t, PetscScalar **u_x)
2615: {

2620:   PetscDSSetUp(prob);
2624:   return(0);
2625: }

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

2633:   PetscDSSetUp(prob);
2640:   return(0);
2641: }

2643: PetscErrorCode PetscDSGetRefCoordArrays(PetscDS prob, PetscReal **x, PetscScalar **refSpaceDer)
2644: {

2649:   PetscDSSetUp(prob);
2652:   return(0);
2653: }

2655: /*@C
2656:   PetscDSAddBoundary - Add a boundary condition to the model

2658:   Input Parameters:
2659: + ds          - The PetscDS object
2660: . type        - The type of condition, e.g. DM_BC_ESSENTIAL/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
2661: . name        - The BC name
2662: . labelname   - The label defining constrained points
2663: . field       - The field to constrain
2664: . numcomps    - The number of constrained field components (0 will constrain all fields)
2665: . comps       - An array of constrained component numbers
2666: . bcFunc      - A pointwise function giving boundary values
2667: . numids      - The number of DMLabel ids for constrained points
2668: . ids         - An array of ids for constrained points
2669: - ctx         - An optional user context for bcFunc

2671:   Options Database Keys:
2672: + -bc_<boundary name> <num> - Overrides the boundary ids
2673: - -bc_<boundary name>_comp <num> - Overrides the boundary components

2675:   Level: developer

2677: .seealso: PetscDSGetBoundary()
2678: @*/
2679: PetscErrorCode PetscDSAddBoundary(PetscDS ds, DMBoundaryConditionType type, const char name[], const char labelname[], PetscInt field, PetscInt numcomps, const PetscInt *comps, void (*bcFunc)(void), PetscInt numids, const PetscInt *ids, void *ctx)
2680: {
2681:   DSBoundary     b;

2686:   PetscNew(&b);
2687:   PetscStrallocpy(name, (char **) &b->name);
2688:   PetscStrallocpy(labelname, (char **) &b->labelname);
2689:   PetscMalloc1(numcomps, &b->comps);
2690:   if (numcomps) {PetscMemcpy(b->comps, comps, numcomps*sizeof(PetscInt));}
2691:   PetscMalloc1(numids, &b->ids);
2692:   if (numids) {PetscMemcpy(b->ids, ids, numids*sizeof(PetscInt));}
2693:   b->type            = type;
2694:   b->field           = field;
2695:   b->numcomps        = numcomps;
2696:   b->func            = bcFunc;
2697:   b->numids          = numids;
2698:   b->ctx             = ctx;
2699:   b->next            = ds->boundary;
2700:   ds->boundary       = b;
2701:   return(0);
2702: }

2704: /*@C
2705:   PetscDSUpdateBoundary - Change a boundary condition for the model

2707:   Input Parameters:
2708: + ds          - The PetscDS object
2709: . bd          - The boundary condition number
2710: . type        - The type of condition, e.g. DM_BC_ESSENTIAL/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
2711: . name        - The BC name
2712: . labelname   - The label defining constrained points
2713: . field       - The field to constrain
2714: . numcomps    - The number of constrained field components
2715: . comps       - An array of constrained component numbers
2716: . bcFunc      - A pointwise function giving boundary values
2717: . numids      - The number of DMLabel ids for constrained points
2718: . ids         - An array of ids for constrained points
2719: - ctx         - An optional user context for bcFunc

2721:   Note: The boundary condition number is the order in which it was registered. The user can get the number of boundary conditions from PetscDSGetNumBoundary().

2723:   Level: developer

2725: .seealso: PetscDSAddBoundary(), PetscDSGetBoundary(), PetscDSGetNumBoundary()
2726: @*/
2727: PetscErrorCode PetscDSUpdateBoundary(PetscDS ds, PetscInt bd, DMBoundaryConditionType type, const char name[], const char labelname[], PetscInt field, PetscInt numcomps, const PetscInt *comps, void (*bcFunc)(void), PetscInt numids, const PetscInt *ids, void *ctx)
2728: {
2729:   DSBoundary     b = ds->boundary;
2730:   PetscInt       n = 0;

2735:   while (b) {
2736:     if (n == bd) break;
2737:     b = b->next;
2738:     ++n;
2739:   }
2740:   if (!b) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %d is not in [0, %d)", bd, n);
2741:   if (name) {
2742:     PetscFree(b->name);
2743:     PetscStrallocpy(name, (char **) &b->name);
2744:   }
2745:   if (labelname) {
2746:     PetscFree(b->labelname);
2747:     PetscStrallocpy(labelname, (char **) &b->labelname);
2748:   }
2749:   if (numcomps >= 0 && numcomps != b->numcomps) {
2750:     b->numcomps = numcomps;
2751:     PetscFree(b->comps);
2752:     PetscMalloc1(numcomps, &b->comps);
2753:     if (numcomps) {PetscMemcpy(b->comps, comps, numcomps*sizeof(PetscInt));}
2754:   }
2755:   if (numids >= 0 && numids != b->numids) {
2756:     b->numids = numids;
2757:     PetscFree(b->ids);
2758:     PetscMalloc1(numids, &b->ids);
2759:     if (numids) {PetscMemcpy(b->ids, ids, numids*sizeof(PetscInt));}
2760:   }
2761:   b->type = type;
2762:   if (field >= 0) {b->field  = field;}
2763:   if (bcFunc)     {b->func   = bcFunc;}
2764:   if (ctx)        {b->ctx    = ctx;}
2765:   return(0);
2766: }

2768: /*@
2769:   PetscDSGetNumBoundary - Get the number of registered BC

2771:   Input Parameters:
2772: . ds - The PetscDS object

2774:   Output Parameters:
2775: . numBd - The number of BC

2777:   Level: intermediate

2779: .seealso: PetscDSAddBoundary(), PetscDSGetBoundary()
2780: @*/
2781: PetscErrorCode PetscDSGetNumBoundary(PetscDS ds, PetscInt *numBd)
2782: {
2783:   DSBoundary b = ds->boundary;

2788:   *numBd = 0;
2789:   while (b) {++(*numBd); b = b->next;}
2790:   return(0);
2791: }

2793: /*@C
2794:   PetscDSGetBoundary - Gets a boundary condition to the model

2796:   Input Parameters:
2797: + ds          - The PetscDS object
2798: - bd          - The BC number

2800:   Output Parameters:
2801: + type        - The type of condition, e.g. DM_BC_ESSENTIAL/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
2802: . name        - The BC name
2803: . labelname   - The label defining constrained points
2804: . field       - The field to constrain
2805: . numcomps    - The number of constrained field components
2806: . comps       - An array of constrained component numbers
2807: . bcFunc      - A pointwise function giving boundary values
2808: . numids      - The number of DMLabel ids for constrained points
2809: . ids         - An array of ids for constrained points
2810: - ctx         - An optional user context for bcFunc

2812:   Options Database Keys:
2813: + -bc_<boundary name> <num> - Overrides the boundary ids
2814: - -bc_<boundary name>_comp <num> - Overrides the boundary components

2816:   Level: developer

2818: .seealso: PetscDSAddBoundary()
2819: @*/
2820: PetscErrorCode PetscDSGetBoundary(PetscDS ds, PetscInt bd, DMBoundaryConditionType *type, const char **name, const char **labelname, PetscInt *field, PetscInt *numcomps, const PetscInt **comps, void (**func)(void), PetscInt *numids, const PetscInt **ids, void **ctx)
2821: {
2822:   DSBoundary b    = ds->boundary;
2823:   PetscInt   n    = 0;

2827:   while (b) {
2828:     if (n == bd) break;
2829:     b = b->next;
2830:     ++n;
2831:   }
2832:   if (!b) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %d is not in [0, %d)", bd, n);
2833:   if (type) {
2835:     *type = b->type;
2836:   }
2837:   if (name) {
2839:     *name = b->name;
2840:   }
2841:   if (labelname) {
2843:     *labelname = b->labelname;
2844:   }
2845:   if (field) {
2847:     *field = b->field;
2848:   }
2849:   if (numcomps) {
2851:     *numcomps = b->numcomps;
2852:   }
2853:   if (comps) {
2855:     *comps = b->comps;
2856:   }
2857:   if (func) {
2859:     *func = b->func;
2860:   }
2861:   if (numids) {
2863:     *numids = b->numids;
2864:   }
2865:   if (ids) {
2867:     *ids = b->ids;
2868:   }
2869:   if (ctx) {
2871:     *ctx = b->ctx;
2872:   }
2873:   return(0);
2874: }

2876: /*@
2877:   PetscDSCopyBoundary - Copy all boundary condition objects to the new problem

2879:   Not collective

2881:   Input Parameter:
2882: . prob - The PetscDS object

2884:   Output Parameter:
2885: . newprob - The PetscDS copy

2887:   Level: intermediate

2889: .seealso: PetscDSCopyEquations(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
2890: @*/
2891: PetscErrorCode PetscDSCopyBoundary(PetscDS probA, PetscDS probB)
2892: {
2893:   DSBoundary     b, next, *lastnext;

2899:   if (probA == probB) return(0);
2900:   next = probB->boundary;
2901:   while (next) {
2902:     DSBoundary b = next;

2904:     next = b->next;
2905:     PetscFree(b->comps);
2906:     PetscFree(b->ids);
2907:     PetscFree(b->name);
2908:     PetscFree(b->labelname);
2909:     PetscFree(b);
2910:   }
2911:   lastnext = &(probB->boundary);
2912:   for (b = probA->boundary; b; b = b->next) {
2913:     DSBoundary bNew;

2915:     PetscNew(&bNew);
2916:     bNew->numcomps = b->numcomps;
2917:     PetscMalloc1(bNew->numcomps, &bNew->comps);
2918:     PetscMemcpy(bNew->comps, b->comps, bNew->numcomps*sizeof(PetscInt));
2919:     bNew->numids = b->numids;
2920:     PetscMalloc1(bNew->numids, &bNew->ids);
2921:     PetscMemcpy(bNew->ids, b->ids, bNew->numids*sizeof(PetscInt));
2922:     PetscStrallocpy(b->labelname,(char **) &(bNew->labelname));
2923:     PetscStrallocpy(b->name,(char **) &(bNew->name));
2924:     bNew->ctx   = b->ctx;
2925:     bNew->type  = b->type;
2926:     bNew->field = b->field;
2927:     bNew->func  = b->func;

2929:     *lastnext = bNew;
2930:     lastnext = &(bNew->next);
2931:   }
2932:   return(0);
2933: }

2935: /*@C
2936:   PetscDSSelectEquations - Copy pointwise function pointers to the new problem with different field layout

2938:   Not collective

2940:   Input Parameter:
2941: + prob - The PetscDS object
2942: . numFields - Number of new fields
2943: - fields - Old field number for each new field

2945:   Output Parameter:
2946: . newprob - The PetscDS copy

2948:   Level: intermediate

2950: .seealso: PetscDSCopyBoundary(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
2951: @*/
2952: PetscErrorCode PetscDSSelectEquations(PetscDS prob, PetscInt numFields, const PetscInt fields[], PetscDS newprob)
2953: {
2954:   PetscInt       Nf, Nfn, fn, gn;

2961:   PetscDSGetNumFields(prob, &Nf);
2962:   PetscDSGetNumFields(newprob, &Nfn);
2963:   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);
2964:   for (fn = 0; fn < numFields; ++fn) {
2965:     const PetscInt   f = fields ? fields[fn] : fn;
2966:     PetscPointFunc   obj;
2967:     PetscPointFunc   f0, f1;
2968:     PetscBdPointFunc f0Bd, f1Bd;
2969:     PetscRiemannFunc r;

2971:     if (f >= Nf) SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_SIZ, "Field %D must be in [0, %D)", f, Nf);
2972:     PetscDSGetObjective(prob, f, &obj);
2973:     PetscDSGetResidual(prob, f, &f0, &f1);
2974:     PetscDSGetBdResidual(prob, f, &f0Bd, &f1Bd);
2975:     PetscDSGetRiemannSolver(prob, f, &r);
2976:     PetscDSSetObjective(newprob, fn, obj);
2977:     PetscDSSetResidual(newprob, fn, f0, f1);
2978:     PetscDSSetBdResidual(newprob, fn, f0Bd, f1Bd);
2979:     PetscDSSetRiemannSolver(newprob, fn, r);
2980:     for (gn = 0; gn < numFields; ++gn) {
2981:       const PetscInt  g = fields ? fields[gn] : gn;
2982:       PetscPointJac   g0, g1, g2, g3;
2983:       PetscPointJac   g0p, g1p, g2p, g3p;
2984:       PetscBdPointJac g0Bd, g1Bd, g2Bd, g3Bd;

2986:       if (g >= Nf) SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_SIZ, "Field %D must be in [0, %D)", g, Nf);
2987:       PetscDSGetJacobian(prob, f, g, &g0, &g1, &g2, &g3);
2988:       PetscDSGetJacobianPreconditioner(prob, f, g, &g0p, &g1p, &g2p, &g3p);
2989:       PetscDSGetBdJacobian(prob, f, g, &g0Bd, &g1Bd, &g2Bd, &g3Bd);
2990:       PetscDSSetJacobian(newprob, fn, gn, g0, g1, g2, g3);
2991:       PetscDSSetJacobianPreconditioner(prob, fn, gn, g0p, g1p, g2p, g3p);
2992:       PetscDSSetBdJacobian(newprob, fn, gn, g0Bd, g1Bd, g2Bd, g3Bd);
2993:     }
2994:   }
2995:   return(0);
2996: }

2998: /*@
2999:   PetscDSCopyEquations - Copy all pointwise function pointers to the new problem

3001:   Not collective

3003:   Input Parameter:
3004: . prob - The PetscDS object

3006:   Output Parameter:
3007: . newprob - The PetscDS copy

3009:   Level: intermediate

3011: .seealso: PetscDSCopyBoundary(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3012: @*/
3013: PetscErrorCode PetscDSCopyEquations(PetscDS prob, PetscDS newprob)
3014: {
3015:   PetscInt       Nf, Ng;

3021:   PetscDSGetNumFields(prob, &Nf);
3022:   PetscDSGetNumFields(newprob, &Ng);
3023:   if (Nf != Ng) SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_SIZ, "Number of fields must match %D != %D", Nf, Ng);
3024:   PetscDSSelectEquations(prob, Nf, NULL, newprob);
3025:   return(0);
3026: }
3027: /*@
3028:   PetscDSCopyConstants - Copy all constants to the new problem

3030:   Not collective

3032:   Input Parameter:
3033: . prob - The PetscDS object

3035:   Output Parameter:
3036: . newprob - The PetscDS copy

3038:   Level: intermediate

3040: .seealso: PetscDSCopyBoundary(), PetscDSCopyEquations(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3041: @*/
3042: PetscErrorCode PetscDSCopyConstants(PetscDS prob, PetscDS newprob)
3043: {
3044:   PetscInt           Nc;
3045:   const PetscScalar *constants;
3046:   PetscErrorCode     ierr;

3051:   PetscDSGetConstants(prob, &Nc, &constants);
3052:   PetscDSSetConstants(newprob, Nc, (PetscScalar *) constants);
3053:   return(0);
3054: }

3056: PetscErrorCode PetscDSGetHeightSubspace(PetscDS prob, PetscInt height, PetscDS *subprob)
3057: {
3058:   PetscInt       dim, Nf, f;

3064:   if (height == 0) {*subprob = prob; return(0);}
3065:   PetscDSGetNumFields(prob, &Nf);
3066:   PetscDSGetSpatialDimension(prob, &dim);
3067:   if (height > dim) SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_OUTOFRANGE, "DS can only handle height in [0, %D], not %D", dim, height);
3068:   if (!prob->subprobs) {PetscCalloc1(dim, &prob->subprobs);}
3069:   if (!prob->subprobs[height-1]) {
3070:     PetscInt cdim;

3072:     PetscDSCreate(PetscObjectComm((PetscObject) prob), &prob->subprobs[height-1]);
3073:     PetscDSGetCoordinateDimension(prob, &cdim);
3074:     PetscDSSetCoordinateDimension(prob->subprobs[height-1], cdim);
3075:     for (f = 0; f < Nf; ++f) {
3076:       PetscFE      subfe;
3077:       PetscObject  obj;
3078:       PetscClassId id;

3080:       PetscDSGetDiscretization(prob, f, &obj);
3081:       PetscObjectGetClassId(obj, &id);
3082:       if (id == PETSCFE_CLASSID) {PetscFEGetHeightSubspace((PetscFE) obj, height, &subfe);}
3083:       else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unsupported discretization type for field %d", f);
3084:       PetscDSSetDiscretization(prob->subprobs[height-1], f, (PetscObject) subfe);
3085:     }
3086:   }
3087:   *subprob = prob->subprobs[height-1];
3088:   return(0);
3089: }

3091: static PetscErrorCode PetscDSDestroy_Basic(PetscDS prob)
3092: {
3093:   PetscErrorCode      ierr;

3096:   PetscFree(prob->data);
3097:   return(0);
3098: }

3100: static PetscErrorCode PetscDSInitialize_Basic(PetscDS prob)
3101: {
3103:   prob->ops->setfromoptions = NULL;
3104:   prob->ops->setup          = NULL;
3105:   prob->ops->view           = NULL;
3106:   prob->ops->destroy        = PetscDSDestroy_Basic;
3107:   return(0);
3108: }

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

3113:   Level: intermediate

3115: .seealso: PetscDSType, PetscDSCreate(), PetscDSSetType()
3116: M*/

3118: PETSC_EXTERN PetscErrorCode PetscDSCreate_Basic(PetscDS prob)
3119: {
3120:   PetscDS_Basic *b;
3121:   PetscErrorCode      ierr;

3125:   PetscNewLog(prob, &b);
3126:   prob->data = b;

3128:   PetscDSInitialize_Basic(prob);
3129:   return(0);
3130: }