Actual source code: dtds.c

petsc-master 2019-11-13
Report Typos and Errors
  1:  #include <petsc/private/petscdsimpl.h>

  3: PetscClassId PETSCDS_CLASSID = 0;

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

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

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

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

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

 24:   Not Collective

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

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

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

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

 48:   Level: advanced

 50:    Not available from Fortran

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

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

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

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

 67:   Collective on prob

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

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

 76:   Level: intermediate

 78:    Not available from Fortran

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

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

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

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

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

109:   Not Collective

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

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

117:   Level: intermediate

119:    Not available from Fortran

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

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

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

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

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

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

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

226: /*@C
227:    PetscDSViewFromOptions - View from Options

229:    Collective on PetscDS

231:    Input Parameters:
232: +  A - the PetscDS object
233: .  obj - Optional object
234: -  name - command line option

236:    Level: intermediate
237: .seealso:  PetscDS, PetscDSView, PetscObjectViewFromOptions(), PetscDSCreate()
238: @*/
239: PetscErrorCode  PetscDSViewFromOptions(PetscDS A,PetscObject obj,const char name[])
240: {

245:   PetscObjectViewFromOptions((PetscObject)A,obj,name);
246:   return(0);
247: }

249: /*@C
250:   PetscDSView - Views a PetscDS

252:   Collective on prob

254:   Input Parameter:
255: + prob - the PetscDS object to view
256: - v  - the viewer

258:   Level: developer

260: .seealso PetscDSDestroy()
261: @*/
262: PetscErrorCode PetscDSView(PetscDS prob, PetscViewer v)
263: {
264:   PetscBool      iascii;

269:   if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) prob), &v);}
271:   PetscObjectTypeCompare((PetscObject) v, PETSCVIEWERASCII, &iascii);
272:   if (iascii) {PetscDSView_Ascii(prob, v);}
273:   if (prob->ops->view) {(*prob->ops->view)(prob, v);}
274:   return(0);
275: }

277: /*@
278:   PetscDSSetFromOptions - sets parameters in a PetscDS from the options database

280:   Collective on prob

282:   Input Parameter:
283: . prob - the PetscDS object to set options for

285:   Options Database:
286: + -petscds_type <type>     : Set the DS type
287: . -petscds_view <view opt> : View the DS
288: . -petscds_jac_pre         : Turn formation of a separate Jacobian preconditioner on and off
289: . -bc_<name> <ids>         : Specify a list of label ids for a boundary condition
290: - -bc_<name>_comp <comps>  : Specify a list of field components to constrain for a boundary condition

292:   Level: developer

294: .seealso PetscDSView()
295: @*/
296: PetscErrorCode PetscDSSetFromOptions(PetscDS prob)
297: {
298:   DSBoundary     b;
299:   const char    *defaultType;
300:   char           name[256];
301:   PetscBool      flg;

306:   if (!((PetscObject) prob)->type_name) {
307:     defaultType = PETSCDSBASIC;
308:   } else {
309:     defaultType = ((PetscObject) prob)->type_name;
310:   }
311:   PetscDSRegisterAll();

313:   PetscObjectOptionsBegin((PetscObject) prob);
314:   for (b = prob->boundary; b; b = b->next) {
315:     char       optname[1024];
316:     PetscInt   ids[1024], len = 1024;
317:     PetscBool  flg;

319:     PetscSNPrintf(optname, sizeof(optname), "-bc_%s", b->name);
320:     PetscMemzero(ids, sizeof(ids));
321:     PetscOptionsIntArray(optname, "List of boundary IDs", "", ids, &len, &flg);
322:     if (flg) {
323:       b->numids = len;
324:       PetscFree(b->ids);
325:       PetscMalloc1(len, &b->ids);
326:       PetscArraycpy(b->ids, ids, len);
327:     }
328:     len = 1024;
329:     PetscSNPrintf(optname, sizeof(optname), "-bc_%s_comp", b->name);
330:     PetscMemzero(ids, sizeof(ids));
331:     PetscOptionsIntArray(optname, "List of boundary field components", "", ids, &len, &flg);
332:     if (flg) {
333:       b->numcomps = len;
334:       PetscFree(b->comps);
335:       PetscMalloc1(len, &b->comps);
336:       PetscArraycpy(b->comps, ids, len);
337:     }
338:   }
339:   PetscOptionsFList("-petscds_type", "Discrete System", "PetscDSSetType", PetscDSList, defaultType, name, 256, &flg);
340:   if (flg) {
341:     PetscDSSetType(prob, name);
342:   } else if (!((PetscObject) prob)->type_name) {
343:     PetscDSSetType(prob, defaultType);
344:   }
345:   PetscOptionsBool("-petscds_jac_pre", "Discrete System", "PetscDSUseJacobianPreconditioner", prob->useJacPre, &prob->useJacPre, &flg);
346:   if (prob->ops->setfromoptions) {(*prob->ops->setfromoptions)(prob);}
347:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
348:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) prob);
349:   PetscOptionsEnd();
350:   if (prob->Nf) {PetscDSViewFromOptions(prob, NULL, "-petscds_view");}
351:   return(0);
352: }

354: /*@C
355:   PetscDSSetUp - Construct data structures for the PetscDS

357:   Collective on prob

359:   Input Parameter:
360: . prob - the PetscDS object to setup

362:   Level: developer

364: .seealso PetscDSView(), PetscDSDestroy()
365: @*/
366: PetscErrorCode PetscDSSetUp(PetscDS prob)
367: {
368:   const PetscInt Nf = prob->Nf;
369:   PetscInt       dim, dimEmbed, NbMax = 0, NcMax = 0, NqMax = 0, NsMax = 1, f;

374:   if (prob->setup) return(0);
375:   /* Calculate sizes */
376:   PetscDSGetSpatialDimension(prob, &dim);
377:   PetscDSGetCoordinateDimension(prob, &dimEmbed);
378:   prob->totDim = prob->totComp = 0;
379:   PetscMalloc2(Nf,&prob->Nc,Nf,&prob->Nb);
380:   PetscCalloc2(Nf+1,&prob->off,Nf+1,&prob->offDer);
381:   PetscMalloc4(Nf,&prob->basis,Nf,&prob->basisDer,Nf,&prob->basisFace,Nf,&prob->basisDerFace);
382:   for (f = 0; f < Nf; ++f) {
383:     PetscObject     obj;
384:     PetscClassId    id;
385:     PetscQuadrature q;
386:     PetscInt        Nq = 0, Nb, Nc;

388:     PetscDSGetDiscretization(prob, f, &obj);
389:     PetscObjectGetClassId(obj, &id);
390:     if (id == PETSCFE_CLASSID)      {
391:       PetscFE fe = (PetscFE) obj;

393:       PetscFEGetQuadrature(fe, &q);
394:       PetscFEGetDimension(fe, &Nb);
395:       PetscFEGetNumComponents(fe, &Nc);
396:       PetscFEGetDefaultTabulation(fe, &prob->basis[f], &prob->basisDer[f], NULL);
397:       PetscFEGetFaceTabulation(fe, &prob->basisFace[f], &prob->basisDerFace[f], NULL);
398:     } else if (id == PETSCFV_CLASSID) {
399:       PetscFV fv = (PetscFV) obj;

401:       PetscFVGetQuadrature(fv, &q);
402:       PetscFVGetNumComponents(fv, &Nc);
403:       Nb   = Nc;
404:       PetscFVGetDefaultTabulation(fv, &prob->basis[f], &prob->basisDer[f], NULL);
405:       /* TODO: should PetscFV also have face tabulation? Otherwise there will be a null pointer in prob->basisFace */
406:     } else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
407:     prob->Nc[f]       = Nc;
408:     prob->Nb[f]       = Nb;
409:     prob->off[f+1]    = Nc     + prob->off[f];
410:     prob->offDer[f+1] = Nc*dim + prob->offDer[f];
411:     if (q) {PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);}
412:     NqMax          = PetscMax(NqMax, Nq);
413:     NbMax          = PetscMax(NbMax, Nb);
414:     NcMax          = PetscMax(NcMax, Nc);
415:     prob->totDim  += Nb;
416:     prob->totComp += Nc;
417:     /* There are two faces for all fields but the cohesive field on a hybrid cell */
418:     if (prob->isHybrid && (f < Nf-1)) prob->totDim += Nb;
419:   }
420:   /* Allocate works space */
421:   if (prob->isHybrid) NsMax = 2;
422:   PetscMalloc3(NsMax*prob->totComp,&prob->u,NsMax*prob->totComp,&prob->u_t,NsMax*prob->totComp*dimEmbed,&prob->u_x);
423:   PetscMalloc5(dimEmbed,&prob->x,NbMax*NcMax,&prob->basisReal,NbMax*NcMax*dimEmbed,&prob->basisDerReal,NbMax*NcMax,&prob->testReal,NbMax*NcMax*dimEmbed,&prob->testDerReal);
424:   PetscMalloc6(NsMax*NqMax*NcMax,&prob->f0,NsMax*NqMax*NcMax*dim,&prob->f1,
425:                       NsMax*NsMax*NqMax*NcMax*NcMax,&prob->g0,NsMax*NsMax*NqMax*NcMax*NcMax*dim,&prob->g1,
426:                       NsMax*NsMax*NqMax*NcMax*NcMax*dim,&prob->g2,NsMax*NsMax*NqMax*NcMax*NcMax*dim*dim,&prob->g3);
427:   if (prob->ops->setup) {(*prob->ops->setup)(prob);}
428:   prob->setup = PETSC_TRUE;
429:   return(0);
430: }

432: static PetscErrorCode PetscDSDestroyStructs_Static(PetscDS prob)
433: {

437:   PetscFree2(prob->Nc,prob->Nb);
438:   PetscFree2(prob->off,prob->offDer);
439:   PetscFree4(prob->basis,prob->basisDer,prob->basisFace,prob->basisDerFace);
440:   PetscFree3(prob->u,prob->u_t,prob->u_x);
441:   PetscFree5(prob->x,prob->basisReal, prob->basisDerReal,prob->testReal,prob->testDerReal);
442:   PetscFree6(prob->f0,prob->f1,prob->g0,prob->g1,prob->g2,prob->g3);
443:   return(0);
444: }

446: static PetscErrorCode PetscDSEnlarge_Static(PetscDS prob, PetscInt NfNew)
447: {
448:   PetscObject      *tmpd;
449:   PetscBool        *tmpi;
450:   PetscPointFunc   *tmpobj, *tmpf, *tmpup;
451:   PetscPointJac    *tmpg, *tmpgp, *tmpgt;
452:   PetscBdPointFunc *tmpfbd;
453:   PetscBdPointJac  *tmpgbd;
454:   PetscRiemannFunc *tmpr;
455:   PetscSimplePointFunc *tmpexactSol;
456:   void                **tmpexactCtx;
457:   void            **tmpctx;
458:   PetscInt          Nf = prob->Nf, f;
459:   PetscErrorCode    ierr;

462:   if (Nf >= NfNew) return(0);
463:   prob->setup = PETSC_FALSE;
464:   PetscDSDestroyStructs_Static(prob);
465:   PetscMalloc2(NfNew, &tmpd, NfNew, &tmpi);
466:   for (f = 0; f < Nf; ++f) {tmpd[f] = prob->disc[f]; tmpi[f] = prob->implicit[f];}
467:   for (f = Nf; f < NfNew; ++f) {tmpd[f] = NULL; tmpi[f] = PETSC_TRUE;}
468:   PetscFree2(prob->disc, prob->implicit);
469:   prob->Nf        = NfNew;
470:   prob->disc      = tmpd;
471:   prob->implicit  = tmpi;
472:   PetscCalloc7(NfNew, &tmpobj, NfNew*2, &tmpf, NfNew*NfNew*4, &tmpg, NfNew*NfNew*4, &tmpgp, NfNew*NfNew*4, &tmpgt, NfNew, &tmpr, NfNew, &tmpctx);
473:   PetscCalloc1(NfNew, &tmpup);
474:   for (f = 0; f < Nf; ++f) tmpobj[f] = prob->obj[f];
475:   for (f = 0; f < Nf*2; ++f) tmpf[f] = prob->f[f];
476:   for (f = 0; f < Nf*Nf*4; ++f) tmpg[f] = prob->g[f];
477:   for (f = 0; f < Nf*Nf*4; ++f) tmpgp[f] = prob->gp[f];
478:   for (f = 0; f < Nf; ++f) tmpr[f] = prob->r[f];
479:   for (f = 0; f < Nf; ++f) tmpup[f] = prob->update[f];
480:   for (f = 0; f < Nf; ++f) tmpctx[f] = prob->ctx[f];
481:   for (f = Nf; f < NfNew; ++f) tmpobj[f] = NULL;
482:   for (f = Nf*2; f < NfNew*2; ++f) tmpf[f] = NULL;
483:   for (f = Nf*Nf*4; f < NfNew*NfNew*4; ++f) tmpg[f] = NULL;
484:   for (f = Nf*Nf*4; f < NfNew*NfNew*4; ++f) tmpgp[f] = NULL;
485:   for (f = Nf*Nf*4; f < NfNew*NfNew*4; ++f) tmpgt[f] = NULL;
486:   for (f = Nf; f < NfNew; ++f) tmpr[f] = NULL;
487:   for (f = Nf; f < NfNew; ++f) tmpup[f] = NULL;
488:   for (f = Nf; f < NfNew; ++f) tmpctx[f] = NULL;
489:   PetscFree7(prob->obj, prob->f, prob->g, prob->gp, prob->gt, prob->r, prob->ctx);
490:   PetscFree(prob->update);
491:   prob->obj = tmpobj;
492:   prob->f   = tmpf;
493:   prob->g   = tmpg;
494:   prob->gp  = tmpgp;
495:   prob->gt  = tmpgt;
496:   prob->r   = tmpr;
497:   prob->update = tmpup;
498:   prob->ctx = tmpctx;
499:   PetscCalloc4(NfNew*2, &tmpfbd, NfNew*NfNew*4, &tmpgbd, NfNew, &tmpexactSol, NfNew, &tmpexactCtx);
500:   for (f = 0; f < Nf*2; ++f) tmpfbd[f] = prob->fBd[f];
501:   for (f = 0; f < Nf*Nf*4; ++f) tmpgbd[f] = prob->gBd[f];
502:   for (f = 0; f < Nf; ++f) tmpexactSol[f] = prob->exactSol[f];
503:   for (f = 0; f < Nf; ++f) tmpexactCtx[f] = prob->exactCtx[f];
504:   for (f = Nf*2; f < NfNew*2; ++f) tmpfbd[f] = NULL;
505:   for (f = Nf*Nf*4; f < NfNew*NfNew*4; ++f) tmpgbd[f] = NULL;
506:   for (f = Nf; f < NfNew; ++f) tmpexactSol[f] = NULL;
507:   for (f = Nf; f < NfNew; ++f) tmpexactCtx[f] = NULL;
508:   PetscFree4(prob->fBd, prob->gBd, prob->exactSol, prob->exactCtx);
509:   prob->fBd = tmpfbd;
510:   prob->gBd = tmpgbd;
511:   prob->exactSol = tmpexactSol;
512:   prob->exactCtx = tmpexactCtx;
513:   return(0);
514: }

516: /*@
517:   PetscDSDestroy - Destroys a PetscDS object

519:   Collective on prob

521:   Input Parameter:
522: . prob - the PetscDS object to destroy

524:   Level: developer

526: .seealso PetscDSView()
527: @*/
528: PetscErrorCode PetscDSDestroy(PetscDS *prob)
529: {
530:   PetscInt       f;
531:   DSBoundary     next;

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

538:   if (--((PetscObject)(*prob))->refct > 0) {*prob = 0; return(0);}
539:   ((PetscObject) (*prob))->refct = 0;
540:   if ((*prob)->subprobs) {
541:     PetscInt dim, d;

543:     PetscDSGetSpatialDimension(*prob, &dim);
544:     for (d = 0; d < dim; ++d) {PetscDSDestroy(&(*prob)->subprobs[d]);}
545:   }
546:   PetscFree((*prob)->subprobs);
547:   PetscDSDestroyStructs_Static(*prob);
548:   for (f = 0; f < (*prob)->Nf; ++f) {
549:     PetscObjectDereference((*prob)->disc[f]);
550:   }
551:   PetscFree2((*prob)->disc, (*prob)->implicit);
552:   PetscFree7((*prob)->obj,(*prob)->f,(*prob)->g,(*prob)->gp,(*prob)->gt,(*prob)->r,(*prob)->ctx);
553:   PetscFree((*prob)->update);
554:   PetscFree4((*prob)->fBd,(*prob)->gBd,(*prob)->exactSol,(*prob)->exactCtx);
555:   if ((*prob)->ops->destroy) {(*(*prob)->ops->destroy)(*prob);}
556:   next = (*prob)->boundary;
557:   while (next) {
558:     DSBoundary b = next;

560:     next = b->next;
561:     PetscFree(b->comps);
562:     PetscFree(b->ids);
563:     PetscFree(b->name);
564:     PetscFree(b->labelname);
565:     PetscFree(b);
566:   }
567:   PetscFree((*prob)->constants);
568:   PetscHeaderDestroy(prob);
569:   return(0);
570: }

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

575:   Collective

577:   Input Parameter:
578: . comm - The communicator for the PetscDS object

580:   Output Parameter:
581: . prob - The PetscDS object

583:   Level: beginner

585: .seealso: PetscDSSetType(), PETSCDSBASIC
586: @*/
587: PetscErrorCode PetscDSCreate(MPI_Comm comm, PetscDS *prob)
588: {
589:   PetscDS   p;

594:   *prob  = NULL;
595:   PetscDSInitializePackage();

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

599:   p->Nf           = 0;
600:   p->setup        = PETSC_FALSE;
601:   p->numConstants = 0;
602:   p->constants    = NULL;
603:   p->dimEmbed     = -1;
604:   p->useJacPre    = PETSC_TRUE;

606:   *prob = p;
607:   return(0);
608: }

610: /*@
611:   PetscDSGetNumFields - Returns the number of fields in the DS

613:   Not collective

615:   Input Parameter:
616: . prob - The PetscDS object

618:   Output Parameter:
619: . Nf - The number of fields

621:   Level: beginner

623: .seealso: PetscDSGetSpatialDimension(), PetscDSCreate()
624: @*/
625: PetscErrorCode PetscDSGetNumFields(PetscDS prob, PetscInt *Nf)
626: {
630:   *Nf = prob->Nf;
631:   return(0);
632: }

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

637:   Not collective

639:   Input Parameter:
640: . prob - The PetscDS object

642:   Output Parameter:
643: . dim - The spatial dimension

645:   Level: beginner

647: .seealso: PetscDSGetCoordinateDimension(), PetscDSGetNumFields(), PetscDSCreate()
648: @*/
649: PetscErrorCode PetscDSGetSpatialDimension(PetscDS prob, PetscInt *dim)
650: {

656:   *dim = 0;
657:   if (prob->Nf) {
658:     PetscObject  obj;
659:     PetscClassId id;

661:     PetscDSGetDiscretization(prob, 0, &obj);
662:     PetscObjectGetClassId(obj, &id);
663:     if (id == PETSCFE_CLASSID)      {PetscFEGetSpatialDimension((PetscFE) obj, dim);}
664:     else if (id == PETSCFV_CLASSID) {PetscFVGetSpatialDimension((PetscFV) obj, dim);}
665:     else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", 0);
666:   }
667:   return(0);
668: }

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

673:   Not collective

675:   Input Parameter:
676: . prob - The PetscDS object

678:   Output Parameter:
679: . dimEmbed - The coordinate dimension

681:   Level: beginner

683: .seealso: PetscDSSetCoordinateDimension(), PetscDSGetSpatialDimension(), PetscDSGetNumFields(), PetscDSCreate()
684: @*/
685: PetscErrorCode PetscDSGetCoordinateDimension(PetscDS prob, PetscInt *dimEmbed)
686: {
690:   if (prob->dimEmbed < 0) SETERRQ(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONGSTATE, "No coordinate dimension set for this DS");
691:   *dimEmbed = prob->dimEmbed;
692:   return(0);
693: }

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

698:   Logically collective on prob

700:   Input Parameters:
701: + prob - The PetscDS object
702: - dimEmbed - The coordinate dimension

704:   Level: beginner

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

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

720:   Not collective

722:   Input Parameter:
723: . prob - The PetscDS object

725:   Output Parameter:
726: . isHybrid - The flag

728:   Level: developer

730: .seealso: PetscDSSetHybrid(), PetscDSCreate()
731: @*/
732: PetscErrorCode PetscDSGetHybrid(PetscDS prob, PetscBool *isHybrid)
733: {
737:   *isHybrid = prob->isHybrid;
738:   return(0);
739: }

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

744:   Not collective

746:   Input Parameters:
747: + prob - The PetscDS object
748: - isHybrid - The flag

750:   Level: developer

752: .seealso: PetscDSGetHybrid(), PetscDSCreate()
753: @*/
754: PetscErrorCode PetscDSSetHybrid(PetscDS prob, PetscBool isHybrid)
755: {
758:   prob->isHybrid = isHybrid;
759:   return(0);
760: }

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

765:   Not collective

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

770:   Output Parameter:
771: . dim - The total problem dimension

773:   Level: beginner

775: .seealso: PetscDSGetNumFields(), PetscDSCreate()
776: @*/
777: PetscErrorCode PetscDSGetTotalDimension(PetscDS prob, PetscInt *dim)
778: {

783:   PetscDSSetUp(prob);
785:   *dim = prob->totDim;
786:   return(0);
787: }

789: /*@
790:   PetscDSGetTotalComponents - Returns the total number of components in this system

792:   Not collective

794:   Input Parameter:
795: . prob - The PetscDS object

797:   Output Parameter:
798: . dim - The total number of components

800:   Level: beginner

802: .seealso: PetscDSGetNumFields(), PetscDSCreate()
803: @*/
804: PetscErrorCode PetscDSGetTotalComponents(PetscDS prob, PetscInt *Nc)
805: {

810:   PetscDSSetUp(prob);
812:   *Nc = prob->totComp;
813:   return(0);
814: }

816: /*@
817:   PetscDSGetDiscretization - Returns the discretization object for the given field

819:   Not collective

821:   Input Parameters:
822: + prob - The PetscDS object
823: - f - The field number

825:   Output Parameter:
826: . disc - The discretization object

828:   Level: beginner

830: .seealso: PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
831: @*/
832: PetscErrorCode PetscDSGetDiscretization(PetscDS prob, PetscInt f, PetscObject *disc)
833: {
837:   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);
838:   *disc = prob->disc[f];
839:   return(0);
840: }

842: /*@
843:   PetscDSSetDiscretization - Sets the discretization object for the given field

845:   Not collective

847:   Input Parameters:
848: + prob - The PetscDS object
849: . f - The field number
850: - disc - The discretization object

852:   Level: beginner

854: .seealso: PetscDSGetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
855: @*/
856: PetscErrorCode PetscDSSetDiscretization(PetscDS prob, PetscInt f, PetscObject disc)
857: {

863:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
864:   PetscDSEnlarge_Static(prob, f+1);
865:   PetscObjectDereference(prob->disc[f]);
866:   prob->disc[f] = disc;
867:   PetscObjectReference(disc);
868:   {
869:     PetscClassId id;

871:     PetscObjectGetClassId(disc, &id);
872:     if (id == PETSCFE_CLASSID) {
873:       PetscDSSetImplicit(prob, f, PETSC_TRUE);
874:     } else if (id == PETSCFV_CLASSID) {
875:       PetscDSSetImplicit(prob, f, PETSC_FALSE);
876:     }
877:   }
878:   return(0);
879: }

881: /*@
882:   PetscDSAddDiscretization - Adds a discretization object

884:   Not collective

886:   Input Parameters:
887: + prob - The PetscDS object
888: - disc - The boundary discretization object

890:   Level: beginner

892: .seealso: PetscDSGetDiscretization(), PetscDSSetDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
893: @*/
894: PetscErrorCode PetscDSAddDiscretization(PetscDS prob, PetscObject disc)
895: {

899:   PetscDSSetDiscretization(prob, prob->Nf, disc);
900:   return(0);
901: }

903: /*@
904:   PetscDSGetImplicit - Returns 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

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

915:   Level: developer

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

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

932:   Not collective

934:   Input Parameters:
935: + prob - The PetscDS object
936: . f - The field number
937: - implicit - The flag indicating what kind of solve to use for this field

939:   Level: developer

941: .seealso: PetscDSGetImplicit(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
942: @*/
943: PetscErrorCode PetscDSSetImplicit(PetscDS prob, PetscInt f, PetscBool implicit)
944: {
947:   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);
948:   prob->implicit[f] = implicit;
949:   return(0);
950: }

952: PetscErrorCode PetscDSGetObjective(PetscDS prob, PetscInt f,
953:                                    void (**obj)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
954:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
955:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
956:                                                 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[]))
957: {
961:   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);
962:   *obj = prob->obj[f];
963:   return(0);
964: }

966: PetscErrorCode PetscDSSetObjective(PetscDS prob, PetscInt f,
967:                                    void (*obj)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
968:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
969:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
970:                                                PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[]))
971: {

977:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
978:   PetscDSEnlarge_Static(prob, f+1);
979:   prob->obj[f] = obj;
980:   return(0);
981: }

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

986:   Not collective

988:   Input Parameters:
989: + prob - The PetscDS
990: - f    - The test field number

992:   Output Parameters:
993: + f0 - integrand for the test function term
994: - f1 - integrand for the test function gradient term

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

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

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

1002: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1003: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1004: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1005: $    PetscReal t, const PetscReal x[], PetscScalar f0[])

1007: + dim - the spatial dimension
1008: . Nf - the number of fields
1009: . uOff - the offset into u[] and u_t[] for each field
1010: . uOff_x - the offset into u_x[] for each field
1011: . u - each field evaluated at the current point
1012: . u_t - the time derivative of each field evaluated at the current point
1013: . u_x - the gradient of each field evaluated at the current point
1014: . aOff - the offset into a[] and a_t[] for each auxiliary field
1015: . aOff_x - the offset into a_x[] for each auxiliary field
1016: . a - each auxiliary field evaluated at the current point
1017: . a_t - the time derivative of each auxiliary field evaluated at the current point
1018: . a_x - the gradient of auxiliary each field evaluated at the current point
1019: . t - current time
1020: . x - coordinates of the current point
1021: . numConstants - number of constant parameters
1022: . constants - constant parameters
1023: - f0 - output values at the current point

1025:   Level: intermediate

1027: .seealso: PetscDSSetResidual()
1028: @*/
1029: PetscErrorCode PetscDSGetResidual(PetscDS prob, PetscInt f,
1030:                                   void (**f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1031:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1032:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1033:                                               PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]),
1034:                                   void (**f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1035:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1036:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1037:                                               PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
1038: {
1041:   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);
1044:   return(0);
1045: }

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

1050:   Not collective

1052:   Input Parameters:
1053: + prob - The PetscDS
1054: . f    - The test field number
1055: . f0 - integrand for the test function term
1056: - f1 - integrand for the test function gradient term

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

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

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

1064: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1065: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1066: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1067: $    PetscReal t, const PetscReal x[], PetscScalar f0[])

1069: + dim - the spatial dimension
1070: . Nf - the number of fields
1071: . uOff - the offset into u[] and u_t[] for each field
1072: . uOff_x - the offset into u_x[] for each field
1073: . u - each field evaluated at the current point
1074: . u_t - the time derivative of each field evaluated at the current point
1075: . u_x - the gradient of each field evaluated at the current point
1076: . aOff - the offset into a[] and a_t[] for each auxiliary field
1077: . aOff_x - the offset into a_x[] for each auxiliary field
1078: . a - each auxiliary field evaluated at the current point
1079: . a_t - the time derivative of each auxiliary field evaluated at the current point
1080: . a_x - the gradient of auxiliary each field evaluated at the current point
1081: . t - current time
1082: . x - coordinates of the current point
1083: . numConstants - number of constant parameters
1084: . constants - constant parameters
1085: - f0 - output values at the current point

1087:   Level: intermediate

1089: .seealso: PetscDSGetResidual()
1090: @*/
1091: PetscErrorCode PetscDSSetResidual(PetscDS prob, PetscInt f,
1092:                                   void (*f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1093:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1094:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1095:                                              PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]),
1096:                                   void (*f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1097:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1098:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1099:                                              PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
1100: {

1107:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1108:   PetscDSEnlarge_Static(prob, f+1);
1109:   prob->f[f*2+0] = f0;
1110:   prob->f[f*2+1] = f1;
1111:   return(0);
1112: }

1114: /*@C
1115:   PetscDSHasJacobian - Signals that Jacobian functions have been set

1117:   Not collective

1119:   Input Parameter:
1120: . prob - The PetscDS

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

1125:   Level: intermediate

1127: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1128: @*/
1129: PetscErrorCode PetscDSHasJacobian(PetscDS prob, PetscBool *hasJac)
1130: {
1131:   PetscInt f, g, h;

1135:   *hasJac = PETSC_FALSE;
1136:   for (f = 0; f < prob->Nf; ++f) {
1137:     for (g = 0; g < prob->Nf; ++g) {
1138:       for (h = 0; h < 4; ++h) {
1139:         if (prob->g[(f*prob->Nf + g)*4+h]) *hasJac = PETSC_TRUE;
1140:       }
1141:     }
1142:   }
1143:   return(0);
1144: }

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

1149:   Not collective

1151:   Input Parameters:
1152: + prob - The PetscDS
1153: . f    - The test field number
1154: - g    - The field number

1156:   Output Parameters:
1157: + g0 - integrand for the test and basis function term
1158: . g1 - integrand for the test function and basis function gradient term
1159: . g2 - integrand for the test function gradient and basis function term
1160: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

1168: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1169: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1170: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1171: $    PetscReal t, const PetscReal u_tShift, const PetscReal x[], PetscScalar g0[])

1173: + dim - the spatial dimension
1174: . Nf - the number of fields
1175: . uOff - the offset into u[] and u_t[] for each field
1176: . uOff_x - the offset into u_x[] for each field
1177: . u - each field evaluated at the current point
1178: . u_t - the time derivative of each field evaluated at the current point
1179: . u_x - the gradient of each field evaluated at the current point
1180: . aOff - the offset into a[] and a_t[] for each auxiliary field
1181: . aOff_x - the offset into a_x[] for each auxiliary field
1182: . a - each auxiliary field evaluated at the current point
1183: . a_t - the time derivative of each auxiliary field evaluated at the current point
1184: . a_x - the gradient of auxiliary each field evaluated at the current point
1185: . t - current time
1186: . u_tShift - the multiplier a for dF/dU_t
1187: . x - coordinates of the current point
1188: . numConstants - number of constant parameters
1189: . constants - constant parameters
1190: - g0 - output values at the current point

1192:   Level: intermediate

1194: .seealso: PetscDSSetJacobian()
1195: @*/
1196: PetscErrorCode PetscDSGetJacobian(PetscDS prob, PetscInt f, PetscInt g,
1197:                                   void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1198:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1199:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1200:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1201:                                   void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1202:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1203:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1204:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1205:                                   void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1206:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1207:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1208:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1209:                                   void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1210:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1211:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1212:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1213: {
1216:   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);
1217:   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);
1222:   return(0);
1223: }

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

1228:   Not collective

1230:   Input Parameters:
1231: + prob - The PetscDS
1232: . f    - The test field number
1233: . g    - The field number
1234: . g0 - integrand for the test and basis function term
1235: . g1 - integrand for the test function and basis function gradient term
1236: . g2 - integrand for the test function gradient and basis function term
1237: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

1245: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1246: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1247: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1248: $    PetscReal t, const PetscReal x[], PetscScalar g0[])

1250: + dim - the spatial dimension
1251: . Nf - the number of fields
1252: . uOff - the offset into u[] and u_t[] for each field
1253: . uOff_x - the offset into u_x[] for each field
1254: . u - each field evaluated at the current point
1255: . u_t - the time derivative of each field evaluated at the current point
1256: . u_x - the gradient of each field evaluated at the current point
1257: . aOff - the offset into a[] and a_t[] for each auxiliary field
1258: . aOff_x - the offset into a_x[] for each auxiliary field
1259: . a - each auxiliary field evaluated at the current point
1260: . a_t - the time derivative of each auxiliary field evaluated at the current point
1261: . a_x - the gradient of auxiliary each field evaluated at the current point
1262: . t - current time
1263: . u_tShift - the multiplier a for dF/dU_t
1264: . x - coordinates of the current point
1265: . numConstants - number of constant parameters
1266: . constants - constant parameters
1267: - g0 - output values at the current point

1269:   Level: intermediate

1271: .seealso: PetscDSGetJacobian()
1272: @*/
1273: PetscErrorCode PetscDSSetJacobian(PetscDS prob, PetscInt f, PetscInt g,
1274:                                   void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1275:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1276:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1277:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1278:                                   void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1279:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1280:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1281:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1282:                                   void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1283:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1284:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1285:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1286:                                   void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1287:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1288:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1289:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1290: {

1299:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1300:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
1301:   PetscDSEnlarge_Static(prob, PetscMax(f, g)+1);
1302:   prob->g[(f*prob->Nf + g)*4+0] = g0;
1303:   prob->g[(f*prob->Nf + g)*4+1] = g1;
1304:   prob->g[(f*prob->Nf + g)*4+2] = g2;
1305:   prob->g[(f*prob->Nf + g)*4+3] = g3;
1306:   return(0);
1307: }

1309: /*@C
1310:   PetscDSUseJacobianPreconditioner - Whether to construct a Jacobian preconditioner

1312:   Not collective

1314:   Input Parameters:
1315: + prob - The PetscDS
1316: - useJacPre - flag that enables construction of a Jacobian preconditioner

1318:   Level: intermediate

1320: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1321: @*/
1322: PetscErrorCode PetscDSUseJacobianPreconditioner(PetscDS prob, PetscBool useJacPre)
1323: {
1326:   prob->useJacPre = useJacPre;
1327:   return(0);
1328: }

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

1333:   Not collective

1335:   Input Parameter:
1336: . prob - The PetscDS

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

1341:   Level: intermediate

1343: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1344: @*/
1345: PetscErrorCode PetscDSHasJacobianPreconditioner(PetscDS prob, PetscBool *hasJacPre)
1346: {
1347:   PetscInt f, g, h;

1351:   *hasJacPre = PETSC_FALSE;
1352:   if (!prob->useJacPre) return(0);
1353:   for (f = 0; f < prob->Nf; ++f) {
1354:     for (g = 0; g < prob->Nf; ++g) {
1355:       for (h = 0; h < 4; ++h) {
1356:         if (prob->gp[(f*prob->Nf + g)*4+h]) *hasJacPre = PETSC_TRUE;
1357:       }
1358:     }
1359:   }
1360:   return(0);
1361: }

1363: /*@C
1364:   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.

1366:   Not collective

1368:   Input Parameters:
1369: + prob - The PetscDS
1370: . f    - The test field number
1371: - g    - The field number

1373:   Output Parameters:
1374: + g0 - integrand for the test and basis function term
1375: . g1 - integrand for the test function and basis function gradient term
1376: . g2 - integrand for the test function gradient and basis function term
1377: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

1385: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1386: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1387: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1388: $    PetscReal t, const PetscReal u_tShift, const PetscReal x[], PetscScalar g0[])

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

1409:   Level: intermediate

1411: .seealso: PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1412: @*/
1413: PetscErrorCode PetscDSGetJacobianPreconditioner(PetscDS prob, PetscInt f, PetscInt g,
1414:                                   void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1415:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1416:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1417:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1418:                                   void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1419:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1420:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1421:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1422:                                   void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1423:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1424:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1425:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1426:                                   void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1427:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1428:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1429:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1430: {
1433:   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);
1434:   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);
1439:   return(0);
1440: }

1442: /*@C
1443:   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.

1445:   Not collective

1447:   Input Parameters:
1448: + prob - The PetscDS
1449: . f    - The test field number
1450: . g    - The field number
1451: . g0 - integrand for the test and basis function term
1452: . g1 - integrand for the test function and basis function gradient term
1453: . g2 - integrand for the test function gradient and basis function term
1454: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

1462: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1463: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1464: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1465: $    PetscReal t, const PetscReal x[], PetscScalar g0[])

1467: + dim - the spatial dimension
1468: . Nf - the number of fields
1469: . uOff - the offset into u[] and u_t[] for each field
1470: . uOff_x - the offset into u_x[] for each field
1471: . u - each field evaluated at the current point
1472: . u_t - the time derivative of each field evaluated at the current point
1473: . u_x - the gradient of each field evaluated at the current point
1474: . aOff - the offset into a[] and a_t[] for each auxiliary field
1475: . aOff_x - the offset into a_x[] for each auxiliary field
1476: . a - each auxiliary field evaluated at the current point
1477: . a_t - the time derivative of each auxiliary field evaluated at the current point
1478: . a_x - the gradient of auxiliary each field evaluated at the current point
1479: . t - current time
1480: . u_tShift - the multiplier a for dF/dU_t
1481: . x - coordinates of the current point
1482: . numConstants - number of constant parameters
1483: . constants - constant parameters
1484: - g0 - output values at the current point

1486:   Level: intermediate

1488: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobian()
1489: @*/
1490: PetscErrorCode PetscDSSetJacobianPreconditioner(PetscDS prob, PetscInt f, PetscInt g,
1491:                                   void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1492:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1493:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1494:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1495:                                   void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1496:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1497:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1498:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1499:                                   void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1500:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1501:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1502:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1503:                                   void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1504:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1505:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1506:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1507: {

1516:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1517:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
1518:   PetscDSEnlarge_Static(prob, PetscMax(f, g)+1);
1519:   prob->gp[(f*prob->Nf + g)*4+0] = g0;
1520:   prob->gp[(f*prob->Nf + g)*4+1] = g1;
1521:   prob->gp[(f*prob->Nf + g)*4+2] = g2;
1522:   prob->gp[(f*prob->Nf + g)*4+3] = g3;
1523:   return(0);
1524: }

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

1529:   Not collective

1531:   Input Parameter:
1532: . prob - The PetscDS

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

1537:   Level: intermediate

1539: .seealso: PetscDSGetDynamicJacobian(), PetscDSSetDynamicJacobian(), PetscDSGetJacobian()
1540: @*/
1541: PetscErrorCode PetscDSHasDynamicJacobian(PetscDS prob, PetscBool *hasDynJac)
1542: {
1543:   PetscInt f, g, h;

1547:   *hasDynJac = PETSC_FALSE;
1548:   for (f = 0; f < prob->Nf; ++f) {
1549:     for (g = 0; g < prob->Nf; ++g) {
1550:       for (h = 0; h < 4; ++h) {
1551:         if (prob->gt[(f*prob->Nf + g)*4+h]) *hasDynJac = PETSC_TRUE;
1552:       }
1553:     }
1554:   }
1555:   return(0);
1556: }

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

1561:   Not collective

1563:   Input Parameters:
1564: + prob - The PetscDS
1565: . f    - The test field number
1566: - g    - The field number

1568:   Output Parameters:
1569: + g0 - integrand for the test and basis function term
1570: . g1 - integrand for the test function and basis function gradient term
1571: . g2 - integrand for the test function gradient and basis function term
1572: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

1580: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1581: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1582: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1583: $    PetscReal t, const PetscReal u_tShift, const PetscReal x[], PetscScalar g0[])

1585: + dim - the spatial dimension
1586: . Nf - the number of fields
1587: . uOff - the offset into u[] and u_t[] for each field
1588: . uOff_x - the offset into u_x[] for each field
1589: . u - each field evaluated at the current point
1590: . u_t - the time derivative of each field evaluated at the current point
1591: . u_x - the gradient of each field evaluated at the current point
1592: . aOff - the offset into a[] and a_t[] for each auxiliary field
1593: . aOff_x - the offset into a_x[] for each auxiliary field
1594: . a - each auxiliary field evaluated at the current point
1595: . a_t - the time derivative of each auxiliary field evaluated at the current point
1596: . a_x - the gradient of auxiliary each field evaluated at the current point
1597: . t - current time
1598: . u_tShift - the multiplier a for dF/dU_t
1599: . x - coordinates of the current point
1600: . numConstants - number of constant parameters
1601: . constants - constant parameters
1602: - g0 - output values at the current point

1604:   Level: intermediate

1606: .seealso: PetscDSSetJacobian()
1607: @*/
1608: PetscErrorCode PetscDSGetDynamicJacobian(PetscDS prob, PetscInt f, PetscInt g,
1609:                                          void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1610:                                                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1611:                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1612:                                                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1613:                                          void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1614:                                                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1615:                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1616:                                                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1617:                                          void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1618:                                                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1619:                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1620:                                                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1621:                                          void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1622:                                                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1623:                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1624:                                                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1625: {
1628:   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);
1629:   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);
1634:   return(0);
1635: }

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

1640:   Not collective

1642:   Input Parameters:
1643: + prob - The PetscDS
1644: . f    - The test field number
1645: . g    - The field number
1646: . g0 - integrand for the test and basis function term
1647: . g1 - integrand for the test function and basis function gradient term
1648: . g2 - integrand for the test function gradient and basis function term
1649: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

1657: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1658: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1659: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1660: $    PetscReal t, const PetscReal x[], PetscScalar g0[])

1662: + dim - the spatial dimension
1663: . Nf - the number of fields
1664: . uOff - the offset into u[] and u_t[] for each field
1665: . uOff_x - the offset into u_x[] for each field
1666: . u - each field evaluated at the current point
1667: . u_t - the time derivative of each field evaluated at the current point
1668: . u_x - the gradient of each field evaluated at the current point
1669: . aOff - the offset into a[] and a_t[] for each auxiliary field
1670: . aOff_x - the offset into a_x[] for each auxiliary field
1671: . a - each auxiliary field evaluated at the current point
1672: . a_t - the time derivative of each auxiliary field evaluated at the current point
1673: . a_x - the gradient of auxiliary each field evaluated at the current point
1674: . t - current time
1675: . u_tShift - the multiplier a for dF/dU_t
1676: . x - coordinates of the current point
1677: . numConstants - number of constant parameters
1678: . constants - constant parameters
1679: - g0 - output values at the current point

1681:   Level: intermediate

1683: .seealso: PetscDSGetJacobian()
1684: @*/
1685: PetscErrorCode PetscDSSetDynamicJacobian(PetscDS prob, PetscInt f, PetscInt g,
1686:                                          void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1687:                                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1688:                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1689:                                                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1690:                                          void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1691:                                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1692:                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1693:                                                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1694:                                          void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1695:                                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1696:                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1697:                                                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1698:                                          void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1699:                                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1700:                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1701:                                                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1702: {

1711:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1712:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
1713:   PetscDSEnlarge_Static(prob, PetscMax(f, g)+1);
1714:   prob->gt[(f*prob->Nf + g)*4+0] = g0;
1715:   prob->gt[(f*prob->Nf + g)*4+1] = g1;
1716:   prob->gt[(f*prob->Nf + g)*4+2] = g2;
1717:   prob->gt[(f*prob->Nf + g)*4+3] = g3;
1718:   return(0);
1719: }

1721: /*@C
1722:   PetscDSGetRiemannSolver - Returns the Riemann solver for the given field

1724:   Not collective

1726:   Input Arguments:
1727: + prob - The PetscDS object
1728: - f    - The field number

1730:   Output Argument:
1731: . r    - Riemann solver

1733:   Calling sequence for r:

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

1737: + dim  - The spatial dimension
1738: . Nf   - The number of fields
1739: . x    - The coordinates at a point on the interface
1740: . n    - The normal vector to the interface
1741: . uL   - The state vector to the left of the interface
1742: . uR   - The state vector to the right of the interface
1743: . flux - output array of flux through the interface
1744: . numConstants - number of constant parameters
1745: . constants - constant parameters
1746: - ctx  - optional user context

1748:   Level: intermediate

1750: .seealso: PetscDSSetRiemannSolver()
1751: @*/
1752: PetscErrorCode PetscDSGetRiemannSolver(PetscDS prob, PetscInt f,
1753:                                        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))
1754: {
1757:   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);
1759:   *r = prob->r[f];
1760:   return(0);
1761: }

1763: /*@C
1764:   PetscDSSetRiemannSolver - Sets the Riemann solver for the given field

1766:   Not collective

1768:   Input Arguments:
1769: + prob - The PetscDS object
1770: . f    - The field number
1771: - r    - Riemann solver

1773:   Calling sequence for r:

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

1777: + dim  - The spatial dimension
1778: . Nf   - The number of fields
1779: . x    - The coordinates at a point on the interface
1780: . n    - The normal vector to the interface
1781: . uL   - The state vector to the left of the interface
1782: . uR   - The state vector to the right of the interface
1783: . flux - output array of flux through the interface
1784: . numConstants - number of constant parameters
1785: . constants - constant parameters
1786: - ctx  - optional user context

1788:   Level: intermediate

1790: .seealso: PetscDSGetRiemannSolver()
1791: @*/
1792: PetscErrorCode PetscDSSetRiemannSolver(PetscDS prob, PetscInt f,
1793:                                        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))
1794: {

1800:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1801:   PetscDSEnlarge_Static(prob, f+1);
1802:   prob->r[f] = r;
1803:   return(0);
1804: }

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

1809:   Not collective

1811:   Input Parameters:
1812: + prob - The PetscDS
1813: - f    - The field number

1815:   Output Parameters:
1816: . update - update function

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

1820: $ 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[], PetscScalar uNew[])

1825: + dim - the spatial dimension
1826: . Nf - the number of fields
1827: . uOff - the offset into u[] and u_t[] for each field
1828: . uOff_x - the offset into u_x[] for each field
1829: . u - each field evaluated at the current point
1830: . u_t - the time derivative of each field evaluated at the current point
1831: . u_x - the gradient of each field evaluated at the current point
1832: . aOff - the offset into a[] and a_t[] for each auxiliary field
1833: . aOff_x - the offset into a_x[] for each auxiliary field
1834: . a - each auxiliary field evaluated at the current point
1835: . a_t - the time derivative of each auxiliary field evaluated at the current point
1836: . a_x - the gradient of auxiliary each field evaluated at the current point
1837: . t - current time
1838: . x - coordinates of the current point
1839: - uNew - new value for field at the current point

1841:   Level: intermediate

1843: .seealso: PetscDSSetUpdate(), PetscDSSetResidual()
1844: @*/
1845: PetscErrorCode PetscDSGetUpdate(PetscDS prob, PetscInt f,
1846:                                   void (**update)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1847:                                                   const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1848:                                                   const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1849:                                                   PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar uNew[]))
1850: {
1853:   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);
1855:   return(0);
1856: }

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

1861:   Not collective

1863:   Input Parameters:
1864: + prob   - The PetscDS
1865: . f      - The field number
1866: - update - update function

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

1870: $ 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[], PetscScalar uNew[])

1875: + dim - the spatial dimension
1876: . Nf - the number of fields
1877: . uOff - the offset into u[] and u_t[] for each field
1878: . uOff_x - the offset into u_x[] for each field
1879: . u - each field evaluated at the current point
1880: . u_t - the time derivative of each field evaluated at the current point
1881: . u_x - the gradient of each field evaluated at the current point
1882: . aOff - the offset into a[] and a_t[] for each auxiliary field
1883: . aOff_x - the offset into a_x[] for each auxiliary field
1884: . a - each auxiliary field evaluated at the current point
1885: . a_t - the time derivative of each auxiliary field evaluated at the current point
1886: . a_x - the gradient of auxiliary each field evaluated at the current point
1887: . t - current time
1888: . x - coordinates of the current point
1889: - uNew - new field values at the current point

1891:   Level: intermediate

1893: .seealso: PetscDSGetResidual()
1894: @*/
1895: PetscErrorCode PetscDSSetUpdate(PetscDS prob, PetscInt f,
1896:                                 void (*update)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1897:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1898:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1899:                                                PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar uNew[]))
1900: {

1906:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1907:   PetscDSEnlarge_Static(prob, f+1);
1908:   prob->update[f] = update;
1909:   return(0);
1910: }

1912: PetscErrorCode PetscDSGetContext(PetscDS prob, PetscInt f, void **ctx)
1913: {
1916:   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);
1918:   *ctx = prob->ctx[f];
1919:   return(0);
1920: }

1922: PetscErrorCode PetscDSSetContext(PetscDS prob, PetscInt f, void *ctx)
1923: {

1928:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1929:   PetscDSEnlarge_Static(prob, f+1);
1930:   prob->ctx[f] = ctx;
1931:   return(0);
1932: }

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

1937:   Not collective

1939:   Input Parameters:
1940: + prob - The PetscDS
1941: - f    - The test field number

1943:   Output Parameters:
1944: + f0 - boundary integrand for the test function term
1945: - f1 - boundary integrand for the test function gradient term

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

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

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

1953: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1954: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1955: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1956: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[])

1958: + dim - the spatial dimension
1959: . Nf - the number of fields
1960: . uOff - the offset into u[] and u_t[] for each field
1961: . uOff_x - the offset into u_x[] for each field
1962: . u - each field evaluated at the current point
1963: . u_t - the time derivative of each field evaluated at the current point
1964: . u_x - the gradient of each field evaluated at the current point
1965: . aOff - the offset into a[] and a_t[] for each auxiliary field
1966: . aOff_x - the offset into a_x[] for each auxiliary field
1967: . a - each auxiliary field evaluated at the current point
1968: . a_t - the time derivative of each auxiliary field evaluated at the current point
1969: . a_x - the gradient of auxiliary each field evaluated at the current point
1970: . t - current time
1971: . x - coordinates of the current point
1972: . n - unit normal at the current point
1973: . numConstants - number of constant parameters
1974: . constants - constant parameters
1975: - f0 - output values at the current point

1977:   Level: intermediate

1979: .seealso: PetscDSSetBdResidual()
1980: @*/
1981: PetscErrorCode PetscDSGetBdResidual(PetscDS prob, PetscInt f,
1982:                                     void (**f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1983:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1984:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1985:                                                 PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]),
1986:                                     void (**f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1987:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1988:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1989:                                                 PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
1990: {
1993:   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
1996:   return(0);
1997: }

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

2002:   Not collective

2004:   Input Parameters:
2005: + prob - The PetscDS
2006: . f    - The test field number
2007: . f0 - boundary integrand for the test function term
2008: - f1 - boundary integrand for the test function gradient term

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

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

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

2016: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2017: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2018: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2019: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[])

2021: + dim - the spatial dimension
2022: . Nf - the number of fields
2023: . uOff - the offset into u[] and u_t[] for each field
2024: . uOff_x - the offset into u_x[] for each field
2025: . u - each field evaluated at the current point
2026: . u_t - the time derivative of each field evaluated at the current point
2027: . u_x - the gradient of each field evaluated at the current point
2028: . aOff - the offset into a[] and a_t[] for each auxiliary field
2029: . aOff_x - the offset into a_x[] for each auxiliary field
2030: . a - each auxiliary field evaluated at the current point
2031: . a_t - the time derivative of each auxiliary field evaluated at the current point
2032: . a_x - the gradient of auxiliary each field evaluated at the current point
2033: . t - current time
2034: . x - coordinates of the current point
2035: . n - unit normal at the current point
2036: . numConstants - number of constant parameters
2037: . constants - constant parameters
2038: - f0 - output values at the current point

2040:   Level: intermediate

2042: .seealso: PetscDSGetBdResidual()
2043: @*/
2044: PetscErrorCode PetscDSSetBdResidual(PetscDS prob, PetscInt f,
2045:                                     void (*f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2046:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2047:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2048:                                                PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]),
2049:                                     void (*f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2050:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2051:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2052:                                                PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
2053: {

2058:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2059:   PetscDSEnlarge_Static(prob, f+1);
2062:   return(0);
2063: }

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

2068:   Not collective

2070:   Input Parameters:
2071: + prob - The PetscDS
2072: . f    - The test field number
2073: - g    - The field number

2075:   Output Parameters:
2076: + g0 - integrand for the test and basis function term
2077: . g1 - integrand for the test function and basis function gradient term
2078: . g2 - integrand for the test function gradient and basis function term
2079: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

2087: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2088: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2089: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2090: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar g0[])

2092: + dim - the spatial dimension
2093: . Nf - the number of fields
2094: . uOff - the offset into u[] and u_t[] for each field
2095: . uOff_x - the offset into u_x[] for each field
2096: . u - each field evaluated at the current point
2097: . u_t - the time derivative of each field evaluated at the current point
2098: . u_x - the gradient of each field evaluated at the current point
2099: . aOff - the offset into a[] and a_t[] for each auxiliary field
2100: . aOff_x - the offset into a_x[] for each auxiliary field
2101: . a - each auxiliary field evaluated at the current point
2102: . a_t - the time derivative of each auxiliary field evaluated at the current point
2103: . a_x - the gradient of auxiliary each field evaluated at the current point
2104: . t - current time
2105: . u_tShift - the multiplier a for dF/dU_t
2106: . x - coordinates of the current point
2107: . n - normal at the current point
2108: . numConstants - number of constant parameters
2109: . constants - constant parameters
2110: - g0 - output values at the current point

2112:   Level: intermediate

2114: .seealso: PetscDSSetBdJacobian()
2115: @*/
2116: PetscErrorCode PetscDSGetBdJacobian(PetscDS prob, PetscInt f, PetscInt g,
2117:                                     void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2118:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2119:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2120:                                                 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
2121:                                     void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2122:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2123:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2124:                                                 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
2125:                                     void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2126:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2127:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2128:                                                 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
2129:                                     void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2130:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2131:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2132:                                                 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
2133: {
2136:   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);
2137:   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);
2142:   return(0);
2143: }

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

2148:   Not collective

2150:   Input Parameters:
2151: + prob - The PetscDS
2152: . f    - The test field number
2153: . g    - The field number
2154: . g0 - integrand for the test and basis function term
2155: . g1 - integrand for the test function and basis function gradient term
2156: . g2 - integrand for the test function gradient and basis function term
2157: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

2165: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2166: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2167: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2168: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar g0[])

2170: + dim - the spatial dimension
2171: . Nf - the number of fields
2172: . uOff - the offset into u[] and u_t[] for each field
2173: . uOff_x - the offset into u_x[] for each field
2174: . u - each field evaluated at the current point
2175: . u_t - the time derivative of each field evaluated at the current point
2176: . u_x - the gradient of each field evaluated at the current point
2177: . aOff - the offset into a[] and a_t[] for each auxiliary field
2178: . aOff_x - the offset into a_x[] for each auxiliary field
2179: . a - each auxiliary field evaluated at the current point
2180: . a_t - the time derivative of each auxiliary field evaluated at the current point
2181: . a_x - the gradient of auxiliary each field evaluated at the current point
2182: . t - current time
2183: . u_tShift - the multiplier a for dF/dU_t
2184: . x - coordinates of the current point
2185: . n - normal at the current point
2186: . numConstants - number of constant parameters
2187: . constants - constant parameters
2188: - g0 - output values at the current point

2190:   Level: intermediate

2192: .seealso: PetscDSGetBdJacobian()
2193: @*/
2194: PetscErrorCode PetscDSSetBdJacobian(PetscDS prob, PetscInt f, PetscInt g,
2195:                                     void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2196:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2197:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2198:                                                PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
2199:                                     void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2200:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2201:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2202:                                                PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
2203:                                     void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2204:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2205:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2206:                                                PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
2207:                                     void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2208:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2209:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2210:                                                PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
2211: {

2220:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2221:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
2222:   PetscDSEnlarge_Static(prob, PetscMax(f, g)+1);
2223:   prob->gBd[(f*prob->Nf + g)*4+0] = g0;
2224:   prob->gBd[(f*prob->Nf + g)*4+1] = g1;
2225:   prob->gBd[(f*prob->Nf + g)*4+2] = g2;
2226:   prob->gBd[(f*prob->Nf + g)*4+3] = g3;
2227:   return(0);
2228: }

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

2233:   Not collective

2235:   Input Parameters:
2236: + prob - The PetscDS
2237: - f    - The test field number

2239:   Output Parameter:
2240: + exactSol - exact solution for the test field
2241: - exactCtx - exact solution context

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

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

2247: + dim - the spatial dimension
2248: . t - current time
2249: . x - coordinates of the current point
2250: . Nc - the number of field components
2251: . u - the solution field evaluated at the current point
2252: - ctx - a user context

2254:   Level: intermediate

2256: .seealso: PetscDSSetExactSolution()
2257: @*/
2258: PetscErrorCode PetscDSGetExactSolution(PetscDS prob, PetscInt f, PetscErrorCode (**sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void **ctx)
2259: {
2262:   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);
2265:   return(0);
2266: }

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

2271:   Not collective

2273:   Input Parameters:
2274: + prob - The PetscDS
2275: . f    - The test field number
2276: . sol  - solution function for the test fields
2277: - ctx  - solution context or NULL

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

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

2283: + dim - the spatial dimension
2284: . t - current time
2285: . x - coordinates of the current point
2286: . Nc - the number of field components
2287: . u - the solution field evaluated at the current point
2288: - ctx - a user context

2290:   Level: intermediate

2292: .seealso: PetscDSGetExactSolution()
2293: @*/
2294: PetscErrorCode PetscDSSetExactSolution(PetscDS prob, PetscInt f, PetscErrorCode (*sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void *ctx)
2295: {

2300:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2301:   PetscDSEnlarge_Static(prob, f+1);
2304:   return(0);
2305: }

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

2310:   Not collective

2312:   Input Parameter:
2313: . prob - The PetscDS object

2315:   Output Parameters:
2316: + numConstants - The number of constants
2317: - constants    - The array of constants, NULL if there are none

2319:   Level: intermediate

2321: .seealso: PetscDSSetConstants(), PetscDSCreate()
2322: @*/
2323: PetscErrorCode PetscDSGetConstants(PetscDS prob, PetscInt *numConstants, const PetscScalar *constants[])
2324: {
2329:   return(0);
2330: }

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

2335:   Not collective

2337:   Input Parameters:
2338: + prob         - The PetscDS object
2339: . numConstants - The number of constants
2340: - constants    - The array of constants, NULL if there are none

2342:   Level: intermediate

2344: .seealso: PetscDSGetConstants(), PetscDSCreate()
2345: @*/
2346: PetscErrorCode PetscDSSetConstants(PetscDS prob, PetscInt numConstants, PetscScalar constants[])
2347: {

2352:   if (numConstants != prob->numConstants) {
2353:     PetscFree(prob->constants);
2354:     prob->numConstants = numConstants;
2355:     if (prob->numConstants) {
2356:       PetscMalloc1(prob->numConstants, &prob->constants);
2357:     } else {
2358:       prob->constants = NULL;
2359:     }
2360:   }
2361:   if (prob->numConstants) {
2363:     PetscArraycpy(prob->constants, constants, prob->numConstants);
2364:   }
2365:   return(0);
2366: }

2368: /*@
2369:   PetscDSGetFieldIndex - Returns the index of the given field

2371:   Not collective

2373:   Input Parameters:
2374: + prob - The PetscDS object
2375: - disc - The discretization object

2377:   Output Parameter:
2378: . f - The field number

2380:   Level: beginner

2382: .seealso: PetscGetDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
2383: @*/
2384: PetscErrorCode PetscDSGetFieldIndex(PetscDS prob, PetscObject disc, PetscInt *f)
2385: {
2386:   PetscInt g;

2391:   *f = -1;
2392:   for (g = 0; g < prob->Nf; ++g) {if (disc == prob->disc[g]) break;}
2393:   if (g == prob->Nf) SETERRQ(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Field not found in PetscDS.");
2394:   *f = g;
2395:   return(0);
2396: }

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

2401:   Not collective

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

2407:   Output Parameter:
2408: . size - The size

2410:   Level: beginner

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

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:   PetscDSSetUp(prob);
2423:   *size = prob->Nb[f];
2424:   return(0);
2425: }

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

2430:   Not collective

2432:   Input Parameters:
2433: + prob - The PetscDS object
2434: - f - The field number

2436:   Output Parameter:
2437: . off - The offset

2439:   Level: beginner

2441: .seealso: PetscDSGetFieldSize(), PetscDSGetNumFields(), PetscDSCreate()
2442: @*/
2443: PetscErrorCode PetscDSGetFieldOffset(PetscDS prob, PetscInt f, PetscInt *off)
2444: {
2445:   PetscInt       size, g;

2451:   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);
2452:   *off = 0;
2453:   for (g = 0; g < f; ++g) {
2454:     PetscDSGetFieldSize(prob, g, &size);
2455:     *off += size;
2456:   }
2457:   return(0);
2458: }

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

2463:   Not collective

2465:   Input Parameter:
2466: . prob - The PetscDS object

2468:   Output Parameter:
2469: . dimensions - The number of dimensions

2471:   Level: beginner

2473: .seealso: PetscDSGetComponentOffsets(), PetscDSGetNumFields(), PetscDSCreate()
2474: @*/
2475: PetscErrorCode PetscDSGetDimensions(PetscDS prob, PetscInt *dimensions[])
2476: {

2481:   PetscDSSetUp(prob);
2483:   *dimensions = prob->Nb;
2484:   return(0);
2485: }

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

2490:   Not collective

2492:   Input Parameter:
2493: . prob - The PetscDS object

2495:   Output Parameter:
2496: . components - The number of components

2498:   Level: beginner

2500: .seealso: PetscDSGetComponentOffsets(), PetscDSGetNumFields(), PetscDSCreate()
2501: @*/
2502: PetscErrorCode PetscDSGetComponents(PetscDS prob, PetscInt *components[])
2503: {

2508:   PetscDSSetUp(prob);
2510:   *components = prob->Nc;
2511:   return(0);
2512: }

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

2517:   Not collective

2519:   Input Parameters:
2520: + prob - The PetscDS object
2521: - f - The field number

2523:   Output Parameter:
2524: . off - The offset

2526:   Level: beginner

2528: .seealso: PetscDSGetNumFields(), PetscDSCreate()
2529: @*/
2530: PetscErrorCode PetscDSGetComponentOffset(PetscDS prob, PetscInt f, PetscInt *off)
2531: {
2535:   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);
2536:   *off = prob->off[f];
2537:   return(0);
2538: }

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

2543:   Not collective

2545:   Input Parameter:
2546: . prob - The PetscDS object

2548:   Output Parameter:
2549: . offsets - The offsets

2551:   Level: beginner

2553: .seealso: PetscDSGetNumFields(), PetscDSCreate()
2554: @*/
2555: PetscErrorCode PetscDSGetComponentOffsets(PetscDS prob, PetscInt *offsets[])
2556: {
2560:   *offsets = prob->off;
2561:   return(0);
2562: }

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

2567:   Not collective

2569:   Input Parameter:
2570: . prob - The PetscDS object

2572:   Output Parameter:
2573: . offsets - The offsets

2575:   Level: beginner

2577: .seealso: PetscDSGetNumFields(), PetscDSCreate()
2578: @*/
2579: PetscErrorCode PetscDSGetComponentDerivativeOffsets(PetscDS prob, PetscInt *offsets[])
2580: {
2584:   *offsets = prob->offDer;
2585:   return(0);
2586: }

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

2591:   Not collective

2593:   Input Parameter:
2594: . prob - The PetscDS object

2596:   Output Parameters:
2597: + basis - The basis function tabulation at quadrature points
2598: - basisDer - The basis function derivative tabulation at quadrature points

2600:   Level: intermediate

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

2610:   PetscDSSetUp(prob);
2613:   return(0);
2614: }

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

2619:   Not collective

2621:   Input Parameter:
2622: . prob - The PetscDS object

2624:   Output Parameters:
2625: + basisFace - The basis function tabulation at quadrature points
2626: - basisDerFace - The basis function derivative tabulation at quadrature points

2628:   Level: intermediate

2630: .seealso: PetscDSGetTabulation(), PetscDSCreate()
2631: @*/
2632: PetscErrorCode PetscDSGetFaceTabulation(PetscDS prob, PetscReal ***basis, PetscReal ***basisDer)
2633: {

2638:   PetscDSSetUp(prob);
2641:   return(0);
2642: }

2644: PetscErrorCode PetscDSGetEvaluationArrays(PetscDS prob, PetscScalar **u, PetscScalar **u_t, PetscScalar **u_x)
2645: {

2650:   PetscDSSetUp(prob);
2654:   return(0);
2655: }

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

2663:   PetscDSSetUp(prob);
2670:   return(0);
2671: }

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

2679:   PetscDSSetUp(prob);
2685:   return(0);
2686: }

2688: /*@C
2689:   PetscDSAddBoundary - Add a boundary condition to the model

2691:   Input Parameters:
2692: + ds          - The PetscDS object
2693: . type        - The type of condition, e.g. DM_BC_ESSENTIAL/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
2694: . name        - The BC name
2695: . labelname   - The label defining constrained points
2696: . field       - The field to constrain
2697: . numcomps    - The number of constrained field components (0 will constrain all fields)
2698: . comps       - An array of constrained component numbers
2699: . bcFunc      - A pointwise function giving boundary values
2700: . numids      - The number of DMLabel ids for constrained points
2701: . ids         - An array of ids for constrained points
2702: - ctx         - An optional user context for bcFunc

2704:   Options Database Keys:
2705: + -bc_<boundary name> <num> - Overrides the boundary ids
2706: - -bc_<boundary name>_comp <num> - Overrides the boundary components

2708:   Level: developer

2710: .seealso: PetscDSGetBoundary()
2711: @*/
2712: 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)
2713: {
2714:   DSBoundary     b;

2719:   PetscNew(&b);
2720:   PetscStrallocpy(name, (char **) &b->name);
2721:   PetscStrallocpy(labelname, (char **) &b->labelname);
2722:   PetscMalloc1(numcomps, &b->comps);
2723:   if (numcomps) {PetscArraycpy(b->comps, comps, numcomps);}
2724:   PetscMalloc1(numids, &b->ids);
2725:   if (numids) {PetscArraycpy(b->ids, ids, numids);}
2726:   b->type            = type;
2727:   b->field           = field;
2728:   b->numcomps        = numcomps;
2729:   b->func            = bcFunc;
2730:   b->numids          = numids;
2731:   b->ctx             = ctx;
2732:   b->next            = ds->boundary;
2733:   ds->boundary       = b;
2734:   return(0);
2735: }

2737: /*@C
2738:   PetscDSUpdateBoundary - Change a boundary condition for the model

2740:   Input Parameters:
2741: + ds          - The PetscDS object
2742: . bd          - The boundary condition number
2743: . type        - The type of condition, e.g. DM_BC_ESSENTIAL/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
2744: . name        - The BC name
2745: . labelname   - The label defining constrained points
2746: . field       - The field to constrain
2747: . numcomps    - The number of constrained field components
2748: . comps       - An array of constrained component numbers
2749: . bcFunc      - A pointwise function giving boundary values
2750: . numids      - The number of DMLabel ids for constrained points
2751: . ids         - An array of ids for constrained points
2752: - ctx         - An optional user context for bcFunc

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

2756:   Level: developer

2758: .seealso: PetscDSAddBoundary(), PetscDSGetBoundary(), PetscDSGetNumBoundary()
2759: @*/
2760: 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)
2761: {
2762:   DSBoundary     b = ds->boundary;
2763:   PetscInt       n = 0;

2768:   while (b) {
2769:     if (n == bd) break;
2770:     b = b->next;
2771:     ++n;
2772:   }
2773:   if (!b) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %d is not in [0, %d)", bd, n);
2774:   if (name) {
2775:     PetscFree(b->name);
2776:     PetscStrallocpy(name, (char **) &b->name);
2777:   }
2778:   if (labelname) {
2779:     PetscFree(b->labelname);
2780:     PetscStrallocpy(labelname, (char **) &b->labelname);
2781:   }
2782:   if (numcomps >= 0 && numcomps != b->numcomps) {
2783:     b->numcomps = numcomps;
2784:     PetscFree(b->comps);
2785:     PetscMalloc1(numcomps, &b->comps);
2786:     if (numcomps) {PetscArraycpy(b->comps, comps, numcomps);}
2787:   }
2788:   if (numids >= 0 && numids != b->numids) {
2789:     b->numids = numids;
2790:     PetscFree(b->ids);
2791:     PetscMalloc1(numids, &b->ids);
2792:     if (numids) {PetscArraycpy(b->ids, ids, numids);}
2793:   }
2794:   b->type = type;
2795:   if (field >= 0) {b->field  = field;}
2796:   if (bcFunc)     {b->func   = bcFunc;}
2797:   if (ctx)        {b->ctx    = ctx;}
2798:   return(0);
2799: }

2801: /*@
2802:   PetscDSGetNumBoundary - Get the number of registered BC

2804:   Input Parameters:
2805: . ds - The PetscDS object

2807:   Output Parameters:
2808: . numBd - The number of BC

2810:   Level: intermediate

2812: .seealso: PetscDSAddBoundary(), PetscDSGetBoundary()
2813: @*/
2814: PetscErrorCode PetscDSGetNumBoundary(PetscDS ds, PetscInt *numBd)
2815: {
2816:   DSBoundary b = ds->boundary;

2821:   *numBd = 0;
2822:   while (b) {++(*numBd); b = b->next;}
2823:   return(0);
2824: }

2826: /*@C
2827:   PetscDSGetBoundary - Gets a boundary condition to the model

2829:   Input Parameters:
2830: + ds          - The PetscDS object
2831: - bd          - The BC number

2833:   Output Parameters:
2834: + type        - The type of condition, e.g. DM_BC_ESSENTIAL/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
2835: . name        - The BC name
2836: . labelname   - The label defining constrained points
2837: . field       - The field to constrain
2838: . numcomps    - The number of constrained field components
2839: . comps       - An array of constrained component numbers
2840: . bcFunc      - A pointwise function giving boundary values
2841: . numids      - The number of DMLabel ids for constrained points
2842: . ids         - An array of ids for constrained points
2843: - ctx         - An optional user context for bcFunc

2845:   Options Database Keys:
2846: + -bc_<boundary name> <num> - Overrides the boundary ids
2847: - -bc_<boundary name>_comp <num> - Overrides the boundary components

2849:   Level: developer

2851: .seealso: PetscDSAddBoundary()
2852: @*/
2853: 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)
2854: {
2855:   DSBoundary b    = ds->boundary;
2856:   PetscInt   n    = 0;

2860:   while (b) {
2861:     if (n == bd) break;
2862:     b = b->next;
2863:     ++n;
2864:   }
2865:   if (!b) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %d is not in [0, %d)", bd, n);
2866:   if (type) {
2868:     *type = b->type;
2869:   }
2870:   if (name) {
2872:     *name = b->name;
2873:   }
2874:   if (labelname) {
2876:     *labelname = b->labelname;
2877:   }
2878:   if (field) {
2880:     *field = b->field;
2881:   }
2882:   if (numcomps) {
2884:     *numcomps = b->numcomps;
2885:   }
2886:   if (comps) {
2888:     *comps = b->comps;
2889:   }
2890:   if (func) {
2892:     *func = b->func;
2893:   }
2894:   if (numids) {
2896:     *numids = b->numids;
2897:   }
2898:   if (ids) {
2900:     *ids = b->ids;
2901:   }
2902:   if (ctx) {
2904:     *ctx = b->ctx;
2905:   }
2906:   return(0);
2907: }

2909: /*@
2910:   PetscDSCopyBoundary - Copy all boundary condition objects to the new problem

2912:   Not collective

2914:   Input Parameter:
2915: . prob - The PetscDS object

2917:   Output Parameter:
2918: . newprob - The PetscDS copy

2920:   Level: intermediate

2922: .seealso: PetscDSCopyEquations(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
2923: @*/
2924: PetscErrorCode PetscDSCopyBoundary(PetscDS probA, PetscDS probB)
2925: {
2926:   DSBoundary     b, next, *lastnext;

2932:   if (probA == probB) return(0);
2933:   next = probB->boundary;
2934:   while (next) {
2935:     DSBoundary b = next;

2937:     next = b->next;
2938:     PetscFree(b->comps);
2939:     PetscFree(b->ids);
2940:     PetscFree(b->name);
2941:     PetscFree(b->labelname);
2942:     PetscFree(b);
2943:   }
2944:   lastnext = &(probB->boundary);
2945:   for (b = probA->boundary; b; b = b->next) {
2946:     DSBoundary bNew;

2948:     PetscNew(&bNew);
2949:     bNew->numcomps = b->numcomps;
2950:     PetscMalloc1(bNew->numcomps, &bNew->comps);
2951:     PetscArraycpy(bNew->comps, b->comps, bNew->numcomps);
2952:     bNew->numids = b->numids;
2953:     PetscMalloc1(bNew->numids, &bNew->ids);
2954:     PetscArraycpy(bNew->ids, b->ids, bNew->numids);
2955:     PetscStrallocpy(b->labelname,(char **) &(bNew->labelname));
2956:     PetscStrallocpy(b->name,(char **) &(bNew->name));
2957:     bNew->ctx   = b->ctx;
2958:     bNew->type  = b->type;
2959:     bNew->field = b->field;
2960:     bNew->func  = b->func;

2962:     *lastnext = bNew;
2963:     lastnext = &(bNew->next);
2964:   }
2965:   return(0);
2966: }

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

2971:   Not collective

2973:   Input Parameter:
2974: + prob - The PetscDS object
2975: . numFields - Number of new fields
2976: - fields - Old field number for each new field

2978:   Output Parameter:
2979: . newprob - The PetscDS copy

2981:   Level: intermediate

2983: .seealso: PetscDSCopyBoundary(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
2984: @*/
2985: PetscErrorCode PetscDSSelectEquations(PetscDS prob, PetscInt numFields, const PetscInt fields[], PetscDS newprob)
2986: {
2987:   PetscInt       Nf, Nfn, fn, gn;

2994:   PetscDSGetNumFields(prob, &Nf);
2995:   PetscDSGetNumFields(newprob, &Nfn);
2996:   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);
2997:   for (fn = 0; fn < numFields; ++fn) {
2998:     const PetscInt   f = fields ? fields[fn] : fn;
2999:     PetscPointFunc   obj;
3000:     PetscPointFunc   f0, f1;
3001:     PetscBdPointFunc f0Bd, f1Bd;
3002:     PetscRiemannFunc r;

3004:     if (f >= Nf) continue;
3005:     PetscDSGetObjective(prob, f, &obj);
3006:     PetscDSGetResidual(prob, f, &f0, &f1);
3007:     PetscDSGetBdResidual(prob, f, &f0Bd, &f1Bd);
3008:     PetscDSGetRiemannSolver(prob, f, &r);
3009:     PetscDSSetObjective(newprob, fn, obj);
3010:     PetscDSSetResidual(newprob, fn, f0, f1);
3011:     PetscDSSetBdResidual(newprob, fn, f0Bd, f1Bd);
3012:     PetscDSSetRiemannSolver(newprob, fn, r);
3013:     for (gn = 0; gn < numFields; ++gn) {
3014:       const PetscInt  g = fields ? fields[gn] : gn;
3015:       PetscPointJac   g0, g1, g2, g3;
3016:       PetscPointJac   g0p, g1p, g2p, g3p;
3017:       PetscBdPointJac g0Bd, g1Bd, g2Bd, g3Bd;

3019:       if (g >= Nf) continue;
3020:       PetscDSGetJacobian(prob, f, g, &g0, &g1, &g2, &g3);
3021:       PetscDSGetJacobianPreconditioner(prob, f, g, &g0p, &g1p, &g2p, &g3p);
3022:       PetscDSGetBdJacobian(prob, f, g, &g0Bd, &g1Bd, &g2Bd, &g3Bd);
3023:       PetscDSSetJacobian(newprob, fn, gn, g0, g1, g2, g3);
3024:       PetscDSSetJacobianPreconditioner(prob, fn, gn, g0p, g1p, g2p, g3p);
3025:       PetscDSSetBdJacobian(newprob, fn, gn, g0Bd, g1Bd, g2Bd, g3Bd);
3026:     }
3027:   }
3028:   return(0);
3029: }

3031: /*@
3032:   PetscDSCopyEquations - Copy all pointwise function pointers to the new problem

3034:   Not collective

3036:   Input Parameter:
3037: . prob - The PetscDS object

3039:   Output Parameter:
3040: . newprob - The PetscDS copy

3042:   Level: intermediate

3044: .seealso: PetscDSCopyBoundary(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3045: @*/
3046: PetscErrorCode PetscDSCopyEquations(PetscDS prob, PetscDS newprob)
3047: {
3048:   PetscInt       Nf, Ng;

3054:   PetscDSGetNumFields(prob, &Nf);
3055:   PetscDSGetNumFields(newprob, &Ng);
3056:   if (Nf != Ng) SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_SIZ, "Number of fields must match %D != %D", Nf, Ng);
3057:   PetscDSSelectEquations(prob, Nf, NULL, newprob);
3058:   return(0);
3059: }
3060: /*@
3061:   PetscDSCopyConstants - Copy all constants to the new problem

3063:   Not collective

3065:   Input Parameter:
3066: . prob - The PetscDS object

3068:   Output Parameter:
3069: . newprob - The PetscDS copy

3071:   Level: intermediate

3073: .seealso: PetscDSCopyBoundary(), PetscDSCopyEquations(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3074: @*/
3075: PetscErrorCode PetscDSCopyConstants(PetscDS prob, PetscDS newprob)
3076: {
3077:   PetscInt           Nc;
3078:   const PetscScalar *constants;
3079:   PetscErrorCode     ierr;

3084:   PetscDSGetConstants(prob, &Nc, &constants);
3085:   PetscDSSetConstants(newprob, Nc, (PetscScalar *) constants);
3086:   return(0);
3087: }

3089: PetscErrorCode PetscDSGetHeightSubspace(PetscDS prob, PetscInt height, PetscDS *subprob)
3090: {
3091:   PetscInt       dim, Nf, f;

3097:   if (height == 0) {*subprob = prob; return(0);}
3098:   PetscDSGetNumFields(prob, &Nf);
3099:   PetscDSGetSpatialDimension(prob, &dim);
3100:   if (height > dim) SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_OUTOFRANGE, "DS can only handle height in [0, %D], not %D", dim, height);
3101:   if (!prob->subprobs) {PetscCalloc1(dim, &prob->subprobs);}
3102:   if (!prob->subprobs[height-1]) {
3103:     PetscInt cdim;

3105:     PetscDSCreate(PetscObjectComm((PetscObject) prob), &prob->subprobs[height-1]);
3106:     PetscDSGetCoordinateDimension(prob, &cdim);
3107:     PetscDSSetCoordinateDimension(prob->subprobs[height-1], cdim);
3108:     for (f = 0; f < Nf; ++f) {
3109:       PetscFE      subfe;
3110:       PetscObject  obj;
3111:       PetscClassId id;

3113:       PetscDSGetDiscretization(prob, f, &obj);
3114:       PetscObjectGetClassId(obj, &id);
3115:       if (id == PETSCFE_CLASSID) {PetscFEGetHeightSubspace((PetscFE) obj, height, &subfe);}
3116:       else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unsupported discretization type for field %d", f);
3117:       PetscDSSetDiscretization(prob->subprobs[height-1], f, (PetscObject) subfe);
3118:     }
3119:   }
3120:   *subprob = prob->subprobs[height-1];
3121:   return(0);
3122: }

3124: static PetscErrorCode PetscDSDestroy_Basic(PetscDS prob)
3125: {
3126:   PetscErrorCode      ierr;

3129:   PetscFree(prob->data);
3130:   return(0);
3131: }

3133: static PetscErrorCode PetscDSInitialize_Basic(PetscDS prob)
3134: {
3136:   prob->ops->setfromoptions = NULL;
3137:   prob->ops->setup          = NULL;
3138:   prob->ops->view           = NULL;
3139:   prob->ops->destroy        = PetscDSDestroy_Basic;
3140:   return(0);
3141: }

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

3146:   Level: intermediate

3148: .seealso: PetscDSType, PetscDSCreate(), PetscDSSetType()
3149: M*/

3151: PETSC_EXTERN PetscErrorCode PetscDSCreate_Basic(PetscDS prob)
3152: {
3153:   PetscDS_Basic *b;
3154:   PetscErrorCode      ierr;

3158:   PetscNewLog(prob, &b);
3159:   prob->data = b;

3161:   PetscDSInitialize_Basic(prob);
3162:   return(0);
3163: }