Actual source code: dtds.c

petsc-master 2020-05-26
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;
462:   void                **tmpexactCtx;
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:   PetscCalloc5(NfNew*2, &tmpfbd, NfNew*NfNew*4, &tmpgbd, NfNew*NfNew*4, &tmpgpbd, NfNew, &tmpexactSol, NfNew, &tmpexactCtx);
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 = Nf*2; f < NfNew*2; ++f) tmpfbd[f] = NULL;
512:   for (f = Nf*Nf*4; f < NfNew*NfNew*4; ++f) tmpgbd[f] = NULL;
513:   for (f = Nf*Nf*4; f < NfNew*NfNew*4; ++f) tmpgpbd[f] = NULL;
514:   for (f = Nf; f < NfNew; ++f) tmpexactSol[f] = NULL;
515:   for (f = Nf; f < NfNew; ++f) tmpexactCtx[f] = NULL;
516:   PetscFree5(prob->fBd, prob->gBd, prob->gpBd, prob->exactSol, prob->exactCtx);
517:   prob->fBd = tmpfbd;
518:   prob->gBd = tmpgbd;
519:   prob->gpBd = tmpgpbd;
520:   prob->exactSol = tmpexactSol;
521:   prob->exactCtx = tmpexactCtx;
522:   return(0);
523: }

525: /*@
526:   PetscDSDestroy - Destroys a PetscDS object

528:   Collective on prob

530:   Input Parameter:
531: . prob - the PetscDS object to destroy

533:   Level: developer

535: .seealso PetscDSView()
536: @*/
537: PetscErrorCode PetscDSDestroy(PetscDS *prob)
538: {
539:   PetscInt       f;
540:   DSBoundary     next;

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

547:   if (--((PetscObject)(*prob))->refct > 0) {*prob = 0; return(0);}
548:   ((PetscObject) (*prob))->refct = 0;
549:   if ((*prob)->subprobs) {
550:     PetscInt dim, d;

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

569:     next = b->next;
570:     PetscFree(b->comps);
571:     PetscFree(b->ids);
572:     PetscFree(b->name);
573:     PetscFree(b->labelname);
574:     PetscFree(b);
575:   }
576:   PetscFree((*prob)->constants);
577:   PetscHeaderDestroy(prob);
578:   return(0);
579: }

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

584:   Collective

586:   Input Parameter:
587: . comm - The communicator for the PetscDS object

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

592:   Level: beginner

594: .seealso: PetscDSSetType(), PETSCDSBASIC
595: @*/
596: PetscErrorCode PetscDSCreate(MPI_Comm comm, PetscDS *prob)
597: {
598:   PetscDS   p;

603:   *prob  = NULL;
604:   PetscDSInitializePackage();

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

608:   p->Nf           = 0;
609:   p->setup        = PETSC_FALSE;
610:   p->numConstants = 0;
611:   p->constants    = NULL;
612:   p->dimEmbed     = -1;
613:   p->useJacPre    = PETSC_TRUE;

615:   *prob = p;
616:   return(0);
617: }

619: /*@
620:   PetscDSGetNumFields - Returns the number of fields in the DS

622:   Not collective

624:   Input Parameter:
625: . prob - The PetscDS object

627:   Output Parameter:
628: . Nf - The number of fields

630:   Level: beginner

632: .seealso: PetscDSGetSpatialDimension(), PetscDSCreate()
633: @*/
634: PetscErrorCode PetscDSGetNumFields(PetscDS prob, PetscInt *Nf)
635: {
639:   *Nf = prob->Nf;
640:   return(0);
641: }

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

646:   Not collective

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

651:   Output Parameter:
652: . dim - The spatial dimension

654:   Level: beginner

656: .seealso: PetscDSGetCoordinateDimension(), PetscDSGetNumFields(), PetscDSCreate()
657: @*/
658: PetscErrorCode PetscDSGetSpatialDimension(PetscDS prob, PetscInt *dim)
659: {

665:   *dim = 0;
666:   if (prob->Nf) {
667:     PetscObject  obj;
668:     PetscClassId id;

670:     PetscDSGetDiscretization(prob, 0, &obj);
671:     if (obj) {
672:       PetscObjectGetClassId(obj, &id);
673:       if (id == PETSCFE_CLASSID)      {PetscFEGetSpatialDimension((PetscFE) obj, dim);}
674:       else if (id == PETSCFV_CLASSID) {PetscFVGetSpatialDimension((PetscFV) obj, dim);}
675:       else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", 0);
676:     }
677:   }
678:   return(0);
679: }

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

684:   Not collective

686:   Input Parameter:
687: . prob - The PetscDS object

689:   Output Parameter:
690: . dimEmbed - The coordinate dimension

692:   Level: beginner

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

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

709:   Logically collective on prob

711:   Input Parameters:
712: + prob - The PetscDS object
713: - dimEmbed - The coordinate dimension

715:   Level: beginner

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

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

731:   Not collective

733:   Input Parameter:
734: . prob - The PetscDS object

736:   Output Parameter:
737: . isHybrid - The flag

739:   Level: developer

741: .seealso: PetscDSSetHybrid(), PetscDSCreate()
742: @*/
743: PetscErrorCode PetscDSGetHybrid(PetscDS prob, PetscBool *isHybrid)
744: {
748:   *isHybrid = prob->isHybrid;
749:   return(0);
750: }

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

755:   Not collective

757:   Input Parameters:
758: + prob - The PetscDS object
759: - isHybrid - The flag

761:   Level: developer

763: .seealso: PetscDSGetHybrid(), PetscDSCreate()
764: @*/
765: PetscErrorCode PetscDSSetHybrid(PetscDS prob, PetscBool isHybrid)
766: {
769:   prob->isHybrid = isHybrid;
770:   return(0);
771: }

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

776:   Not collective

778:   Input Parameter:
779: . prob - The PetscDS object

781:   Output Parameter:
782: . dim - The total problem dimension

784:   Level: beginner

786: .seealso: PetscDSGetNumFields(), PetscDSCreate()
787: @*/
788: PetscErrorCode PetscDSGetTotalDimension(PetscDS prob, PetscInt *dim)
789: {

794:   PetscDSSetUp(prob);
796:   *dim = prob->totDim;
797:   return(0);
798: }

800: /*@
801:   PetscDSGetTotalComponents - Returns the total number of components in this system

803:   Not collective

805:   Input Parameter:
806: . prob - The PetscDS object

808:   Output Parameter:
809: . dim - The total number of components

811:   Level: beginner

813: .seealso: PetscDSGetNumFields(), PetscDSCreate()
814: @*/
815: PetscErrorCode PetscDSGetTotalComponents(PetscDS prob, PetscInt *Nc)
816: {

821:   PetscDSSetUp(prob);
823:   *Nc = prob->totComp;
824:   return(0);
825: }

827: /*@
828:   PetscDSGetDiscretization - Returns the discretization object for the given field

830:   Not collective

832:   Input Parameters:
833: + prob - The PetscDS object
834: - f - The field number

836:   Output Parameter:
837: . disc - The discretization object

839:   Level: beginner

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

853: /*@
854:   PetscDSSetDiscretization - Sets the discretization object for the given field

856:   Not collective

858:   Input Parameters:
859: + prob - The PetscDS object
860: . f - The field number
861: - disc - The discretization object

863:   Level: beginner

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

874:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
875:   PetscDSEnlarge_Static(prob, f+1);
876:   PetscObjectDereference(prob->disc[f]);
877:   prob->disc[f] = disc;
878:   PetscObjectReference(disc);
879:   if (disc) {
880:     PetscClassId id;

882:     PetscObjectGetClassId(disc, &id);
883:     if (id == PETSCFE_CLASSID) {
884:       PetscDSSetImplicit(prob, f, PETSC_TRUE);
885:     } else if (id == PETSCFV_CLASSID) {
886:       PetscDSSetImplicit(prob, f, PETSC_FALSE);
887:     }
888:   }
889:   return(0);
890: }

892: /*@
893:   PetscDSAddDiscretization - Adds a discretization object

895:   Not collective

897:   Input Parameters:
898: + prob - The PetscDS object
899: - disc - The boundary discretization object

901:   Level: beginner

903: .seealso: PetscDSGetDiscretization(), PetscDSSetDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
904: @*/
905: PetscErrorCode PetscDSAddDiscretization(PetscDS prob, PetscObject disc)
906: {

910:   PetscDSSetDiscretization(prob, prob->Nf, disc);
911:   return(0);
912: }

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

917:   Not collective

919:   Input Parameter:
920: . prob - The PetscDS object

922:   Output Parameter:
923: . q - The quadrature object

925: Level: intermediate

927: .seealso: PetscDSSetImplicit(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
928: @*/
929: PetscErrorCode PetscDSGetQuadrature(PetscDS prob, PetscQuadrature *q)
930: {
931:   PetscObject    obj;
932:   PetscClassId   id;

936:   *q = NULL;
937:   if (!prob->Nf) return(0);
938:   PetscDSGetDiscretization(prob, 0, &obj);
939:   PetscObjectGetClassId(obj, &id);
940:   if      (id == PETSCFE_CLASSID) {PetscFEGetQuadrature((PetscFE) obj, q);}
941:   else if (id == PETSCFV_CLASSID) {PetscFVGetQuadrature((PetscFV) obj, q);}
942:   else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", 0);
943:   return(0);
944: }

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

949:   Not collective

951:   Input Parameters:
952: + prob - The PetscDS object
953: - f - The field number

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

958:   Level: developer

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

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

975:   Not collective

977:   Input Parameters:
978: + prob - The PetscDS object
979: . f - The field number
980: - implicit - The flag indicating what kind of solve to use for this field

982:   Level: developer

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

995: PetscErrorCode PetscDSGetObjective(PetscDS prob, PetscInt f,
996:                                    void (**obj)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
997:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
998:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
999:                                                 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[]))
1000: {
1004:   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);
1005:   *obj = prob->obj[f];
1006:   return(0);
1007: }

1009: PetscErrorCode PetscDSSetObjective(PetscDS prob, PetscInt f,
1010:                                    void (*obj)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1011:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1012:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1013:                                                PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[]))
1014: {

1020:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1021:   PetscDSEnlarge_Static(prob, f+1);
1022:   prob->obj[f] = obj;
1023:   return(0);
1024: }

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

1029:   Not collective

1031:   Input Parameters:
1032: + prob - The PetscDS
1033: - f    - The test field number

1035:   Output Parameters:
1036: + f0 - integrand for the test function term
1037: - f1 - integrand for the test function gradient term

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

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

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

1045: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1046: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1047: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1048: $    PetscReal t, const PetscReal x[], PetscScalar f0[])

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

1068:   Level: intermediate

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

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

1093:   Not collective

1095:   Input Parameters:
1096: + prob - The PetscDS
1097: . f    - The test field number
1098: . f0 - integrand for the test function term
1099: - f1 - integrand for the test function gradient term

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

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

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

1107: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1108: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1109: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1110: $    PetscReal t, const PetscReal x[], PetscScalar f0[])

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

1130:   Level: intermediate

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

1150:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1151:   PetscDSEnlarge_Static(prob, f+1);
1152:   prob->f[f*2+0] = f0;
1153:   prob->f[f*2+1] = f1;
1154:   return(0);
1155: }

1157: /*@C
1158:   PetscDSHasJacobian - Signals that Jacobian functions have been set

1160:   Not collective

1162:   Input Parameter:
1163: . prob - The PetscDS

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

1168:   Level: intermediate

1170: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1171: @*/
1172: PetscErrorCode PetscDSHasJacobian(PetscDS prob, PetscBool *hasJac)
1173: {
1174:   PetscInt f, g, h;

1178:   *hasJac = PETSC_FALSE;
1179:   for (f = 0; f < prob->Nf; ++f) {
1180:     for (g = 0; g < prob->Nf; ++g) {
1181:       for (h = 0; h < 4; ++h) {
1182:         if (prob->g[(f*prob->Nf + g)*4+h]) *hasJac = PETSC_TRUE;
1183:       }
1184:     }
1185:   }
1186:   return(0);
1187: }

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

1192:   Not collective

1194:   Input Parameters:
1195: + prob - The PetscDS
1196: . f    - The test field number
1197: - g    - The field number

1199:   Output Parameters:
1200: + g0 - integrand for the test and basis function term
1201: . g1 - integrand for the test function and basis function gradient term
1202: . g2 - integrand for the test function gradient and basis function term
1203: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

1211: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1212: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1213: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1214: $    PetscReal t, const PetscReal u_tShift, const PetscReal x[], PetscScalar g0[])

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

1235:   Level: intermediate

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

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

1271:   Not collective

1273:   Input Parameters:
1274: + prob - The PetscDS
1275: . f    - The test field number
1276: . g    - The field number
1277: . g0 - integrand for the test and basis function term
1278: . g1 - integrand for the test function and basis function gradient term
1279: . g2 - integrand for the test function gradient and basis function term
1280: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

1288: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1289: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1290: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1291: $    PetscReal t, const PetscReal x[], PetscScalar g0[])

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

1312:   Level: intermediate

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

1342:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1343:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
1344:   PetscDSEnlarge_Static(prob, PetscMax(f, g)+1);
1345:   prob->g[(f*prob->Nf + g)*4+0] = g0;
1346:   prob->g[(f*prob->Nf + g)*4+1] = g1;
1347:   prob->g[(f*prob->Nf + g)*4+2] = g2;
1348:   prob->g[(f*prob->Nf + g)*4+3] = g3;
1349:   return(0);
1350: }

1352: /*@C
1353:   PetscDSUseJacobianPreconditioner - Whether to construct a Jacobian preconditioner

1355:   Not collective

1357:   Input Parameters:
1358: + prob - The PetscDS
1359: - useJacPre - flag that enables construction of a Jacobian preconditioner

1361:   Level: intermediate

1363: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1364: @*/
1365: PetscErrorCode PetscDSUseJacobianPreconditioner(PetscDS prob, PetscBool useJacPre)
1366: {
1369:   prob->useJacPre = useJacPre;
1370:   return(0);
1371: }

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

1376:   Not collective

1378:   Input Parameter:
1379: . prob - The PetscDS

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

1384:   Level: intermediate

1386: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1387: @*/
1388: PetscErrorCode PetscDSHasJacobianPreconditioner(PetscDS prob, PetscBool *hasJacPre)
1389: {
1390:   PetscInt f, g, h;

1394:   *hasJacPre = PETSC_FALSE;
1395:   if (!prob->useJacPre) return(0);
1396:   for (f = 0; f < prob->Nf; ++f) {
1397:     for (g = 0; g < prob->Nf; ++g) {
1398:       for (h = 0; h < 4; ++h) {
1399:         if (prob->gp[(f*prob->Nf + g)*4+h]) *hasJacPre = PETSC_TRUE;
1400:       }
1401:     }
1402:   }
1403:   return(0);
1404: }

1406: /*@C
1407:   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.

1409:   Not collective

1411:   Input Parameters:
1412: + prob - The PetscDS
1413: . f    - The test field number
1414: - g    - The field number

1416:   Output Parameters:
1417: + g0 - integrand for the test and basis function term
1418: . g1 - integrand for the test function and basis function gradient term
1419: . g2 - integrand for the test function gradient and basis function term
1420: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

1428: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1429: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1430: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1431: $    PetscReal t, const PetscReal u_tShift, const PetscReal x[], PetscScalar g0[])

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

1452:   Level: intermediate

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

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

1488:   Not collective

1490:   Input Parameters:
1491: + prob - The PetscDS
1492: . f    - The test field number
1493: . g    - The field number
1494: . g0 - integrand for the test and basis function term
1495: . g1 - integrand for the test function and basis function gradient term
1496: . g2 - integrand for the test function gradient and basis function term
1497: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

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

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

1529:   Level: intermediate

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

1559:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1560:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
1561:   PetscDSEnlarge_Static(prob, PetscMax(f, g)+1);
1562:   prob->gp[(f*prob->Nf + g)*4+0] = g0;
1563:   prob->gp[(f*prob->Nf + g)*4+1] = g1;
1564:   prob->gp[(f*prob->Nf + g)*4+2] = g2;
1565:   prob->gp[(f*prob->Nf + g)*4+3] = g3;
1566:   return(0);
1567: }

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

1572:   Not collective

1574:   Input Parameter:
1575: . prob - The PetscDS

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

1580:   Level: intermediate

1582: .seealso: PetscDSGetDynamicJacobian(), PetscDSSetDynamicJacobian(), PetscDSGetJacobian()
1583: @*/
1584: PetscErrorCode PetscDSHasDynamicJacobian(PetscDS prob, PetscBool *hasDynJac)
1585: {
1586:   PetscInt f, g, h;

1590:   *hasDynJac = PETSC_FALSE;
1591:   for (f = 0; f < prob->Nf; ++f) {
1592:     for (g = 0; g < prob->Nf; ++g) {
1593:       for (h = 0; h < 4; ++h) {
1594:         if (prob->gt[(f*prob->Nf + g)*4+h]) *hasDynJac = PETSC_TRUE;
1595:       }
1596:     }
1597:   }
1598:   return(0);
1599: }

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

1604:   Not collective

1606:   Input Parameters:
1607: + prob - The PetscDS
1608: . f    - The test field number
1609: - g    - The field number

1611:   Output Parameters:
1612: + g0 - integrand for the test and basis function term
1613: . g1 - integrand for the test function and basis function gradient term
1614: . g2 - integrand for the test function gradient and basis function term
1615: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

1623: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1624: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1625: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1626: $    PetscReal t, const PetscReal u_tShift, const PetscReal x[], PetscScalar g0[])

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

1647:   Level: intermediate

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

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

1683:   Not collective

1685:   Input Parameters:
1686: + prob - The PetscDS
1687: . f    - The test field number
1688: . g    - The field number
1689: . g0 - integrand for the test and basis function term
1690: . g1 - integrand for the test function and basis function gradient term
1691: . g2 - integrand for the test function gradient and basis function term
1692: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

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

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

1724:   Level: intermediate

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

1754:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1755:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
1756:   PetscDSEnlarge_Static(prob, PetscMax(f, g)+1);
1757:   prob->gt[(f*prob->Nf + g)*4+0] = g0;
1758:   prob->gt[(f*prob->Nf + g)*4+1] = g1;
1759:   prob->gt[(f*prob->Nf + g)*4+2] = g2;
1760:   prob->gt[(f*prob->Nf + g)*4+3] = g3;
1761:   return(0);
1762: }

1764: /*@C
1765:   PetscDSGetRiemannSolver - Returns the Riemann solver for the given field

1767:   Not collective

1769:   Input Arguments:
1770: + prob - The PetscDS object
1771: - f    - The field number

1773:   Output Argument:
1774: . r    - Riemann solver

1776:   Calling sequence for r:

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

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

1791:   Level: intermediate

1793: .seealso: PetscDSSetRiemannSolver()
1794: @*/
1795: PetscErrorCode PetscDSGetRiemannSolver(PetscDS prob, PetscInt f,
1796:                                        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))
1797: {
1800:   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);
1802:   *r = prob->r[f];
1803:   return(0);
1804: }

1806: /*@C
1807:   PetscDSSetRiemannSolver - Sets the Riemann solver for the given field

1809:   Not collective

1811:   Input Arguments:
1812: + prob - The PetscDS object
1813: . f    - The field number
1814: - r    - Riemann solver

1816:   Calling sequence for r:

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

1820: + dim  - The spatial dimension
1821: . Nf   - The number of fields
1822: . x    - The coordinates at a point on the interface
1823: . n    - The normal vector to the interface
1824: . uL   - The state vector to the left of the interface
1825: . uR   - The state vector to the right of the interface
1826: . flux - output array of flux through the interface
1827: . numConstants - number of constant parameters
1828: . constants - constant parameters
1829: - ctx  - optional user context

1831:   Level: intermediate

1833: .seealso: PetscDSGetRiemannSolver()
1834: @*/
1835: PetscErrorCode PetscDSSetRiemannSolver(PetscDS prob, PetscInt f,
1836:                                        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))
1837: {

1843:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1844:   PetscDSEnlarge_Static(prob, f+1);
1845:   prob->r[f] = r;
1846:   return(0);
1847: }

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

1852:   Not collective

1854:   Input Parameters:
1855: + prob - The PetscDS
1856: - f    - The field number

1858:   Output Parameters:
1859: . update - update function

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

1863: $ update(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1864: $        const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1865: $        const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1866: $        PetscReal t, const PetscReal x[], PetscScalar uNew[])

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

1884:   Level: intermediate

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

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

1904:   Not collective

1906:   Input Parameters:
1907: + prob   - The PetscDS
1908: . f      - The field number
1909: - update - update function

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

1913: $ update(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1914: $        const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1915: $        const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1916: $        PetscReal t, const PetscReal x[], PetscScalar uNew[])

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

1934:   Level: intermediate

1936: .seealso: PetscDSGetResidual()
1937: @*/
1938: PetscErrorCode PetscDSSetUpdate(PetscDS prob, PetscInt f,
1939:                                 void (*update)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1940:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1941:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1942:                                                PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar uNew[]))
1943: {

1949:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1950:   PetscDSEnlarge_Static(prob, f+1);
1951:   prob->update[f] = update;
1952:   return(0);
1953: }

1955: PetscErrorCode PetscDSGetContext(PetscDS prob, PetscInt f, void **ctx)
1956: {
1959:   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);
1961:   *ctx = prob->ctx[f];
1962:   return(0);
1963: }

1965: PetscErrorCode PetscDSSetContext(PetscDS prob, PetscInt f, void *ctx)
1966: {

1971:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1972:   PetscDSEnlarge_Static(prob, f+1);
1973:   prob->ctx[f] = ctx;
1974:   return(0);
1975: }

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

1980:   Not collective

1982:   Input Parameters:
1983: + prob - The PetscDS
1984: - f    - The test field number

1986:   Output Parameters:
1987: + f0 - boundary integrand for the test function term
1988: - f1 - boundary integrand for the test function gradient term

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

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

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

1996: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1997: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1998: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1999: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[])

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

2020:   Level: intermediate

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

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

2045:   Not collective

2047:   Input Parameters:
2048: + prob - The PetscDS
2049: . f    - The test field number
2050: . f0 - boundary integrand for the test function term
2051: - f1 - boundary integrand for the test function gradient term

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

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

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

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

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

2083:   Level: intermediate

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

2101:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2102:   PetscDSEnlarge_Static(prob, f+1);
2105:   return(0);
2106: }

2108: /*@
2109:   PetscDSHasBdJacobian - Signals that boundary Jacobian functions have been set

2111:   Not collective

2113:   Input Parameter:
2114: . prob - The PetscDS

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

2119:   Level: intermediate

2121: .seealso: PetscDSHasJacobian(), PetscDSSetBdJacobian(), PetscDSGetBdJacobian()
2122: @*/
2123: PetscErrorCode PetscDSHasBdJacobian(PetscDS prob, PetscBool *hasBdJac)
2124: {
2125:   PetscInt f, g, h;

2129:   *hasBdJac = PETSC_FALSE;
2130:   for (f = 0; f < prob->Nf; ++f) {
2131:     for (g = 0; g < prob->Nf; ++g) {
2132:       for (h = 0; h < 4; ++h) {
2133:         if (prob->gBd[(f*prob->Nf + g)*4+h]) *hasBdJac = PETSC_TRUE;
2134:       }
2135:     }
2136:   }
2137:   return(0);
2138: }

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

2143:   Not collective

2145:   Input Parameters:
2146: + prob - The PetscDS
2147: . f    - The test field number
2148: - g    - The field number

2150:   Output Parameters:
2151: + g0 - integrand for the test and basis function term
2152: . g1 - integrand for the test function and basis function gradient term
2153: . g2 - integrand for the test function gradient and basis function term
2154: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

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

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

2187:   Level: intermediate

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

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

2223:   Not collective

2225:   Input Parameters:
2226: + prob - The PetscDS
2227: . f    - The test field number
2228: . g    - The field number
2229: . g0 - integrand for the test and basis function term
2230: . g1 - integrand for the test function and basis function gradient term
2231: . g2 - integrand for the test function gradient and basis function term
2232: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

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

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

2265:   Level: intermediate

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

2295:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2296:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
2297:   PetscDSEnlarge_Static(prob, PetscMax(f, g)+1);
2298:   prob->gBd[(f*prob->Nf + g)*4+0] = g0;
2299:   prob->gBd[(f*prob->Nf + g)*4+1] = g1;
2300:   prob->gBd[(f*prob->Nf + g)*4+2] = g2;
2301:   prob->gBd[(f*prob->Nf + g)*4+3] = g3;
2302:   return(0);
2303: }

2305: /*@
2306:   PetscDSHasBdJacobianPreconditioner - Signals that boundary Jacobian preconditioner functions have been set

2308:   Not collective

2310:   Input Parameter:
2311: . prob - The PetscDS

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

2316:   Level: intermediate

2318: .seealso: PetscDSHasJacobian(), PetscDSSetBdJacobian(), PetscDSGetBdJacobian()
2319: @*/
2320: PetscErrorCode PetscDSHasBdJacobianPreconditioner(PetscDS prob, PetscBool *hasBdJacPre)
2321: {
2322:   PetscInt f, g, h;

2326:   *hasBdJacPre = PETSC_FALSE;
2327:   for (f = 0; f < prob->Nf; ++f) {
2328:     for (g = 0; g < prob->Nf; ++g) {
2329:       for (h = 0; h < 4; ++h) {
2330:         if (prob->gpBd[(f*prob->Nf + g)*4+h]) *hasBdJacPre = PETSC_TRUE;
2331:       }
2332:     }
2333:   }
2334:   return(0);
2335: }

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

2340:   Not collective

2342:   Input Parameters:
2343: + prob - The PetscDS
2344: . f    - The test field number
2345: - g    - The field number

2347:   Output Parameters:
2348: + g0 - integrand for the test and basis function term
2349: . g1 - integrand for the test function and basis function gradient term
2350: . g2 - integrand for the test function gradient and basis function term
2351: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

2359: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2360: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2361: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2362: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])

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

2385:   This is not yet available in Fortran.

2387:   Level: intermediate

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

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

2423:   Not collective

2425:   Input Parameters:
2426: + prob - The PetscDS
2427: . f    - The test field number
2428: . g    - The field number
2429: . g0 - integrand for the test and basis function term
2430: . g1 - integrand for the test function and basis function gradient term
2431: . g2 - integrand for the test function gradient and basis function term
2432: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

2440: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2441: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2442: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2443: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])

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

2466:   This is not yet available in Fortran.

2468:   Level: intermediate

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

2498:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2499:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
2500:   PetscDSEnlarge_Static(prob, PetscMax(f, g)+1);
2501:   prob->gpBd[(f*prob->Nf + g)*4+0] = g0;
2502:   prob->gpBd[(f*prob->Nf + g)*4+1] = g1;
2503:   prob->gpBd[(f*prob->Nf + g)*4+2] = g2;
2504:   prob->gpBd[(f*prob->Nf + g)*4+3] = g3;
2505:   return(0);
2506: }

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

2511:   Not collective

2513:   Input Parameters:
2514: + prob - The PetscDS
2515: - f    - The test field number

2517:   Output Parameter:
2518: + exactSol - exact solution for the test field
2519: - exactCtx - exact solution context

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

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

2525: + dim - the spatial dimension
2526: . t - current time
2527: . x - coordinates of the current point
2528: . Nc - the number of field components
2529: . u - the solution field evaluated at the current point
2530: - ctx - a user context

2532:   Level: intermediate

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

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

2549:   Not collective

2551:   Input Parameters:
2552: + prob - The PetscDS
2553: . f    - The test field number
2554: . sol  - solution function for the test fields
2555: - ctx  - solution context or NULL

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

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

2561: + dim - the spatial dimension
2562: . t - current time
2563: . x - coordinates of the current point
2564: . Nc - the number of field components
2565: . u - the solution field evaluated at the current point
2566: - ctx - a user context

2568:   Level: intermediate

2570: .seealso: PetscDSGetExactSolution()
2571: @*/
2572: PetscErrorCode PetscDSSetExactSolution(PetscDS prob, PetscInt f, PetscErrorCode (*sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void *ctx)
2573: {

2578:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2579:   PetscDSEnlarge_Static(prob, f+1);
2582:   return(0);
2583: }

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

2588:   Not collective

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

2593:   Output Parameters:
2594: + numConstants - The number of constants
2595: - constants    - The array of constants, NULL if there are none

2597:   Level: intermediate

2599: .seealso: PetscDSSetConstants(), PetscDSCreate()
2600: @*/
2601: PetscErrorCode PetscDSGetConstants(PetscDS prob, PetscInt *numConstants, const PetscScalar *constants[])
2602: {
2607:   return(0);
2608: }

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

2613:   Not collective

2615:   Input Parameters:
2616: + prob         - The PetscDS object
2617: . numConstants - The number of constants
2618: - constants    - The array of constants, NULL if there are none

2620:   Level: intermediate

2622: .seealso: PetscDSGetConstants(), PetscDSCreate()
2623: @*/
2624: PetscErrorCode PetscDSSetConstants(PetscDS prob, PetscInt numConstants, PetscScalar constants[])
2625: {

2630:   if (numConstants != prob->numConstants) {
2631:     PetscFree(prob->constants);
2632:     prob->numConstants = numConstants;
2633:     if (prob->numConstants) {
2634:       PetscMalloc1(prob->numConstants, &prob->constants);
2635:     } else {
2636:       prob->constants = NULL;
2637:     }
2638:   }
2639:   if (prob->numConstants) {
2641:     PetscArraycpy(prob->constants, constants, prob->numConstants);
2642:   }
2643:   return(0);
2644: }

2646: /*@
2647:   PetscDSGetFieldIndex - Returns the index of the given field

2649:   Not collective

2651:   Input Parameters:
2652: + prob - The PetscDS object
2653: - disc - The discretization object

2655:   Output Parameter:
2656: . f - The field number

2658:   Level: beginner

2660: .seealso: PetscGetDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
2661: @*/
2662: PetscErrorCode PetscDSGetFieldIndex(PetscDS prob, PetscObject disc, PetscInt *f)
2663: {
2664:   PetscInt g;

2669:   *f = -1;
2670:   for (g = 0; g < prob->Nf; ++g) {if (disc == prob->disc[g]) break;}
2671:   if (g == prob->Nf) SETERRQ(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Field not found in PetscDS.");
2672:   *f = g;
2673:   return(0);
2674: }

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

2679:   Not collective

2681:   Input Parameters:
2682: + prob - The PetscDS object
2683: - f - The field number

2685:   Output Parameter:
2686: . size - The size

2688:   Level: beginner

2690: .seealso: PetscDSGetFieldOffset(), PetscDSGetNumFields(), PetscDSCreate()
2691: @*/
2692: PetscErrorCode PetscDSGetFieldSize(PetscDS prob, PetscInt f, PetscInt *size)
2693: {

2699:   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);
2700:   PetscDSSetUp(prob);
2701:   *size = prob->Nb[f];
2702:   return(0);
2703: }

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

2708:   Not collective

2710:   Input Parameters:
2711: + prob - The PetscDS object
2712: - f - The field number

2714:   Output Parameter:
2715: . off - The offset

2717:   Level: beginner

2719: .seealso: PetscDSGetFieldSize(), PetscDSGetNumFields(), PetscDSCreate()
2720: @*/
2721: PetscErrorCode PetscDSGetFieldOffset(PetscDS prob, PetscInt f, PetscInt *off)
2722: {
2723:   PetscInt       size, g;

2729:   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);
2730:   *off = 0;
2731:   for (g = 0; g < f; ++g) {
2732:     PetscDSGetFieldSize(prob, g, &size);
2733:     *off += size;
2734:   }
2735:   return(0);
2736: }

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

2741:   Not collective

2743:   Input Parameter:
2744: . prob - The PetscDS object

2746:   Output Parameter:
2747: . dimensions - The number of dimensions

2749:   Level: beginner

2751: .seealso: PetscDSGetComponentOffsets(), PetscDSGetNumFields(), PetscDSCreate()
2752: @*/
2753: PetscErrorCode PetscDSGetDimensions(PetscDS prob, PetscInt *dimensions[])
2754: {

2759:   PetscDSSetUp(prob);
2761:   *dimensions = prob->Nb;
2762:   return(0);
2763: }

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

2768:   Not collective

2770:   Input Parameter:
2771: . prob - The PetscDS object

2773:   Output Parameter:
2774: . components - The number of components

2776:   Level: beginner

2778: .seealso: PetscDSGetComponentOffsets(), PetscDSGetNumFields(), PetscDSCreate()
2779: @*/
2780: PetscErrorCode PetscDSGetComponents(PetscDS prob, PetscInt *components[])
2781: {

2786:   PetscDSSetUp(prob);
2788:   *components = prob->Nc;
2789:   return(0);
2790: }

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

2795:   Not collective

2797:   Input Parameters:
2798: + prob - The PetscDS object
2799: - f - The field number

2801:   Output Parameter:
2802: . off - The offset

2804:   Level: beginner

2806: .seealso: PetscDSGetNumFields(), PetscDSCreate()
2807: @*/
2808: PetscErrorCode PetscDSGetComponentOffset(PetscDS prob, PetscInt f, PetscInt *off)
2809: {
2813:   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);
2814:   *off = prob->off[f];
2815:   return(0);
2816: }

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

2821:   Not collective

2823:   Input Parameter:
2824: . prob - The PetscDS object

2826:   Output Parameter:
2827: . offsets - The offsets

2829:   Level: beginner

2831: .seealso: PetscDSGetNumFields(), PetscDSCreate()
2832: @*/
2833: PetscErrorCode PetscDSGetComponentOffsets(PetscDS prob, PetscInt *offsets[])
2834: {
2838:   *offsets = prob->off;
2839:   return(0);
2840: }

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

2845:   Not collective

2847:   Input Parameter:
2848: . prob - The PetscDS object

2850:   Output Parameter:
2851: . offsets - The offsets

2853:   Level: beginner

2855: .seealso: PetscDSGetNumFields(), PetscDSCreate()
2856: @*/
2857: PetscErrorCode PetscDSGetComponentDerivativeOffsets(PetscDS prob, PetscInt *offsets[])
2858: {
2862:   *offsets = prob->offDer;
2863:   return(0);
2864: }

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

2869:   Not collective

2871:   Input Parameter:
2872: . prob - The PetscDS object

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

2877:   Level: intermediate

2879: .seealso: PetscDSCreate()
2880: @*/
2881: PetscErrorCode PetscDSGetTabulation(PetscDS prob, PetscTabulation *T[])
2882: {

2888:   PetscDSSetUp(prob);
2889:   *T = prob->T;
2890:   return(0);
2891: }

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

2896:   Not collective

2898:   Input Parameter:
2899: . prob - The PetscDS object

2901:   Output Parameter:
2902: . Tf - The basis function and derviative tabulation on each lcoal face at quadrature points for each and field

2904:   Level: intermediate

2906: .seealso: PetscDSGetTabulation(), PetscDSCreate()
2907: @*/
2908: PetscErrorCode PetscDSGetFaceTabulation(PetscDS prob, PetscTabulation *Tf[])
2909: {

2915:   PetscDSSetUp(prob);
2916:   *Tf = prob->Tf;
2917:   return(0);
2918: }

2920: PetscErrorCode PetscDSGetEvaluationArrays(PetscDS prob, PetscScalar **u, PetscScalar **u_t, PetscScalar **u_x)
2921: {

2926:   PetscDSSetUp(prob);
2930:   return(0);
2931: }

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

2939:   PetscDSSetUp(prob);
2946:   return(0);
2947: }

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

2955:   PetscDSSetUp(prob);
2961:   return(0);
2962: }

2964: /*@C
2965:   PetscDSAddBoundary - Add a boundary condition to the model

2967:   Collective on ds

2969:   Input Parameters:
2970: + ds          - The PetscDS object
2971: . type        - The type of condition, e.g. DM_BC_ESSENTIAL/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
2972: . name        - The BC name
2973: . labelname   - The label defining constrained points
2974: . field       - The field to constrain
2975: . numcomps    - The number of constrained field components (0 will constrain all fields)
2976: . comps       - An array of constrained component numbers
2977: . bcFunc      - A pointwise function giving boundary values
2978: . numids      - The number of DMLabel ids for constrained points
2979: . ids         - An array of ids for constrained points
2980: - ctx         - An optional user context for bcFunc

2982:   Options Database Keys:
2983: + -bc_<boundary name> <num> - Overrides the boundary ids
2984: - -bc_<boundary name>_comp <num> - Overrides the boundary components

2986:   Level: developer

2988: .seealso: PetscDSGetBoundary()
2989: @*/
2990: PetscErrorCode PetscDSAddBoundary(PetscDS ds, DMBoundaryConditionType type, const char name[], const char labelname[], PetscInt field, PetscInt numcomps, const PetscInt *comps, void (*bcFunc)(void), PetscInt numids, const PetscInt *ids, void *ctx)
2991: {
2992:   DSBoundary     b;

3001:   PetscNew(&b);
3002:   PetscStrallocpy(name, (char **) &b->name);
3003:   PetscStrallocpy(labelname, (char **) &b->labelname);
3004:   PetscMalloc1(numcomps, &b->comps);
3005:   if (numcomps) {PetscArraycpy(b->comps, comps, numcomps);}
3006:   PetscMalloc1(numids, &b->ids);
3007:   if (numids) {PetscArraycpy(b->ids, ids, numids);}
3008:   b->type            = type;
3009:   b->field           = field;
3010:   b->numcomps        = numcomps;
3011:   b->func            = bcFunc;
3012:   b->numids          = numids;
3013:   b->ctx             = ctx;
3014:   b->next            = ds->boundary;
3015:   ds->boundary       = b;
3016:   return(0);
3017: }

3019: /*@C
3020:   PetscDSUpdateBoundary - Change a boundary condition for the model

3022:   Input Parameters:
3023: + ds          - The PetscDS object
3024: . bd          - The boundary condition number
3025: . type        - The type of condition, e.g. DM_BC_ESSENTIAL/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
3026: . name        - The BC name
3027: . labelname   - The label defining constrained points
3028: . field       - The field to constrain
3029: . numcomps    - The number of constrained field components
3030: . comps       - An array of constrained component numbers
3031: . bcFunc      - A pointwise function giving boundary values
3032: . numids      - The number of DMLabel ids for constrained points
3033: . ids         - An array of ids for constrained points
3034: - ctx         - An optional user context for bcFunc

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

3038:   Level: developer

3040: .seealso: PetscDSAddBoundary(), PetscDSGetBoundary(), PetscDSGetNumBoundary()
3041: @*/
3042: PetscErrorCode PetscDSUpdateBoundary(PetscDS ds, PetscInt bd, DMBoundaryConditionType type, const char name[], const char labelname[], PetscInt field, PetscInt numcomps, const PetscInt *comps, void (*bcFunc)(void), PetscInt numids, const PetscInt *ids, void *ctx)
3043: {
3044:   DSBoundary     b = ds->boundary;
3045:   PetscInt       n = 0;

3050:   while (b) {
3051:     if (n == bd) break;
3052:     b = b->next;
3053:     ++n;
3054:   }
3055:   if (!b) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %d is not in [0, %d)", bd, n);
3056:   if (name) {
3057:     PetscFree(b->name);
3058:     PetscStrallocpy(name, (char **) &b->name);
3059:   }
3060:   if (labelname) {
3061:     PetscFree(b->labelname);
3062:     PetscStrallocpy(labelname, (char **) &b->labelname);
3063:   }
3064:   if (numcomps >= 0 && numcomps != b->numcomps) {
3065:     b->numcomps = numcomps;
3066:     PetscFree(b->comps);
3067:     PetscMalloc1(numcomps, &b->comps);
3068:     if (numcomps) {PetscArraycpy(b->comps, comps, numcomps);}
3069:   }
3070:   if (numids >= 0 && numids != b->numids) {
3071:     b->numids = numids;
3072:     PetscFree(b->ids);
3073:     PetscMalloc1(numids, &b->ids);
3074:     if (numids) {PetscArraycpy(b->ids, ids, numids);}
3075:   }
3076:   b->type = type;
3077:   if (field >= 0) {b->field  = field;}
3078:   if (bcFunc)     {b->func   = bcFunc;}
3079:   if (ctx)        {b->ctx    = ctx;}
3080:   return(0);
3081: }

3083: /*@
3084:   PetscDSGetNumBoundary - Get the number of registered BC

3086:   Input Parameters:
3087: . ds - The PetscDS object

3089:   Output Parameters:
3090: . numBd - The number of BC

3092:   Level: intermediate

3094: .seealso: PetscDSAddBoundary(), PetscDSGetBoundary()
3095: @*/
3096: PetscErrorCode PetscDSGetNumBoundary(PetscDS ds, PetscInt *numBd)
3097: {
3098:   DSBoundary b = ds->boundary;

3103:   *numBd = 0;
3104:   while (b) {++(*numBd); b = b->next;}
3105:   return(0);
3106: }

3108: /*@C
3109:   PetscDSGetBoundary - Gets a boundary condition to the model

3111:   Input Parameters:
3112: + ds          - The PetscDS object
3113: - bd          - The BC number

3115:   Output Parameters:
3116: + type        - The type of condition, e.g. DM_BC_ESSENTIAL/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
3117: . name        - The BC name
3118: . labelname   - The label defining constrained points
3119: . field       - The field to constrain
3120: . numcomps    - The number of constrained field components
3121: . comps       - An array of constrained component numbers
3122: . bcFunc      - A pointwise function giving boundary values
3123: . numids      - The number of DMLabel ids for constrained points
3124: . ids         - An array of ids for constrained points
3125: - ctx         - An optional user context for bcFunc

3127:   Options Database Keys:
3128: + -bc_<boundary name> <num> - Overrides the boundary ids
3129: - -bc_<boundary name>_comp <num> - Overrides the boundary components

3131:   Level: developer

3133: .seealso: PetscDSAddBoundary()
3134: @*/
3135: PetscErrorCode PetscDSGetBoundary(PetscDS ds, PetscInt bd, DMBoundaryConditionType *type, const char **name, const char **labelname, PetscInt *field, PetscInt *numcomps, const PetscInt **comps, void (**func)(void), PetscInt *numids, const PetscInt **ids, void **ctx)
3136: {
3137:   DSBoundary b    = ds->boundary;
3138:   PetscInt   n    = 0;

3142:   while (b) {
3143:     if (n == bd) break;
3144:     b = b->next;
3145:     ++n;
3146:   }
3147:   if (!b) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %d is not in [0, %d)", bd, n);
3148:   if (type) {
3150:     *type = b->type;
3151:   }
3152:   if (name) {
3154:     *name = b->name;
3155:   }
3156:   if (labelname) {
3158:     *labelname = b->labelname;
3159:   }
3160:   if (field) {
3162:     *field = b->field;
3163:   }
3164:   if (numcomps) {
3166:     *numcomps = b->numcomps;
3167:   }
3168:   if (comps) {
3170:     *comps = b->comps;
3171:   }
3172:   if (func) {
3174:     *func = b->func;
3175:   }
3176:   if (numids) {
3178:     *numids = b->numids;
3179:   }
3180:   if (ids) {
3182:     *ids = b->ids;
3183:   }
3184:   if (ctx) {
3186:     *ctx = b->ctx;
3187:   }
3188:   return(0);
3189: }

3191: /*@
3192:   PetscDSCopyBoundary - Copy all boundary condition objects to the new problem

3194:   Not collective

3196:   Input Parameter:
3197: . prob - The PetscDS object

3199:   Output Parameter:
3200: . newprob - The PetscDS copy

3202:   Level: intermediate

3204: .seealso: PetscDSCopyEquations(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3205: @*/
3206: PetscErrorCode PetscDSCopyBoundary(PetscDS probA, PetscDS probB)
3207: {
3208:   DSBoundary     b, next, *lastnext;

3214:   if (probA == probB) return(0);
3215:   next = probB->boundary;
3216:   while (next) {
3217:     DSBoundary b = next;

3219:     next = b->next;
3220:     PetscFree(b->comps);
3221:     PetscFree(b->ids);
3222:     PetscFree(b->name);
3223:     PetscFree(b->labelname);
3224:     PetscFree(b);
3225:   }
3226:   lastnext = &(probB->boundary);
3227:   for (b = probA->boundary; b; b = b->next) {
3228:     DSBoundary bNew;

3230:     PetscNew(&bNew);
3231:     bNew->numcomps = b->numcomps;
3232:     PetscMalloc1(bNew->numcomps, &bNew->comps);
3233:     PetscArraycpy(bNew->comps, b->comps, bNew->numcomps);
3234:     bNew->numids = b->numids;
3235:     PetscMalloc1(bNew->numids, &bNew->ids);
3236:     PetscArraycpy(bNew->ids, b->ids, bNew->numids);
3237:     PetscStrallocpy(b->labelname,(char **) &(bNew->labelname));
3238:     PetscStrallocpy(b->name,(char **) &(bNew->name));
3239:     bNew->ctx   = b->ctx;
3240:     bNew->type  = b->type;
3241:     bNew->field = b->field;
3242:     bNew->func  = b->func;

3244:     *lastnext = bNew;
3245:     lastnext = &(bNew->next);
3246:   }
3247:   return(0);
3248: }

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

3253:   Not collective

3255:   Input Parameter:
3256: + prob - The PetscDS object
3257: . numFields - Number of new fields
3258: - fields - Old field number for each new field

3260:   Output Parameter:
3261: . newprob - The PetscDS copy

3263:   Level: intermediate

3265: .seealso: PetscDSCopyBoundary(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3266: @*/
3267: PetscErrorCode PetscDSSelectEquations(PetscDS prob, PetscInt numFields, const PetscInt fields[], PetscDS newprob)
3268: {
3269:   PetscInt       Nf, Nfn, fn, gn;

3276:   PetscDSGetNumFields(prob, &Nf);
3277:   PetscDSGetNumFields(newprob, &Nfn);
3278:   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);
3279:   for (fn = 0; fn < numFields; ++fn) {
3280:     const PetscInt   f = fields ? fields[fn] : fn;
3281:     PetscPointFunc   obj;
3282:     PetscPointFunc   f0, f1;
3283:     PetscBdPointFunc f0Bd, f1Bd;
3284:     PetscRiemannFunc r;

3286:     if (f >= Nf) continue;
3287:     PetscDSGetObjective(prob, f, &obj);
3288:     PetscDSGetResidual(prob, f, &f0, &f1);
3289:     PetscDSGetBdResidual(prob, f, &f0Bd, &f1Bd);
3290:     PetscDSGetRiemannSolver(prob, f, &r);
3291:     PetscDSSetObjective(newprob, fn, obj);
3292:     PetscDSSetResidual(newprob, fn, f0, f1);
3293:     PetscDSSetBdResidual(newprob, fn, f0Bd, f1Bd);
3294:     PetscDSSetRiemannSolver(newprob, fn, r);
3295:     for (gn = 0; gn < numFields; ++gn) {
3296:       const PetscInt  g = fields ? fields[gn] : gn;
3297:       PetscPointJac   g0, g1, g2, g3;
3298:       PetscPointJac   g0p, g1p, g2p, g3p;
3299:       PetscBdPointJac g0Bd, g1Bd, g2Bd, g3Bd;

3301:       if (g >= Nf) continue;
3302:       PetscDSGetJacobian(prob, f, g, &g0, &g1, &g2, &g3);
3303:       PetscDSGetJacobianPreconditioner(prob, f, g, &g0p, &g1p, &g2p, &g3p);
3304:       PetscDSGetBdJacobian(prob, f, g, &g0Bd, &g1Bd, &g2Bd, &g3Bd);
3305:       PetscDSSetJacobian(newprob, fn, gn, g0, g1, g2, g3);
3306:       PetscDSSetJacobianPreconditioner(prob, fn, gn, g0p, g1p, g2p, g3p);
3307:       PetscDSSetBdJacobian(newprob, fn, gn, g0Bd, g1Bd, g2Bd, g3Bd);
3308:     }
3309:   }
3310:   return(0);
3311: }

3313: /*@
3314:   PetscDSCopyEquations - Copy all pointwise function pointers to the new problem

3316:   Not collective

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

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

3324:   Level: intermediate

3326: .seealso: PetscDSCopyBoundary(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3327: @*/
3328: PetscErrorCode PetscDSCopyEquations(PetscDS prob, PetscDS newprob)
3329: {
3330:   PetscInt       Nf, Ng;

3336:   PetscDSGetNumFields(prob, &Nf);
3337:   PetscDSGetNumFields(newprob, &Ng);
3338:   if (Nf != Ng) SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_SIZ, "Number of fields must match %D != %D", Nf, Ng);
3339:   PetscDSSelectEquations(prob, Nf, NULL, newprob);
3340:   return(0);
3341: }
3342: /*@
3343:   PetscDSCopyConstants - Copy all constants to the new problem

3345:   Not collective

3347:   Input Parameter:
3348: . prob - The PetscDS object

3350:   Output Parameter:
3351: . newprob - The PetscDS copy

3353:   Level: intermediate

3355: .seealso: PetscDSCopyBoundary(), PetscDSCopyEquations(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3356: @*/
3357: PetscErrorCode PetscDSCopyConstants(PetscDS prob, PetscDS newprob)
3358: {
3359:   PetscInt           Nc;
3360:   const PetscScalar *constants;
3361:   PetscErrorCode     ierr;

3366:   PetscDSGetConstants(prob, &Nc, &constants);
3367:   PetscDSSetConstants(newprob, Nc, (PetscScalar *) constants);
3368:   return(0);
3369: }

3371: PetscErrorCode PetscDSGetHeightSubspace(PetscDS prob, PetscInt height, PetscDS *subprob)
3372: {
3373:   PetscInt       dim, Nf, f;

3379:   if (height == 0) {*subprob = prob; return(0);}
3380:   PetscDSGetNumFields(prob, &Nf);
3381:   PetscDSGetSpatialDimension(prob, &dim);
3382:   if (height > dim) SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_OUTOFRANGE, "DS can only handle height in [0, %D], not %D", dim, height);
3383:   if (!prob->subprobs) {PetscCalloc1(dim, &prob->subprobs);}
3384:   if (!prob->subprobs[height-1]) {
3385:     PetscInt cdim;

3387:     PetscDSCreate(PetscObjectComm((PetscObject) prob), &prob->subprobs[height-1]);
3388:     PetscDSGetCoordinateDimension(prob, &cdim);
3389:     PetscDSSetCoordinateDimension(prob->subprobs[height-1], cdim);
3390:     for (f = 0; f < Nf; ++f) {
3391:       PetscFE      subfe;
3392:       PetscObject  obj;
3393:       PetscClassId id;

3395:       PetscDSGetDiscretization(prob, f, &obj);
3396:       PetscObjectGetClassId(obj, &id);
3397:       if (id == PETSCFE_CLASSID) {PetscFEGetHeightSubspace((PetscFE) obj, height, &subfe);}
3398:       else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unsupported discretization type for field %d", f);
3399:       PetscDSSetDiscretization(prob->subprobs[height-1], f, (PetscObject) subfe);
3400:     }
3401:   }
3402:   *subprob = prob->subprobs[height-1];
3403:   return(0);
3404: }

3406: PetscErrorCode PetscDSGetDiscType_Internal(PetscDS ds, PetscInt f, PetscDiscType *disctype)
3407: {
3408:   PetscObject    obj;
3409:   PetscClassId   id;
3410:   PetscInt       Nf;

3416:   *disctype = PETSC_DISC_NONE;
3417:   PetscDSGetNumFields(ds, &Nf);
3418:   if (f >= Nf) SETERRQ2(PetscObjectComm((PetscObject) ds), PETSC_ERR_ARG_SIZ, "Field %D must be in [0, %D)", f, Nf);
3419:   PetscDSGetDiscretization(ds, f, &obj);
3420:   if (obj) {
3421:     PetscObjectGetClassId(obj, &id);
3422:     if (id == PETSCFE_CLASSID) *disctype = PETSC_DISC_FE;
3423:     else                       *disctype = PETSC_DISC_FV;
3424:   }
3425:   return(0);
3426: }

3428: static PetscErrorCode PetscDSDestroy_Basic(PetscDS prob)
3429: {
3430:   PetscErrorCode      ierr;

3433:   PetscFree(prob->data);
3434:   return(0);
3435: }

3437: static PetscErrorCode PetscDSInitialize_Basic(PetscDS prob)
3438: {
3440:   prob->ops->setfromoptions = NULL;
3441:   prob->ops->setup          = NULL;
3442:   prob->ops->view           = NULL;
3443:   prob->ops->destroy        = PetscDSDestroy_Basic;
3444:   return(0);
3445: }

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

3450:   Level: intermediate

3452: .seealso: PetscDSType, PetscDSCreate(), PetscDSSetType()
3453: M*/

3455: PETSC_EXTERN PetscErrorCode PetscDSCreate_Basic(PetscDS prob)
3456: {
3457:   PetscDS_Basic *b;
3458:   PetscErrorCode      ierr;

3462:   PetscNewLog(prob, &b);
3463:   prob->data = b;

3465:   PetscDSInitialize_Basic(prob);
3466:   return(0);
3467: }