Actual source code: dtds.c

petsc-master 2020-11-24
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:   PetscMalloc2(Nf,&prob->T,Nf,&prob->Tf);
382:   for (f = 0; f < Nf; ++f) {
383:     PetscObject     obj;
384:     PetscClassId    id;
385:     PetscQuadrature q = NULL;
386:     PetscInt        Nq = 0, Nb, Nc;

388:     PetscDSGetDiscretization(prob, f, &obj);
389:     if (!obj) {
390:       /* Empty mesh */
391:       Nb = Nc = 0;
392:       prob->T[f] = prob->Tf[f] = NULL;
393:     } else {
394:       PetscObjectGetClassId(obj, &id);
395:       if (id == PETSCFE_CLASSID)      {
396:         PetscFE fe = (PetscFE) obj;

398:         PetscFEGetQuadrature(fe, &q);
399:         PetscFEGetDimension(fe, &Nb);
400:         PetscFEGetNumComponents(fe, &Nc);
401:         PetscFEGetCellTabulation(fe, &prob->T[f]);
402:         PetscFEGetFaceTabulation(fe, &prob->Tf[f]);
403:       } else if (id == PETSCFV_CLASSID) {
404:         PetscFV fv = (PetscFV) obj;

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

438: static PetscErrorCode PetscDSDestroyStructs_Static(PetscDS prob)
439: {

443:   PetscFree2(prob->Nc,prob->Nb);
444:   PetscFree2(prob->off,prob->offDer);
445:   PetscFree2(prob->T,prob->Tf);
446:   PetscFree3(prob->u,prob->u_t,prob->u_x);
447:   PetscFree5(prob->x,prob->basisReal, prob->basisDerReal,prob->testReal,prob->testDerReal);
448:   PetscFree6(prob->f0,prob->f1,prob->g0,prob->g1,prob->g2,prob->g3);
449:   return(0);
450: }

452: static PetscErrorCode PetscDSEnlarge_Static(PetscDS prob, PetscInt NfNew)
453: {
454:   PetscObject      *tmpd;
455:   PetscBool        *tmpi;
456:   PetscPointFunc   *tmpobj, *tmpf, *tmpup;
457:   PetscPointJac    *tmpg, *tmpgp, *tmpgt;
458:   PetscBdPointFunc *tmpfbd;
459:   PetscBdPointJac  *tmpgbd, *tmpgpbd;
460:   PetscRiemannFunc *tmpr;
461:   PetscSimplePointFunc *tmpexactSol,  *tmpexactSol_t;
462:   void                **tmpexactCtx, **tmpexactCtx_t;
463:   void            **tmpctx;
464:   PetscInt          Nf = prob->Nf, f;
465:   PetscErrorCode    ierr;

468:   if (Nf >= NfNew) return(0);
469:   prob->setup = PETSC_FALSE;
470:   PetscDSDestroyStructs_Static(prob);
471:   PetscMalloc2(NfNew, &tmpd, NfNew, &tmpi);
472:   for (f = 0; f < Nf; ++f) {tmpd[f] = prob->disc[f]; tmpi[f] = prob->implicit[f];}
473:   for (f = Nf; f < NfNew; ++f) {tmpd[f] = NULL; tmpi[f] = PETSC_TRUE;}
474:   PetscFree2(prob->disc, prob->implicit);
475:   prob->Nf        = NfNew;
476:   prob->disc      = tmpd;
477:   prob->implicit  = tmpi;
478:   PetscCalloc7(NfNew, &tmpobj, NfNew*2, &tmpf, NfNew*NfNew*4, &tmpg, NfNew*NfNew*4, &tmpgp, NfNew*NfNew*4, &tmpgt, NfNew, &tmpr, NfNew, &tmpctx);
479:   PetscCalloc1(NfNew, &tmpup);
480:   for (f = 0; f < Nf; ++f) tmpobj[f] = prob->obj[f];
481:   for (f = 0; f < Nf*2; ++f) tmpf[f] = prob->f[f];
482:   for (f = 0; f < Nf*Nf*4; ++f) tmpg[f] = prob->g[f];
483:   for (f = 0; f < Nf*Nf*4; ++f) tmpgp[f] = prob->gp[f];
484:   for (f = 0; f < Nf; ++f) tmpr[f] = prob->r[f];
485:   for (f = 0; f < Nf; ++f) tmpup[f] = prob->update[f];
486:   for (f = 0; f < Nf; ++f) tmpctx[f] = prob->ctx[f];
487:   for (f = Nf; f < NfNew; ++f) tmpobj[f] = NULL;
488:   for (f = Nf*2; f < NfNew*2; ++f) tmpf[f] = NULL;
489:   for (f = Nf*Nf*4; f < NfNew*NfNew*4; ++f) tmpg[f] = NULL;
490:   for (f = Nf*Nf*4; f < NfNew*NfNew*4; ++f) tmpgp[f] = NULL;
491:   for (f = Nf*Nf*4; f < NfNew*NfNew*4; ++f) tmpgt[f] = NULL;
492:   for (f = Nf; f < NfNew; ++f) tmpr[f] = NULL;
493:   for (f = Nf; f < NfNew; ++f) tmpup[f] = NULL;
494:   for (f = Nf; f < NfNew; ++f) tmpctx[f] = NULL;
495:   PetscFree7(prob->obj, prob->f, prob->g, prob->gp, prob->gt, prob->r, prob->ctx);
496:   PetscFree(prob->update);
497:   prob->obj = tmpobj;
498:   prob->f   = tmpf;
499:   prob->g   = tmpg;
500:   prob->gp  = tmpgp;
501:   prob->gt  = tmpgt;
502:   prob->r   = tmpr;
503:   prob->update = tmpup;
504:   prob->ctx = tmpctx;
505:   PetscCalloc7(NfNew*2, &tmpfbd, NfNew*NfNew*4, &tmpgbd, NfNew*NfNew*4, &tmpgpbd, NfNew, &tmpexactSol, NfNew, &tmpexactCtx, NfNew, &tmpexactSol_t, NfNew, &tmpexactCtx_t);
506:   for (f = 0; f < Nf*2; ++f) tmpfbd[f] = prob->fBd[f];
507:   for (f = 0; f < Nf*Nf*4; ++f) tmpgbd[f] = prob->gBd[f];
508:   for (f = 0; f < Nf*Nf*4; ++f) tmpgpbd[f] = prob->gpBd[f];
509:   for (f = 0; f < Nf; ++f) tmpexactSol[f] = prob->exactSol[f];
510:   for (f = 0; f < Nf; ++f) tmpexactCtx[f] = prob->exactCtx[f];
511:   for (f = 0; f < Nf; ++f) tmpexactSol_t[f] = prob->exactSol_t[f];
512:   for (f = 0; f < Nf; ++f) tmpexactCtx_t[f] = prob->exactCtx_t[f];
513:   for (f = Nf*2; f < NfNew*2; ++f) tmpfbd[f] = NULL;
514:   for (f = Nf*Nf*4; f < NfNew*NfNew*4; ++f) tmpgbd[f] = NULL;
515:   for (f = Nf*Nf*4; f < NfNew*NfNew*4; ++f) tmpgpbd[f] = NULL;
516:   for (f = Nf; f < NfNew; ++f) tmpexactSol[f] = NULL;
517:   for (f = Nf; f < NfNew; ++f) tmpexactCtx[f] = NULL;
518:   for (f = Nf; f < NfNew; ++f) tmpexactSol_t[f] = NULL;
519:   for (f = Nf; f < NfNew; ++f) tmpexactCtx_t[f] = NULL;
520:   PetscFree7(prob->fBd, prob->gBd, prob->gpBd, prob->exactSol, prob->exactCtx, prob->exactSol_t, prob->exactCtx_t);
521:   prob->fBd = tmpfbd;
522:   prob->gBd = tmpgbd;
523:   prob->gpBd = tmpgpbd;
524:   prob->exactSol = tmpexactSol;
525:   prob->exactCtx = tmpexactCtx;
526:   prob->exactSol_t = tmpexactSol_t;
527:   prob->exactCtx_t = tmpexactCtx_t;
528:   return(0);
529: }

531: /*@
532:   PetscDSDestroy - Destroys a PetscDS object

534:   Collective on prob

536:   Input Parameter:
537: . prob - the PetscDS object to destroy

539:   Level: developer

541: .seealso PetscDSView()
542: @*/
543: PetscErrorCode PetscDSDestroy(PetscDS *prob)
544: {
545:   PetscInt       f;
546:   DSBoundary     next;

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

553:   if (--((PetscObject)(*prob))->refct > 0) {*prob = NULL; return(0);}
554:   ((PetscObject) (*prob))->refct = 0;
555:   if ((*prob)->subprobs) {
556:     PetscInt dim, d;

558:     PetscDSGetSpatialDimension(*prob, &dim);
559:     for (d = 0; d < dim; ++d) {PetscDSDestroy(&(*prob)->subprobs[d]);}
560:   }
561:   PetscFree((*prob)->subprobs);
562:   PetscDSDestroyStructs_Static(*prob);
563:   for (f = 0; f < (*prob)->Nf; ++f) {
564:     PetscObjectDereference((*prob)->disc[f]);
565:   }
566:   PetscFree2((*prob)->disc, (*prob)->implicit);
567:   PetscFree7((*prob)->obj,(*prob)->f,(*prob)->g,(*prob)->gp,(*prob)->gt,(*prob)->r,(*prob)->ctx);
568:   PetscFree((*prob)->update);
569:   PetscFree7((*prob)->fBd,(*prob)->gBd,(*prob)->gpBd,(*prob)->exactSol,(*prob)->exactCtx,(*prob)->exactSol_t,(*prob)->exactCtx_t);
570:   if ((*prob)->ops->destroy) {(*(*prob)->ops->destroy)(*prob);}
571:   next = (*prob)->boundary;
572:   while (next) {
573:     DSBoundary b = next;

575:     next = b->next;
576:     PetscFree(b->comps);
577:     PetscFree(b->ids);
578:     PetscFree(b->name);
579:     PetscFree(b->labelname);
580:     PetscFree(b);
581:   }
582:   PetscFree((*prob)->constants);
583:   PetscHeaderDestroy(prob);
584:   return(0);
585: }

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

590:   Collective

592:   Input Parameter:
593: . comm - The communicator for the PetscDS object

595:   Output Parameter:
596: . prob - The PetscDS object

598:   Level: beginner

600: .seealso: PetscDSSetType(), PETSCDSBASIC
601: @*/
602: PetscErrorCode PetscDSCreate(MPI_Comm comm, PetscDS *prob)
603: {
604:   PetscDS   p;

609:   *prob  = NULL;
610:   PetscDSInitializePackage();

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

614:   p->Nf           = 0;
615:   p->setup        = PETSC_FALSE;
616:   p->numConstants = 0;
617:   p->constants    = NULL;
618:   p->dimEmbed     = -1;
619:   p->useJacPre    = PETSC_TRUE;

621:   *prob = p;
622:   return(0);
623: }

625: /*@
626:   PetscDSGetNumFields - Returns the number of fields in the DS

628:   Not collective

630:   Input Parameter:
631: . prob - The PetscDS object

633:   Output Parameter:
634: . Nf - The number of fields

636:   Level: beginner

638: .seealso: PetscDSGetSpatialDimension(), PetscDSCreate()
639: @*/
640: PetscErrorCode PetscDSGetNumFields(PetscDS prob, PetscInt *Nf)
641: {
645:   *Nf = prob->Nf;
646:   return(0);
647: }

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

652:   Not collective

654:   Input Parameter:
655: . prob - The PetscDS object

657:   Output Parameter:
658: . dim - The spatial dimension

660:   Level: beginner

662: .seealso: PetscDSGetCoordinateDimension(), PetscDSGetNumFields(), PetscDSCreate()
663: @*/
664: PetscErrorCode PetscDSGetSpatialDimension(PetscDS prob, PetscInt *dim)
665: {

671:   *dim = 0;
672:   if (prob->Nf) {
673:     PetscObject  obj;
674:     PetscClassId id;

676:     PetscDSGetDiscretization(prob, 0, &obj);
677:     if (obj) {
678:       PetscObjectGetClassId(obj, &id);
679:       if (id == PETSCFE_CLASSID)      {PetscFEGetSpatialDimension((PetscFE) obj, dim);}
680:       else if (id == PETSCFV_CLASSID) {PetscFVGetSpatialDimension((PetscFV) obj, dim);}
681:       else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", 0);
682:     }
683:   }
684:   return(0);
685: }

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

690:   Not collective

692:   Input Parameter:
693: . prob - The PetscDS object

695:   Output Parameter:
696: . dimEmbed - The coordinate dimension

698:   Level: beginner

700: .seealso: PetscDSSetCoordinateDimension(), PetscDSGetSpatialDimension(), PetscDSGetNumFields(), PetscDSCreate()
701: @*/
702: PetscErrorCode PetscDSGetCoordinateDimension(PetscDS prob, PetscInt *dimEmbed)
703: {
707:   if (prob->dimEmbed < 0) SETERRQ(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONGSTATE, "No coordinate dimension set for this DS");
708:   *dimEmbed = prob->dimEmbed;
709:   return(0);
710: }

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

715:   Logically collective on prob

717:   Input Parameters:
718: + prob - The PetscDS object
719: - dimEmbed - The coordinate dimension

721:   Level: beginner

723: .seealso: PetscDSGetCoordinateDimension(), PetscDSGetSpatialDimension(), PetscDSGetNumFields(), PetscDSCreate()
724: @*/
725: PetscErrorCode PetscDSSetCoordinateDimension(PetscDS prob, PetscInt dimEmbed)
726: {
729:   if (dimEmbed < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Coordinate dimension must be non-negative, not %D", dimEmbed);
730:   prob->dimEmbed = dimEmbed;
731:   return(0);
732: }

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

737:   Not collective

739:   Input Parameter:
740: . prob - The PetscDS object

742:   Output Parameter:
743: . isHybrid - The flag

745:   Level: developer

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

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

761:   Not collective

763:   Input Parameters:
764: + prob - The PetscDS object
765: - isHybrid - The flag

767:   Level: developer

769: .seealso: PetscDSGetHybrid(), PetscDSCreate()
770: @*/
771: PetscErrorCode PetscDSSetHybrid(PetscDS prob, PetscBool isHybrid)
772: {
775:   prob->isHybrid = isHybrid;
776:   return(0);
777: }

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

782:   Not collective

784:   Input Parameter:
785: . prob - The PetscDS object

787:   Output Parameter:
788: . dim - The total problem dimension

790:   Level: beginner

792: .seealso: PetscDSGetNumFields(), PetscDSCreate()
793: @*/
794: PetscErrorCode PetscDSGetTotalDimension(PetscDS prob, PetscInt *dim)
795: {

800:   PetscDSSetUp(prob);
802:   *dim = prob->totDim;
803:   return(0);
804: }

806: /*@
807:   PetscDSGetTotalComponents - Returns the total number of components in this system

809:   Not collective

811:   Input Parameter:
812: . prob - The PetscDS object

814:   Output Parameter:
815: . dim - The total number of components

817:   Level: beginner

819: .seealso: PetscDSGetNumFields(), PetscDSCreate()
820: @*/
821: PetscErrorCode PetscDSGetTotalComponents(PetscDS prob, PetscInt *Nc)
822: {

827:   PetscDSSetUp(prob);
829:   *Nc = prob->totComp;
830:   return(0);
831: }

833: /*@
834:   PetscDSGetDiscretization - Returns the discretization object for the given field

836:   Not collective

838:   Input Parameters:
839: + prob - The PetscDS object
840: - f - The field number

842:   Output Parameter:
843: . disc - The discretization object

845:   Level: beginner

847: .seealso: PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
848: @*/
849: PetscErrorCode PetscDSGetDiscretization(PetscDS prob, PetscInt f, PetscObject *disc)
850: {
854:   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);
855:   *disc = prob->disc[f];
856:   return(0);
857: }

859: /*@
860:   PetscDSSetDiscretization - Sets the discretization object for the given field

862:   Not collective

864:   Input Parameters:
865: + prob - The PetscDS object
866: . f - The field number
867: - disc - The discretization object

869:   Level: beginner

871: .seealso: PetscDSGetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
872: @*/
873: PetscErrorCode PetscDSSetDiscretization(PetscDS prob, PetscInt f, PetscObject disc)
874: {

880:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
881:   PetscDSEnlarge_Static(prob, f+1);
882:   PetscObjectDereference(prob->disc[f]);
883:   prob->disc[f] = disc;
884:   PetscObjectReference(disc);
885:   if (disc) {
886:     PetscClassId id;

888:     PetscObjectGetClassId(disc, &id);
889:     if (id == PETSCFE_CLASSID) {
890:       PetscDSSetImplicit(prob, f, PETSC_TRUE);
891:     } else if (id == PETSCFV_CLASSID) {
892:       PetscDSSetImplicit(prob, f, PETSC_FALSE);
893:     }
894:   }
895:   return(0);
896: }

898: /*@
899:   PetscDSAddDiscretization - Adds a discretization object

901:   Not collective

903:   Input Parameters:
904: + prob - The PetscDS object
905: - disc - The boundary discretization object

907:   Level: beginner

909: .seealso: PetscDSGetDiscretization(), PetscDSSetDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
910: @*/
911: PetscErrorCode PetscDSAddDiscretization(PetscDS prob, PetscObject disc)
912: {

916:   PetscDSSetDiscretization(prob, prob->Nf, disc);
917:   return(0);
918: }

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

923:   Not collective

925:   Input Parameter:
926: . prob - The PetscDS object

928:   Output Parameter:
929: . q - The quadrature object

931: Level: intermediate

933: .seealso: PetscDSSetImplicit(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
934: @*/
935: PetscErrorCode PetscDSGetQuadrature(PetscDS prob, PetscQuadrature *q)
936: {
937:   PetscObject    obj;
938:   PetscClassId   id;

942:   *q = NULL;
943:   if (!prob->Nf) return(0);
944:   PetscDSGetDiscretization(prob, 0, &obj);
945:   PetscObjectGetClassId(obj, &id);
946:   if      (id == PETSCFE_CLASSID) {PetscFEGetQuadrature((PetscFE) obj, q);}
947:   else if (id == PETSCFV_CLASSID) {PetscFVGetQuadrature((PetscFV) obj, q);}
948:   else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", 0);
949:   return(0);
950: }

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

955:   Not collective

957:   Input Parameters:
958: + prob - The PetscDS object
959: - f - The field number

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

964:   Level: developer

966: .seealso: PetscDSSetImplicit(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
967: @*/
968: PetscErrorCode PetscDSGetImplicit(PetscDS prob, PetscInt f, PetscBool *implicit)
969: {
973:   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);
974:   *implicit = prob->implicit[f];
975:   return(0);
976: }

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

981:   Not collective

983:   Input Parameters:
984: + prob - The PetscDS object
985: . f - The field number
986: - implicit - The flag indicating what kind of solve to use for this field

988:   Level: developer

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

1001: PetscErrorCode PetscDSGetObjective(PetscDS prob, PetscInt f,
1002:                                    void (**obj)(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[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[]))
1006: {
1010:   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);
1011:   *obj = prob->obj[f];
1012:   return(0);
1013: }

1015: PetscErrorCode PetscDSSetObjective(PetscDS prob, PetscInt f,
1016:                                    void (*obj)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1017:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1018:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1019:                                                PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[]))
1020: {

1026:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1027:   PetscDSEnlarge_Static(prob, f+1);
1028:   prob->obj[f] = obj;
1029:   return(0);
1030: }

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

1035:   Not collective

1037:   Input Parameters:
1038: + prob - The PetscDS
1039: - f    - The test field number

1041:   Output Parameters:
1042: + f0 - integrand for the test function term
1043: - f1 - integrand for the test function gradient term

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

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

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

1051: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1052: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1053: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1054: $    PetscReal t, const PetscReal x[], PetscScalar f0[])

1056: + dim - the spatial dimension
1057: . Nf - the number of fields
1058: . uOff - the offset into u[] and u_t[] for each field
1059: . uOff_x - the offset into u_x[] for each field
1060: . u - each field evaluated at the current point
1061: . u_t - the time derivative of each field evaluated at the current point
1062: . u_x - the gradient of each field evaluated at the current point
1063: . aOff - the offset into a[] and a_t[] for each auxiliary field
1064: . aOff_x - the offset into a_x[] for each auxiliary field
1065: . a - each auxiliary field evaluated at the current point
1066: . a_t - the time derivative of each auxiliary field evaluated at the current point
1067: . a_x - the gradient of auxiliary each field evaluated at the current point
1068: . t - current time
1069: . x - coordinates of the current point
1070: . numConstants - number of constant parameters
1071: . constants - constant parameters
1072: - f0 - output values at the current point

1074:   Level: intermediate

1076: .seealso: PetscDSSetResidual()
1077: @*/
1078: PetscErrorCode PetscDSGetResidual(PetscDS prob, PetscInt f,
1079:                                   void (**f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1080:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1081:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1082:                                               PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]),
1083:                                   void (**f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1084:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1085:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1086:                                               PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
1087: {
1090:   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);
1093:   return(0);
1094: }

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

1099:   Not collective

1101:   Input Parameters:
1102: + prob - The PetscDS
1103: . f    - The test field number
1104: . f0 - integrand for the test function term
1105: - f1 - integrand for the test function gradient term

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

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

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

1113: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1114: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1115: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1116: $    PetscReal t, const PetscReal x[], PetscScalar f0[])

1118: + dim - the spatial dimension
1119: . Nf - the number of fields
1120: . uOff - the offset into u[] and u_t[] for each field
1121: . uOff_x - the offset into u_x[] for each field
1122: . u - each field evaluated at the current point
1123: . u_t - the time derivative of each field evaluated at the current point
1124: . u_x - the gradient of each field evaluated at the current point
1125: . aOff - the offset into a[] and a_t[] for each auxiliary field
1126: . aOff_x - the offset into a_x[] for each auxiliary field
1127: . a - each auxiliary field evaluated at the current point
1128: . a_t - the time derivative of each auxiliary field evaluated at the current point
1129: . a_x - the gradient of auxiliary each field evaluated at the current point
1130: . t - current time
1131: . x - coordinates of the current point
1132: . numConstants - number of constant parameters
1133: . constants - constant parameters
1134: - f0 - output values at the current point

1136:   Level: intermediate

1138: .seealso: PetscDSGetResidual()
1139: @*/
1140: PetscErrorCode PetscDSSetResidual(PetscDS prob, PetscInt f,
1141:                                   void (*f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1142:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1143:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1144:                                              PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]),
1145:                                   void (*f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1146:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1147:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1148:                                              PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
1149: {

1156:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1157:   PetscDSEnlarge_Static(prob, f+1);
1158:   prob->f[f*2+0] = f0;
1159:   prob->f[f*2+1] = f1;
1160:   return(0);
1161: }

1163: /*@C
1164:   PetscDSHasJacobian - Signals that Jacobian functions have been set

1166:   Not collective

1168:   Input Parameter:
1169: . prob - The PetscDS

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

1174:   Level: intermediate

1176: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1177: @*/
1178: PetscErrorCode PetscDSHasJacobian(PetscDS prob, PetscBool *hasJac)
1179: {
1180:   PetscInt f, g, h;

1184:   *hasJac = PETSC_FALSE;
1185:   for (f = 0; f < prob->Nf; ++f) {
1186:     for (g = 0; g < prob->Nf; ++g) {
1187:       for (h = 0; h < 4; ++h) {
1188:         if (prob->g[(f*prob->Nf + g)*4+h]) *hasJac = PETSC_TRUE;
1189:       }
1190:     }
1191:   }
1192:   return(0);
1193: }

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

1198:   Not collective

1200:   Input Parameters:
1201: + prob - The PetscDS
1202: . f    - The test field number
1203: - g    - The field number

1205:   Output Parameters:
1206: + g0 - integrand for the test and basis function term
1207: . g1 - integrand for the test function and basis function gradient term
1208: . g2 - integrand for the test function gradient and basis function term
1209: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

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

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

1241:   Level: intermediate

1243: .seealso: PetscDSSetJacobian()
1244: @*/
1245: PetscErrorCode PetscDSGetJacobian(PetscDS prob, PetscInt f, PetscInt g,
1246:                                   void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1247:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1248:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1249:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1250:                                   void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1251:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1252:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1253:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1254:                                   void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1255:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1256:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1257:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1258:                                   void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1259:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1260:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1261:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1262: {
1265:   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);
1266:   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);
1271:   return(0);
1272: }

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

1277:   Not collective

1279:   Input Parameters:
1280: + prob - The PetscDS
1281: . f    - The test field number
1282: . g    - The field number
1283: . g0 - integrand for the test and basis function term
1284: . g1 - integrand for the test function and basis function gradient term
1285: . g2 - integrand for the test function gradient and basis function term
1286: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

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

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

1318:   Level: intermediate

1320: .seealso: PetscDSGetJacobian()
1321: @*/
1322: PetscErrorCode PetscDSSetJacobian(PetscDS prob, PetscInt f, PetscInt g,
1323:                                   void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1324:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1325:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1326:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1327:                                   void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1328:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1329:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1330:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1331:                                   void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1332:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1333:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1334:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1335:                                   void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1336:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1337:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1338:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1339: {

1348:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1349:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
1350:   PetscDSEnlarge_Static(prob, PetscMax(f, g)+1);
1351:   prob->g[(f*prob->Nf + g)*4+0] = g0;
1352:   prob->g[(f*prob->Nf + g)*4+1] = g1;
1353:   prob->g[(f*prob->Nf + g)*4+2] = g2;
1354:   prob->g[(f*prob->Nf + g)*4+3] = g3;
1355:   return(0);
1356: }

1358: /*@C
1359:   PetscDSUseJacobianPreconditioner - Whether to construct a Jacobian preconditioner

1361:   Not collective

1363:   Input Parameters:
1364: + prob - The PetscDS
1365: - useJacPre - flag that enables construction of a Jacobian preconditioner

1367:   Level: intermediate

1369: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1370: @*/
1371: PetscErrorCode PetscDSUseJacobianPreconditioner(PetscDS prob, PetscBool useJacPre)
1372: {
1375:   prob->useJacPre = useJacPre;
1376:   return(0);
1377: }

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

1382:   Not collective

1384:   Input Parameter:
1385: . prob - The PetscDS

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

1390:   Level: intermediate

1392: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1393: @*/
1394: PetscErrorCode PetscDSHasJacobianPreconditioner(PetscDS prob, PetscBool *hasJacPre)
1395: {
1396:   PetscInt f, g, h;

1400:   *hasJacPre = PETSC_FALSE;
1401:   if (!prob->useJacPre) return(0);
1402:   for (f = 0; f < prob->Nf; ++f) {
1403:     for (g = 0; g < prob->Nf; ++g) {
1404:       for (h = 0; h < 4; ++h) {
1405:         if (prob->gp[(f*prob->Nf + g)*4+h]) *hasJacPre = PETSC_TRUE;
1406:       }
1407:     }
1408:   }
1409:   return(0);
1410: }

1412: /*@C
1413:   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.

1415:   Not collective

1417:   Input Parameters:
1418: + prob - The PetscDS
1419: . f    - The test field number
1420: - g    - The field number

1422:   Output Parameters:
1423: + g0 - integrand for the test and basis function term
1424: . g1 - integrand for the test function and basis function gradient term
1425: . g2 - integrand for the test function gradient and basis function term
1426: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

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

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

1458:   Level: intermediate

1460: .seealso: PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1461: @*/
1462: PetscErrorCode PetscDSGetJacobianPreconditioner(PetscDS prob, PetscInt f, PetscInt g,
1463:                                   void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1464:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1465:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1466:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1467:                                   void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1468:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1469:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1470:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1471:                                   void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1472:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1473:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1474:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1475:                                   void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1476:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1477:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1478:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1479: {
1482:   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);
1483:   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);
1488:   return(0);
1489: }

1491: /*@C
1492:   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.

1494:   Not collective

1496:   Input Parameters:
1497: + prob - The PetscDS
1498: . f    - The test field number
1499: . g    - The field number
1500: . g0 - integrand for the test and basis function term
1501: . g1 - integrand for the test function and basis function gradient term
1502: . g2 - integrand for the test function gradient and basis function term
1503: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

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

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

1535:   Level: intermediate

1537: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobian()
1538: @*/
1539: PetscErrorCode PetscDSSetJacobianPreconditioner(PetscDS prob, PetscInt f, PetscInt g,
1540:                                   void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1541:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1542:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1543:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1544:                                   void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1545:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1546:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1547:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1548:                                   void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1549:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1550:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1551:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1552:                                   void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1553:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1554:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1555:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1556: {

1565:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1566:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
1567:   PetscDSEnlarge_Static(prob, PetscMax(f, g)+1);
1568:   prob->gp[(f*prob->Nf + g)*4+0] = g0;
1569:   prob->gp[(f*prob->Nf + g)*4+1] = g1;
1570:   prob->gp[(f*prob->Nf + g)*4+2] = g2;
1571:   prob->gp[(f*prob->Nf + g)*4+3] = g3;
1572:   return(0);
1573: }

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

1578:   Not collective

1580:   Input Parameter:
1581: . prob - The PetscDS

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

1586:   Level: intermediate

1588: .seealso: PetscDSGetDynamicJacobian(), PetscDSSetDynamicJacobian(), PetscDSGetJacobian()
1589: @*/
1590: PetscErrorCode PetscDSHasDynamicJacobian(PetscDS prob, PetscBool *hasDynJac)
1591: {
1592:   PetscInt f, g, h;

1596:   *hasDynJac = PETSC_FALSE;
1597:   for (f = 0; f < prob->Nf; ++f) {
1598:     for (g = 0; g < prob->Nf; ++g) {
1599:       for (h = 0; h < 4; ++h) {
1600:         if (prob->gt[(f*prob->Nf + g)*4+h]) *hasDynJac = PETSC_TRUE;
1601:       }
1602:     }
1603:   }
1604:   return(0);
1605: }

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

1610:   Not collective

1612:   Input Parameters:
1613: + prob - The PetscDS
1614: . f    - The test field number
1615: - g    - The field number

1617:   Output Parameters:
1618: + g0 - integrand for the test and basis function term
1619: . g1 - integrand for the test function and basis function gradient term
1620: . g2 - integrand for the test function gradient and basis function term
1621: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

1629: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1630: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1631: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1632: $    PetscReal t, const PetscReal u_tShift, const PetscReal x[], PetscScalar g0[])

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

1653:   Level: intermediate

1655: .seealso: PetscDSSetJacobian()
1656: @*/
1657: PetscErrorCode PetscDSGetDynamicJacobian(PetscDS prob, PetscInt f, PetscInt g,
1658:                                          void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1659:                                                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1660:                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1661:                                                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1662:                                          void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1663:                                                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1664:                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1665:                                                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1666:                                          void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1667:                                                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1668:                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1669:                                                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1670:                                          void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1671:                                                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1672:                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1673:                                                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1674: {
1677:   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);
1678:   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);
1683:   return(0);
1684: }

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

1689:   Not collective

1691:   Input Parameters:
1692: + prob - The PetscDS
1693: . f    - The test field number
1694: . g    - The field number
1695: . g0 - integrand for the test and basis function term
1696: . g1 - integrand for the test function and basis function gradient term
1697: . g2 - integrand for the test function gradient and basis function term
1698: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

1706: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1707: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1708: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1709: $    PetscReal t, const PetscReal x[], PetscScalar g0[])

1711: + dim - the spatial dimension
1712: . Nf - the number of fields
1713: . uOff - the offset into u[] and u_t[] for each field
1714: . uOff_x - the offset into u_x[] for each field
1715: . u - each field evaluated at the current point
1716: . u_t - the time derivative of each field evaluated at the current point
1717: . u_x - the gradient of each field evaluated at the current point
1718: . aOff - the offset into a[] and a_t[] for each auxiliary field
1719: . aOff_x - the offset into a_x[] for each auxiliary field
1720: . a - each auxiliary field evaluated at the current point
1721: . a_t - the time derivative of each auxiliary field evaluated at the current point
1722: . a_x - the gradient of auxiliary each field evaluated at the current point
1723: . t - current time
1724: . u_tShift - the multiplier a for dF/dU_t
1725: . x - coordinates of the current point
1726: . numConstants - number of constant parameters
1727: . constants - constant parameters
1728: - g0 - output values at the current point

1730:   Level: intermediate

1732: .seealso: PetscDSGetJacobian()
1733: @*/
1734: PetscErrorCode PetscDSSetDynamicJacobian(PetscDS prob, PetscInt f, PetscInt g,
1735:                                          void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1736:                                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1737:                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1738:                                                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1739:                                          void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1740:                                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1741:                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1742:                                                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1743:                                          void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1744:                                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1745:                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1746:                                                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1747:                                          void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1748:                                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1749:                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1750:                                                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1751: {

1760:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1761:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
1762:   PetscDSEnlarge_Static(prob, PetscMax(f, g)+1);
1763:   prob->gt[(f*prob->Nf + g)*4+0] = g0;
1764:   prob->gt[(f*prob->Nf + g)*4+1] = g1;
1765:   prob->gt[(f*prob->Nf + g)*4+2] = g2;
1766:   prob->gt[(f*prob->Nf + g)*4+3] = g3;
1767:   return(0);
1768: }

1770: /*@C
1771:   PetscDSGetRiemannSolver - Returns the Riemann solver for the given field

1773:   Not collective

1775:   Input Arguments:
1776: + prob - The PetscDS object
1777: - f    - The field number

1779:   Output Argument:
1780: . r    - Riemann solver

1782:   Calling sequence for r:

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

1786: + dim  - The spatial dimension
1787: . Nf   - The number of fields
1788: . x    - The coordinates at a point on the interface
1789: . n    - The normal vector to the interface
1790: . uL   - The state vector to the left of the interface
1791: . uR   - The state vector to the right of the interface
1792: . flux - output array of flux through the interface
1793: . numConstants - number of constant parameters
1794: . constants - constant parameters
1795: - ctx  - optional user context

1797:   Level: intermediate

1799: .seealso: PetscDSSetRiemannSolver()
1800: @*/
1801: PetscErrorCode PetscDSGetRiemannSolver(PetscDS prob, PetscInt f,
1802:                                        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))
1803: {
1806:   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);
1808:   *r = prob->r[f];
1809:   return(0);
1810: }

1812: /*@C
1813:   PetscDSSetRiemannSolver - Sets the Riemann solver for the given field

1815:   Not collective

1817:   Input Arguments:
1818: + prob - The PetscDS object
1819: . f    - The field number
1820: - r    - Riemann solver

1822:   Calling sequence for r:

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

1826: + dim  - The spatial dimension
1827: . Nf   - The number of fields
1828: . x    - The coordinates at a point on the interface
1829: . n    - The normal vector to the interface
1830: . uL   - The state vector to the left of the interface
1831: . uR   - The state vector to the right of the interface
1832: . flux - output array of flux through the interface
1833: . numConstants - number of constant parameters
1834: . constants - constant parameters
1835: - ctx  - optional user context

1837:   Level: intermediate

1839: .seealso: PetscDSGetRiemannSolver()
1840: @*/
1841: PetscErrorCode PetscDSSetRiemannSolver(PetscDS prob, PetscInt f,
1842:                                        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))
1843: {

1849:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1850:   PetscDSEnlarge_Static(prob, f+1);
1851:   prob->r[f] = r;
1852:   return(0);
1853: }

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

1858:   Not collective

1860:   Input Parameters:
1861: + prob - The PetscDS
1862: - f    - The field number

1864:   Output Parameters:
1865: . update - update function

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

1869: $ update(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1870: $        const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1871: $        const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1872: $        PetscReal t, const PetscReal x[], PetscScalar uNew[])

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

1890:   Level: intermediate

1892: .seealso: PetscDSSetUpdate(), PetscDSSetResidual()
1893: @*/
1894: PetscErrorCode PetscDSGetUpdate(PetscDS prob, PetscInt f,
1895:                                   void (**update)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1896:                                                   const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1897:                                                   const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1898:                                                   PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar uNew[]))
1899: {
1902:   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);
1904:   return(0);
1905: }

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

1910:   Not collective

1912:   Input Parameters:
1913: + prob   - The PetscDS
1914: . f      - The field number
1915: - update - update function

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

1919: $ update(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1920: $        const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1921: $        const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1922: $        PetscReal t, const PetscReal x[], PetscScalar uNew[])

1924: + dim - the spatial dimension
1925: . Nf - the number of fields
1926: . uOff - the offset into u[] and u_t[] for each field
1927: . uOff_x - the offset into u_x[] for each field
1928: . u - each field evaluated at the current point
1929: . u_t - the time derivative of each field evaluated at the current point
1930: . u_x - the gradient of each field evaluated at the current point
1931: . aOff - the offset into a[] and a_t[] for each auxiliary field
1932: . aOff_x - the offset into a_x[] for each auxiliary field
1933: . a - each auxiliary field evaluated at the current point
1934: . a_t - the time derivative of each auxiliary field evaluated at the current point
1935: . a_x - the gradient of auxiliary each field evaluated at the current point
1936: . t - current time
1937: . x - coordinates of the current point
1938: - uNew - new field values at the current point

1940:   Level: intermediate

1942: .seealso: PetscDSGetResidual()
1943: @*/
1944: PetscErrorCode PetscDSSetUpdate(PetscDS prob, PetscInt f,
1945:                                 void (*update)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1946:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1947:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1948:                                                PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar uNew[]))
1949: {

1955:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1956:   PetscDSEnlarge_Static(prob, f+1);
1957:   prob->update[f] = update;
1958:   return(0);
1959: }

1961: PetscErrorCode PetscDSGetContext(PetscDS prob, PetscInt f, void **ctx)
1962: {
1965:   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);
1967:   *ctx = prob->ctx[f];
1968:   return(0);
1969: }

1971: PetscErrorCode PetscDSSetContext(PetscDS prob, PetscInt f, void *ctx)
1972: {

1977:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1978:   PetscDSEnlarge_Static(prob, f+1);
1979:   prob->ctx[f] = ctx;
1980:   return(0);
1981: }

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

1986:   Not collective

1988:   Input Parameters:
1989: + prob - The PetscDS
1990: - f    - The test field number

1992:   Output Parameters:
1993: + f0 - boundary integrand for the test function term
1994: - f1 - boundary integrand for the test function gradient term

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

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

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

2002: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2003: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2004: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2005: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[])

2007: + dim - the spatial dimension
2008: . Nf - the number of fields
2009: . uOff - the offset into u[] and u_t[] for each field
2010: . uOff_x - the offset into u_x[] for each field
2011: . u - each field evaluated at the current point
2012: . u_t - the time derivative of each field evaluated at the current point
2013: . u_x - the gradient of each field evaluated at the current point
2014: . aOff - the offset into a[] and a_t[] for each auxiliary field
2015: . aOff_x - the offset into a_x[] for each auxiliary field
2016: . a - each auxiliary field evaluated at the current point
2017: . a_t - the time derivative of each auxiliary field evaluated at the current point
2018: . a_x - the gradient of auxiliary each field evaluated at the current point
2019: . t - current time
2020: . x - coordinates of the current point
2021: . n - unit normal at the current point
2022: . numConstants - number of constant parameters
2023: . constants - constant parameters
2024: - f0 - output values at the current point

2026:   Level: intermediate

2028: .seealso: PetscDSSetBdResidual()
2029: @*/
2030: PetscErrorCode PetscDSGetBdResidual(PetscDS prob, PetscInt f,
2031:                                     void (**f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2032:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2033:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2034:                                                 PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]),
2035:                                     void (**f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2036:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2037:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2038:                                                 PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
2039: {
2042:   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);
2045:   return(0);
2046: }

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

2051:   Not collective

2053:   Input Parameters:
2054: + prob - The PetscDS
2055: . f    - The test field number
2056: . f0 - boundary integrand for the test function term
2057: - f1 - boundary integrand for the test function gradient term

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

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

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

2065: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2066: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2067: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2068: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[])

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

2089:   Level: intermediate

2091: .seealso: PetscDSGetBdResidual()
2092: @*/
2093: PetscErrorCode PetscDSSetBdResidual(PetscDS prob, PetscInt f,
2094:                                     void (*f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2095:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2096:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2097:                                                PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]),
2098:                                     void (*f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2099:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2100:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2101:                                                PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
2102: {

2107:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2108:   PetscDSEnlarge_Static(prob, f+1);
2111:   return(0);
2112: }

2114: /*@
2115:   PetscDSHasBdJacobian - Signals that boundary Jacobian functions have been set

2117:   Not collective

2119:   Input Parameter:
2120: . prob - The PetscDS

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

2125:   Level: intermediate

2127: .seealso: PetscDSHasJacobian(), PetscDSSetBdJacobian(), PetscDSGetBdJacobian()
2128: @*/
2129: PetscErrorCode PetscDSHasBdJacobian(PetscDS prob, PetscBool *hasBdJac)
2130: {
2131:   PetscInt f, g, h;

2135:   *hasBdJac = PETSC_FALSE;
2136:   for (f = 0; f < prob->Nf; ++f) {
2137:     for (g = 0; g < prob->Nf; ++g) {
2138:       for (h = 0; h < 4; ++h) {
2139:         if (prob->gBd[(f*prob->Nf + g)*4+h]) *hasBdJac = PETSC_TRUE;
2140:       }
2141:     }
2142:   }
2143:   return(0);
2144: }

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

2149:   Not collective

2151:   Input Parameters:
2152: + prob - The PetscDS
2153: . f    - The test field number
2154: - g    - The field number

2156:   Output Parameters:
2157: + g0 - integrand for the test and basis function term
2158: . g1 - integrand for the test function and basis function gradient term
2159: . g2 - integrand for the test function gradient and basis function term
2160: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

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

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

2193:   Level: intermediate

2195: .seealso: PetscDSSetBdJacobian()
2196: @*/
2197: PetscErrorCode PetscDSGetBdJacobian(PetscDS prob, PetscInt f, PetscInt g,
2198:                                     void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2199:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2200:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2201:                                                 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
2202:                                     void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2203:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2204:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2205:                                                 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
2206:                                     void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2207:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2208:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2209:                                                 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
2210:                                     void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2211:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2212:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2213:                                                 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
2214: {
2217:   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);
2218:   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);
2223:   return(0);
2224: }

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

2229:   Not collective

2231:   Input Parameters:
2232: + prob - The PetscDS
2233: . f    - The test field number
2234: . g    - The field number
2235: . g0 - integrand for the test and basis function term
2236: . g1 - integrand for the test function and basis function gradient term
2237: . g2 - integrand for the test function gradient and basis function term
2238: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

2246: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2247: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2248: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2249: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar g0[])

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

2271:   Level: intermediate

2273: .seealso: PetscDSGetBdJacobian()
2274: @*/
2275: PetscErrorCode PetscDSSetBdJacobian(PetscDS prob, PetscInt f, PetscInt g,
2276:                                     void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2277:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2278:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2279:                                                PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
2280:                                     void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2281:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2282:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2283:                                                PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
2284:                                     void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2285:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2286:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2287:                                                PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
2288:                                     void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2289:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2290:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2291:                                                PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
2292: {

2301:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2302:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
2303:   PetscDSEnlarge_Static(prob, PetscMax(f, g)+1);
2304:   prob->gBd[(f*prob->Nf + g)*4+0] = g0;
2305:   prob->gBd[(f*prob->Nf + g)*4+1] = g1;
2306:   prob->gBd[(f*prob->Nf + g)*4+2] = g2;
2307:   prob->gBd[(f*prob->Nf + g)*4+3] = g3;
2308:   return(0);
2309: }

2311: /*@
2312:   PetscDSHasBdJacobianPreconditioner - Signals that boundary Jacobian preconditioner functions have been set

2314:   Not collective

2316:   Input Parameter:
2317: . prob - The PetscDS

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

2322:   Level: intermediate

2324: .seealso: PetscDSHasJacobian(), PetscDSSetBdJacobian(), PetscDSGetBdJacobian()
2325: @*/
2326: PetscErrorCode PetscDSHasBdJacobianPreconditioner(PetscDS prob, PetscBool *hasBdJacPre)
2327: {
2328:   PetscInt f, g, h;

2332:   *hasBdJacPre = PETSC_FALSE;
2333:   for (f = 0; f < prob->Nf; ++f) {
2334:     for (g = 0; g < prob->Nf; ++g) {
2335:       for (h = 0; h < 4; ++h) {
2336:         if (prob->gpBd[(f*prob->Nf + g)*4+h]) *hasBdJacPre = PETSC_TRUE;
2337:       }
2338:     }
2339:   }
2340:   return(0);
2341: }

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

2346:   Not collective

2348:   Input Parameters:
2349: + prob - The PetscDS
2350: . f    - The test field number
2351: - g    - The field number

2353:   Output Parameters:
2354: + g0 - integrand for the test and basis function term
2355: . g1 - integrand for the test function and basis function gradient term
2356: . g2 - integrand for the test function gradient and basis function term
2357: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

2365: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2366: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2367: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2368: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])

2370: + dim - the spatial dimension
2371: . Nf - the number of fields
2372: . NfAux - the number of auxiliary fields
2373: . uOff - the offset into u[] and u_t[] for each field
2374: . uOff_x - the offset into u_x[] for each field
2375: . u - each field evaluated at the current point
2376: . u_t - the time derivative of each field evaluated at the current point
2377: . u_x - the gradient of each field evaluated at the current point
2378: . aOff - the offset into a[] and a_t[] for each auxiliary field
2379: . aOff_x - the offset into a_x[] for each auxiliary field
2380: . a - each auxiliary field evaluated at the current point
2381: . a_t - the time derivative of each auxiliary field evaluated at the current point
2382: . a_x - the gradient of auxiliary each field evaluated at the current point
2383: . t - current time
2384: . u_tShift - the multiplier a for dF/dU_t
2385: . x - coordinates of the current point
2386: . n - normal at the current point
2387: . numConstants - number of constant parameters
2388: . constants - constant parameters
2389: - g0 - output values at the current point

2391:   This is not yet available in Fortran.

2393:   Level: intermediate

2395: .seealso: PetscDSSetBdJacobianPreconditioner()
2396: @*/
2397: PetscErrorCode PetscDSGetBdJacobianPreconditioner(PetscDS prob, PetscInt f, PetscInt g,
2398:                                                   void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2399:                                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2400:                                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2401:                                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
2402:                                                   void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2403:                                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2404:                                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2405:                                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
2406:                                                   void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2407:                                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2408:                                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2409:                                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
2410:                                                   void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2411:                                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2412:                                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2413:                                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
2414: {
2417:   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);
2418:   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);
2423:   return(0);
2424: }

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

2429:   Not collective

2431:   Input Parameters:
2432: + prob - The PetscDS
2433: . f    - The test field number
2434: . g    - The field number
2435: . g0 - integrand for the test and basis function term
2436: . g1 - integrand for the test function and basis function gradient term
2437: . g2 - integrand for the test function gradient and basis function term
2438: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

2446: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2447: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2448: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2449: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])

2451: + dim - the spatial dimension
2452: . Nf - the number of fields
2453: . NfAux - the number of auxiliary fields
2454: . uOff - the offset into u[] and u_t[] for each field
2455: . uOff_x - the offset into u_x[] for each field
2456: . u - each field evaluated at the current point
2457: . u_t - the time derivative of each field evaluated at the current point
2458: . u_x - the gradient of each field evaluated at the current point
2459: . aOff - the offset into a[] and a_t[] for each auxiliary field
2460: . aOff_x - the offset into a_x[] for each auxiliary field
2461: . a - each auxiliary field evaluated at the current point
2462: . a_t - the time derivative of each auxiliary field evaluated at the current point
2463: . a_x - the gradient of auxiliary each field evaluated at the current point
2464: . t - current time
2465: . u_tShift - the multiplier a for dF/dU_t
2466: . x - coordinates of the current point
2467: . n - normal at the current point
2468: . numConstants - number of constant parameters
2469: . constants - constant parameters
2470: - g0 - output values at the current point

2472:   This is not yet available in Fortran.

2474:   Level: intermediate

2476: .seealso: PetscDSGetBdJacobianPreconditioner()
2477: @*/
2478: PetscErrorCode PetscDSSetBdJacobianPreconditioner(PetscDS prob, PetscInt f, PetscInt g,
2479:                                                   void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2480:                                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2481:                                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2482:                                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
2483:                                                   void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2484:                                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2485:                                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2486:                                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
2487:                                                   void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2488:                                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2489:                                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2490:                                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
2491:                                                   void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2492:                                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2493:                                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2494:                                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
2495: {

2504:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2505:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
2506:   PetscDSEnlarge_Static(prob, PetscMax(f, g)+1);
2507:   prob->gpBd[(f*prob->Nf + g)*4+0] = g0;
2508:   prob->gpBd[(f*prob->Nf + g)*4+1] = g1;
2509:   prob->gpBd[(f*prob->Nf + g)*4+2] = g2;
2510:   prob->gpBd[(f*prob->Nf + g)*4+3] = g3;
2511:   return(0);
2512: }

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

2517:   Not collective

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

2523:   Output Parameter:
2524: + exactSol - exact solution for the test field
2525: - exactCtx - exact solution context

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

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

2531: + dim - the spatial dimension
2532: . t - current time
2533: . x - coordinates of the current point
2534: . Nc - the number of field components
2535: . u - the solution field evaluated at the current point
2536: - ctx - a user context

2538:   Level: intermediate

2540: .seealso: PetscDSSetExactSolution(), PetscDSGetExactSolutionTimeDerivative()
2541: @*/
2542: PetscErrorCode PetscDSGetExactSolution(PetscDS prob, PetscInt f, PetscErrorCode (**sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void **ctx)
2543: {
2546:   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);
2549:   return(0);
2550: }

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

2555:   Not collective

2557:   Input Parameters:
2558: + prob - The PetscDS
2559: . f    - The test field number
2560: . sol  - solution function for the test fields
2561: - ctx  - solution context or NULL

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

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

2567: + dim - the spatial dimension
2568: . t - current time
2569: . x - coordinates of the current point
2570: . Nc - the number of field components
2571: . u - the solution field evaluated at the current point
2572: - ctx - a user context

2574:   Level: intermediate

2576: .seealso: PetscDSGetExactSolution()
2577: @*/
2578: PetscErrorCode PetscDSSetExactSolution(PetscDS prob, PetscInt f, PetscErrorCode (*sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void *ctx)
2579: {

2584:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2585:   PetscDSEnlarge_Static(prob, f+1);
2588:   return(0);
2589: }

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

2594:   Not collective

2596:   Input Parameters:
2597: + prob - The PetscDS
2598: - f    - The test field number

2600:   Output Parameter:
2601: + exactSol - time derivative of the exact solution for the test field
2602: - exactCtx - time derivative of the exact solution context

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

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

2608: + dim - the spatial dimension
2609: . t - current time
2610: . x - coordinates of the current point
2611: . Nc - the number of field components
2612: . u - the solution field evaluated at the current point
2613: - ctx - a user context

2615:   Level: intermediate

2617: .seealso: PetscDSSetExactSolutionTimeDerivative(), PetscDSGetExactSolution()
2618: @*/
2619: PetscErrorCode PetscDSGetExactSolutionTimeDerivative(PetscDS prob, PetscInt f, PetscErrorCode (**sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void **ctx)
2620: {
2623:   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);
2626:   return(0);
2627: }

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

2632:   Not collective

2634:   Input Parameters:
2635: + prob - The PetscDS
2636: . f    - The test field number
2637: . sol  - time derivative of the solution function for the test fields
2638: - ctx  - time derivative of the solution context or NULL

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

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

2644: + dim - the spatial dimension
2645: . t - current time
2646: . x - coordinates of the current point
2647: . Nc - the number of field components
2648: . u - the solution field evaluated at the current point
2649: - ctx - a user context

2651:   Level: intermediate

2653: .seealso: PetscDSGetExactSolutionTimeDerivative(), PetscDSSetExactSolution()
2654: @*/
2655: PetscErrorCode PetscDSSetExactSolutionTimeDerivative(PetscDS prob, PetscInt f, PetscErrorCode (*sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void *ctx)
2656: {

2661:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2662:   PetscDSEnlarge_Static(prob, f+1);
2665:   return(0);
2666: }

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

2671:   Not collective

2673:   Input Parameter:
2674: . prob - The PetscDS object

2676:   Output Parameters:
2677: + numConstants - The number of constants
2678: - constants    - The array of constants, NULL if there are none

2680:   Level: intermediate

2682: .seealso: PetscDSSetConstants(), PetscDSCreate()
2683: @*/
2684: PetscErrorCode PetscDSGetConstants(PetscDS prob, PetscInt *numConstants, const PetscScalar *constants[])
2685: {
2690:   return(0);
2691: }

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

2696:   Not collective

2698:   Input Parameters:
2699: + prob         - The PetscDS object
2700: . numConstants - The number of constants
2701: - constants    - The array of constants, NULL if there are none

2703:   Level: intermediate

2705: .seealso: PetscDSGetConstants(), PetscDSCreate()
2706: @*/
2707: PetscErrorCode PetscDSSetConstants(PetscDS prob, PetscInt numConstants, PetscScalar constants[])
2708: {

2713:   if (numConstants != prob->numConstants) {
2714:     PetscFree(prob->constants);
2715:     prob->numConstants = numConstants;
2716:     if (prob->numConstants) {
2717:       PetscMalloc1(prob->numConstants, &prob->constants);
2718:     } else {
2719:       prob->constants = NULL;
2720:     }
2721:   }
2722:   if (prob->numConstants) {
2724:     PetscArraycpy(prob->constants, constants, prob->numConstants);
2725:   }
2726:   return(0);
2727: }

2729: /*@
2730:   PetscDSGetFieldIndex - Returns the index of the given field

2732:   Not collective

2734:   Input Parameters:
2735: + prob - The PetscDS object
2736: - disc - The discretization object

2738:   Output Parameter:
2739: . f - The field number

2741:   Level: beginner

2743: .seealso: PetscGetDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
2744: @*/
2745: PetscErrorCode PetscDSGetFieldIndex(PetscDS prob, PetscObject disc, PetscInt *f)
2746: {
2747:   PetscInt g;

2752:   *f = -1;
2753:   for (g = 0; g < prob->Nf; ++g) {if (disc == prob->disc[g]) break;}
2754:   if (g == prob->Nf) SETERRQ(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Field not found in PetscDS.");
2755:   *f = g;
2756:   return(0);
2757: }

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

2762:   Not collective

2764:   Input Parameters:
2765: + prob - The PetscDS object
2766: - f - The field number

2768:   Output Parameter:
2769: . size - The size

2771:   Level: beginner

2773: .seealso: PetscDSGetFieldOffset(), PetscDSGetNumFields(), PetscDSCreate()
2774: @*/
2775: PetscErrorCode PetscDSGetFieldSize(PetscDS prob, PetscInt f, PetscInt *size)
2776: {

2782:   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);
2783:   PetscDSSetUp(prob);
2784:   *size = prob->Nb[f];
2785:   return(0);
2786: }

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

2791:   Not collective

2793:   Input Parameters:
2794: + prob - The PetscDS object
2795: - f - The field number

2797:   Output Parameter:
2798: . off - The offset

2800:   Level: beginner

2802: .seealso: PetscDSGetFieldSize(), PetscDSGetNumFields(), PetscDSCreate()
2803: @*/
2804: PetscErrorCode PetscDSGetFieldOffset(PetscDS prob, PetscInt f, PetscInt *off)
2805: {
2806:   PetscInt       size, g;

2812:   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);
2813:   *off = 0;
2814:   for (g = 0; g < f; ++g) {
2815:     PetscDSGetFieldSize(prob, g, &size);
2816:     *off += size;
2817:   }
2818:   return(0);
2819: }

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

2824:   Not collective

2826:   Input Parameter:
2827: . prob - The PetscDS object

2829:   Output Parameter:
2830: . dimensions - The number of dimensions

2832:   Level: beginner

2834: .seealso: PetscDSGetComponentOffsets(), PetscDSGetNumFields(), PetscDSCreate()
2835: @*/
2836: PetscErrorCode PetscDSGetDimensions(PetscDS prob, PetscInt *dimensions[])
2837: {

2842:   PetscDSSetUp(prob);
2844:   *dimensions = prob->Nb;
2845:   return(0);
2846: }

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

2851:   Not collective

2853:   Input Parameter:
2854: . prob - The PetscDS object

2856:   Output Parameter:
2857: . components - The number of components

2859:   Level: beginner

2861: .seealso: PetscDSGetComponentOffsets(), PetscDSGetNumFields(), PetscDSCreate()
2862: @*/
2863: PetscErrorCode PetscDSGetComponents(PetscDS prob, PetscInt *components[])
2864: {

2869:   PetscDSSetUp(prob);
2871:   *components = prob->Nc;
2872:   return(0);
2873: }

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

2878:   Not collective

2880:   Input Parameters:
2881: + prob - The PetscDS object
2882: - f - The field number

2884:   Output Parameter:
2885: . off - The offset

2887:   Level: beginner

2889: .seealso: PetscDSGetNumFields(), PetscDSCreate()
2890: @*/
2891: PetscErrorCode PetscDSGetComponentOffset(PetscDS prob, PetscInt f, PetscInt *off)
2892: {
2896:   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);
2897:   *off = prob->off[f];
2898:   return(0);
2899: }

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

2904:   Not collective

2906:   Input Parameter:
2907: . prob - The PetscDS object

2909:   Output Parameter:
2910: . offsets - The offsets

2912:   Level: beginner

2914: .seealso: PetscDSGetNumFields(), PetscDSCreate()
2915: @*/
2916: PetscErrorCode PetscDSGetComponentOffsets(PetscDS prob, PetscInt *offsets[])
2917: {
2921:   *offsets = prob->off;
2922:   return(0);
2923: }

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

2928:   Not collective

2930:   Input Parameter:
2931: . prob - The PetscDS object

2933:   Output Parameter:
2934: . offsets - The offsets

2936:   Level: beginner

2938: .seealso: PetscDSGetNumFields(), PetscDSCreate()
2939: @*/
2940: PetscErrorCode PetscDSGetComponentDerivativeOffsets(PetscDS prob, PetscInt *offsets[])
2941: {
2945:   *offsets = prob->offDer;
2946:   return(0);
2947: }

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

2952:   Not collective

2954:   Input Parameter:
2955: . prob - The PetscDS object

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

2960:   Level: intermediate

2962: .seealso: PetscDSCreate()
2963: @*/
2964: PetscErrorCode PetscDSGetTabulation(PetscDS prob, PetscTabulation *T[])
2965: {

2971:   PetscDSSetUp(prob);
2972:   *T = prob->T;
2973:   return(0);
2974: }

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

2979:   Not collective

2981:   Input Parameter:
2982: . prob - The PetscDS object

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

2987:   Level: intermediate

2989: .seealso: PetscDSGetTabulation(), PetscDSCreate()
2990: @*/
2991: PetscErrorCode PetscDSGetFaceTabulation(PetscDS prob, PetscTabulation *Tf[])
2992: {

2998:   PetscDSSetUp(prob);
2999:   *Tf = prob->Tf;
3000:   return(0);
3001: }

3003: PetscErrorCode PetscDSGetEvaluationArrays(PetscDS prob, PetscScalar **u, PetscScalar **u_t, PetscScalar **u_x)
3004: {

3009:   PetscDSSetUp(prob);
3013:   return(0);
3014: }

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

3022:   PetscDSSetUp(prob);
3029:   return(0);
3030: }

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

3038:   PetscDSSetUp(prob);
3044:   return(0);
3045: }

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

3050:   Collective on ds

3052:   Input Parameters:
3053: + ds          - The PetscDS object
3054: . type        - The type of condition, e.g. DM_BC_ESSENTIAL/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
3055: . name        - The BC name
3056: . labelname   - The label defining constrained points
3057: . field       - The field to constrain
3058: . numcomps    - The number of constrained field components (0 will constrain all fields)
3059: . comps       - An array of constrained component numbers
3060: . bcFunc      - A pointwise function giving boundary values
3061: . bcFunc_t    - A pointwise function giving the time derviative of the boundary values, or NULL
3062: . numids      - The number of DMLabel ids for constrained points
3063: . ids         - An array of ids for constrained points
3064: - ctx         - An optional user context for bcFunc

3066:   Options Database Keys:
3067: + -bc_<boundary name> <num> - Overrides the boundary ids
3068: - -bc_<boundary name>_comp <num> - Overrides the boundary components

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

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

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

3077: $ bcFunc(PetscInt dim, PetscInt Nf, PetscInt NfAux,
3078: $        const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
3079: $        const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
3080: $        PetscReal time, const PetscReal x[], PetscScalar bcval[])

3082: + dim - the spatial dimension
3083: . Nf - the number of fields
3084: . uOff - the offset into u[] and u_t[] for each field
3085: . uOff_x - the offset into u_x[] for each field
3086: . u - each field evaluated at the current point
3087: . u_t - the time derivative of each field evaluated at the current point
3088: . u_x - the gradient of each field evaluated at the current point
3089: . aOff - the offset into a[] and a_t[] for each auxiliary field
3090: . aOff_x - the offset into a_x[] for each auxiliary field
3091: . a - each auxiliary field evaluated at the current point
3092: . a_t - the time derivative of each auxiliary field evaluated at the current point
3093: . a_x - the gradient of auxiliary each field evaluated at the current point
3094: . t - current time
3095: . x - coordinates of the current point
3096: . numConstants - number of constant parameters
3097: . constants - constant parameters
3098: - bcval - output values at the current point

3100:   Level: developer

3102: .seealso: PetscDSGetBoundary(), PetscDSSetResidual(), PetscDSSetBdResidual()
3103: @*/
3104: PetscErrorCode PetscDSAddBoundary(PetscDS ds, DMBoundaryConditionType type, const char name[], const char labelname[], PetscInt field, PetscInt numcomps, const PetscInt *comps, void (*bcFunc)(void), void (*bcFunc_t)(void), PetscInt numids, const PetscInt *ids, void *ctx)
3105: {
3106:   DSBoundary     b;

3115:   PetscNew(&b);
3116:   PetscStrallocpy(name, (char **) &b->name);
3117:   PetscStrallocpy(labelname, (char **) &b->labelname);
3118:   PetscMalloc1(numcomps, &b->comps);
3119:   if (numcomps) {PetscArraycpy(b->comps, comps, numcomps);}
3120:   PetscMalloc1(numids, &b->ids);
3121:   if (numids) {PetscArraycpy(b->ids, ids, numids);}
3122:   b->type            = type;
3123:   b->field           = field;
3124:   b->numcomps        = numcomps;
3125:   b->func            = bcFunc;
3126:   b->func_t          = bcFunc_t;
3127:   b->numids          = numids;
3128:   b->ctx             = ctx;
3129:   b->next            = ds->boundary;
3130:   ds->boundary       = b;
3131:   return(0);
3132: }

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

3137:   Input Parameters:
3138: + ds          - The PetscDS object
3139: . bd          - The boundary condition number
3140: . type        - The type of condition, e.g. DM_BC_ESSENTIAL/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
3141: . name        - The BC name
3142: . labelname   - The label defining constrained points
3143: . field       - The field to constrain
3144: . numcomps    - The number of constrained field components
3145: . comps       - An array of constrained component numbers
3146: . bcFunc      - A pointwise function giving boundary values
3147: . bcFunc_t    - A pointwise function giving the time derviative of the boundary values, or NULL
3148: . numids      - The number of DMLabel ids for constrained points
3149: . ids         - An array of ids for constrained points
3150: - ctx         - An optional user context for bcFunc

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

3155:   Level: developer

3157: .seealso: PetscDSAddBoundary(), PetscDSGetBoundary(), PetscDSGetNumBoundary()
3158: @*/
3159: PetscErrorCode PetscDSUpdateBoundary(PetscDS ds, PetscInt bd, DMBoundaryConditionType type, const char name[], const char labelname[], PetscInt field, PetscInt numcomps, const PetscInt *comps, void (*bcFunc)(void), void (*bcFunc_t)(void), PetscInt numids, const PetscInt *ids, void *ctx)
3160: {
3161:   DSBoundary     b = ds->boundary;
3162:   PetscInt       n = 0;

3167:   while (b) {
3168:     if (n == bd) break;
3169:     b = b->next;
3170:     ++n;
3171:   }
3172:   if (!b) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %d is not in [0, %d)", bd, n);
3173:   if (name) {
3174:     PetscFree(b->name);
3175:     PetscStrallocpy(name, (char **) &b->name);
3176:   }
3177:   if (labelname) {
3178:     PetscFree(b->labelname);
3179:     PetscStrallocpy(labelname, (char **) &b->labelname);
3180:   }
3181:   if (numcomps >= 0 && numcomps != b->numcomps) {
3182:     b->numcomps = numcomps;
3183:     PetscFree(b->comps);
3184:     PetscMalloc1(numcomps, &b->comps);
3185:     if (numcomps) {PetscArraycpy(b->comps, comps, numcomps);}
3186:   }
3187:   if (numids >= 0 && numids != b->numids) {
3188:     b->numids = numids;
3189:     PetscFree(b->ids);
3190:     PetscMalloc1(numids, &b->ids);
3191:     if (numids) {PetscArraycpy(b->ids, ids, numids);}
3192:   }
3193:   b->type = type;
3194:   if (field >= 0) {b->field  = field;}
3195:   if (bcFunc)     {b->func   = bcFunc;}
3196:   if (bcFunc_t)   {b->func_t = bcFunc_t;}
3197:   if (ctx)        {b->ctx    = ctx;}
3198:   return(0);
3199: }

3201: /*@
3202:   PetscDSGetNumBoundary - Get the number of registered BC

3204:   Input Parameters:
3205: . ds - The PetscDS object

3207:   Output Parameters:
3208: . numBd - The number of BC

3210:   Level: intermediate

3212: .seealso: PetscDSAddBoundary(), PetscDSGetBoundary()
3213: @*/
3214: PetscErrorCode PetscDSGetNumBoundary(PetscDS ds, PetscInt *numBd)
3215: {
3216:   DSBoundary b = ds->boundary;

3221:   *numBd = 0;
3222:   while (b) {++(*numBd); b = b->next;}
3223:   return(0);
3224: }

3226: /*@C
3227:   PetscDSGetBoundary - Gets a boundary condition to the model

3229:   Input Parameters:
3230: + ds          - The PetscDS object
3231: - bd          - The BC number

3233:   Output Parameters:
3234: + type        - The type of condition, e.g. DM_BC_ESSENTIAL/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
3235: . name        - The BC name
3236: . labelname   - The label defining constrained points
3237: . field       - The field to constrain
3238: . numcomps    - The number of constrained field components
3239: . comps       - An array of constrained component numbers
3240: . bcFunc      - A pointwise function giving boundary values
3241: . bcFunc_t    - A pointwise function giving the time derviative of the boundary values
3242: . numids      - The number of DMLabel ids for constrained points
3243: . ids         - An array of ids for constrained points
3244: - ctx         - An optional user context for bcFunc

3246:   Options Database Keys:
3247: + -bc_<boundary name> <num> - Overrides the boundary ids
3248: - -bc_<boundary name>_comp <num> - Overrides the boundary components

3250:   Level: developer

3252: .seealso: PetscDSAddBoundary()
3253: @*/
3254: PetscErrorCode PetscDSGetBoundary(PetscDS ds, PetscInt bd, DMBoundaryConditionType *type, const char **name, const char **labelname, PetscInt *field, PetscInt *numcomps, const PetscInt **comps, void (**func)(void), void (**func_t)(void), PetscInt *numids, const PetscInt **ids, void **ctx)
3255: {
3256:   DSBoundary b    = ds->boundary;
3257:   PetscInt   n    = 0;

3261:   while (b) {
3262:     if (n == bd) break;
3263:     b = b->next;
3264:     ++n;
3265:   }
3266:   if (!b) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %d is not in [0, %d)", bd, n);
3267:   if (type) {
3269:     *type = b->type;
3270:   }
3271:   if (name) {
3273:     *name = b->name;
3274:   }
3275:   if (labelname) {
3277:     *labelname = b->labelname;
3278:   }
3279:   if (field) {
3281:     *field = b->field;
3282:   }
3283:   if (numcomps) {
3285:     *numcomps = b->numcomps;
3286:   }
3287:   if (comps) {
3289:     *comps = b->comps;
3290:   }
3291:   if (func) {
3293:     *func = b->func;
3294:   }
3295:   if (func_t) {
3297:     *func_t = b->func_t;
3298:   }
3299:   if (numids) {
3301:     *numids = b->numids;
3302:   }
3303:   if (ids) {
3305:     *ids = b->ids;
3306:   }
3307:   if (ctx) {
3309:     *ctx = b->ctx;
3310:   }
3311:   return(0);
3312: }

3314: /*@
3315:   PetscDSCopyBoundary - Copy all boundary condition objects to the new problem

3317:   Not collective

3319:   Input Parameter:
3320: . prob - The PetscDS object

3322:   Output Parameter:
3323: . newprob - The PetscDS copy

3325:   Level: intermediate

3327: .seealso: PetscDSCopyEquations(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3328: @*/
3329: PetscErrorCode PetscDSCopyBoundary(PetscDS probA, PetscDS probB)
3330: {
3331:   DSBoundary     b, next, *lastnext;

3337:   if (probA == probB) return(0);
3338:   next = probB->boundary;
3339:   while (next) {
3340:     DSBoundary b = next;

3342:     next = b->next;
3343:     PetscFree(b->comps);
3344:     PetscFree(b->ids);
3345:     PetscFree(b->name);
3346:     PetscFree(b->labelname);
3347:     PetscFree(b);
3348:   }
3349:   lastnext = &(probB->boundary);
3350:   for (b = probA->boundary; b; b = b->next) {
3351:     DSBoundary bNew;

3353:     PetscNew(&bNew);
3354:     bNew->numcomps = b->numcomps;
3355:     PetscMalloc1(bNew->numcomps, &bNew->comps);
3356:     PetscArraycpy(bNew->comps, b->comps, bNew->numcomps);
3357:     bNew->numids = b->numids;
3358:     PetscMalloc1(bNew->numids, &bNew->ids);
3359:     PetscArraycpy(bNew->ids, b->ids, bNew->numids);
3360:     PetscStrallocpy(b->labelname,(char **) &(bNew->labelname));
3361:     PetscStrallocpy(b->name,(char **) &(bNew->name));
3362:     bNew->ctx   = b->ctx;
3363:     bNew->type  = b->type;
3364:     bNew->field = b->field;
3365:     bNew->func  = b->func;

3367:     *lastnext = bNew;
3368:     lastnext = &(bNew->next);
3369:   }
3370:   return(0);
3371: }

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

3376:   Not collective

3378:   Input Parameter:
3379: + prob - The PetscDS object
3380: . numFields - Number of new fields
3381: - fields - Old field number for each new field

3383:   Output Parameter:
3384: . newprob - The PetscDS copy

3386:   Level: intermediate

3388: .seealso: PetscDSSelectEquations(), PetscDSCopyBoundary(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3389: @*/
3390: PetscErrorCode PetscDSSelectDiscretizations(PetscDS prob, PetscInt numFields, const PetscInt fields[], PetscDS newprob)
3391: {
3392:   PetscInt       Nf, Nfn, fn;

3399:   PetscDSGetNumFields(prob, &Nf);
3400:   PetscDSGetNumFields(newprob, &Nfn);
3401:   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);
3402:   for (fn = 0; fn < numFields; ++fn) {
3403:     const PetscInt f = fields ? fields[fn] : fn;
3404:     PetscObject    disc;

3406:     if (f >= Nf) continue;
3407:     PetscDSGetDiscretization(prob, f, &disc);
3408:     PetscDSSetDiscretization(newprob, fn, disc);
3409:   }
3410:   return(0);
3411: }

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

3416:   Not collective

3418:   Input Parameter:
3419: + prob - The PetscDS object
3420: . numFields - Number of new fields
3421: - fields - Old field number for each new field

3423:   Output Parameter:
3424: . newprob - The PetscDS copy

3426:   Level: intermediate

3428: .seealso: PetscDSSelectDiscretizations(), PetscDSCopyBoundary(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3429: @*/
3430: PetscErrorCode PetscDSSelectEquations(PetscDS prob, PetscInt numFields, const PetscInt fields[], PetscDS newprob)
3431: {
3432:   PetscInt       Nf, Nfn, fn, gn;

3439:   PetscDSGetNumFields(prob, &Nf);
3440:   PetscDSGetNumFields(newprob, &Nfn);
3441:   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);
3442:   for (fn = 0; fn < numFields; ++fn) {
3443:     const PetscInt   f = fields ? fields[fn] : fn;
3444:     PetscPointFunc   obj;
3445:     PetscPointFunc   f0, f1;
3446:     PetscBdPointFunc f0Bd, f1Bd;
3447:     PetscRiemannFunc r;

3449:     if (f >= Nf) continue;
3450:     PetscDSGetObjective(prob, f, &obj);
3451:     PetscDSGetResidual(prob, f, &f0, &f1);
3452:     PetscDSGetBdResidual(prob, f, &f0Bd, &f1Bd);
3453:     PetscDSGetRiemannSolver(prob, f, &r);
3454:     PetscDSSetObjective(newprob, fn, obj);
3455:     PetscDSSetResidual(newprob, fn, f0, f1);
3456:     PetscDSSetBdResidual(newprob, fn, f0Bd, f1Bd);
3457:     PetscDSSetRiemannSolver(newprob, fn, r);
3458:     for (gn = 0; gn < numFields; ++gn) {
3459:       const PetscInt  g = fields ? fields[gn] : gn;
3460:       PetscPointJac   g0, g1, g2, g3;
3461:       PetscPointJac   g0p, g1p, g2p, g3p;
3462:       PetscBdPointJac g0Bd, g1Bd, g2Bd, g3Bd;

3464:       if (g >= Nf) continue;
3465:       PetscDSGetJacobian(prob, f, g, &g0, &g1, &g2, &g3);
3466:       PetscDSGetJacobianPreconditioner(prob, f, g, &g0p, &g1p, &g2p, &g3p);
3467:       PetscDSGetBdJacobian(prob, f, g, &g0Bd, &g1Bd, &g2Bd, &g3Bd);
3468:       PetscDSSetJacobian(newprob, fn, gn, g0, g1, g2, g3);
3469:       PetscDSSetJacobianPreconditioner(newprob, fn, gn, g0p, g1p, g2p, g3p);
3470:       PetscDSSetBdJacobian(newprob, fn, gn, g0Bd, g1Bd, g2Bd, g3Bd);
3471:     }
3472:   }
3473:   return(0);
3474: }

3476: /*@
3477:   PetscDSCopyEquations - Copy all pointwise function pointers to the new problem

3479:   Not collective

3481:   Input Parameter:
3482: . prob - The PetscDS object

3484:   Output Parameter:
3485: . newprob - The PetscDS copy

3487:   Level: intermediate

3489: .seealso: PetscDSCopyBoundary(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3490: @*/
3491: PetscErrorCode PetscDSCopyEquations(PetscDS prob, PetscDS newprob)
3492: {
3493:   PetscInt       Nf, Ng;

3499:   PetscDSGetNumFields(prob, &Nf);
3500:   PetscDSGetNumFields(newprob, &Ng);
3501:   if (Nf != Ng) SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_SIZ, "Number of fields must match %D != %D", Nf, Ng);
3502:   PetscDSSelectEquations(prob, Nf, NULL, newprob);
3503:   return(0);
3504: }
3505: /*@
3506:   PetscDSCopyConstants - Copy all constants to the new problem

3508:   Not collective

3510:   Input Parameter:
3511: . prob - The PetscDS object

3513:   Output Parameter:
3514: . newprob - The PetscDS copy

3516:   Level: intermediate

3518: .seealso: PetscDSCopyBoundary(), PetscDSCopyEquations(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3519: @*/
3520: PetscErrorCode PetscDSCopyConstants(PetscDS prob, PetscDS newprob)
3521: {
3522:   PetscInt           Nc;
3523:   const PetscScalar *constants;
3524:   PetscErrorCode     ierr;

3529:   PetscDSGetConstants(prob, &Nc, &constants);
3530:   PetscDSSetConstants(newprob, Nc, (PetscScalar *) constants);
3531:   return(0);
3532: }

3534: PetscErrorCode PetscDSGetHeightSubspace(PetscDS prob, PetscInt height, PetscDS *subprob)
3535: {
3536:   PetscInt       dim, Nf, f;

3542:   if (height == 0) {*subprob = prob; return(0);}
3543:   PetscDSGetNumFields(prob, &Nf);
3544:   PetscDSGetSpatialDimension(prob, &dim);
3545:   if (height > dim) SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_OUTOFRANGE, "DS can only handle height in [0, %D], not %D", dim, height);
3546:   if (!prob->subprobs) {PetscCalloc1(dim, &prob->subprobs);}
3547:   if (!prob->subprobs[height-1]) {
3548:     PetscInt cdim;

3550:     PetscDSCreate(PetscObjectComm((PetscObject) prob), &prob->subprobs[height-1]);
3551:     PetscDSGetCoordinateDimension(prob, &cdim);
3552:     PetscDSSetCoordinateDimension(prob->subprobs[height-1], cdim);
3553:     for (f = 0; f < Nf; ++f) {
3554:       PetscFE      subfe;
3555:       PetscObject  obj;
3556:       PetscClassId id;

3558:       PetscDSGetDiscretization(prob, f, &obj);
3559:       PetscObjectGetClassId(obj, &id);
3560:       if (id == PETSCFE_CLASSID) {PetscFEGetHeightSubspace((PetscFE) obj, height, &subfe);}
3561:       else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unsupported discretization type for field %d", f);
3562:       PetscDSSetDiscretization(prob->subprobs[height-1], f, (PetscObject) subfe);
3563:     }
3564:   }
3565:   *subprob = prob->subprobs[height-1];
3566:   return(0);
3567: }

3569: PetscErrorCode PetscDSGetDiscType_Internal(PetscDS ds, PetscInt f, PetscDiscType *disctype)
3570: {
3571:   PetscObject    obj;
3572:   PetscClassId   id;
3573:   PetscInt       Nf;

3579:   *disctype = PETSC_DISC_NONE;
3580:   PetscDSGetNumFields(ds, &Nf);
3581:   if (f >= Nf) SETERRQ2(PetscObjectComm((PetscObject) ds), PETSC_ERR_ARG_SIZ, "Field %D must be in [0, %D)", f, Nf);
3582:   PetscDSGetDiscretization(ds, f, &obj);
3583:   if (obj) {
3584:     PetscObjectGetClassId(obj, &id);
3585:     if (id == PETSCFE_CLASSID) *disctype = PETSC_DISC_FE;
3586:     else                       *disctype = PETSC_DISC_FV;
3587:   }
3588:   return(0);
3589: }

3591: static PetscErrorCode PetscDSDestroy_Basic(PetscDS prob)
3592: {
3593:   PetscErrorCode      ierr;

3596:   PetscFree(prob->data);
3597:   return(0);
3598: }

3600: static PetscErrorCode PetscDSInitialize_Basic(PetscDS prob)
3601: {
3603:   prob->ops->setfromoptions = NULL;
3604:   prob->ops->setup          = NULL;
3605:   prob->ops->view           = NULL;
3606:   prob->ops->destroy        = PetscDSDestroy_Basic;
3607:   return(0);
3608: }

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

3613:   Level: intermediate

3615: .seealso: PetscDSType, PetscDSCreate(), PetscDSSetType()
3616: M*/

3618: PETSC_EXTERN PetscErrorCode PetscDSCreate_Basic(PetscDS prob)
3619: {
3620:   PetscDS_Basic *b;
3621:   PetscErrorCode      ierr;

3625:   PetscNewLog(prob, &b);
3626:   prob->data = b;

3628:   PetscDSInitialize_Basic(prob);
3629:   return(0);
3630: }