Actual source code: dtds.c

petsc-3.10.5 2019-03-28
Report Typos and Errors
  1:  #include <petsc/private/petscdsimpl.h>

  3: PetscClassId PETSCDS_CLASSID = 0;

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

  8: /*@C
  9:   PetscDSRegister - Adds a new PetscDS implementation

 11:   Not Collective

 13:   Input Parameters:
 14: + name        - The name of a new user-defined creation routine
 15: - create_func - The creation routine itself

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

 20:   Sample usage:
 21: .vb
 22:     PetscDSRegister("my_ds", MyPetscDSCreate);
 23: .ve

 25:   Then, your PetscDS type can be chosen with the procedural interface via
 26: .vb
 27:     PetscDSCreate(MPI_Comm, PetscDS *);
 28:     PetscDSSetType(PetscDS, "my_ds");
 29: .ve
 30:    or at runtime via the option
 31: .vb
 32:     -petscds_type my_ds
 33: .ve

 35:   Level: advanced

 37:    Not available from Fortran

 39: .keywords: PetscDS, register
 40: .seealso: PetscDSRegisterAll(), PetscDSRegisterDestroy()

 42: @*/
 43: PetscErrorCode PetscDSRegister(const char sname[], PetscErrorCode (*function)(PetscDS))
 44: {

 48:   PetscFunctionListAdd(&PetscDSList, sname, function);
 49:   return(0);
 50: }

 52: /*@C
 53:   PetscDSSetType - Builds a particular PetscDS

 55:   Collective on PetscDS

 57:   Input Parameters:
 58: + prob - The PetscDS object
 59: - name - The kind of system

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

 64:   Level: intermediate

 66:    Not available from Fortran

 68: .keywords: PetscDS, set, type
 69: .seealso: PetscDSGetType(), PetscDSCreate()
 70: @*/
 71: PetscErrorCode PetscDSSetType(PetscDS prob, PetscDSType name)
 72: {
 73:   PetscErrorCode (*r)(PetscDS);
 74:   PetscBool      match;

 79:   PetscObjectTypeCompare((PetscObject) prob, name, &match);
 80:   if (match) return(0);

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

 86:   if (prob->ops->destroy) {
 87:     (*prob->ops->destroy)(prob);
 88:     prob->ops->destroy = NULL;
 89:   }
 90:   (*r)(prob);
 91:   PetscObjectChangeTypeName((PetscObject) prob, name);
 92:   return(0);
 93: }

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

 98:   Not Collective

100:   Input Parameter:
101: . prob  - The PetscDS

103:   Output Parameter:
104: . name - The PetscDS type name

106:   Level: intermediate

108:    Not available from Fortran

110: .keywords: PetscDS, get, type, name
111: .seealso: PetscDSSetType(), PetscDSCreate()
112: @*/
113: PetscErrorCode PetscDSGetType(PetscDS prob, PetscDSType *name)
114: {

120:   PetscDSRegisterAll();
121:   *name = ((PetscObject) prob)->type_name;
122:   return(0);
123: }

125: static PetscErrorCode PetscDSView_Ascii(PetscDS prob, PetscViewer viewer)
126: {
127:   PetscViewerFormat  format;
128:   const PetscScalar *constants;
129:   PetscInt           numConstants, f;
130:   PetscErrorCode     ierr;

133:   PetscViewerGetFormat(viewer, &format);
134:   PetscViewerASCIIPrintf(viewer, "Discrete System with %d fields\n", prob->Nf);
135:   PetscViewerASCIIPushTab(viewer);
136:   for (f = 0; f < prob->Nf; ++f) {
137:     PetscObject     obj;
138:     PetscClassId    id;
139:     PetscQuadrature q;
140:     const char     *name;
141:     PetscInt        Nc, Nq, Nqc;

143:     PetscDSGetDiscretization(prob, f, &obj);
144:     PetscObjectGetClassId(obj, &id);
145:     PetscObjectGetName(obj, &name);
146:     PetscViewerASCIIPrintf(viewer, "Field %s", name ? name : "<unknown>");
147:     if (id == PETSCFE_CLASSID)      {
148:       PetscFEGetNumComponents((PetscFE) obj, &Nc);
149:       PetscFEGetQuadrature((PetscFE) obj, &q);
150:       PetscViewerASCIIPrintf(viewer, " FEM");
151:     } else if (id == PETSCFV_CLASSID) {
152:       PetscFVGetNumComponents((PetscFV) obj, &Nc);
153:       PetscFVGetQuadrature((PetscFV) obj, &q);
154:       PetscViewerASCIIPrintf(viewer, " FVM");
155:     }
156:     else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", f);
157:     if (Nc > 1) {PetscViewerASCIIPrintf(viewer, "%D components", Nc);}
158:     else        {PetscViewerASCIIPrintf(viewer, "%D component ", Nc);}
159:     if (prob->implicit[f]) {PetscViewerASCIIPrintf(viewer, " (implicit)");}
160:     else                   {PetscViewerASCIIPrintf(viewer, " (explicit)");}
161:     if (prob->adjacency[f*2+0]) {
162:       if (prob->adjacency[f*2+1]) {PetscViewerASCIIPrintf(viewer, " (adj FVM++)");}
163:       else                        {PetscViewerASCIIPrintf(viewer, " (adj FVM)");}
164:     } else {
165:       if (prob->adjacency[f*2+1]) {PetscViewerASCIIPrintf(viewer, " (adj FEM)");}
166:       else                        {PetscViewerASCIIPrintf(viewer, " (adj FUNKY)");}
167:     }
168:     if (q) {
169:       PetscQuadratureGetData(q, NULL, &Nqc, &Nq, NULL, NULL);
170:       PetscViewerASCIIPrintf(viewer, " (Nq %D Nqc %D)", Nq, Nqc);
171:     }
172:     PetscViewerASCIIPrintf(viewer, "\n");
173:     PetscViewerASCIIPushTab(viewer);
174:     if (id == PETSCFE_CLASSID)      {PetscFEView((PetscFE) obj, viewer);}
175:     else if (id == PETSCFV_CLASSID) {PetscFVView((PetscFV) obj, viewer);}
176:     PetscViewerASCIIPopTab(viewer);
177:   }
178:   PetscDSGetConstants(prob, &numConstants, &constants);
179:   if (numConstants) {
180:     PetscViewerASCIIPrintf(viewer, "%D constants\n", numConstants);
181:     PetscViewerASCIIPushTab(viewer);
182:     for (f = 0; f < numConstants; ++f) {PetscViewerASCIIPrintf(viewer, "%g\n", (double) PetscRealPart(constants[f]));}
183:     PetscViewerASCIIPopTab(viewer);
184:   }
185:   PetscViewerASCIIPopTab(viewer);
186:   return(0);
187: }

189: /*@C
190:   PetscDSView - Views a PetscDS

192:   Collective on PetscDS

194:   Input Parameter:
195: + prob - the PetscDS object to view
196: - v  - the viewer

198:   Level: developer

200: .seealso PetscDSDestroy()
201: @*/
202: PetscErrorCode PetscDSView(PetscDS prob, PetscViewer v)
203: {
204:   PetscBool      iascii;

209:   if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) prob), &v);}
211:   PetscObjectTypeCompare((PetscObject) v, PETSCVIEWERASCII, &iascii);
212:   if (iascii) {PetscDSView_Ascii(prob, v);}
213:   if (prob->ops->view) {(*prob->ops->view)(prob, v);}
214:   return(0);
215: }

217: /*@
218:   PetscDSSetFromOptions - sets parameters in a PetscDS from the options database

220:   Collective on PetscDS

222:   Input Parameter:
223: . prob - the PetscDS object to set options for

225:   Options Database:
226: + -petscds_type <type>     : Set the DS type
227: . -petscds_view <view opt> : View the DS
228: . -petscds_jac_pre         : Turn formation of a separate Jacobian preconditioner on and off
229: . -bc_<name> <ids>         : Specify a list of label ids for a boundary condition
230: - -bc_<name>_comp <comps>  : Specify a list of field components to constrain for a boundary condition

232:   Level: developer

234: .seealso PetscDSView()
235: @*/
236: PetscErrorCode PetscDSSetFromOptions(PetscDS prob)
237: {
238:   DSBoundary     b;
239:   const char    *defaultType;
240:   char           name[256];
241:   PetscBool      flg;

246:   if (!((PetscObject) prob)->type_name) {
247:     defaultType = PETSCDSBASIC;
248:   } else {
249:     defaultType = ((PetscObject) prob)->type_name;
250:   }
251:   PetscDSRegisterAll();

253:   PetscObjectOptionsBegin((PetscObject) prob);
254:   for (b = prob->boundary; b; b = b->next) {
255:     char       optname[1024];
256:     PetscInt   ids[1024], len = 1024;
257:     PetscBool  flg;

259:     PetscSNPrintf(optname, sizeof(optname), "-bc_%s", b->name);
260:     PetscMemzero(ids, sizeof(ids));
261:     PetscOptionsIntArray(optname, "List of boundary IDs", "", ids, &len, &flg);
262:     if (flg) {
263:       b->numids = len;
264:       PetscFree(b->ids);
265:       PetscMalloc1(len, &b->ids);
266:       PetscMemcpy(b->ids, ids, len*sizeof(PetscInt));
267:     }
268:     len = 1024;
269:     PetscSNPrintf(optname, sizeof(optname), "-bc_%s_comp", b->name);
270:     PetscMemzero(ids, sizeof(ids));
271:     PetscOptionsIntArray(optname, "List of boundary field components", "", ids, &len, &flg);
272:     if (flg) {
273:       b->numcomps = len;
274:       PetscFree(b->comps);
275:       PetscMalloc1(len, &b->comps);
276:       PetscMemcpy(b->comps, ids, len*sizeof(PetscInt));
277:     }
278:   }
279:   PetscOptionsFList("-petscds_type", "Discrete System", "PetscDSSetType", PetscDSList, defaultType, name, 256, &flg);
280:   if (flg) {
281:     PetscDSSetType(prob, name);
282:   } else if (!((PetscObject) prob)->type_name) {
283:     PetscDSSetType(prob, defaultType);
284:   }
285:   PetscOptionsBool("-petscds_jac_pre", "Discrete System", "PetscDSUseJacobianPreconditioner", prob->useJacPre, &prob->useJacPre, &flg);
286:   if (prob->ops->setfromoptions) {(*prob->ops->setfromoptions)(prob);}
287:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
288:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) prob);
289:   PetscOptionsEnd();
290:   if (prob->Nf) {PetscDSViewFromOptions(prob, NULL, "-petscds_view");}
291:   return(0);
292: }

294: /*@C
295:   PetscDSSetUp - Construct data structures for the PetscDS

297:   Collective on PetscDS

299:   Input Parameter:
300: . prob - the PetscDS object to setup

302:   Level: developer

304: .seealso PetscDSView(), PetscDSDestroy()
305: @*/
306: PetscErrorCode PetscDSSetUp(PetscDS prob)
307: {
308:   const PetscInt Nf = prob->Nf;
309:   PetscInt       dim, dimEmbed, work, NcMax = 0, NqMax = 0, f;

314:   if (prob->setup) return(0);
315:   /* Calculate sizes */
316:   PetscDSGetSpatialDimension(prob, &dim);
317:   PetscDSGetCoordinateDimension(prob, &dimEmbed);
318:   prob->totDim = prob->totComp = 0;
319:   PetscMalloc2(Nf,&prob->Nc,Nf,&prob->Nb);
320:   PetscCalloc2(Nf+1,&prob->off,Nf+1,&prob->offDer);
321:   PetscMalloc4(Nf,&prob->basis,Nf,&prob->basisDer,Nf,&prob->basisFace,Nf,&prob->basisDerFace);
322:   for (f = 0; f < Nf; ++f) {
323:     PetscObject     obj;
324:     PetscClassId    id;
325:     PetscQuadrature q;
326:     PetscInt        Nq = 0, Nb, Nc;

328:     PetscDSGetDiscretization(prob, f, &obj);
329:     PetscObjectGetClassId(obj, &id);
330:     if (id == PETSCFE_CLASSID)      {
331:       PetscFE fe = (PetscFE) obj;

333:       PetscFEGetQuadrature(fe, &q);
334:       PetscFEGetDimension(fe, &Nb);
335:       PetscFEGetNumComponents(fe, &Nc);
336:       PetscFEGetDefaultTabulation(fe, &prob->basis[f], &prob->basisDer[f], NULL);
337:       PetscFEGetFaceTabulation(fe, &prob->basisFace[f], &prob->basisDerFace[f], NULL);
338:     } else if (id == PETSCFV_CLASSID) {
339:       PetscFV fv = (PetscFV) obj;

341:       PetscFVGetQuadrature(fv, &q);
342:       PetscFVGetNumComponents(fv, &Nc);
343:       Nb   = Nc;
344:       PetscFVGetDefaultTabulation(fv, &prob->basis[f], &prob->basisDer[f], NULL);
345:       /* TODO: should PetscFV also have face tabulation? Otherwise there will be a null pointer in prob->basisFace */
346:     } else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
347:     prob->Nc[f]       = Nc;
348:     prob->Nb[f]       = Nb;
349:     prob->off[f+1]    = Nc     + prob->off[f];
350:     prob->offDer[f+1] = Nc*dim + prob->offDer[f];
351:     if (q) {PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);}
352:     NqMax          = PetscMax(NqMax, Nq);
353:     NcMax          = PetscMax(NcMax, Nc);
354:     prob->totDim  += Nb;
355:     prob->totComp += Nc;
356:   }
357:   work = PetscMax(prob->totComp*dim, PetscSqr(NcMax*dim));
358:   /* Allocate works space */
359:   PetscMalloc5(prob->totComp,&prob->u,prob->totComp,&prob->u_t,prob->totComp*dimEmbed,&prob->u_x,dimEmbed,&prob->x,work,&prob->refSpaceDer);
360:   PetscMalloc6(NqMax*NcMax,&prob->f0,NqMax*NcMax*dim,&prob->f1,NqMax*NcMax*NcMax,&prob->g0,NqMax*NcMax*NcMax*dim,&prob->g1,NqMax*NcMax*NcMax*dim,&prob->g2,NqMax*NcMax*NcMax*dim*dim,&prob->g3);
361:   if (prob->ops->setup) {(*prob->ops->setup)(prob);}
362:   prob->setup = PETSC_TRUE;
363:   return(0);
364: }

366: static PetscErrorCode PetscDSDestroyStructs_Static(PetscDS prob)
367: {

371:   PetscFree2(prob->Nc,prob->Nb);
372:   PetscFree2(prob->off,prob->offDer);
373:   PetscFree4(prob->basis,prob->basisDer,prob->basisFace,prob->basisDerFace);
374:   PetscFree5(prob->u,prob->u_t,prob->u_x,prob->x,prob->refSpaceDer);
375:   PetscFree6(prob->f0,prob->f1,prob->g0,prob->g1,prob->g2,prob->g3);
376:   return(0);
377: }

379: static PetscErrorCode PetscDSEnlarge_Static(PetscDS prob, PetscInt NfNew)
380: {
381:   PetscObject      *tmpd;
382:   PetscBool        *tmpi, *tmpa;
383:   PetscPointFunc   *tmpobj, *tmpf, *tmpup;
384:   PetscPointJac    *tmpg, *tmpgp, *tmpgt;
385:   PetscBdPointFunc *tmpfbd;
386:   PetscBdPointJac  *tmpgbd;
387:   PetscRiemannFunc *tmpr;
388:   PetscSimplePointFunc *tmpexactSol;
389:   void            **tmpctx;
390:   PetscInt          Nf = prob->Nf, f, i;
391:   PetscErrorCode    ierr;

394:   if (Nf >= NfNew) return(0);
395:   prob->setup = PETSC_FALSE;
396:   PetscDSDestroyStructs_Static(prob);
397:   PetscMalloc3(NfNew, &tmpd, NfNew, &tmpi, NfNew*2, &tmpa);
398:   for (f = 0; f < Nf; ++f) {tmpd[f] = prob->disc[f]; tmpi[f] = prob->implicit[f]; for (i = 0; i < 2; ++i) tmpa[f*2+i] = prob->adjacency[f*2+i];}
399:   for (f = Nf; f < NfNew; ++f) {tmpd[f] = NULL; tmpi[f] = PETSC_TRUE; tmpa[f*2+0] = PETSC_FALSE; tmpa[f*2+1] = PETSC_TRUE;}
400:   PetscFree3(prob->disc, prob->implicit, prob->adjacency);
401:   prob->Nf        = NfNew;
402:   prob->disc      = tmpd;
403:   prob->implicit  = tmpi;
404:   prob->adjacency = tmpa;
405:   PetscCalloc7(NfNew, &tmpobj, NfNew*2, &tmpf, NfNew*NfNew*4, &tmpg, NfNew*NfNew*4, &tmpgp, NfNew*NfNew*4, &tmpgt, NfNew, &tmpr, NfNew, &tmpctx);
406:   PetscCalloc1(NfNew, &tmpup);
407:   for (f = 0; f < Nf; ++f) tmpobj[f] = prob->obj[f];
408:   for (f = 0; f < Nf*2; ++f) tmpf[f] = prob->f[f];
409:   for (f = 0; f < Nf*Nf*4; ++f) tmpg[f] = prob->g[f];
410:   for (f = 0; f < Nf*Nf*4; ++f) tmpgp[f] = prob->gp[f];
411:   for (f = 0; f < Nf; ++f) tmpr[f] = prob->r[f];
412:   for (f = 0; f < Nf; ++f) tmpup[f] = prob->update[f];
413:   for (f = 0; f < Nf; ++f) tmpctx[f] = prob->ctx[f];
414:   for (f = Nf; f < NfNew; ++f) tmpobj[f] = NULL;
415:   for (f = Nf*2; f < NfNew*2; ++f) tmpf[f] = NULL;
416:   for (f = Nf*Nf*4; f < NfNew*NfNew*4; ++f) tmpg[f] = NULL;
417:   for (f = Nf*Nf*4; f < NfNew*NfNew*4; ++f) tmpgp[f] = NULL;
418:   for (f = Nf*Nf*4; f < NfNew*NfNew*4; ++f) tmpgt[f] = NULL;
419:   for (f = Nf; f < NfNew; ++f) tmpr[f] = NULL;
420:   for (f = Nf; f < NfNew; ++f) tmpup[f] = NULL;
421:   for (f = Nf; f < NfNew; ++f) tmpctx[f] = NULL;
422:   PetscFree7(prob->obj, prob->f, prob->g, prob->gp, prob->gt, prob->r, prob->ctx);
423:   PetscFree(prob->update);
424:   prob->obj = tmpobj;
425:   prob->f   = tmpf;
426:   prob->g   = tmpg;
427:   prob->gp  = tmpgp;
428:   prob->gt  = tmpgt;
429:   prob->r   = tmpr;
430:   prob->update = tmpup;
431:   prob->ctx = tmpctx;
432:   PetscCalloc3(NfNew*2, &tmpfbd, NfNew*NfNew*4, &tmpgbd, NfNew, &tmpexactSol);
433:   for (f = 0; f < Nf*2; ++f) tmpfbd[f] = prob->fBd[f];
434:   for (f = 0; f < Nf*Nf*4; ++f) tmpgbd[f] = prob->gBd[f];
435:   for (f = 0; f < Nf; ++f) tmpexactSol[f] = prob->exactSol[f];
436:   for (f = Nf*2; f < NfNew*2; ++f) tmpfbd[f] = NULL;
437:   for (f = Nf*Nf*4; f < NfNew*NfNew*4; ++f) tmpgbd[f] = NULL;
438:   for (f = Nf; f < NfNew; ++f) tmpexactSol[f] = NULL;
439:   PetscFree3(prob->fBd, prob->gBd, prob->exactSol);
440:   prob->fBd = tmpfbd;
441:   prob->gBd = tmpgbd;
442:   prob->exactSol = tmpexactSol;
443:   return(0);
444: }

446: /*@
447:   PetscDSDestroy - Destroys a PetscDS object

449:   Collective on PetscDS

451:   Input Parameter:
452: . prob - the PetscDS object to destroy

454:   Level: developer

456: .seealso PetscDSView()
457: @*/
458: PetscErrorCode PetscDSDestroy(PetscDS *prob)
459: {
460:   PetscInt       f;
461:   DSBoundary     next;

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

468:   if (--((PetscObject)(*prob))->refct > 0) {*prob = 0; return(0);}
469:   ((PetscObject) (*prob))->refct = 0;
470:   if ((*prob)->subprobs) {
471:     PetscInt dim, d;

473:     PetscDSGetSpatialDimension(*prob, &dim);
474:     for (d = 0; d < dim; ++d) {PetscDSDestroy(&(*prob)->subprobs[d]);}
475:   }
476:   PetscFree((*prob)->subprobs);
477:   PetscDSDestroyStructs_Static(*prob);
478:   for (f = 0; f < (*prob)->Nf; ++f) {
479:     PetscObjectDereference((*prob)->disc[f]);
480:   }
481:   PetscFree3((*prob)->disc, (*prob)->implicit, (*prob)->adjacency);
482:   PetscFree7((*prob)->obj,(*prob)->f,(*prob)->g,(*prob)->gp,(*prob)->gt,(*prob)->r,(*prob)->ctx);
483:   PetscFree((*prob)->update);
484:   PetscFree3((*prob)->fBd,(*prob)->gBd,(*prob)->exactSol);
485:   if ((*prob)->ops->destroy) {(*(*prob)->ops->destroy)(*prob);}
486:   next = (*prob)->boundary;
487:   while (next) {
488:     DSBoundary b = next;

490:     next = b->next;
491:     PetscFree(b->comps);
492:     PetscFree(b->ids);
493:     PetscFree(b->name);
494:     PetscFree(b->labelname);
495:     PetscFree(b);
496:   }
497:   PetscFree((*prob)->constants);
498:   PetscHeaderDestroy(prob);
499:   return(0);
500: }

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

505:   Collective on MPI_Comm

507:   Input Parameter:
508: . comm - The communicator for the PetscDS object

510:   Output Parameter:
511: . prob - The PetscDS object

513:   Level: beginner

515: .seealso: PetscDSSetType(), PETSCDSBASIC
516: @*/
517: PetscErrorCode PetscDSCreate(MPI_Comm comm, PetscDS *prob)
518: {
519:   PetscDS   p;

524:   *prob  = NULL;
525:   PetscDSInitializePackage();

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

529:   p->Nf    = 0;
530:   p->setup = PETSC_FALSE;
531:   p->numConstants  = 0;
532:   p->constants     = NULL;
533:   p->dimEmbed      = -1;
534:   p->defaultAdj[0] = PETSC_FALSE;
535:   p->defaultAdj[1] = PETSC_TRUE;
536:   p->useJacPre     = PETSC_TRUE;

538:   *prob = p;
539:   return(0);
540: }

542: /*@
543:   PetscDSGetNumFields - Returns the number of fields in the DS

545:   Not collective

547:   Input Parameter:
548: . prob - The PetscDS object

550:   Output Parameter:
551: . Nf - The number of fields

553:   Level: beginner

555: .seealso: PetscDSGetSpatialDimension(), PetscDSCreate()
556: @*/
557: PetscErrorCode PetscDSGetNumFields(PetscDS prob, PetscInt *Nf)
558: {
562:   *Nf = prob->Nf;
563:   return(0);
564: }

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

569:   Not collective

571:   Input Parameter:
572: . prob - The PetscDS object

574:   Output Parameter:
575: . dim - The spatial dimension

577:   Level: beginner

579: .seealso: PetscDSGetCoordinateDimension(), PetscDSGetNumFields(), PetscDSCreate()
580: @*/
581: PetscErrorCode PetscDSGetSpatialDimension(PetscDS prob, PetscInt *dim)
582: {

588:   *dim = 0;
589:   if (prob->Nf) {
590:     PetscObject  obj;
591:     PetscClassId id;

593:     PetscDSGetDiscretization(prob, 0, &obj);
594:     PetscObjectGetClassId(obj, &id);
595:     if (id == PETSCFE_CLASSID)      {PetscFEGetSpatialDimension((PetscFE) obj, dim);}
596:     else if (id == PETSCFV_CLASSID) {PetscFVGetSpatialDimension((PetscFV) obj, dim);}
597:     else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", 0);
598:   }
599:   return(0);
600: }

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

605:   Not collective

607:   Input Parameter:
608: . prob - The PetscDS object

610:   Output Parameter:
611: . dimEmbed - The coordinate dimension

613:   Level: beginner

615: .seealso: PetscDSSetCoordinateDimension(), PetscDSGetSpatialDimension(), PetscDSGetNumFields(), PetscDSCreate()
616: @*/
617: PetscErrorCode PetscDSGetCoordinateDimension(PetscDS prob, PetscInt *dimEmbed)
618: {
622:   if (prob->dimEmbed < 0) SETERRQ(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONGSTATE, "No coordinate dimension set for this DS");
623:   *dimEmbed = prob->dimEmbed;
624:   return(0);
625: }

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

630:   Not collective

632:   Input Parameters:
633: + prob - The PetscDS object
634: - dimEmbed - The coordinate dimension

636:   Level: beginner

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

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

651:   Not collective

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

656:   Output Parameter:
657: . dim - The total problem dimension

659:   Level: beginner

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

669:   PetscDSSetUp(prob);
671:   *dim = prob->totDim;
672:   return(0);
673: }

675: /*@
676:   PetscDSGetTotalComponents - Returns the total number of components in this system

678:   Not collective

680:   Input Parameter:
681: . prob - The PetscDS object

683:   Output Parameter:
684: . dim - The total number of components

686:   Level: beginner

688: .seealso: PetscDSGetNumFields(), PetscDSCreate()
689: @*/
690: PetscErrorCode PetscDSGetTotalComponents(PetscDS prob, PetscInt *Nc)
691: {

696:   PetscDSSetUp(prob);
698:   *Nc = prob->totComp;
699:   return(0);
700: }

702: /*@
703:   PetscDSGetDiscretization - Returns the discretization object for the given field

705:   Not collective

707:   Input Parameters:
708: + prob - The PetscDS object
709: - f - The field number

711:   Output Parameter:
712: . disc - The discretization object

714:   Level: beginner

716: .seealso: PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
717: @*/
718: PetscErrorCode PetscDSGetDiscretization(PetscDS prob, PetscInt f, PetscObject *disc)
719: {
723:   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);
724:   *disc = prob->disc[f];
725:   return(0);
726: }

728: /*@
729:   PetscDSSetDiscretization - Sets the discretization object for the given field

731:   Not collective

733:   Input Parameters:
734: + prob - The PetscDS object
735: . f - The field number
736: - disc - The discretization object

738:   Level: beginner

740: .seealso: PetscDSGetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
741: @*/
742: PetscErrorCode PetscDSSetDiscretization(PetscDS prob, PetscInt f, PetscObject disc)
743: {

749:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
750:   PetscDSEnlarge_Static(prob, f+1);
751:   if (prob->disc[f]) {PetscObjectDereference(prob->disc[f]);}
752:   prob->disc[f] = disc;
753:   PetscObjectReference(disc);
754:   {
755:     PetscClassId id;

757:     PetscObjectGetClassId(disc, &id);
758:     if (id == PETSCFE_CLASSID) {
759:       PetscDSSetImplicit(prob, f, PETSC_TRUE);
760:       PetscDSSetAdjacency(prob, f, PETSC_FALSE, PETSC_TRUE);
761:     } else if (id == PETSCFV_CLASSID) {
762:       PetscDSSetImplicit(prob, f, PETSC_FALSE);
763:       PetscDSSetAdjacency(prob, f, PETSC_TRUE, PETSC_FALSE);
764:     }
765:   }
766:   return(0);
767: }

769: /*@
770:   PetscDSAddDiscretization - Adds a discretization object

772:   Not collective

774:   Input Parameters:
775: + prob - The PetscDS object
776: - disc - The boundary discretization object

778:   Level: beginner

780: .seealso: PetscDSGetDiscretization(), PetscDSSetDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
781: @*/
782: PetscErrorCode PetscDSAddDiscretization(PetscDS prob, PetscObject disc)
783: {

787:   PetscDSSetDiscretization(prob, prob->Nf, disc);
788:   return(0);
789: }

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

794:   Not collective

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

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

803:   Level: developer

805: .seealso: PetscDSSetImplicit(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
806: @*/
807: PetscErrorCode PetscDSGetImplicit(PetscDS prob, PetscInt f, PetscBool *implicit)
808: {
812:   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);
813:   *implicit = prob->implicit[f];
814:   return(0);
815: }

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

820:   Not collective

822:   Input Parameters:
823: + prob - The PetscDS object
824: . f - The field number
825: - implicit - The flag indicating what kind of solve to use for this field

827:   Level: developer

829: .seealso: PetscDSGetImplicit(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
830: @*/
831: PetscErrorCode PetscDSSetImplicit(PetscDS prob, PetscInt f, PetscBool implicit)
832: {
835:   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);
836:   prob->implicit[f] = implicit;
837:   return(0);
838: }

840: /*@
841:   PetscDSGetAdjacency - Returns the flags for determining variable influence

843:   Not collective

845:   Input Parameters:
846: + prob - The PetscDS object
847: - f - The field number

849:   Output Parameter:
850: + useCone    - Flag for variable influence starting with the cone operation
851: - useClosure - Flag for variable influence using transitive closure

853:   Note: See the discussion in DMPlexGetAdjacencyUseCone() and DMPlexGetAdjacencyUseClosure()

855:   Level: developer

857: .seealso: PetscDSSetAdjacency(), DMPlexGetAdjacencyUseCone(), DMPlexGetAdjacencyUseClosure(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
858: @*/
859: PetscErrorCode PetscDSGetAdjacency(PetscDS prob, PetscInt f, PetscBool *useCone, PetscBool *useClosure)
860: {
865:   if (f < 0) {
866:     if (useCone)    *useCone    = prob->defaultAdj[0];
867:     if (useClosure) *useClosure = prob->defaultAdj[1];
868:   } else {
869:     if (f >= prob->Nf) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
870:     if (useCone)    *useCone    = prob->adjacency[f*2+0];
871:     if (useClosure) *useClosure = prob->adjacency[f*2+1];
872:   }
873:   return(0);
874: }

876: /*@
877:   PetscDSSetAdjacency - Set the flags for determining variable influence

879:   Not collective

881:   Input Parameters:
882: + prob - The PetscDS object
883: . f - The field number
884: . useCone    - Flag for variable influence starting with the cone operation
885: - useClosure - Flag for variable influence using transitive closure

887:   Note: See the discussion in DMPlexGetAdjacencyUseCone() and DMPlexGetAdjacencyUseClosure()

889:   Level: developer

891: .seealso: PetscDSGetAdjacency(), DMPlexGetAdjacencyUseCone(), DMPlexGetAdjacencyUseClosure(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
892: @*/
893: PetscErrorCode PetscDSSetAdjacency(PetscDS prob, PetscInt f, PetscBool useCone, PetscBool useClosure)
894: {
897:   if (f < 0) {
898:     prob->defaultAdj[0] = useCone;
899:     prob->defaultAdj[1] = useClosure;
900:   } else {
901:     if (f >= prob->Nf) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
902:     prob->adjacency[f*2+0] = useCone;
903:     prob->adjacency[f*2+1] = useClosure;
904:   }
905:   return(0);
906: }

908: PetscErrorCode PetscDSGetObjective(PetscDS prob, PetscInt f,
909:                                    void (**obj)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
910:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
911:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
912:                                                 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[]))
913: {
917:   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);
918:   *obj = prob->obj[f];
919:   return(0);
920: }

922: PetscErrorCode PetscDSSetObjective(PetscDS prob, PetscInt f,
923:                                    void (*obj)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
924:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
925:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
926:                                                PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[]))
927: {

933:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
934:   PetscDSEnlarge_Static(prob, f+1);
935:   prob->obj[f] = obj;
936:   return(0);
937: }

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

942:   Not collective

944:   Input Parameters:
945: + prob - The PetscDS
946: - f    - The test field number

948:   Output Parameters:
949: + f0 - integrand for the test function term
950: - f1 - integrand for the test function gradient term

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

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

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

958: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
959: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
960: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
961: $    PetscReal t, const PetscReal x[], PetscScalar f0[])

963: + dim - the spatial dimension
964: . Nf - the number of fields
965: . uOff - the offset into u[] and u_t[] for each field
966: . uOff_x - the offset into u_x[] for each field
967: . u - each field evaluated at the current point
968: . u_t - the time derivative of each field evaluated at the current point
969: . u_x - the gradient of each field evaluated at the current point
970: . aOff - the offset into a[] and a_t[] for each auxiliary field
971: . aOff_x - the offset into a_x[] for each auxiliary field
972: . a - each auxiliary field evaluated at the current point
973: . a_t - the time derivative of each auxiliary field evaluated at the current point
974: . a_x - the gradient of auxiliary each field evaluated at the current point
975: . t - current time
976: . x - coordinates of the current point
977: . numConstants - number of constant parameters
978: . constants - constant parameters
979: - f0 - output values at the current point

981:   Level: intermediate

983: .seealso: PetscDSSetResidual()
984: @*/
985: PetscErrorCode PetscDSGetResidual(PetscDS prob, PetscInt f,
986:                                   void (**f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
987:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
988:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
989:                                               PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]),
990:                                   void (**f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
991:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
992:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
993:                                               PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
994: {
997:   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);
1000:   return(0);
1001: }

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

1006:   Not collective

1008:   Input Parameters:
1009: + prob - The PetscDS
1010: . f    - The test field number
1011: . f0 - integrand for the test function term
1012: - f1 - integrand for the test function gradient term

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

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

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

1020: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1021: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1022: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1023: $    PetscReal t, const PetscReal x[], PetscScalar f0[])

1025: + dim - the spatial dimension
1026: . Nf - the number of fields
1027: . uOff - the offset into u[] and u_t[] for each field
1028: . uOff_x - the offset into u_x[] for each field
1029: . u - each field evaluated at the current point
1030: . u_t - the time derivative of each field evaluated at the current point
1031: . u_x - the gradient of each field evaluated at the current point
1032: . aOff - the offset into a[] and a_t[] for each auxiliary field
1033: . aOff_x - the offset into a_x[] for each auxiliary field
1034: . a - each auxiliary field evaluated at the current point
1035: . a_t - the time derivative of each auxiliary field evaluated at the current point
1036: . a_x - the gradient of auxiliary each field evaluated at the current point
1037: . t - current time
1038: . x - coordinates of the current point
1039: . numConstants - number of constant parameters
1040: . constants - constant parameters
1041: - f0 - output values at the current point

1043:   Level: intermediate

1045: .seealso: PetscDSGetResidual()
1046: @*/
1047: PetscErrorCode PetscDSSetResidual(PetscDS prob, PetscInt f,
1048:                                   void (*f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1049:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1050:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1051:                                              PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]),
1052:                                   void (*f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1053:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1054:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1055:                                              PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
1056: {

1063:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1064:   PetscDSEnlarge_Static(prob, f+1);
1065:   prob->f[f*2+0] = f0;
1066:   prob->f[f*2+1] = f1;
1067:   return(0);
1068: }

1070: /*@C
1071:   PetscDSHasJacobian - Signals that Jacobian functions have been set

1073:   Not collective

1075:   Input Parameter:
1076: . prob - The PetscDS

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

1081:   Level: intermediate

1083: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1084: @*/
1085: PetscErrorCode PetscDSHasJacobian(PetscDS prob, PetscBool *hasJac)
1086: {
1087:   PetscInt f, g, h;

1091:   *hasJac = PETSC_FALSE;
1092:   for (f = 0; f < prob->Nf; ++f) {
1093:     for (g = 0; g < prob->Nf; ++g) {
1094:       for (h = 0; h < 4; ++h) {
1095:         if (prob->g[(f*prob->Nf + g)*4+h]) *hasJac = PETSC_TRUE;
1096:       }
1097:     }
1098:   }
1099:   return(0);
1100: }

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

1105:   Not collective

1107:   Input Parameters:
1108: + prob - The PetscDS
1109: . f    - The test field number
1110: - g    - The field number

1112:   Output Parameters:
1113: + g0 - integrand for the test and basis function term
1114: . g1 - integrand for the test function and basis function gradient term
1115: . g2 - integrand for the test function gradient and basis function term
1116: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

1124: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1125: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1126: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1127: $    PetscReal t, const PetscReal u_tShift, const PetscReal x[], PetscScalar g0[])

1129: + dim - the spatial dimension
1130: . Nf - the number of fields
1131: . uOff - the offset into u[] and u_t[] for each field
1132: . uOff_x - the offset into u_x[] for each field
1133: . u - each field evaluated at the current point
1134: . u_t - the time derivative of each field evaluated at the current point
1135: . u_x - the gradient of each field evaluated at the current point
1136: . aOff - the offset into a[] and a_t[] for each auxiliary field
1137: . aOff_x - the offset into a_x[] for each auxiliary field
1138: . a - each auxiliary field evaluated at the current point
1139: . a_t - the time derivative of each auxiliary field evaluated at the current point
1140: . a_x - the gradient of auxiliary each field evaluated at the current point
1141: . t - current time
1142: . u_tShift - the multiplier a for dF/dU_t
1143: . x - coordinates of the current point
1144: . numConstants - number of constant parameters
1145: . constants - constant parameters
1146: - g0 - output values at the current point

1148:   Level: intermediate

1150: .seealso: PetscDSSetJacobian()
1151: @*/
1152: PetscErrorCode PetscDSGetJacobian(PetscDS prob, PetscInt f, PetscInt g,
1153:                                   void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1154:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1155:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1156:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1157:                                   void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1158:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1159:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1160:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1161:                                   void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1162:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1163:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1164:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1165:                                   void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1166:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1167:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1168:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1169: {
1172:   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);
1173:   if ((g < 0) || (g >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, prob->Nf);
1178:   return(0);
1179: }

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

1184:   Not collective

1186:   Input Parameters:
1187: + prob - The PetscDS
1188: . f    - The test field number
1189: . g    - The field number
1190: . g0 - integrand for the test and basis function term
1191: . g1 - integrand for the test function and basis function gradient term
1192: . g2 - integrand for the test function gradient and basis function term
1193: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

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

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

1225:   Level: intermediate

1227: .seealso: PetscDSGetJacobian()
1228: @*/
1229: PetscErrorCode PetscDSSetJacobian(PetscDS prob, PetscInt f, PetscInt g,
1230:                                   void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1231:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1232:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1233:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1234:                                   void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1235:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1236:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1237:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1238:                                   void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1239:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1240:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1241:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1242:                                   void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1243:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1244:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1245:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1246: {

1255:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1256:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
1257:   PetscDSEnlarge_Static(prob, PetscMax(f, g)+1);
1258:   prob->g[(f*prob->Nf + g)*4+0] = g0;
1259:   prob->g[(f*prob->Nf + g)*4+1] = g1;
1260:   prob->g[(f*prob->Nf + g)*4+2] = g2;
1261:   prob->g[(f*prob->Nf + g)*4+3] = g3;
1262:   return(0);
1263: }

1265: /*@C
1266:   PetscDSUseJacobianPreconditioner - Whether to construct a Jacobian preconditioner

1268:   Not collective

1270:   Input Parameters:
1271: + prob - The PetscDS
1272: - useJacPre - flag that enables construction of a Jacobian preconditioner

1274:   Level: intermediate

1276: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1277: @*/
1278: PetscErrorCode PetscDSUseJacobianPreconditioner(PetscDS prob, PetscBool useJacPre)
1279: {
1282:   prob->useJacPre = useJacPre;
1283:   return(0);
1284: }

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

1289:   Not collective

1291:   Input Parameter:
1292: . prob - The PetscDS

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

1297:   Level: intermediate

1299: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1300: @*/
1301: PetscErrorCode PetscDSHasJacobianPreconditioner(PetscDS prob, PetscBool *hasJacPre)
1302: {
1303:   PetscInt f, g, h;

1307:   *hasJacPre = PETSC_FALSE;
1308:   if (!prob->useJacPre) return(0);
1309:   for (f = 0; f < prob->Nf; ++f) {
1310:     for (g = 0; g < prob->Nf; ++g) {
1311:       for (h = 0; h < 4; ++h) {
1312:         if (prob->gp[(f*prob->Nf + g)*4+h]) *hasJacPre = PETSC_TRUE;
1313:       }
1314:     }
1315:   }
1316:   return(0);
1317: }

1319: /*@C
1320:   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.

1322:   Not collective

1324:   Input Parameters:
1325: + prob - The PetscDS
1326: . f    - The test field number
1327: - g    - The field number

1329:   Output Parameters:
1330: + g0 - integrand for the test and basis function term
1331: . g1 - integrand for the test function and basis function gradient term
1332: . g2 - integrand for the test function gradient and basis function term
1333: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

1341: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1342: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1343: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1344: $    PetscReal t, const PetscReal u_tShift, const PetscReal x[], PetscScalar g0[])

1346: + dim - the spatial dimension
1347: . Nf - the number of fields
1348: . uOff - the offset into u[] and u_t[] for each field
1349: . uOff_x - the offset into u_x[] for each field
1350: . u - each field evaluated at the current point
1351: . u_t - the time derivative of each field evaluated at the current point
1352: . u_x - the gradient of each field evaluated at the current point
1353: . aOff - the offset into a[] and a_t[] for each auxiliary field
1354: . aOff_x - the offset into a_x[] for each auxiliary field
1355: . a - each auxiliary field evaluated at the current point
1356: . a_t - the time derivative of each auxiliary field evaluated at the current point
1357: . a_x - the gradient of auxiliary each field evaluated at the current point
1358: . t - current time
1359: . u_tShift - the multiplier a for dF/dU_t
1360: . x - coordinates of the current point
1361: . numConstants - number of constant parameters
1362: . constants - constant parameters
1363: - g0 - output values at the current point

1365:   Level: intermediate

1367: .seealso: PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1368: @*/
1369: PetscErrorCode PetscDSGetJacobianPreconditioner(PetscDS prob, PetscInt f, PetscInt g,
1370:                                   void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1371:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1372:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1373:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1374:                                   void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1375:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1376:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1377:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1378:                                   void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1379:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1380:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1381:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1382:                                   void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1383:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1384:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1385:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1386: {
1389:   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);
1390:   if ((g < 0) || (g >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, prob->Nf);
1395:   return(0);
1396: }

1398: /*@C
1399:   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.

1401:   Not collective

1403:   Input Parameters:
1404: + prob - The PetscDS
1405: . f    - The test field number
1406: . g    - The field number
1407: . g0 - integrand for the test and basis function term
1408: . g1 - integrand for the test function and basis function gradient term
1409: . g2 - integrand for the test function gradient and basis function term
1410: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

1418: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1419: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1420: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1421: $    PetscReal t, const PetscReal x[], PetscScalar g0[])

1423: + dim - the spatial dimension
1424: . Nf - the number of fields
1425: . uOff - the offset into u[] and u_t[] for each field
1426: . uOff_x - the offset into u_x[] for each field
1427: . u - each field evaluated at the current point
1428: . u_t - the time derivative of each field evaluated at the current point
1429: . u_x - the gradient of each field evaluated at the current point
1430: . aOff - the offset into a[] and a_t[] for each auxiliary field
1431: . aOff_x - the offset into a_x[] for each auxiliary field
1432: . a - each auxiliary field evaluated at the current point
1433: . a_t - the time derivative of each auxiliary field evaluated at the current point
1434: . a_x - the gradient of auxiliary each field evaluated at the current point
1435: . t - current time
1436: . u_tShift - the multiplier a for dF/dU_t
1437: . x - coordinates of the current point
1438: . numConstants - number of constant parameters
1439: . constants - constant parameters
1440: - g0 - output values at the current point

1442:   Level: intermediate

1444: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobian()
1445: @*/
1446: PetscErrorCode PetscDSSetJacobianPreconditioner(PetscDS prob, PetscInt f, PetscInt g,
1447:                                   void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1448:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1449:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1450:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1451:                                   void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1452:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1453:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1454:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1455:                                   void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1456:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1457:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1458:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1459:                                   void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1460:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1461:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1462:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1463: {

1472:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1473:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
1474:   PetscDSEnlarge_Static(prob, PetscMax(f, g)+1);
1475:   prob->gp[(f*prob->Nf + g)*4+0] = g0;
1476:   prob->gp[(f*prob->Nf + g)*4+1] = g1;
1477:   prob->gp[(f*prob->Nf + g)*4+2] = g2;
1478:   prob->gp[(f*prob->Nf + g)*4+3] = g3;
1479:   return(0);
1480: }

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

1485:   Not collective

1487:   Input Parameter:
1488: . prob - The PetscDS

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

1493:   Level: intermediate

1495: .seealso: PetscDSGetDynamicJacobian(), PetscDSSetDynamicJacobian(), PetscDSGetJacobian()
1496: @*/
1497: PetscErrorCode PetscDSHasDynamicJacobian(PetscDS prob, PetscBool *hasDynJac)
1498: {
1499:   PetscInt f, g, h;

1503:   *hasDynJac = PETSC_FALSE;
1504:   for (f = 0; f < prob->Nf; ++f) {
1505:     for (g = 0; g < prob->Nf; ++g) {
1506:       for (h = 0; h < 4; ++h) {
1507:         if (prob->gt[(f*prob->Nf + g)*4+h]) *hasDynJac = PETSC_TRUE;
1508:       }
1509:     }
1510:   }
1511:   return(0);
1512: }

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

1517:   Not collective

1519:   Input Parameters:
1520: + prob - The PetscDS
1521: . f    - The test field number
1522: - g    - The field number

1524:   Output Parameters:
1525: + g0 - integrand for the test and basis function term
1526: . g1 - integrand for the test function and basis function gradient term
1527: . g2 - integrand for the test function gradient and basis function term
1528: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

1536: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1537: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1538: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1539: $    PetscReal t, const PetscReal u_tShift, const PetscReal x[], PetscScalar g0[])

1541: + dim - the spatial dimension
1542: . Nf - the number of fields
1543: . uOff - the offset into u[] and u_t[] for each field
1544: . uOff_x - the offset into u_x[] for each field
1545: . u - each field evaluated at the current point
1546: . u_t - the time derivative of each field evaluated at the current point
1547: . u_x - the gradient of each field evaluated at the current point
1548: . aOff - the offset into a[] and a_t[] for each auxiliary field
1549: . aOff_x - the offset into a_x[] for each auxiliary field
1550: . a - each auxiliary field evaluated at the current point
1551: . a_t - the time derivative of each auxiliary field evaluated at the current point
1552: . a_x - the gradient of auxiliary each field evaluated at the current point
1553: . t - current time
1554: . u_tShift - the multiplier a for dF/dU_t
1555: . x - coordinates of the current point
1556: . numConstants - number of constant parameters
1557: . constants - constant parameters
1558: - g0 - output values at the current point

1560:   Level: intermediate

1562: .seealso: PetscDSSetJacobian()
1563: @*/
1564: PetscErrorCode PetscDSGetDynamicJacobian(PetscDS prob, PetscInt f, PetscInt g,
1565:                                          void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1566:                                                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1567:                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1568:                                                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1569:                                          void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1570:                                                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1571:                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1572:                                                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1573:                                          void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1574:                                                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1575:                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1576:                                                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1577:                                          void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1578:                                                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1579:                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1580:                                                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1581: {
1584:   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);
1585:   if ((g < 0) || (g >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, prob->Nf);
1590:   return(0);
1591: }

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

1596:   Not collective

1598:   Input Parameters:
1599: + prob - The PetscDS
1600: . f    - The test field number
1601: . g    - The field number
1602: . g0 - integrand for the test and basis function term
1603: . g1 - integrand for the test function and basis function gradient term
1604: . g2 - integrand for the test function gradient and basis function term
1605: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

1613: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1614: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1615: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1616: $    PetscReal t, const PetscReal x[], PetscScalar g0[])

1618: + dim - the spatial dimension
1619: . Nf - the number of fields
1620: . uOff - the offset into u[] and u_t[] for each field
1621: . uOff_x - the offset into u_x[] for each field
1622: . u - each field evaluated at the current point
1623: . u_t - the time derivative of each field evaluated at the current point
1624: . u_x - the gradient of each field evaluated at the current point
1625: . aOff - the offset into a[] and a_t[] for each auxiliary field
1626: . aOff_x - the offset into a_x[] for each auxiliary field
1627: . a - each auxiliary field evaluated at the current point
1628: . a_t - the time derivative of each auxiliary field evaluated at the current point
1629: . a_x - the gradient of auxiliary each field evaluated at the current point
1630: . t - current time
1631: . u_tShift - the multiplier a for dF/dU_t
1632: . x - coordinates of the current point
1633: . numConstants - number of constant parameters
1634: . constants - constant parameters
1635: - g0 - output values at the current point

1637:   Level: intermediate

1639: .seealso: PetscDSGetJacobian()
1640: @*/
1641: PetscErrorCode PetscDSSetDynamicJacobian(PetscDS prob, PetscInt f, PetscInt g,
1642:                                          void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1643:                                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1644:                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1645:                                                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1646:                                          void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1647:                                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1648:                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1649:                                                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1650:                                          void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1651:                                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1652:                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1653:                                                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1654:                                          void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1655:                                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1656:                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1657:                                                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1658: {

1667:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1668:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
1669:   PetscDSEnlarge_Static(prob, PetscMax(f, g)+1);
1670:   prob->gt[(f*prob->Nf + g)*4+0] = g0;
1671:   prob->gt[(f*prob->Nf + g)*4+1] = g1;
1672:   prob->gt[(f*prob->Nf + g)*4+2] = g2;
1673:   prob->gt[(f*prob->Nf + g)*4+3] = g3;
1674:   return(0);
1675: }

1677: /*@C
1678:   PetscDSGetRiemannSolver - Returns the Riemann solver for the given field

1680:   Not collective

1682:   Input Arguments:
1683: + prob - The PetscDS object
1684: - f    - The field number

1686:   Output Argument:
1687: . r    - Riemann solver

1689:   Calling sequence for r:

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

1693: + dim  - The spatial dimension
1694: . Nf   - The number of fields
1695: . x    - The coordinates at a point on the interface
1696: . n    - The normal vector to the interface
1697: . uL   - The state vector to the left of the interface
1698: . uR   - The state vector to the right of the interface
1699: . flux - output array of flux through the interface
1700: . numConstants - number of constant parameters
1701: . constants - constant parameters
1702: - ctx  - optional user context

1704:   Level: intermediate

1706: .seealso: PetscDSSetRiemannSolver()
1707: @*/
1708: PetscErrorCode PetscDSGetRiemannSolver(PetscDS prob, PetscInt f,
1709:                                        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))
1710: {
1713:   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);
1715:   *r = prob->r[f];
1716:   return(0);
1717: }

1719: /*@C
1720:   PetscDSSetRiemannSolver - Sets the Riemann solver for the given field

1722:   Not collective

1724:   Input Arguments:
1725: + prob - The PetscDS object
1726: . f    - The field number
1727: - r    - Riemann solver

1729:   Calling sequence for r:

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

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

1744:   Level: intermediate

1746: .seealso: PetscDSGetRiemannSolver()
1747: @*/
1748: PetscErrorCode PetscDSSetRiemannSolver(PetscDS prob, PetscInt f,
1749:                                        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))
1750: {

1756:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1757:   PetscDSEnlarge_Static(prob, f+1);
1758:   prob->r[f] = r;
1759:   return(0);
1760: }

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

1765:   Not collective

1767:   Input Parameters:
1768: + prob - The PetscDS
1769: - f    - The field number

1771:   Output Parameters:
1772: + update - update function

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

1776: $ update(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1777: $        const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1778: $        const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1779: $        PetscReal t, const PetscReal x[], PetscScalar uNew[])

1781: + dim - the spatial dimension
1782: . Nf - the number of fields
1783: . uOff - the offset into u[] and u_t[] for each field
1784: . uOff_x - the offset into u_x[] for each field
1785: . u - each field evaluated at the current point
1786: . u_t - the time derivative of each field evaluated at the current point
1787: . u_x - the gradient of each field evaluated at the current point
1788: . aOff - the offset into a[] and a_t[] for each auxiliary field
1789: . aOff_x - the offset into a_x[] for each auxiliary field
1790: . a - each auxiliary field evaluated at the current point
1791: . a_t - the time derivative of each auxiliary field evaluated at the current point
1792: . a_x - the gradient of auxiliary each field evaluated at the current point
1793: . t - current time
1794: . x - coordinates of the current point
1795: - uNew - new value for field at the current point

1797:   Level: intermediate

1799: .seealso: PetscDSSetUpdate(), PetscDSSetResidual()
1800: @*/
1801: PetscErrorCode PetscDSGetUpdate(PetscDS prob, PetscInt f,
1802:                                   void (**update)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1803:                                                   const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1804:                                                   const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1805:                                                   PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar uNew[]))
1806: {
1809:   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);
1811:   return(0);
1812: }

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

1817:   Not collective

1819:   Input Parameters:
1820: + prob   - The PetscDS
1821: . f      - The field number
1822: - update - update function

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

1826: $ update(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1827: $        const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1828: $        const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1829: $        PetscReal t, const PetscReal x[], PetscScalar uNew[])

1831: + dim - the spatial dimension
1832: . Nf - the number of fields
1833: . uOff - the offset into u[] and u_t[] for each field
1834: . uOff_x - the offset into u_x[] for each field
1835: . u - each field evaluated at the current point
1836: . u_t - the time derivative of each field evaluated at the current point
1837: . u_x - the gradient of each field evaluated at the current point
1838: . aOff - the offset into a[] and a_t[] for each auxiliary field
1839: . aOff_x - the offset into a_x[] for each auxiliary field
1840: . a - each auxiliary field evaluated at the current point
1841: . a_t - the time derivative of each auxiliary field evaluated at the current point
1842: . a_x - the gradient of auxiliary each field evaluated at the current point
1843: . t - current time
1844: . x - coordinates of the current point
1845: - uNew - new field values at the current point

1847:   Level: intermediate

1849: .seealso: PetscDSGetResidual()
1850: @*/
1851: PetscErrorCode PetscDSSetUpdate(PetscDS prob, PetscInt f,
1852:                                 void (*update)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1853:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1854:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1855:                                                PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar uNew[]))
1856: {

1862:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1863:   PetscDSEnlarge_Static(prob, f+1);
1864:   prob->update[f] = update;
1865:   return(0);
1866: }

1868: PetscErrorCode PetscDSGetContext(PetscDS prob, PetscInt f, void **ctx)
1869: {
1872:   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);
1874:   *ctx = prob->ctx[f];
1875:   return(0);
1876: }

1878: PetscErrorCode PetscDSSetContext(PetscDS prob, PetscInt f, void *ctx)
1879: {

1884:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1885:   PetscDSEnlarge_Static(prob, f+1);
1886:   prob->ctx[f] = ctx;
1887:   return(0);
1888: }

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

1893:   Not collective

1895:   Input Parameters:
1896: + prob - The PetscDS
1897: - f    - The test field number

1899:   Output Parameters:
1900: + f0 - boundary integrand for the test function term
1901: - f1 - boundary integrand for the test function gradient term

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

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

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

1909: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1910: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1911: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1912: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[])

1914: + dim - the spatial dimension
1915: . Nf - the number of fields
1916: . uOff - the offset into u[] and u_t[] for each field
1917: . uOff_x - the offset into u_x[] for each field
1918: . u - each field evaluated at the current point
1919: . u_t - the time derivative of each field evaluated at the current point
1920: . u_x - the gradient of each field evaluated at the current point
1921: . aOff - the offset into a[] and a_t[] for each auxiliary field
1922: . aOff_x - the offset into a_x[] for each auxiliary field
1923: . a - each auxiliary field evaluated at the current point
1924: . a_t - the time derivative of each auxiliary field evaluated at the current point
1925: . a_x - the gradient of auxiliary each field evaluated at the current point
1926: . t - current time
1927: . x - coordinates of the current point
1928: . n - unit normal at the current point
1929: . numConstants - number of constant parameters
1930: . constants - constant parameters
1931: - f0 - output values at the current point

1933:   Level: intermediate

1935: .seealso: PetscDSSetBdResidual()
1936: @*/
1937: PetscErrorCode PetscDSGetBdResidual(PetscDS prob, PetscInt f,
1938:                                     void (**f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1939:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1940:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1941:                                                 PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]),
1942:                                     void (**f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1943:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1944:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1945:                                                 PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
1946: {
1949:   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);
1952:   return(0);
1953: }

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

1958:   Not collective

1960:   Input Parameters:
1961: + prob - The PetscDS
1962: . f    - The test field number
1963: . f0 - boundary integrand for the test function term
1964: - f1 - boundary integrand for the test function gradient term

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

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

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

1972: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1973: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1974: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1975: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[])

1977: + dim - the spatial dimension
1978: . Nf - the number of fields
1979: . uOff - the offset into u[] and u_t[] for each field
1980: . uOff_x - the offset into u_x[] for each field
1981: . u - each field evaluated at the current point
1982: . u_t - the time derivative of each field evaluated at the current point
1983: . u_x - the gradient of each field evaluated at the current point
1984: . aOff - the offset into a[] and a_t[] for each auxiliary field
1985: . aOff_x - the offset into a_x[] for each auxiliary field
1986: . a - each auxiliary field evaluated at the current point
1987: . a_t - the time derivative of each auxiliary field evaluated at the current point
1988: . a_x - the gradient of auxiliary each field evaluated at the current point
1989: . t - current time
1990: . x - coordinates of the current point
1991: . n - unit normal at the current point
1992: . numConstants - number of constant parameters
1993: . constants - constant parameters
1994: - f0 - output values at the current point

1996:   Level: intermediate

1998: .seealso: PetscDSGetBdResidual()
1999: @*/
2000: PetscErrorCode PetscDSSetBdResidual(PetscDS prob, PetscInt f,
2001:                                     void (*f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2002:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2003:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2004:                                                PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]),
2005:                                     void (*f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2006:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2007:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2008:                                                PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
2009: {

2014:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2015:   PetscDSEnlarge_Static(prob, f+1);
2018:   return(0);
2019: }

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

2024:   Not collective

2026:   Input Parameters:
2027: + prob - The PetscDS
2028: . f    - The test field number
2029: - g    - The field number

2031:   Output Parameters:
2032: + g0 - integrand for the test and basis function term
2033: . g1 - integrand for the test function and basis function gradient term
2034: . g2 - integrand for the test function gradient and basis function term
2035: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

2043: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2044: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2045: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2046: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar g0[])

2048: + dim - the spatial dimension
2049: . Nf - the number of fields
2050: . uOff - the offset into u[] and u_t[] for each field
2051: . uOff_x - the offset into u_x[] for each field
2052: . u - each field evaluated at the current point
2053: . u_t - the time derivative of each field evaluated at the current point
2054: . u_x - the gradient of each field evaluated at the current point
2055: . aOff - the offset into a[] and a_t[] for each auxiliary field
2056: . aOff_x - the offset into a_x[] for each auxiliary field
2057: . a - each auxiliary field evaluated at the current point
2058: . a_t - the time derivative of each auxiliary field evaluated at the current point
2059: . a_x - the gradient of auxiliary each field evaluated at the current point
2060: . t - current time
2061: . u_tShift - the multiplier a for dF/dU_t
2062: . x - coordinates of the current point
2063: . n - normal at the current point
2064: . numConstants - number of constant parameters
2065: . constants - constant parameters
2066: - g0 - output values at the current point

2068:   Level: intermediate

2070: .seealso: PetscDSSetBdJacobian()
2071: @*/
2072: PetscErrorCode PetscDSGetBdJacobian(PetscDS prob, PetscInt f, PetscInt g,
2073:                                     void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2074:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2075:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2076:                                                 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
2077:                                     void (**g1)(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, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
2081:                                     void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2082:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2083:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2084:                                                 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
2085:                                     void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2086:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2087:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2088:                                                 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
2089: {
2092:   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);
2093:   if ((g < 0) || (g >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, prob->Nf);
2098:   return(0);
2099: }

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

2104:   Not collective

2106:   Input Parameters:
2107: + prob - The PetscDS
2108: . f    - The test field number
2109: . g    - The field number
2110: . g0 - integrand for the test and basis function term
2111: . g1 - integrand for the test function and basis function gradient term
2112: . g2 - integrand for the test function gradient and basis function term
2113: - g3 - integrand for the test function gradient and basis function gradient term

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

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

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

2121: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2122: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2123: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2124: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar g0[])

2126: + dim - the spatial dimension
2127: . Nf - the number of fields
2128: . uOff - the offset into u[] and u_t[] for each field
2129: . uOff_x - the offset into u_x[] for each field
2130: . u - each field evaluated at the current point
2131: . u_t - the time derivative of each field evaluated at the current point
2132: . u_x - the gradient of each field evaluated at the current point
2133: . aOff - the offset into a[] and a_t[] for each auxiliary field
2134: . aOff_x - the offset into a_x[] for each auxiliary field
2135: . a - each auxiliary field evaluated at the current point
2136: . a_t - the time derivative of each auxiliary field evaluated at the current point
2137: . a_x - the gradient of auxiliary each field evaluated at the current point
2138: . t - current time
2139: . u_tShift - the multiplier a for dF/dU_t
2140: . x - coordinates of the current point
2141: . n - normal at the current point
2142: . numConstants - number of constant parameters
2143: . constants - constant parameters
2144: - g0 - output values at the current point

2146:   Level: intermediate

2148: .seealso: PetscDSGetBdJacobian()
2149: @*/
2150: PetscErrorCode PetscDSSetBdJacobian(PetscDS prob, PetscInt f, PetscInt g,
2151:                                     void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2152:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2153:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2154:                                                PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
2155:                                     void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2156:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2157:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2158:                                                PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
2159:                                     void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2160:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2161:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2162:                                                PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
2163:                                     void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2164:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2165:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2166:                                                PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
2167: {

2176:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2177:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
2178:   PetscDSEnlarge_Static(prob, PetscMax(f, g)+1);
2179:   prob->gBd[(f*prob->Nf + g)*4+0] = g0;
2180:   prob->gBd[(f*prob->Nf + g)*4+1] = g1;
2181:   prob->gBd[(f*prob->Nf + g)*4+2] = g2;
2182:   prob->gBd[(f*prob->Nf + g)*4+3] = g3;
2183:   return(0);
2184: }

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

2189:   Not collective

2191:   Input Parameters:
2192: + prob - The PetscDS
2193: - f    - The test field number

2195:   Output Parameter:
2196: . exactSol - exact solution for the test field

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

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

2202: + dim - the spatial dimension
2203: . t - current time
2204: . x - coordinates of the current point
2205: . Nc - the number of field components
2206: . u - the solution field evaluated at the current point
2207: - ctx - a user context

2209:   Level: intermediate

2211: .seealso: PetscDSSetExactSolution()
2212: @*/
2213: PetscErrorCode PetscDSGetExactSolution(PetscDS prob, PetscInt f, PetscErrorCode (**sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx))
2214: {
2217:   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
2219:   return(0);
2220: }

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

2225:   Not collective

2227:   Input Parameters:
2228: + prob - The PetscDS
2229: . f    - The test field number
2230: - sol  - solution function for the test fields

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

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

2236: + dim - the spatial dimension
2237: . t - current time
2238: . x - coordinates of the current point
2239: . Nc - the number of field components
2240: . u - the solution field evaluated at the current point
2241: - ctx - a user context

2243:   Level: intermediate

2245: .seealso: PetscDSGetExactSolution()
2246: @*/
2247: PetscErrorCode PetscDSSetExactSolution(PetscDS prob, PetscInt f, PetscErrorCode (*sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx))
2248: {

2253:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2254:   PetscDSEnlarge_Static(prob, f+1);
2256:   return(0);
2257: }

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

2262:   Not collective

2264:   Input Parameter:
2265: . prob - The PetscDS object

2267:   Output Parameters:
2268: + numConstants - The number of constants
2269: - constants    - The array of constants, NULL if there are none

2271:   Level: intermediate

2273: .seealso: PetscDSSetConstants(), PetscDSCreate()
2274: @*/
2275: PetscErrorCode PetscDSGetConstants(PetscDS prob, PetscInt *numConstants, const PetscScalar *constants[])
2276: {
2281:   return(0);
2282: }

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

2287:   Not collective

2289:   Input Parameters:
2290: + prob         - The PetscDS object
2291: . numConstants - The number of constants
2292: - constants    - The array of constants, NULL if there are none

2294:   Level: intermediate

2296: .seealso: PetscDSGetConstants(), PetscDSCreate()
2297: @*/
2298: PetscErrorCode PetscDSSetConstants(PetscDS prob, PetscInt numConstants, PetscScalar constants[])
2299: {

2304:   if (numConstants != prob->numConstants) {
2305:     PetscFree(prob->constants);
2306:     prob->numConstants = numConstants;
2307:     if (prob->numConstants) {
2308:       PetscMalloc1(prob->numConstants, &prob->constants);
2309:     } else {
2310:       prob->constants = NULL;
2311:     }
2312:   }
2313:   if (prob->numConstants) {
2315:     PetscMemcpy(prob->constants, constants, prob->numConstants * sizeof(PetscScalar));
2316:   }
2317:   return(0);
2318: }

2320: /*@
2321:   PetscDSGetFieldIndex - Returns the index of the given field

2323:   Not collective

2325:   Input Parameters:
2326: + prob - The PetscDS object
2327: - disc - The discretization object

2329:   Output Parameter:
2330: . f - The field number

2332:   Level: beginner

2334: .seealso: PetscGetDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
2335: @*/
2336: PetscErrorCode PetscDSGetFieldIndex(PetscDS prob, PetscObject disc, PetscInt *f)
2337: {
2338:   PetscInt g;

2343:   *f = -1;
2344:   for (g = 0; g < prob->Nf; ++g) {if (disc == prob->disc[g]) break;}
2345:   if (g == prob->Nf) SETERRQ(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Field not found in PetscDS.");
2346:   *f = g;
2347:   return(0);
2348: }

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

2353:   Not collective

2355:   Input Parameters:
2356: + prob - The PetscDS object
2357: - f - The field number

2359:   Output Parameter:
2360: . size - The size

2362:   Level: beginner

2364: .seealso: PetscDSGetFieldOffset(), PetscDSGetNumFields(), PetscDSCreate()
2365: @*/
2366: PetscErrorCode PetscDSGetFieldSize(PetscDS prob, PetscInt f, PetscInt *size)
2367: {

2373:   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);
2374:   PetscDSSetUp(prob);
2375:   *size = prob->Nb[f];
2376:   return(0);
2377: }

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

2382:   Not collective

2384:   Input Parameters:
2385: + prob - The PetscDS object
2386: - f - The field number

2388:   Output Parameter:
2389: . off - The offset

2391:   Level: beginner

2393: .seealso: PetscDSGetFieldSize(), PetscDSGetNumFields(), PetscDSCreate()
2394: @*/
2395: PetscErrorCode PetscDSGetFieldOffset(PetscDS prob, PetscInt f, PetscInt *off)
2396: {
2397:   PetscInt       size, g;

2403:   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);
2404:   *off = 0;
2405:   for (g = 0; g < f; ++g) {
2406:     PetscDSGetFieldSize(prob, g, &size);
2407:     *off += size;
2408:   }
2409:   return(0);
2410: }

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

2415:   Not collective

2417:   Input Parameter:
2418: . prob - The PetscDS object

2420:   Output Parameter:
2421: . dimensions - The number of dimensions

2423:   Level: beginner

2425: .seealso: PetscDSGetComponentOffsets(), PetscDSGetNumFields(), PetscDSCreate()
2426: @*/
2427: PetscErrorCode PetscDSGetDimensions(PetscDS prob, PetscInt *dimensions[])
2428: {

2433:   PetscDSSetUp(prob);
2435:   *dimensions = prob->Nb;
2436:   return(0);
2437: }

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

2442:   Not collective

2444:   Input Parameter:
2445: . prob - The PetscDS object

2447:   Output Parameter:
2448: . components - The number of components

2450:   Level: beginner

2452: .seealso: PetscDSGetComponentOffsets(), PetscDSGetNumFields(), PetscDSCreate()
2453: @*/
2454: PetscErrorCode PetscDSGetComponents(PetscDS prob, PetscInt *components[])
2455: {

2460:   PetscDSSetUp(prob);
2462:   *components = prob->Nc;
2463:   return(0);
2464: }

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

2469:   Not collective

2471:   Input Parameters:
2472: + prob - The PetscDS object
2473: - f - The field number

2475:   Output Parameter:
2476: . off - The offset

2478:   Level: beginner

2480: .seealso: PetscDSGetNumFields(), PetscDSCreate()
2481: @*/
2482: PetscErrorCode PetscDSGetComponentOffset(PetscDS prob, PetscInt f, PetscInt *off)
2483: {
2487:   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);
2488:   *off = prob->off[f];
2489:   return(0);
2490: }

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

2495:   Not collective

2497:   Input Parameter:
2498: . prob - The PetscDS object

2500:   Output Parameter:
2501: . offsets - The offsets

2503:   Level: beginner

2505: .seealso: PetscDSGetNumFields(), PetscDSCreate()
2506: @*/
2507: PetscErrorCode PetscDSGetComponentOffsets(PetscDS prob, PetscInt *offsets[])
2508: {
2512:   *offsets = prob->off;
2513:   return(0);
2514: }

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

2519:   Not collective

2521:   Input Parameter:
2522: . prob - The PetscDS object

2524:   Output Parameter:
2525: . offsets - The offsets

2527:   Level: beginner

2529: .seealso: PetscDSGetNumFields(), PetscDSCreate()
2530: @*/
2531: PetscErrorCode PetscDSGetComponentDerivativeOffsets(PetscDS prob, PetscInt *offsets[])
2532: {
2536:   *offsets = prob->offDer;
2537:   return(0);
2538: }

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

2543:   Not collective

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

2548:   Output Parameters:
2549: + basis - The basis function tabulation at quadrature points
2550: - basisDer - The basis function derivative tabulation at quadrature points

2552:   Level: intermediate

2554: .seealso: PetscDSCreate()
2555: @*/
2556: PetscErrorCode PetscDSGetTabulation(PetscDS prob, PetscReal ***basis, PetscReal ***basisDer)
2557: {

2562:   PetscDSSetUp(prob);
2565:   return(0);
2566: }

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

2571:   Not collective

2573:   Input Parameter:
2574: . prob - The PetscDS object

2576:   Output Parameters:
2577: + basisFace - The basis function tabulation at quadrature points
2578: - basisDerFace - The basis function derivative tabulation at quadrature points

2580:   Level: intermediate

2582: .seealso: PetscDSGetTabulation(), PetscDSCreate()
2583: @*/
2584: PetscErrorCode PetscDSGetFaceTabulation(PetscDS prob, PetscReal ***basis, PetscReal ***basisDer)
2585: {

2590:   PetscDSSetUp(prob);
2593:   return(0);
2594: }

2596: PetscErrorCode PetscDSGetEvaluationArrays(PetscDS prob, PetscScalar **u, PetscScalar **u_t, PetscScalar **u_x)
2597: {

2602:   PetscDSSetUp(prob);
2606:   return(0);
2607: }

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

2615:   PetscDSSetUp(prob);
2622:   return(0);
2623: }

2625: PetscErrorCode PetscDSGetRefCoordArrays(PetscDS prob, PetscReal **x, PetscScalar **refSpaceDer)
2626: {

2631:   PetscDSSetUp(prob);
2634:   return(0);
2635: }

2637: /*@C
2638:   PetscDSAddBoundary - Add a boundary condition to the model

2640:   Input Parameters:
2641: + ds          - The PetscDS object
2642: . type        - The type of condition, e.g. DM_BC_ESSENTIAL/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
2643: . name        - The BC name
2644: . labelname   - The label defining constrained points
2645: . field       - The field to constrain
2646: . numcomps    - The number of constrained field components
2647: . comps       - An array of constrained component numbers
2648: . bcFunc      - A pointwise function giving boundary values
2649: . numids      - The number of DMLabel ids for constrained points
2650: . ids         - An array of ids for constrained points
2651: - ctx         - An optional user context for bcFunc

2653:   Options Database Keys:
2654: + -bc_<boundary name> <num> - Overrides the boundary ids
2655: - -bc_<boundary name>_comp <num> - Overrides the boundary components

2657:   Level: developer

2659: .seealso: PetscDSGetBoundary()
2660: @*/
2661: PetscErrorCode PetscDSAddBoundary(PetscDS ds, DMBoundaryConditionType type, const char name[], const char labelname[], PetscInt field, PetscInt numcomps, const PetscInt *comps, void (*bcFunc)(void), PetscInt numids, const PetscInt *ids, void *ctx)
2662: {
2663:   DSBoundary     b;

2668:   PetscNew(&b);
2669:   PetscStrallocpy(name, (char **) &b->name);
2670:   PetscStrallocpy(labelname, (char **) &b->labelname);
2671:   PetscMalloc1(numcomps, &b->comps);
2672:   if (numcomps) {PetscMemcpy(b->comps, comps, numcomps*sizeof(PetscInt));}
2673:   PetscMalloc1(numids, &b->ids);
2674:   if (numids) {PetscMemcpy(b->ids, ids, numids*sizeof(PetscInt));}
2675:   b->type            = type;
2676:   b->field           = field;
2677:   b->numcomps        = numcomps;
2678:   b->func            = bcFunc;
2679:   b->numids          = numids;
2680:   b->ctx             = ctx;
2681:   b->next            = ds->boundary;
2682:   ds->boundary       = b;
2683:   return(0);
2684: }

2686: /*@
2687:   PetscDSGetNumBoundary - Get the number of registered BC

2689:   Input Parameters:
2690: . ds - The PetscDS object

2692:   Output Parameters:
2693: . numBd - The number of BC

2695:   Level: intermediate

2697: .seealso: PetscDSAddBoundary(), PetscDSGetBoundary()
2698: @*/
2699: PetscErrorCode PetscDSGetNumBoundary(PetscDS ds, PetscInt *numBd)
2700: {
2701:   DSBoundary b = ds->boundary;

2706:   *numBd = 0;
2707:   while (b) {++(*numBd); b = b->next;}
2708:   return(0);
2709: }

2711: /*@C
2712:   PetscDSGetBoundary - Add a boundary condition to the model

2714:   Input Parameters:
2715: + ds          - The PetscDS object
2716: - bd          - The BC number

2718:   Output Parameters:
2719: + type        - The type of condition, e.g. DM_BC_ESSENTIAL/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
2720: . name        - The BC name
2721: . labelname   - The label defining constrained points
2722: . field       - The field to constrain
2723: . numcomps    - The number of constrained field components
2724: . comps       - An array of constrained component numbers
2725: . bcFunc      - A pointwise function giving boundary values
2726: . numids      - The number of DMLabel ids for constrained points
2727: . ids         - An array of ids for constrained points
2728: - ctx         - An optional user context for bcFunc

2730:   Options Database Keys:
2731: + -bc_<boundary name> <num> - Overrides the boundary ids
2732: - -bc_<boundary name>_comp <num> - Overrides the boundary components

2734:   Level: developer

2736: .seealso: PetscDSAddBoundary()
2737: @*/
2738: PetscErrorCode PetscDSGetBoundary(PetscDS ds, PetscInt bd, DMBoundaryConditionType *type, const char **name, const char **labelname, PetscInt *field, PetscInt *numcomps, const PetscInt **comps, void (**func)(void), PetscInt *numids, const PetscInt **ids, void **ctx)
2739: {
2740:   DSBoundary b    = ds->boundary;
2741:   PetscInt   n    = 0;

2745:   while (b) {
2746:     if (n == bd) break;
2747:     b = b->next;
2748:     ++n;
2749:   }
2750:   if (!b) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %d is not in [0, %d)", bd, n);
2751:   if (type) {
2753:     *type = b->type;
2754:   }
2755:   if (name) {
2757:     *name = b->name;
2758:   }
2759:   if (labelname) {
2761:     *labelname = b->labelname;
2762:   }
2763:   if (field) {
2765:     *field = b->field;
2766:   }
2767:   if (numcomps) {
2769:     *numcomps = b->numcomps;
2770:   }
2771:   if (comps) {
2773:     *comps = b->comps;
2774:   }
2775:   if (func) {
2777:     *func = b->func;
2778:   }
2779:   if (numids) {
2781:     *numids = b->numids;
2782:   }
2783:   if (ids) {
2785:     *ids = b->ids;
2786:   }
2787:   if (ctx) {
2789:     *ctx = b->ctx;
2790:   }
2791:   return(0);
2792: }

2794: /*@
2795:   PetscDSCopyBoundary - Copy all boundary condition objects to the new problem

2797:   Not collective

2799:   Input Parameter:
2800: . prob - The PetscDS object

2802:   Output Parameter:
2803: . newprob - The PetscDS copy

2805:   Level: intermediate

2807: .seealso: PetscDSCopyEquations(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
2808: @*/
2809: PetscErrorCode PetscDSCopyBoundary(PetscDS probA, PetscDS probB)
2810: {
2811:   DSBoundary     b, next, *lastnext;

2817:   if (probA == probB) return(0);
2818:   next = probB->boundary;
2819:   while (next) {
2820:     DSBoundary b = next;

2822:     next = b->next;
2823:     PetscFree(b->comps);
2824:     PetscFree(b->ids);
2825:     PetscFree(b->name);
2826:     PetscFree(b->labelname);
2827:     PetscFree(b);
2828:   }
2829:   lastnext = &(probB->boundary);
2830:   for (b = probA->boundary; b; b = b->next) {
2831:     DSBoundary bNew;

2833:     PetscNew(&bNew);
2834:     bNew->numcomps = b->numcomps;
2835:     PetscMalloc1(bNew->numcomps, &bNew->comps);
2836:     PetscMemcpy(bNew->comps, b->comps, bNew->numcomps*sizeof(PetscInt));
2837:     bNew->numids = b->numids;
2838:     PetscMalloc1(bNew->numids, &bNew->ids);
2839:     PetscMemcpy(bNew->ids, b->ids, bNew->numids*sizeof(PetscInt));
2840:     PetscStrallocpy(b->labelname,(char **) &(bNew->labelname));
2841:     PetscStrallocpy(b->name,(char **) &(bNew->name));
2842:     bNew->ctx   = b->ctx;
2843:     bNew->type  = b->type;
2844:     bNew->field = b->field;
2845:     bNew->func  = b->func;

2847:     *lastnext = bNew;
2848:     lastnext = &(bNew->next);
2849:   }
2850:   return(0);
2851: }

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

2856:   Not collective

2858:   Input Parameter:
2859: + prob - The PetscDS object
2860: . numFields - Number of new fields
2861: - fields - Old field number for each new field

2863:   Output Parameter:
2864: . newprob - The PetscDS copy

2866:   Level: intermediate

2868: .seealso: PetscDSCopyBoundary(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
2869: @*/
2870: PetscErrorCode PetscDSSelectEquations(PetscDS prob, PetscInt numFields, const PetscInt fields[], PetscDS newprob)
2871: {
2872:   PetscInt       Nf, Nfn, fn, gn;

2879:   PetscDSGetNumFields(prob, &Nf);
2880:   PetscDSGetNumFields(newprob, &Nfn);
2881:   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);
2882:   for (fn = 0; fn < numFields; ++fn) {
2883:     const PetscInt   f = fields ? fields[fn] : fn;
2884:     PetscPointFunc   obj;
2885:     PetscPointFunc   f0, f1;
2886:     PetscBdPointFunc f0Bd, f1Bd;
2887:     PetscRiemannFunc r;

2889:     if (f >= Nf) SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_SIZ, "Field %D must be in [0, %D)", f, Nf);
2890:     PetscDSGetObjective(prob, f, &obj);
2891:     PetscDSGetResidual(prob, f, &f0, &f1);
2892:     PetscDSGetBdResidual(prob, f, &f0Bd, &f1Bd);
2893:     PetscDSGetRiemannSolver(prob, f, &r);
2894:     PetscDSSetObjective(newprob, fn, obj);
2895:     PetscDSSetResidual(newprob, fn, f0, f1);
2896:     PetscDSSetBdResidual(newprob, fn, f0Bd, f1Bd);
2897:     PetscDSSetRiemannSolver(newprob, fn, r);
2898:     for (gn = 0; gn < numFields; ++gn) {
2899:       const PetscInt  g = fields ? fields[gn] : gn;
2900:       PetscPointJac   g0, g1, g2, g3;
2901:       PetscPointJac   g0p, g1p, g2p, g3p;
2902:       PetscBdPointJac g0Bd, g1Bd, g2Bd, g3Bd;

2904:       if (g >= Nf) SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_SIZ, "Field %D must be in [0, %D)", g, Nf);
2905:       PetscDSGetJacobian(prob, f, g, &g0, &g1, &g2, &g3);
2906:       PetscDSGetJacobianPreconditioner(prob, f, g, &g0p, &g1p, &g2p, &g3p);
2907:       PetscDSGetBdJacobian(prob, f, g, &g0Bd, &g1Bd, &g2Bd, &g3Bd);
2908:       PetscDSSetJacobian(newprob, fn, gn, g0, g1, g2, g3);
2909:       PetscDSSetJacobianPreconditioner(prob, fn, gn, g0p, g1p, g2p, g3p);
2910:       PetscDSSetBdJacobian(newprob, fn, gn, g0Bd, g1Bd, g2Bd, g3Bd);
2911:     }
2912:   }
2913:   return(0);
2914: }

2916: /*@
2917:   PetscDSCopyEquations - Copy all pointwise function pointers to the new problem

2919:   Not collective

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

2924:   Output Parameter:
2925: . newprob - The PetscDS copy

2927:   Level: intermediate

2929: .seealso: PetscDSCopyBoundary(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
2930: @*/
2931: PetscErrorCode PetscDSCopyEquations(PetscDS prob, PetscDS newprob)
2932: {
2933:   PetscInt       Nf, Ng;

2939:   PetscDSGetNumFields(prob, &Nf);
2940:   PetscDSGetNumFields(newprob, &Ng);
2941:   if (Nf != Ng) SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_SIZ, "Number of fields must match %D != %D", Nf, Ng);
2942:   PetscDSSelectEquations(prob, Nf, NULL, newprob);
2943:   return(0);
2944: }
2945: /*@
2946:   PetscDSCopyConstants - Copy all constants to the new problem

2948:   Not collective

2950:   Input Parameter:
2951: . prob - The PetscDS object

2953:   Output Parameter:
2954: . newprob - The PetscDS copy

2956:   Level: intermediate

2958: .seealso: PetscDSCopyBoundary(), PetscDSCopyEquations(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
2959: @*/
2960: PetscErrorCode PetscDSCopyConstants(PetscDS prob, PetscDS newprob)
2961: {
2962:   PetscInt           Nc;
2963:   const PetscScalar *constants;
2964:   PetscErrorCode     ierr;

2969:   PetscDSGetConstants(prob, &Nc, &constants);
2970:   PetscDSSetConstants(newprob, Nc, (PetscScalar *) constants);
2971:   return(0);
2972: }

2974: PetscErrorCode PetscDSGetHeightSubspace(PetscDS prob, PetscInt height, PetscDS *subprob)
2975: {
2976:   PetscInt       dim, Nf, f;

2982:   if (height == 0) {*subprob = prob; return(0);}
2983:   PetscDSGetNumFields(prob, &Nf);
2984:   PetscDSGetSpatialDimension(prob, &dim);
2985:   if (height > dim) SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_OUTOFRANGE, "DS can only handle height in [0, %D], not %D", dim, height);
2986:   if (!prob->subprobs) {PetscCalloc1(dim, &prob->subprobs);}
2987:   if (!prob->subprobs[height-1]) {
2988:     PetscInt cdim;

2990:     PetscDSCreate(PetscObjectComm((PetscObject) prob), &prob->subprobs[height-1]);
2991:     PetscDSGetCoordinateDimension(prob, &cdim);
2992:     PetscDSSetCoordinateDimension(prob->subprobs[height-1], cdim);
2993:     for (f = 0; f < Nf; ++f) {
2994:       PetscFE      subfe;
2995:       PetscObject  obj;
2996:       PetscClassId id;

2998:       PetscDSGetDiscretization(prob, f, &obj);
2999:       PetscObjectGetClassId(obj, &id);
3000:       if (id == PETSCFE_CLASSID) {PetscFEGetHeightSubspace((PetscFE) obj, height, &subfe);}
3001:       else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unsupported discretization type for field %d", f);
3002:       PetscDSSetDiscretization(prob->subprobs[height-1], f, (PetscObject) subfe);
3003:     }
3004:   }
3005:   *subprob = prob->subprobs[height-1];
3006:   return(0);
3007: }

3009: static PetscErrorCode PetscDSDestroy_Basic(PetscDS prob)
3010: {
3011:   PetscErrorCode      ierr;

3014:   PetscFree(prob->data);
3015:   return(0);
3016: }

3018: static PetscErrorCode PetscDSInitialize_Basic(PetscDS prob)
3019: {
3021:   prob->ops->setfromoptions = NULL;
3022:   prob->ops->setup          = NULL;
3023:   prob->ops->view           = NULL;
3024:   prob->ops->destroy        = PetscDSDestroy_Basic;
3025:   return(0);
3026: }

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

3031:   Level: intermediate

3033: .seealso: PetscDSType, PetscDSCreate(), PetscDSSetType()
3034: M*/

3036: PETSC_EXTERN PetscErrorCode PetscDSCreate_Basic(PetscDS prob)
3037: {
3038:   PetscDS_Basic *b;
3039:   PetscErrorCode      ierr;

3043:   PetscNewLog(prob, &b);
3044:   prob->data = b;

3046:   PetscDSInitialize_Basic(prob);
3047:   return(0);
3048: }