Actual source code: dtds.c

  1: #include <petsc/private/petscdsimpl.h>

  3: PetscClassId PETSCDS_CLASSID = 0;

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

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

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

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

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

 24:   Not Collective

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

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

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

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

 48:   Level: advanced

 50:    Not available from Fortran

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

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

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

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

 67:   Collective on prob

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

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

 76:   Level: intermediate

 78:    Not available from Fortran

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

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

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

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

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

109:   Not Collective

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

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

117:   Level: intermediate

119:    Not available from Fortran

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

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

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

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

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

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

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

242: /*@C
243:    PetscDSViewFromOptions - View from Options

245:    Collective on PetscDS

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

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

261:   PetscObjectViewFromOptions((PetscObject)A,obj,name);
262:   return(0);
263: }

265: /*@C
266:   PetscDSView - Views a PetscDS

268:   Collective on prob

270:   Input Parameter:
271: + prob - the PetscDS object to view
272: - v  - the viewer

274:   Level: developer

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

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

293: /*@
294:   PetscDSSetFromOptions - sets parameters in a PetscDS from the options database

296:   Collective on prob

298:   Input Parameter:
299: . prob - the PetscDS object to set options for

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

308:   Level: developer

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

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

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

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

371: /*@C
372:   PetscDSSetUp - Construct data structures for the PetscDS

374:   Collective on prob

376:   Input Parameter:
377: . prob - the PetscDS object to setup

379:   Level: developer

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

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

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

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

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

457: static PetscErrorCode PetscDSDestroyStructs_Static(PetscDS prob)
458: {

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

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

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

521: /*@
522:   PetscDSDestroy - Destroys a PetscDS object

524:   Collective on prob

526:   Input Parameter:
527: . prob - the PetscDS object to destroy

529:   Level: developer

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

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

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

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

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

569:   Collective

571:   Input Parameter:
572: . comm - The communicator for the PetscDS object

574:   Output Parameter:
575: . ds   - The PetscDS object

577:   Level: beginner

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

588:   *ds  = NULL;
589:   PetscDSInitializePackage();

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

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

601:   *ds = p;
602:   return(0);
603: }

605: /*@
606:   PetscDSGetNumFields - Returns the number of fields in the DS

608:   Not collective

610:   Input Parameter:
611: . prob - The PetscDS object

613:   Output Parameter:
614: . Nf - The number of fields

616:   Level: beginner

618: .seealso: PetscDSGetSpatialDimension(), PetscDSCreate()
619: @*/
620: PetscErrorCode PetscDSGetNumFields(PetscDS prob, PetscInt *Nf)
621: {
625:   *Nf = prob->Nf;
626:   return(0);
627: }

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

632:   Not collective

634:   Input Parameter:
635: . prob - The PetscDS object

637:   Output Parameter:
638: . dim - The spatial dimension

640:   Level: beginner

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

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

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

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

670:   Not collective

672:   Input Parameter:
673: . prob - The PetscDS object

675:   Output Parameter:
676: . dimEmbed - The coordinate dimension

678:   Level: beginner

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

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

695:   Logically collective on prob

697:   Input Parameters:
698: + prob - The PetscDS object
699: - dimEmbed - The coordinate dimension

701:   Level: beginner

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

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

717:   Not collective

719:   Input Parameter:
720: . prob - The PetscDS object

722:   Output Parameter:
723: . isHybrid - The flag

725:   Level: developer

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

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

741:   Not collective

743:   Input Parameters:
744: + prob - The PetscDS object
745: - isHybrid - The flag

747:   Level: developer

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

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

762:   Not collective

764:   Input Parameter:
765: . prob - The PetscDS object

767:   Output Parameter:
768: . dim - The total problem dimension

770:   Level: beginner

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

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

786: /*@
787:   PetscDSGetTotalComponents - Returns the total number of components in this system

789:   Not collective

791:   Input Parameter:
792: . prob - The PetscDS object

794:   Output Parameter:
795: . dim - The total number of components

797:   Level: beginner

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

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

813: /*@
814:   PetscDSGetDiscretization - Returns the discretization object for the given field

816:   Not collective

818:   Input Parameters:
819: + prob - The PetscDS object
820: - f - The field number

822:   Output Parameter:
823: . disc - The discretization object

825:   Level: beginner

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

839: /*@
840:   PetscDSSetDiscretization - Sets the discretization object for the given field

842:   Not collective

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

849:   Level: beginner

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

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

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

879: /*@
880:   PetscDSGetWeakForm - Returns the weak form object

882:   Not collective

884:   Input Parameter:
885: . ds - The PetscDS object

887:   Output Parameter:
888: . wf - The weak form object

890:   Level: beginner

892: .seealso: PetscDSSetWeakForm(), PetscDSGetNumFields(), PetscDSCreate()
893: @*/
894: PetscErrorCode PetscDSGetWeakForm(PetscDS ds, PetscWeakForm *wf)
895: {
899:   *wf = ds->wf;
900:   return(0);
901: }

903: /*@
904:   PetscDSSetWeakForm - Sets the weak form object

906:   Not collective

908:   Input Parameters:
909: + ds - The PetscDS object
910: - wf - The weak form object

912:   Level: beginner

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

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

930: /*@
931:   PetscDSAddDiscretization - Adds a discretization object

933:   Not collective

935:   Input Parameters:
936: + prob - The PetscDS object
937: - disc - The boundary discretization object

939:   Level: beginner

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

948:   PetscDSSetDiscretization(prob, prob->Nf, disc);
949:   return(0);
950: }

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

955:   Not collective

957:   Input Parameter:
958: . prob - The PetscDS object

960:   Output Parameter:
961: . q - The quadrature object

963: Level: intermediate

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

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

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

987:   Not collective

989:   Input Parameters:
990: + prob - The PetscDS object
991: - f - The field number

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

996:   Level: developer

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

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

1013:   Not collective

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

1020:   Level: developer

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

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

1036:   Not collective

1038:   Input Parameters:
1039: + ds - The PetscDS object
1040: - f  - The field number

1042:   Output Parameter:
1043: . k  - The highest derivative we need to tabulate

1045:   Level: developer

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

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

1062:   Not collective

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

1069:   Level: developer

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

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

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

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

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

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

1120:   Not collective

1122:   Input Parameters:
1123: + ds - The PetscDS
1124: - f  - The test field number

1126:   Output Parameters:
1127: + f0 - integrand for the test function term
1128: - f1 - integrand for the test function gradient term

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

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

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

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

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

1159:   Level: intermediate

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

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

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

1189:   Not collective

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

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

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

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

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

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

1226:   Level: intermediate

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

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

1251: /*@C
1252:   PetscDSHasJacobian - Signals that Jacobian functions have been set

1254:   Not collective

1256:   Input Parameter:
1257: . prob - The PetscDS

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

1262:   Level: intermediate

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

1272:   PetscWeakFormHasJacobian(ds->wf, hasJac);
1273:   return(0);
1274: }

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

1279:   Not collective

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

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

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

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

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

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

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

1322:   Level: intermediate

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

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

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

1363:   Not collective

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

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

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

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

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

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

1404:   Level: intermediate

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

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

1440: /*@C
1441:   PetscDSUseJacobianPreconditioner - Whether to construct a Jacobian preconditioner

1443:   Not collective

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

1449:   Level: intermediate

1451: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1452: @*/
1453: PetscErrorCode PetscDSUseJacobianPreconditioner(PetscDS prob, PetscBool useJacPre)
1454: {
1457:   prob->useJacPre = useJacPre;
1458:   return(0);
1459: }

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

1464:   Not collective

1466:   Input Parameter:
1467: . prob - The PetscDS

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

1472:   Level: intermediate

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

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

1488: /*@C
1489:   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.

1491:   Not collective

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

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

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

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

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

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

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

1534:   Level: intermediate

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

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

1572: /*@C
1573:   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.

1575:   Not collective

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

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

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

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

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

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

1616:   Level: intermediate

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

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

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

1655:   Not collective

1657:   Input Parameter:
1658: . ds - The PetscDS

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

1663:   Level: intermediate

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

1673:   PetscWeakFormHasDynamicJacobian(ds->wf, hasDynJac);
1674:   return(0);
1675: }

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

1680:   Not collective

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

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

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

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

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

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

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

1723:   Level: intermediate

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

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

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

1764:   Not collective

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

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

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

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

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

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

1805:   Level: intermediate

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

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

1841: /*@C
1842:   PetscDSGetRiemannSolver - Returns the Riemann solver for the given field

1844:   Not collective

1846:   Input Arguments:
1847: + ds - The PetscDS object
1848: - f  - The field number

1850:   Output Argument:
1851: . r    - Riemann solver

1853:   Calling sequence for r:

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

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

1868:   Level: intermediate

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

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

1888: /*@C
1889:   PetscDSSetRiemannSolver - Sets the Riemann solver for the given field

1891:   Not collective

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

1898:   Calling sequence for r:

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

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

1913:   Level: intermediate

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

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

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

1933:   Not collective

1935:   Input Parameters:
1936: + ds - The PetscDS
1937: - f  - The field number

1939:   Output Parameters:
1940: . update - update function

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

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

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

1965:   Level: intermediate

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

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

1985:   Not collective

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

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

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

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

2015:   Level: intermediate

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

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

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

2046: PetscErrorCode PetscDSSetContext(PetscDS ds, PetscInt f, void *ctx)
2047: {

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

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

2061:   Not collective

2063:   Input Parameters:
2064: + ds - The PetscDS
2065: - f  - The test field number

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

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

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

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

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

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

2101:   Level: intermediate

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

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

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

2131:   Not collective

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

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

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

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

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

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

2169:   Level: intermediate

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

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

2192: /*@
2193:   PetscDSHasBdJacobian - Signals that boundary Jacobian functions have been set

2195:   Not collective

2197:   Input Parameter:
2198: . ds - The PetscDS

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

2203:   Level: intermediate

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

2214:   PetscWeakFormHasBdJacobian(ds->wf, hasBdJac);
2215:   return(0);
2216: }

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

2221:   Not collective

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

2228:   Output Parameters:
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: PetscDSSetBdJacobian()
2268: @*/
2269: PetscErrorCode PetscDSGetBdJacobian(PetscDS ds, 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: {
2287:   PetscBdPointJac *tmp0, *tmp1, *tmp2, *tmp3;
2288:   PetscInt         n0, n1, n2, n3;
2289:   PetscErrorCode   ierr;

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

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

2306:   Not collective

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

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

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

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

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

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

2348:   Level: intermediate

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

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

2384: /*@
2385:   PetscDSHasBdJacobianPreconditioner - Signals that boundary Jacobian preconditioner functions have been set

2387:   Not collective

2389:   Input Parameter:
2390: . ds - The PetscDS

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

2395:   Level: intermediate

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

2406:   PetscWeakFormHasBdJacobianPreconditioner(ds->wf, hasBdJacPre);
2407:   return(0);
2408: }

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

2413:   Not collective

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

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

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

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

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

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

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

2458:   This is not yet available in Fortran.

2460:   Level: intermediate

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

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

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

2501:   Not collective

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

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

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

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

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

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

2544:   This is not yet available in Fortran.

2546:   Level: intermediate

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

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

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

2585:   Not collective

2587:   Input Parameters:
2588: + prob - The PetscDS
2589: - f    - The test field number

2591:   Output Parameter:
2592: + exactSol - exact solution for the test field
2593: - exactCtx - exact solution context

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

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

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

2606:   Level: intermediate

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

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

2623:   Not collective

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

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

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

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

2642:   Level: intermediate

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

2652:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2653:   PetscDSEnlarge_Static(prob, f+1);
2656:   return(0);
2657: }

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

2662:   Not collective

2664:   Input Parameters:
2665: + prob - The PetscDS
2666: - f    - The test field number

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

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

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

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

2683:   Level: intermediate

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

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

2700:   Not collective

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

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

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

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

2719:   Level: intermediate

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

2729:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2730:   PetscDSEnlarge_Static(prob, f+1);
2733:   return(0);
2734: }

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

2739:   Not collective

2741:   Input Parameter:
2742: . prob - The PetscDS object

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

2748:   Level: intermediate

2750: .seealso: PetscDSSetConstants(), PetscDSCreate()
2751: @*/
2752: PetscErrorCode PetscDSGetConstants(PetscDS prob, PetscInt *numConstants, const PetscScalar *constants[])
2753: {
2758:   return(0);
2759: }

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

2764:   Not collective

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

2771:   Level: intermediate

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

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

2797: /*@
2798:   PetscDSGetFieldIndex - Returns the index of the given field

2800:   Not collective

2802:   Input Parameters:
2803: + prob - The PetscDS object
2804: - disc - The discretization object

2806:   Output Parameter:
2807: . f - The field number

2809:   Level: beginner

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

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

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

2830:   Not collective

2832:   Input Parameters:
2833: + prob - The PetscDS object
2834: - f - The field number

2836:   Output Parameter:
2837: . size - The size

2839:   Level: beginner

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

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

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

2859:   Not collective

2861:   Input Parameters:
2862: + prob - The PetscDS object
2863: - f - The field number

2865:   Output Parameter:
2866: . off - The offset

2868:   Level: beginner

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

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

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

2892:   Not collective

2894:   Input Parameter:
2895: . prob - The PetscDS object

2897:   Output Parameter:
2898: . dimensions - The number of dimensions

2900:   Level: beginner

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

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

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

2919:   Not collective

2921:   Input Parameter:
2922: . prob - The PetscDS object

2924:   Output Parameter:
2925: . components - The number of components

2927:   Level: beginner

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

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

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

2946:   Not collective

2948:   Input Parameters:
2949: + prob - The PetscDS object
2950: - f - The field number

2952:   Output Parameter:
2953: . off - The offset

2955:   Level: beginner

2957: .seealso: PetscDSGetNumFields(), PetscDSCreate()
2958: @*/
2959: PetscErrorCode PetscDSGetComponentOffset(PetscDS prob, PetscInt f, PetscInt *off)
2960: {

2966:   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);
2967:   PetscDSSetUp(prob);
2968:   *off = prob->off[f];
2969:   return(0);
2970: }

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

2975:   Not collective

2977:   Input Parameter:
2978: . prob - The PetscDS object

2980:   Output Parameter:
2981: . offsets - The offsets

2983:   Level: beginner

2985: .seealso: PetscDSGetNumFields(), PetscDSCreate()
2986: @*/
2987: PetscErrorCode PetscDSGetComponentOffsets(PetscDS prob, PetscInt *offsets[])
2988: {

2994:   PetscDSSetUp(prob);
2995:   *offsets = prob->off;
2996:   return(0);
2997: }

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

3002:   Not collective

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

3007:   Output Parameter:
3008: . offsets - The offsets

3010:   Level: beginner

3012: .seealso: PetscDSGetNumFields(), PetscDSCreate()
3013: @*/
3014: PetscErrorCode PetscDSGetComponentDerivativeOffsets(PetscDS prob, PetscInt *offsets[])
3015: {

3021:   PetscDSSetUp(prob);
3022:   *offsets = prob->offDer;
3023:   return(0);
3024: }

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

3029:   Not collective

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

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

3037:   Level: intermediate

3039: .seealso: PetscDSCreate()
3040: @*/
3041: PetscErrorCode PetscDSGetTabulation(PetscDS prob, PetscTabulation *T[])
3042: {

3048:   PetscDSSetUp(prob);
3049:   *T = prob->T;
3050:   return(0);
3051: }

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

3056:   Not collective

3058:   Input Parameter:
3059: . prob - The PetscDS object

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

3064:   Level: intermediate

3066: .seealso: PetscDSGetTabulation(), PetscDSCreate()
3067: @*/
3068: PetscErrorCode PetscDSGetFaceTabulation(PetscDS prob, PetscTabulation *Tf[])
3069: {

3075:   PetscDSSetUp(prob);
3076:   *Tf = prob->Tf;
3077:   return(0);
3078: }

3080: PetscErrorCode PetscDSGetEvaluationArrays(PetscDS prob, PetscScalar **u, PetscScalar **u_t, PetscScalar **u_x)
3081: {

3086:   PetscDSSetUp(prob);
3090:   return(0);
3091: }

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

3099:   PetscDSSetUp(prob);
3106:   return(0);
3107: }

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

3115:   PetscDSSetUp(prob);
3121:   return(0);
3122: }

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

3127:   Collective on ds

3129:   Input Parameters:
3130: + ds       - The PetscDS object
3131: . type     - The type of condition, e.g. DM_BC_ESSENTIAL/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
3132: . name     - The BC name
3133: . label    - The label defining constrained points
3134: . Nv       - The number of DMLabel values for constrained points
3135: . values   - An array of label values for constrained points
3136: . field    - The field to constrain
3137: . Nc       - The number of constrained field components (0 will constrain all fields)
3138: . comps    - An array of constrained component numbers
3139: . bcFunc   - A pointwise function giving boundary values
3140: . bcFunc_t - A pointwise function giving the time derviative of the boundary values, or NULL
3141: - ctx      - An optional user context for bcFunc

3143:   Output Parameters:
3144: - bd       - The boundary number

3146:   Options Database Keys:
3147: + -bc_<boundary name> <num> - Overrides the boundary ids
3148: - -bc_<boundary name>_comp <num> - Overrides the boundary components

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

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

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

3157: $ bcFunc(PetscInt dim, PetscInt Nf, PetscInt NfAux,
3158: $        const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
3159: $        const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
3160: $        PetscReal time, const PetscReal x[], PetscScalar bcval[])

3162: + dim - the spatial dimension
3163: . Nf - the number of fields
3164: . uOff - the offset into u[] and u_t[] for each field
3165: . uOff_x - the offset into u_x[] for each field
3166: . u - each field evaluated at the current point
3167: . u_t - the time derivative of each field evaluated at the current point
3168: . u_x - the gradient of each field evaluated at the current point
3169: . aOff - the offset into a[] and a_t[] for each auxiliary field
3170: . aOff_x - the offset into a_x[] for each auxiliary field
3171: . a - each auxiliary field evaluated at the current point
3172: . a_t - the time derivative of each auxiliary field evaluated at the current point
3173: . a_x - the gradient of auxiliary each field evaluated at the current point
3174: . t - current time
3175: . x - coordinates of the current point
3176: . numConstants - number of constant parameters
3177: . constants - constant parameters
3178: - bcval - output values at the current point

3180:   Level: developer

3182: .seealso: PetscDSAddBoundaryByName(), PetscDSGetBoundary(), PetscDSSetResidual(), PetscDSSetBdResidual()
3183: @*/
3184: PetscErrorCode PetscDSAddBoundary(PetscDS ds, DMBoundaryConditionType type, const char name[], DMLabel label, PetscInt Nv, const PetscInt values[], PetscInt field, PetscInt Nc, const PetscInt comps[], void (*bcFunc)(void), void (*bcFunc_t)(void), void *ctx, PetscInt *bd)
3185: {
3186:   DSBoundary     head = ds->boundary, b;
3187:   PetscInt       n    = 0;
3188:   const char    *lname;

3199:   if (Nc > 0) {
3200:     PetscInt *fcomps;
3201:     PetscInt  c;

3203:     PetscDSGetComponents(ds, &fcomps);
3204:     if (Nc > fcomps[field]) SETERRQ3(PetscObjectComm((PetscObject) ds), PETSC_ERR_ARG_OUTOFRANGE, "Number of constrained components %D > %D components for field %D", Nc, fcomps[field], field);
3205:     for (c = 0; c < Nc; ++c) {
3206:       if (comps[c] < 0 || comps[c] >= fcomps[field]) SETERRQ4(PetscObjectComm((PetscObject) ds), PETSC_ERR_ARG_OUTOFRANGE, "Constrained component[%D] %D not in [0, %D) components for field %D", c, comps[c], fcomps[field], field);
3207:     }
3208:   }
3209:   PetscNew(&b);
3210:   PetscStrallocpy(name, (char **) &b->name);
3211:   PetscWeakFormCreate(PETSC_COMM_SELF, &b->wf);
3212:   PetscWeakFormSetNumFields(b->wf, ds->Nf);
3213:   PetscMalloc1(Nv, &b->values);
3214:   if (Nv) {PetscArraycpy(b->values, values, Nv);}
3215:   PetscMalloc1(Nc, &b->comps);
3216:   if (Nc) {PetscArraycpy(b->comps, comps, Nc);}
3217:   PetscObjectGetName((PetscObject) label, &lname);
3218:   PetscStrallocpy(lname, (char **) &b->lname);
3219:   b->type   = type;
3220:   b->label  = label;
3221:   b->Nv     = Nv;
3222:   b->field  = field;
3223:   b->Nc     = Nc;
3224:   b->func   = bcFunc;
3225:   b->func_t = bcFunc_t;
3226:   b->ctx    = ctx;
3227:   b->next   = NULL;
3228:   /* Append to linked list so that we can preserve the order */
3229:   if (!head) ds->boundary = b;
3230:   while (head) {
3231:     if (!head->next) {
3232:       head->next = b;
3233:       head       = b;
3234:     }
3235:     head = head->next;
3236:     ++n;
3237:   }
3239:   return(0);
3240: }

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

3245:   Collective on ds

3247:   Input Parameters:
3248: + ds       - The PetscDS object
3249: . type     - The type of condition, e.g. DM_BC_ESSENTIAL/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
3250: . name     - The BC name
3251: . lname    - The naem of the label defining constrained points
3252: . Nv       - The number of DMLabel values for constrained points
3253: . values   - An array of label values for constrained points
3254: . field    - The field to constrain
3255: . Nc       - The number of constrained field components (0 will constrain all fields)
3256: . comps    - An array of constrained component numbers
3257: . bcFunc   - A pointwise function giving boundary values
3258: . bcFunc_t - A pointwise function giving the time derviative of the boundary values, or NULL
3259: - ctx      - An optional user context for bcFunc

3261:   Output Parameters:
3262: - bd       - The boundary number

3264:   Options Database Keys:
3265: + -bc_<boundary name> <num> - Overrides the boundary ids
3266: - -bc_<boundary name>_comp <num> - Overrides the boundary components

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

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

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

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

3277: $ bcFunc(PetscInt dim, PetscInt Nf, PetscInt NfAux,
3278: $        const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
3279: $        const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
3280: $        PetscReal time, const PetscReal x[], PetscScalar bcval[])

3282: + dim - the spatial dimension
3283: . Nf - the number of fields
3284: . uOff - the offset into u[] and u_t[] for each field
3285: . uOff_x - the offset into u_x[] for each field
3286: . u - each field evaluated at the current point
3287: . u_t - the time derivative of each field evaluated at the current point
3288: . u_x - the gradient of each field evaluated at the current point
3289: . aOff - the offset into a[] and a_t[] for each auxiliary field
3290: . aOff_x - the offset into a_x[] for each auxiliary field
3291: . a - each auxiliary field evaluated at the current point
3292: . a_t - the time derivative of each auxiliary field evaluated at the current point
3293: . a_x - the gradient of auxiliary each field evaluated at the current point
3294: . t - current time
3295: . x - coordinates of the current point
3296: . numConstants - number of constant parameters
3297: . constants - constant parameters
3298: - bcval - output values at the current point

3300:   Level: developer

3302: .seealso: PetscDSAddBoundary(), PetscDSGetBoundary(), PetscDSSetResidual(), PetscDSSetBdResidual()
3303: @*/
3304: PetscErrorCode PetscDSAddBoundaryByName(PetscDS ds, DMBoundaryConditionType type, const char name[], const char lname[], PetscInt Nv, const PetscInt values[], PetscInt field, PetscInt Nc, const PetscInt comps[], void (*bcFunc)(void), void (*bcFunc_t)(void), void *ctx, PetscInt *bd)
3305: {
3306:   DSBoundary     head = ds->boundary, b;
3307:   PetscInt       n    = 0;

3318:   PetscNew(&b);
3319:   PetscStrallocpy(name, (char **) &b->name);
3320:   PetscWeakFormCreate(PETSC_COMM_SELF, &b->wf);
3321:   PetscWeakFormSetNumFields(b->wf, ds->Nf);
3322:   PetscMalloc1(Nv, &b->values);
3323:   if (Nv) {PetscArraycpy(b->values, values, Nv);}
3324:   PetscMalloc1(Nc, &b->comps);
3325:   if (Nc) {PetscArraycpy(b->comps, comps, Nc);}
3326:   PetscStrallocpy(lname, (char **) &b->lname);
3327:   b->type   = type;
3328:   b->label  = NULL;
3329:   b->Nv     = Nv;
3330:   b->field  = field;
3331:   b->Nc     = Nc;
3332:   b->func   = bcFunc;
3333:   b->func_t = bcFunc_t;
3334:   b->ctx    = ctx;
3335:   b->next   = NULL;
3336:   /* Append to linked list so that we can preserve the order */
3337:   if (!head) ds->boundary = b;
3338:   while (head) {
3339:     if (!head->next) {
3340:       head->next = b;
3341:       head       = b;
3342:     }
3343:     head = head->next;
3344:     ++n;
3345:   }
3347:   return(0);
3348: }

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

3353:   Input Parameters:
3354: + ds       - The PetscDS object
3355: . bd       - The boundary condition number
3356: . type     - The type of condition, e.g. DM_BC_ESSENTIAL/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
3357: . name     - The BC name
3358: . label    - The label defining constrained points
3359: . Nv       - The number of DMLabel ids for constrained points
3360: . values   - An array of ids for constrained points
3361: . field    - The field to constrain
3362: . Nc       - The number of constrained field components
3363: . comps    - An array of constrained component numbers
3364: . bcFunc   - A pointwise function giving boundary values
3365: . bcFunc_t - A pointwise function giving the time derviative of the boundary values, or NULL
3366: - ctx      - An optional user context for bcFunc

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

3371:   Level: developer

3373: .seealso: PetscDSAddBoundary(), PetscDSGetBoundary(), PetscDSGetNumBoundary()
3374: @*/
3375: PetscErrorCode PetscDSUpdateBoundary(PetscDS ds, PetscInt bd, DMBoundaryConditionType type, const char name[], DMLabel label, PetscInt Nv, const PetscInt values[], PetscInt field, PetscInt Nc, const PetscInt comps[], void (*bcFunc)(void), void (*bcFunc_t)(void), void *ctx)
3376: {
3377:   DSBoundary     b = ds->boundary;
3378:   PetscInt       n = 0;

3383:   while (b) {
3384:     if (n == bd) break;
3385:     b = b->next;
3386:     ++n;
3387:   }
3388:   if (!b) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %d is not in [0, %d)", bd, n);
3389:   if (name) {
3390:     PetscFree(b->name);
3391:     PetscStrallocpy(name, (char **) &b->name);
3392:   }
3393:   b->type = type;
3394:   if (label) {
3395:     const char *name;

3397:     b->label = label;
3398:     PetscFree(b->lname);
3399:     PetscObjectGetName((PetscObject) label, &name);
3400:     PetscStrallocpy(name, (char **) &b->lname);
3401:   }
3402:   if (Nv >= 0) {
3403:     b->Nv = Nv;
3404:     PetscFree(b->values);
3405:     PetscMalloc1(Nv, &b->values);
3406:     if (Nv) {PetscArraycpy(b->values, values, Nv);}
3407:   }
3408:   if (field >= 0) b->field = field;
3409:   if (Nc >= 0) {
3410:     b->Nc = Nc;
3411:     PetscFree(b->comps);
3412:     PetscMalloc1(Nc, &b->comps);
3413:     if (Nc) {PetscArraycpy(b->comps, comps, Nc);}
3414:   }
3415:   if (bcFunc)   b->func   = bcFunc;
3416:   if (bcFunc_t) b->func_t = bcFunc_t;
3417:   if (ctx)      b->ctx    = ctx;
3418:   return(0);
3419: }

3421: /*@
3422:   PetscDSGetNumBoundary - Get the number of registered BC

3424:   Input Parameters:
3425: . ds - The PetscDS object

3427:   Output Parameters:
3428: . numBd - The number of BC

3430:   Level: intermediate

3432: .seealso: PetscDSAddBoundary(), PetscDSGetBoundary()
3433: @*/
3434: PetscErrorCode PetscDSGetNumBoundary(PetscDS ds, PetscInt *numBd)
3435: {
3436:   DSBoundary b = ds->boundary;

3441:   *numBd = 0;
3442:   while (b) {++(*numBd); b = b->next;}
3443:   return(0);
3444: }

3446: /*@C
3447:   PetscDSGetBoundary - Gets a boundary condition to the model

3449:   Input Parameters:
3450: + ds          - The PetscDS object
3451: - bd          - The BC number

3453:   Output Parameters:
3454: + wf       - The PetscWeakForm holding the pointwise functions
3455: . type     - The type of condition, e.g. DM_BC_ESSENTIAL/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
3456: . name     - The BC name
3457: . label    - The label defining constrained points
3458: . Nv       - The number of DMLabel ids for constrained points
3459: . values   - An array of ids for constrained points
3460: . field    - The field to constrain
3461: . Nc       - The number of constrained field components
3462: . comps    - An array of constrained component numbers
3463: . bcFunc   - A pointwise function giving boundary values
3464: . bcFunc_t - A pointwise function giving the time derviative of the boundary values
3465: - ctx      - An optional user context for bcFunc

3467:   Options Database Keys:
3468: + -bc_<boundary name> <num> - Overrides the boundary ids
3469: - -bc_<boundary name>_comp <num> - Overrides the boundary components

3471:   Level: developer

3473: .seealso: PetscDSAddBoundary()
3474: @*/
3475: PetscErrorCode PetscDSGetBoundary(PetscDS ds, PetscInt bd, PetscWeakForm *wf, DMBoundaryConditionType *type, const char *name[], DMLabel *label, PetscInt *Nv, const PetscInt *values[], PetscInt *field, PetscInt *Nc, const PetscInt *comps[], void (**func)(void), void (**func_t)(void), void **ctx)
3476: {
3477:   DSBoundary b = ds->boundary;
3478:   PetscInt   n = 0;

3482:   while (b) {
3483:     if (n == bd) break;
3484:     b = b->next;
3485:     ++n;
3486:   }
3487:   if (!b) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %d is not in [0, %d)", bd, n);
3488:   if (wf) {
3490:     *wf = b->wf;
3491:   }
3492:   if (type) {
3494:     *type = b->type;
3495:   }
3496:   if (name) {
3498:     *name = b->name;
3499:   }
3500:   if (label) {
3502:     *label = b->label;
3503:   }
3504:   if (Nv) {
3506:     *Nv = b->Nv;
3507:   }
3508:   if (values) {
3510:     *values = b->values;
3511:   }
3512:   if (field) {
3514:     *field = b->field;
3515:   }
3516:   if (Nc) {
3518:     *Nc = b->Nc;
3519:   }
3520:   if (comps) {
3522:     *comps = b->comps;
3523:   }
3524:   if (func) {
3526:     *func = b->func;
3527:   }
3528:   if (func_t) {
3530:     *func_t = b->func_t;
3531:   }
3532:   if (ctx) {
3534:     *ctx = b->ctx;
3535:   }
3536:   return(0);
3537: }

3539: static PetscErrorCode DSBoundaryDuplicate_Internal(DSBoundary b, DSBoundary *bNew)
3540: {

3544:   PetscNew(bNew);
3545:   PetscWeakFormCreate(PETSC_COMM_SELF, &(*bNew)->wf);
3546:   PetscWeakFormCopy(b->wf, (*bNew)->wf);
3547:   PetscStrallocpy(b->name,(char **) &((*bNew)->name));
3548:   PetscStrallocpy(b->lname,(char **) &((*bNew)->lname));
3549:   (*bNew)->type   = b->type;
3550:   (*bNew)->label  = b->label;
3551:   (*bNew)->Nv     = b->Nv;
3552:   PetscMalloc1(b->Nv, &(*bNew)->values);
3553:   PetscArraycpy((*bNew)->values, b->values, b->Nv);
3554:   (*bNew)->field  = b->field;
3555:   (*bNew)->Nc     = b->Nc;
3556:   PetscMalloc1(b->Nc, &(*bNew)->comps);
3557:   PetscArraycpy((*bNew)->comps, b->comps, b->Nc);
3558:   (*bNew)->func   = b->func;
3559:   (*bNew)->func_t = b->func_t;
3560:   (*bNew)->ctx    = b->ctx;
3561:   return(0);
3562: }

3564: /*@
3565:   PetscDSCopyBoundary - Copy all boundary condition objects to the new problem

3567:   Not collective

3569:   Input Parameters:
3570: + ds        - The source PetscDS object
3571: . numFields - The number of selected fields, or PETSC_DEFAULT for all fields
3572: - fields    - The selected fields, or NULL for all fields

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

3577:   Level: intermediate

3579: .seealso: PetscDSCopyEquations(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3580: @*/
3581: PetscErrorCode PetscDSCopyBoundary(PetscDS ds, PetscInt numFields, const PetscInt fields[], PetscDS newds)
3582: {
3583:   DSBoundary     b, *lastnext;

3589:   if (ds == newds) return(0);
3590:   PetscDSDestroyBoundary(newds);
3591:   lastnext = &(newds->boundary);
3592:   for (b = ds->boundary; b; b = b->next) {
3593:     DSBoundary bNew;
3594:     PetscInt   fieldNew = -1;

3596:     if (numFields > 0 && fields) {
3597:       PetscInt f;

3599:       for (f = 0; f < numFields; ++f) if (b->field == fields[f]) break;
3600:       if (f == numFields) continue;
3601:       fieldNew = f;
3602:     }
3603:     DSBoundaryDuplicate_Internal(b, &bNew);
3604:     bNew->field = fieldNew < 0 ? b->field : fieldNew;
3605:     *lastnext = bNew;
3606:     lastnext  = &(bNew->next);
3607:   }
3608:   return(0);
3609: }

3611: /*@
3612:   PetscDSDestroyBoundary - Remove all DMBoundary objects from the PetscDS

3614:   Not collective

3616:   Input Parameter:
3617: . ds - The PetscDS object

3619:   Level: intermediate

3621: .seealso: PetscDSCopyBoundary(), PetscDSCopyEquations()
3622: @*/
3623: PetscErrorCode PetscDSDestroyBoundary(PetscDS ds)
3624: {
3625:   DSBoundary     next = ds->boundary;

3629:   while (next) {
3630:     DSBoundary b = next;

3632:     next = b->next;
3633:     PetscWeakFormDestroy(&b->wf);
3634:     PetscFree(b->name);
3635:     PetscFree(b->lname);
3636:     PetscFree(b->values);
3637:     PetscFree(b->comps);
3638:     PetscFree(b);
3639:   }
3640:   return(0);
3641: }

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

3646:   Not collective

3648:   Input Parameter:
3649: + prob - The PetscDS object
3650: . numFields - Number of new fields
3651: - fields - Old field number for each new field

3653:   Output Parameter:
3654: . newprob - The PetscDS copy

3656:   Level: intermediate

3658: .seealso: PetscDSSelectEquations(), PetscDSCopyBoundary(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3659: @*/
3660: PetscErrorCode PetscDSSelectDiscretizations(PetscDS prob, PetscInt numFields, const PetscInt fields[], PetscDS newprob)
3661: {
3662:   PetscInt       Nf, Nfn, fn;

3669:   PetscDSGetNumFields(prob, &Nf);
3670:   PetscDSGetNumFields(newprob, &Nfn);
3671:   numFields = numFields < 0 ? Nf : numFields;
3672:   for (fn = 0; fn < numFields; ++fn) {
3673:     const PetscInt f = fields ? fields[fn] : fn;
3674:     PetscObject    disc;

3676:     if (f >= Nf) continue;
3677:     PetscDSGetDiscretization(prob, f, &disc);
3678:     PetscDSSetDiscretization(newprob, fn, disc);
3679:   }
3680:   return(0);
3681: }

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

3686:   Not collective

3688:   Input Parameter:
3689: + prob - The PetscDS object
3690: . numFields - Number of new fields
3691: - fields - Old field number for each new field

3693:   Output Parameter:
3694: . newprob - The PetscDS copy

3696:   Level: intermediate

3698: .seealso: PetscDSSelectDiscretizations(), PetscDSCopyBoundary(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3699: @*/
3700: PetscErrorCode PetscDSSelectEquations(PetscDS prob, PetscInt numFields, const PetscInt fields[], PetscDS newprob)
3701: {
3702:   PetscInt       Nf, Nfn, fn, gn;

3709:   PetscDSGetNumFields(prob, &Nf);
3710:   PetscDSGetNumFields(newprob, &Nfn);
3711:   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);
3712:   for (fn = 0; fn < numFields; ++fn) {
3713:     const PetscInt   f = fields ? fields[fn] : fn;
3714:     PetscPointFunc   obj;
3715:     PetscPointFunc   f0, f1;
3716:     PetscBdPointFunc f0Bd, f1Bd;
3717:     PetscRiemannFunc r;

3719:     if (f >= Nf) continue;
3720:     PetscDSGetObjective(prob, f, &obj);
3721:     PetscDSGetResidual(prob, f, &f0, &f1);
3722:     PetscDSGetBdResidual(prob, f, &f0Bd, &f1Bd);
3723:     PetscDSGetRiemannSolver(prob, f, &r);
3724:     PetscDSSetObjective(newprob, fn, obj);
3725:     PetscDSSetResidual(newprob, fn, f0, f1);
3726:     PetscDSSetBdResidual(newprob, fn, f0Bd, f1Bd);
3727:     PetscDSSetRiemannSolver(newprob, fn, r);
3728:     for (gn = 0; gn < numFields; ++gn) {
3729:       const PetscInt  g = fields ? fields[gn] : gn;
3730:       PetscPointJac   g0, g1, g2, g3;
3731:       PetscPointJac   g0p, g1p, g2p, g3p;
3732:       PetscBdPointJac g0Bd, g1Bd, g2Bd, g3Bd;

3734:       if (g >= Nf) continue;
3735:       PetscDSGetJacobian(prob, f, g, &g0, &g1, &g2, &g3);
3736:       PetscDSGetJacobianPreconditioner(prob, f, g, &g0p, &g1p, &g2p, &g3p);
3737:       PetscDSGetBdJacobian(prob, f, g, &g0Bd, &g1Bd, &g2Bd, &g3Bd);
3738:       PetscDSSetJacobian(newprob, fn, gn, g0, g1, g2, g3);
3739:       PetscDSSetJacobianPreconditioner(newprob, fn, gn, g0p, g1p, g2p, g3p);
3740:       PetscDSSetBdJacobian(newprob, fn, gn, g0Bd, g1Bd, g2Bd, g3Bd);
3741:     }
3742:   }
3743:   return(0);
3744: }

3746: /*@
3747:   PetscDSCopyEquations - Copy all pointwise function pointers to the new problem

3749:   Not collective

3751:   Input Parameter:
3752: . prob - The PetscDS object

3754:   Output Parameter:
3755: . newprob - The PetscDS copy

3757:   Level: intermediate

3759: .seealso: PetscDSCopyBoundary(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3760: @*/
3761: PetscErrorCode PetscDSCopyEquations(PetscDS prob, PetscDS newprob)
3762: {
3763:   PetscInt       Nf, Ng;

3769:   PetscDSGetNumFields(prob, &Nf);
3770:   PetscDSGetNumFields(newprob, &Ng);
3771:   if (Nf != Ng) SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_SIZ, "Number of fields must match %D != %D", Nf, Ng);
3772:   PetscDSSelectEquations(prob, Nf, NULL, newprob);
3773:   return(0);
3774: }

3776: /*@
3777:   PetscDSCopyConstants - Copy all constants to the new problem

3779:   Not collective

3781:   Input Parameter:
3782: . prob - The PetscDS object

3784:   Output Parameter:
3785: . newprob - The PetscDS copy

3787:   Level: intermediate

3789: .seealso: PetscDSCopyBoundary(), PetscDSCopyEquations(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3790: @*/
3791: PetscErrorCode PetscDSCopyConstants(PetscDS prob, PetscDS newprob)
3792: {
3793:   PetscInt           Nc;
3794:   const PetscScalar *constants;
3795:   PetscErrorCode     ierr;

3800:   PetscDSGetConstants(prob, &Nc, &constants);
3801:   PetscDSSetConstants(newprob, Nc, (PetscScalar *) constants);
3802:   return(0);
3803: }

3805: /*@
3806:   PetscDSCopyExactSolutions - Copy all exact solutions to the new problem

3808:   Not collective

3810:   Input Parameter:
3811: . ds - The PetscDS object

3813:   Output Parameter:
3814: . newds - The PetscDS copy

3816:   Level: intermediate

3818: .seealso: PetscDSCopyBoundary(), PetscDSCopyEquations(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3819: @*/
3820: PetscErrorCode PetscDSCopyExactSolutions(PetscDS ds, PetscDS newds)
3821: {
3822:   PetscSimplePointFunc sol;
3823:   void                *ctx;
3824:   PetscInt             Nf, f;
3825:   PetscErrorCode       ierr;

3830:   PetscDSGetNumFields(ds, &Nf);
3831:   for (f = 0; f < Nf; ++f) {
3832:     PetscDSGetExactSolution(ds,    f, &sol, &ctx);
3833:     PetscDSSetExactSolution(newds, f,  sol,  ctx);
3834:     PetscDSGetExactSolutionTimeDerivative(ds,    f, &sol, &ctx);
3835:     PetscDSSetExactSolutionTimeDerivative(newds, f,  sol,  ctx);
3836:   }
3837:   return(0);
3838: }

3840: PetscErrorCode PetscDSGetHeightSubspace(PetscDS prob, PetscInt height, PetscDS *subprob)
3841: {
3842:   PetscInt       dim, Nf, f;

3848:   if (height == 0) {*subprob = prob; return(0);}
3849:   PetscDSGetNumFields(prob, &Nf);
3850:   PetscDSGetSpatialDimension(prob, &dim);
3851:   if (height > dim) SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_OUTOFRANGE, "DS can only handle height in [0, %D], not %D", dim, height);
3852:   if (!prob->subprobs) {PetscCalloc1(dim, &prob->subprobs);}
3853:   if (!prob->subprobs[height-1]) {
3854:     PetscInt cdim;

3856:     PetscDSCreate(PetscObjectComm((PetscObject) prob), &prob->subprobs[height-1]);
3857:     PetscDSGetCoordinateDimension(prob, &cdim);
3858:     PetscDSSetCoordinateDimension(prob->subprobs[height-1], cdim);
3859:     for (f = 0; f < Nf; ++f) {
3860:       PetscFE      subfe;
3861:       PetscObject  obj;
3862:       PetscClassId id;

3864:       PetscDSGetDiscretization(prob, f, &obj);
3865:       PetscObjectGetClassId(obj, &id);
3866:       if (id == PETSCFE_CLASSID) {PetscFEGetHeightSubspace((PetscFE) obj, height, &subfe);}
3867:       else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unsupported discretization type for field %d", f);
3868:       PetscDSSetDiscretization(prob->subprobs[height-1], f, (PetscObject) subfe);
3869:     }
3870:   }
3871:   *subprob = prob->subprobs[height-1];
3872:   return(0);
3873: }

3875: PetscErrorCode PetscDSGetDiscType_Internal(PetscDS ds, PetscInt f, PetscDiscType *disctype)
3876: {
3877:   PetscObject    obj;
3878:   PetscClassId   id;
3879:   PetscInt       Nf;

3885:   *disctype = PETSC_DISC_NONE;
3886:   PetscDSGetNumFields(ds, &Nf);
3887:   if (f >= Nf) SETERRQ2(PetscObjectComm((PetscObject) ds), PETSC_ERR_ARG_SIZ, "Field %D must be in [0, %D)", f, Nf);
3888:   PetscDSGetDiscretization(ds, f, &obj);
3889:   if (obj) {
3890:     PetscObjectGetClassId(obj, &id);
3891:     if (id == PETSCFE_CLASSID) *disctype = PETSC_DISC_FE;
3892:     else                       *disctype = PETSC_DISC_FV;
3893:   }
3894:   return(0);
3895: }

3897: static PetscErrorCode PetscDSDestroy_Basic(PetscDS ds)
3898: {

3902:   PetscFree(ds->data);
3903:   return(0);
3904: }

3906: static PetscErrorCode PetscDSInitialize_Basic(PetscDS ds)
3907: {
3909:   ds->ops->setfromoptions = NULL;
3910:   ds->ops->setup          = NULL;
3911:   ds->ops->view           = NULL;
3912:   ds->ops->destroy        = PetscDSDestroy_Basic;
3913:   return(0);
3914: }

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

3919:   Level: intermediate

3921: .seealso: PetscDSType, PetscDSCreate(), PetscDSSetType()
3922: M*/

3924: PETSC_EXTERN PetscErrorCode PetscDSCreate_Basic(PetscDS ds)
3925: {
3926:   PetscDS_Basic *b;

3931:   PetscNewLog(ds, &b);
3932:   ds->data = b;

3934:   PetscDSInitialize_Basic(ds);
3935:   return(0);
3936: }