Actual source code: dtds.c
1: #include <petsc/private/petscdsimpl.h>
3: PetscClassId PETSCDS_CLASSID = 0;
5: PetscFunctionList PetscDSList = NULL;
6: PetscBool PetscDSRegisterAllCalled = PETSC_FALSE;
8: /* A PetscDS (Discrete System) encodes a set of equations posed in a discrete space, which represents a set of
9: nonlinear continuum equations. The equations can have multiple fields, each field having a different
10: discretization. In addition, different pieces of the domain can have different field combinations and equations.
12: The DS provides the user a description of the approximation space on any given cell. It also gives pointwise
13: functions representing the equations.
15: Each field is associated with a label, marking the cells on which it is supported. Note that a field can be
16: supported on the closure of a cell not in the label due to overlap of the boundary of neighboring cells. The DM
17: then creates a DS for each set of cells with identical approximation spaces. When assembling, the user asks for
18: the space associated with a given cell. DMPlex uses the labels associated with each DS in the default integration loop.
19: */
21: /*@C
22: PetscDSRegister - Adds a new PetscDS implementation
24: Not Collective
26: Input Parameters:
27: + name - The name of a new user-defined creation routine
28: - create_func - The creation routine itself
30: Notes:
31: PetscDSRegister() may be called multiple times to add several user-defined PetscDSs
33: Sample usage:
34: .vb
35: PetscDSRegister("my_ds", MyPetscDSCreate);
36: .ve
38: Then, your PetscDS type can be chosen with the procedural interface via
39: .vb
40: PetscDSCreate(MPI_Comm, PetscDS *);
41: PetscDSSetType(PetscDS, "my_ds");
42: .ve
43: or at runtime via the option
44: .vb
45: -petscds_type my_ds
46: .ve
48: Level: advanced
50: Not available from Fortran
52: .seealso: PetscDSRegisterAll(), PetscDSRegisterDestroy()
54: @*/
55: PetscErrorCode PetscDSRegister(const char sname[], PetscErrorCode (*function)(PetscDS))
56: {
60: PetscFunctionListAdd(&PetscDSList, sname, function);
61: return(0);
62: }
64: /*@C
65: PetscDSSetType - Builds a particular PetscDS
67: Collective on prob
69: Input Parameters:
70: + prob - The PetscDS object
71: - name - The kind of system
73: Options Database Key:
74: . -petscds_type <type> - Sets the PetscDS type; use -help for a list of available types
76: Level: intermediate
78: Not available from Fortran
80: .seealso: PetscDSGetType(), PetscDSCreate()
81: @*/
82: PetscErrorCode PetscDSSetType(PetscDS prob, PetscDSType name)
83: {
84: PetscErrorCode (*r)(PetscDS);
85: PetscBool match;
90: PetscObjectTypeCompare((PetscObject) prob, name, &match);
91: if (match) return(0);
93: PetscDSRegisterAll();
94: PetscFunctionListFind(PetscDSList, name, &r);
95: if (!r) SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscDS type: %s", name);
97: if (prob->ops->destroy) {
98: (*prob->ops->destroy)(prob);
99: prob->ops->destroy = NULL;
100: }
101: (*r)(prob);
102: PetscObjectChangeTypeName((PetscObject) prob, name);
103: return(0);
104: }
106: /*@C
107: PetscDSGetType - Gets the PetscDS type name (as a string) from the object.
109: Not Collective
111: Input Parameter:
112: . prob - The PetscDS
114: Output Parameter:
115: . name - The PetscDS type name
117: Level: intermediate
119: Not available from Fortran
121: .seealso: PetscDSSetType(), PetscDSCreate()
122: @*/
123: PetscErrorCode PetscDSGetType(PetscDS prob, PetscDSType *name)
124: {
130: PetscDSRegisterAll();
131: *name = ((PetscObject) prob)->type_name;
132: return(0);
133: }
135: static PetscErrorCode PetscDSView_Ascii(PetscDS prob, PetscViewer viewer)
136: {
137: PetscViewerFormat format;
138: const PetscScalar *constants;
139: PetscInt numConstants, f;
140: PetscErrorCode ierr;
143: PetscViewerGetFormat(viewer, &format);
144: PetscViewerASCIIPrintf(viewer, "Discrete System with %d fields\n", prob->Nf);
145: PetscViewerASCIIPushTab(viewer);
146: PetscViewerASCIIPrintf(viewer, " cell total dim %D total comp %D\n", prob->totDim, prob->totComp);
147: if (prob->isHybrid) {PetscViewerASCIIPrintf(viewer, " hybrid cell\n");}
148: for (f = 0; f < prob->Nf; ++f) {
149: DSBoundary b;
150: PetscObject obj;
151: PetscClassId id;
152: PetscQuadrature q;
153: const char *name;
154: PetscInt Nc, Nq, Nqc;
156: PetscDSGetDiscretization(prob, f, &obj);
157: PetscObjectGetClassId(obj, &id);
158: PetscObjectGetName(obj, &name);
159: PetscViewerASCIIPrintf(viewer, "Field %s", name ? name : "<unknown>");
160: PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
161: if (id == PETSCFE_CLASSID) {
162: PetscFEGetNumComponents((PetscFE) obj, &Nc);
163: PetscFEGetQuadrature((PetscFE) obj, &q);
164: PetscViewerASCIIPrintf(viewer, " FEM");
165: } else if (id == PETSCFV_CLASSID) {
166: PetscFVGetNumComponents((PetscFV) obj, &Nc);
167: PetscFVGetQuadrature((PetscFV) obj, &q);
168: PetscViewerASCIIPrintf(viewer, " FVM");
169: }
170: else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", f);
171: if (Nc > 1) {PetscViewerASCIIPrintf(viewer, " %D components", Nc);}
172: else {PetscViewerASCIIPrintf(viewer, " %D component ", Nc);}
173: if (prob->implicit[f]) {PetscViewerASCIIPrintf(viewer, " (implicit)");}
174: else {PetscViewerASCIIPrintf(viewer, " (explicit)");}
175: if (q) {
176: PetscQuadratureGetData(q, NULL, &Nqc, &Nq, NULL, NULL);
177: PetscViewerASCIIPrintf(viewer, " (Nq %D Nqc %D)", Nq, Nqc);
178: }
179: PetscViewerASCIIPrintf(viewer, " %D-jet", prob->jetDegree[f]);
180: PetscViewerASCIIPrintf(viewer, "\n");
181: PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
182: PetscViewerASCIIPushTab(viewer);
183: if (id == PETSCFE_CLASSID) {PetscFEView((PetscFE) obj, viewer);}
184: else if (id == PETSCFV_CLASSID) {PetscFVView((PetscFV) obj, viewer);}
185: PetscViewerASCIIPopTab(viewer);
187: for (b = prob->boundary; b; b = b->next) {
188: char *name;
189: PetscInt c, i;
191: if (b->field != f) continue;
192: PetscViewerASCIIPushTab(viewer);
193: PetscViewerASCIIPrintf(viewer, "Boundary %s (%s) %s\n", b->name, b->lname, DMBoundaryConditionTypes[b->type]);
194: if (!b->Nc) {
195: PetscViewerASCIIPrintf(viewer, " all components\n");
196: } else {
197: PetscViewerASCIIPrintf(viewer, " components: ");
198: PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
199: for (c = 0; c < b->Nc; ++c) {
200: if (c > 0) {PetscViewerASCIIPrintf(viewer, ", ");}
201: PetscViewerASCIIPrintf(viewer, "%D", b->comps[c]);
202: }
203: PetscViewerASCIIPrintf(viewer, "\n");
204: PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
205: }
206: PetscViewerASCIIPrintf(viewer, " values: ");
207: PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
208: for (i = 0; i < b->Nv; ++i) {
209: if (i > 0) {PetscViewerASCIIPrintf(viewer, ", ");}
210: PetscViewerASCIIPrintf(viewer, "%D", b->values[i]);
211: }
212: PetscViewerASCIIPrintf(viewer, "\n");
213: PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
214: if (b->func) {
215: PetscDLAddr(b->func, &name);
216: if (name) {PetscViewerASCIIPrintf(viewer, " func: %s\n", name);}
217: else {PetscViewerASCIIPrintf(viewer, " func: %p\n", b->func);}
218: PetscFree(name);
219: }
220: if (b->func_t) {
221: PetscDLAddr(b->func_t, &name);
222: if (name) {PetscViewerASCIIPrintf(viewer, " func_t: %s\n", name);}
223: else {PetscViewerASCIIPrintf(viewer, " func_t: %p\n", b->func_t);}
224: PetscFree(name);
225: }
226: PetscWeakFormView(b->wf, viewer);
227: PetscViewerASCIIPopTab(viewer);
228: }
229: }
230: PetscDSGetConstants(prob, &numConstants, &constants);
231: if (numConstants) {
232: PetscViewerASCIIPrintf(viewer, "%D constants\n", numConstants);
233: PetscViewerASCIIPushTab(viewer);
234: for (f = 0; f < numConstants; ++f) {PetscViewerASCIIPrintf(viewer, "%g\n", (double) PetscRealPart(constants[f]));}
235: PetscViewerASCIIPopTab(viewer);
236: }
237: PetscWeakFormView(prob->wf, viewer);
238: PetscViewerASCIIPopTab(viewer);
239: return(0);
240: }
242: /*@C
243: PetscDSViewFromOptions - View from Options
245: Collective on PetscDS
247: Input Parameters:
248: + A - the PetscDS object
249: . obj - Optional object
250: - name - command line option
252: Level: intermediate
253: .seealso: PetscDS, PetscDSView, PetscObjectViewFromOptions(), PetscDSCreate()
254: @*/
255: PetscErrorCode PetscDSViewFromOptions(PetscDS A,PetscObject obj,const char name[])
256: {
261: PetscObjectViewFromOptions((PetscObject)A,obj,name);
262: return(0);
263: }
265: /*@C
266: PetscDSView - Views a PetscDS
268: Collective on prob
270: Input Parameters:
271: + prob - the PetscDS object to view
272: - v - the viewer
274: Level: developer
276: .seealso PetscDSDestroy()
277: @*/
278: PetscErrorCode PetscDSView(PetscDS prob, PetscViewer v)
279: {
280: PetscBool iascii;
285: if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) prob), &v);}
287: PetscObjectTypeCompare((PetscObject) v, PETSCVIEWERASCII, &iascii);
288: if (iascii) {PetscDSView_Ascii(prob, v);}
289: if (prob->ops->view) {(*prob->ops->view)(prob, v);}
290: return(0);
291: }
293: /*@
294: PetscDSSetFromOptions - sets parameters in a PetscDS from the options database
296: Collective on prob
298: Input Parameter:
299: . prob - the PetscDS object to set options for
301: Options Database:
302: + -petscds_type <type> : Set the DS type
303: . -petscds_view <view opt> : View the DS
304: . -petscds_jac_pre : Turn formation of a separate Jacobian preconditioner on and off
305: . -bc_<name> <ids> : Specify a list of label ids for a boundary condition
306: - -bc_<name>_comp <comps> : Specify a list of field components to constrain for a boundary condition
308: Level: developer
310: .seealso PetscDSView()
311: @*/
312: PetscErrorCode PetscDSSetFromOptions(PetscDS prob)
313: {
314: DSBoundary b;
315: const char *defaultType;
316: char name[256];
317: PetscBool flg;
322: if (!((PetscObject) prob)->type_name) {
323: defaultType = PETSCDSBASIC;
324: } else {
325: defaultType = ((PetscObject) prob)->type_name;
326: }
327: PetscDSRegisterAll();
329: PetscObjectOptionsBegin((PetscObject) prob);
330: for (b = prob->boundary; b; b = b->next) {
331: char optname[1024];
332: PetscInt ids[1024], len = 1024;
333: PetscBool flg;
335: PetscSNPrintf(optname, sizeof(optname), "-bc_%s", b->name);
336: PetscMemzero(ids, sizeof(ids));
337: PetscOptionsIntArray(optname, "List of boundary IDs", "", ids, &len, &flg);
338: if (flg) {
339: b->Nv = len;
340: PetscFree(b->values);
341: PetscMalloc1(len, &b->values);
342: PetscArraycpy(b->values, ids, len);
343: PetscWeakFormRewriteKeys(b->wf, b->label, len, b->values);
344: }
345: len = 1024;
346: PetscSNPrintf(optname, sizeof(optname), "-bc_%s_comp", b->name);
347: PetscMemzero(ids, sizeof(ids));
348: PetscOptionsIntArray(optname, "List of boundary field components", "", ids, &len, &flg);
349: if (flg) {
350: b->Nc = len;
351: PetscFree(b->comps);
352: PetscMalloc1(len, &b->comps);
353: PetscArraycpy(b->comps, ids, len);
354: }
355: }
356: PetscOptionsFList("-petscds_type", "Discrete System", "PetscDSSetType", PetscDSList, defaultType, name, 256, &flg);
357: if (flg) {
358: PetscDSSetType(prob, name);
359: } else if (!((PetscObject) prob)->type_name) {
360: PetscDSSetType(prob, defaultType);
361: }
362: PetscOptionsBool("-petscds_jac_pre", "Discrete System", "PetscDSUseJacobianPreconditioner", prob->useJacPre, &prob->useJacPre, &flg);
363: if (prob->ops->setfromoptions) {(*prob->ops->setfromoptions)(prob);}
364: /* process any options handlers added with PetscObjectAddOptionsHandler() */
365: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) prob);
366: PetscOptionsEnd();
367: if (prob->Nf) {PetscDSViewFromOptions(prob, NULL, "-petscds_view");}
368: return(0);
369: }
371: /*@C
372: PetscDSSetUp - Construct data structures for the PetscDS
374: Collective on prob
376: Input Parameter:
377: . prob - the PetscDS object to setup
379: Level: developer
381: .seealso PetscDSView(), PetscDSDestroy()
382: @*/
383: PetscErrorCode PetscDSSetUp(PetscDS prob)
384: {
385: const PetscInt Nf = prob->Nf;
386: PetscBool hasH = PETSC_FALSE;
387: PetscInt dim, dimEmbed, NbMax = 0, NcMax = 0, NqMax = 0, NsMax = 1, f;
392: if (prob->setup) return(0);
393: /* Calculate sizes */
394: PetscDSGetSpatialDimension(prob, &dim);
395: PetscDSGetCoordinateDimension(prob, &dimEmbed);
396: prob->totDim = prob->totComp = 0;
397: PetscMalloc2(Nf,&prob->Nc,Nf,&prob->Nb);
398: PetscCalloc2(Nf+1,&prob->off,Nf+1,&prob->offDer);
399: PetscMalloc2(Nf,&prob->T,Nf,&prob->Tf);
400: for (f = 0; f < Nf; ++f) {
401: PetscObject obj;
402: PetscClassId id;
403: PetscQuadrature q = NULL;
404: PetscInt Nq = 0, Nb, Nc;
406: PetscDSGetDiscretization(prob, f, &obj);
407: if (prob->jetDegree[f] > 1) hasH = PETSC_TRUE;
408: if (!obj) {
409: /* Empty mesh */
410: Nb = Nc = 0;
411: prob->T[f] = prob->Tf[f] = NULL;
412: } else {
413: PetscObjectGetClassId(obj, &id);
414: if (id == PETSCFE_CLASSID) {
415: PetscFE fe = (PetscFE) obj;
417: PetscFEGetQuadrature(fe, &q);
418: PetscFEGetDimension(fe, &Nb);
419: PetscFEGetNumComponents(fe, &Nc);
420: PetscFEGetCellTabulation(fe, prob->jetDegree[f], &prob->T[f]);
421: PetscFEGetFaceTabulation(fe, prob->jetDegree[f], &prob->Tf[f]);
422: } else if (id == PETSCFV_CLASSID) {
423: PetscFV fv = (PetscFV) obj;
425: PetscFVGetQuadrature(fv, &q);
426: PetscFVGetNumComponents(fv, &Nc);
427: Nb = Nc;
428: PetscFVGetCellTabulation(fv, &prob->T[f]);
429: /* TODO: should PetscFV also have face tabulation? Otherwise there will be a null pointer in prob->basisFace */
430: } else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
431: }
432: prob->Nc[f] = Nc;
433: prob->Nb[f] = Nb;
434: prob->off[f+1] = Nc + prob->off[f];
435: prob->offDer[f+1] = Nc*dim + prob->offDer[f];
436: if (q) {PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);}
437: NqMax = PetscMax(NqMax, Nq);
438: NbMax = PetscMax(NbMax, Nb);
439: NcMax = PetscMax(NcMax, Nc);
440: prob->totDim += Nb;
441: prob->totComp += Nc;
442: /* There are two faces for all fields but the cohesive field on a hybrid cell */
443: if (prob->isHybrid && (f < Nf-1)) prob->totDim += Nb;
444: }
445: /* Allocate works space */
446: NsMax = 2; /* Even non-hybrid discretizations can be used in a hybrid integration, so we need this extra workspace */
447: PetscMalloc3(NsMax*prob->totComp,&prob->u,NsMax*prob->totComp,&prob->u_t,NsMax*prob->totComp*dimEmbed + (hasH ? NsMax*prob->totComp*dimEmbed*dimEmbed : 0),&prob->u_x);
448: PetscMalloc5(dimEmbed,&prob->x,NbMax*NcMax,&prob->basisReal,NbMax*NcMax*dimEmbed,&prob->basisDerReal,NbMax*NcMax,&prob->testReal,NbMax*NcMax*dimEmbed,&prob->testDerReal);
449: PetscMalloc6(NsMax*NqMax*NcMax,&prob->f0,NsMax*NqMax*NcMax*dimEmbed,&prob->f1,
450: NsMax*NsMax*NqMax*NcMax*NcMax,&prob->g0,NsMax*NsMax*NqMax*NcMax*NcMax*dimEmbed,&prob->g1,
451: NsMax*NsMax*NqMax*NcMax*NcMax*dimEmbed,&prob->g2,NsMax*NsMax*NqMax*NcMax*NcMax*dimEmbed*dimEmbed,&prob->g3);
452: if (prob->ops->setup) {(*prob->ops->setup)(prob);}
453: prob->setup = PETSC_TRUE;
454: return(0);
455: }
457: static PetscErrorCode PetscDSDestroyStructs_Static(PetscDS prob)
458: {
462: PetscFree2(prob->Nc,prob->Nb);
463: PetscFree2(prob->off,prob->offDer);
464: PetscFree2(prob->T,prob->Tf);
465: PetscFree3(prob->u,prob->u_t,prob->u_x);
466: PetscFree5(prob->x,prob->basisReal, prob->basisDerReal,prob->testReal,prob->testDerReal);
467: PetscFree6(prob->f0,prob->f1,prob->g0,prob->g1,prob->g2,prob->g3);
468: return(0);
469: }
471: static PetscErrorCode PetscDSEnlarge_Static(PetscDS prob, PetscInt NfNew)
472: {
473: PetscObject *tmpd;
474: PetscBool *tmpi;
475: PetscInt *tmpk;
476: PetscPointFunc *tmpup;
477: PetscSimplePointFunc *tmpexactSol, *tmpexactSol_t;
478: void **tmpexactCtx, **tmpexactCtx_t;
479: void **tmpctx;
480: PetscInt Nf = prob->Nf, f;
481: PetscErrorCode ierr;
484: if (Nf >= NfNew) return(0);
485: prob->setup = PETSC_FALSE;
486: PetscDSDestroyStructs_Static(prob);
487: PetscMalloc3(NfNew, &tmpd, NfNew, &tmpi, NfNew, &tmpk);
488: for (f = 0; f < Nf; ++f) {tmpd[f] = prob->disc[f]; tmpi[f] = prob->implicit[f]; tmpk[f] = prob->jetDegree[f];}
489: for (f = Nf; f < NfNew; ++f) {tmpd[f] = NULL; tmpi[f] = PETSC_TRUE; tmpk[f] = 1;}
490: PetscFree3(prob->disc, prob->implicit, prob->jetDegree);
491: PetscWeakFormSetNumFields(prob->wf, NfNew);
492: prob->Nf = NfNew;
493: prob->disc = tmpd;
494: prob->implicit = tmpi;
495: prob->jetDegree = tmpk;
496: PetscCalloc2(NfNew, &tmpup, NfNew, &tmpctx);
497: for (f = 0; f < Nf; ++f) tmpup[f] = prob->update[f];
498: for (f = 0; f < Nf; ++f) tmpctx[f] = prob->ctx[f];
499: for (f = Nf; f < NfNew; ++f) tmpup[f] = NULL;
500: for (f = Nf; f < NfNew; ++f) tmpctx[f] = NULL;
501: PetscFree2(prob->update, prob->ctx);
502: prob->update = tmpup;
503: prob->ctx = tmpctx;
504: PetscCalloc4(NfNew, &tmpexactSol, NfNew, &tmpexactCtx, NfNew, &tmpexactSol_t, NfNew, &tmpexactCtx_t);
505: for (f = 0; f < Nf; ++f) tmpexactSol[f] = prob->exactSol[f];
506: for (f = 0; f < Nf; ++f) tmpexactCtx[f] = prob->exactCtx[f];
507: for (f = 0; f < Nf; ++f) tmpexactSol_t[f] = prob->exactSol_t[f];
508: for (f = 0; f < Nf; ++f) tmpexactCtx_t[f] = prob->exactCtx_t[f];
509: for (f = Nf; f < NfNew; ++f) tmpexactSol[f] = NULL;
510: for (f = Nf; f < NfNew; ++f) tmpexactCtx[f] = NULL;
511: for (f = Nf; f < NfNew; ++f) tmpexactSol_t[f] = NULL;
512: for (f = Nf; f < NfNew; ++f) tmpexactCtx_t[f] = NULL;
513: PetscFree4(prob->exactSol, prob->exactCtx, prob->exactSol_t, prob->exactCtx_t);
514: prob->exactSol = tmpexactSol;
515: prob->exactCtx = tmpexactCtx;
516: prob->exactSol_t = tmpexactSol_t;
517: prob->exactCtx_t = tmpexactCtx_t;
518: return(0);
519: }
521: /*@
522: PetscDSDestroy - Destroys a PetscDS object
524: Collective on prob
526: Input Parameter:
527: . prob - the PetscDS object to destroy
529: Level: developer
531: .seealso PetscDSView()
532: @*/
533: PetscErrorCode PetscDSDestroy(PetscDS *ds)
534: {
535: PetscInt f;
539: if (!*ds) return(0);
542: if (--((PetscObject)(*ds))->refct > 0) {*ds = NULL; return(0);}
543: ((PetscObject) (*ds))->refct = 0;
544: if ((*ds)->subprobs) {
545: PetscInt dim, d;
547: PetscDSGetSpatialDimension(*ds, &dim);
548: for (d = 0; d < dim; ++d) {PetscDSDestroy(&(*ds)->subprobs[d]);}
549: }
550: PetscFree((*ds)->subprobs);
551: PetscDSDestroyStructs_Static(*ds);
552: for (f = 0; f < (*ds)->Nf; ++f) {
553: PetscObjectDereference((*ds)->disc[f]);
554: }
555: PetscFree3((*ds)->disc, (*ds)->implicit, (*ds)->jetDegree);
556: PetscWeakFormDestroy(&(*ds)->wf);
557: PetscFree2((*ds)->update,(*ds)->ctx);
558: PetscFree4((*ds)->exactSol,(*ds)->exactCtx,(*ds)->exactSol_t,(*ds)->exactCtx_t);
559: if ((*ds)->ops->destroy) {(*(*ds)->ops->destroy)(*ds);}
560: PetscDSDestroyBoundary(*ds);
561: PetscFree((*ds)->constants);
562: PetscHeaderDestroy(ds);
563: return(0);
564: }
566: /*@
567: PetscDSCreate - Creates an empty PetscDS object. The type can then be set with PetscDSSetType().
569: Collective
571: Input Parameter:
572: . comm - The communicator for the PetscDS object
574: Output Parameter:
575: . ds - The PetscDS object
577: Level: beginner
579: .seealso: PetscDSSetType(), PETSCDSBASIC
580: @*/
581: PetscErrorCode PetscDSCreate(MPI_Comm comm, PetscDS *ds)
582: {
583: PetscDS p;
588: *ds = NULL;
589: PetscDSInitializePackage();
591: PetscHeaderCreate(p, PETSCDS_CLASSID, "PetscDS", "Discrete System", "PetscDS", comm, PetscDSDestroy, PetscDSView);
593: p->Nf = 0;
594: p->setup = PETSC_FALSE;
595: p->numConstants = 0;
596: p->constants = NULL;
597: p->dimEmbed = -1;
598: p->useJacPre = PETSC_TRUE;
599: PetscWeakFormCreate(comm, &p->wf);
601: *ds = p;
602: return(0);
603: }
605: /*@
606: PetscDSGetNumFields - Returns the number of fields in the DS
608: Not collective
610: Input Parameter:
611: . prob - The PetscDS object
613: Output Parameter:
614: . Nf - The number of fields
616: Level: beginner
618: .seealso: PetscDSGetSpatialDimension(), PetscDSCreate()
619: @*/
620: PetscErrorCode PetscDSGetNumFields(PetscDS prob, PetscInt *Nf)
621: {
625: *Nf = prob->Nf;
626: return(0);
627: }
629: /*@
630: PetscDSGetSpatialDimension - Returns the spatial dimension of the DS, meaning the topological dimension of the discretizations
632: Not collective
634: Input Parameter:
635: . prob - The PetscDS object
637: Output Parameter:
638: . dim - The spatial dimension
640: Level: beginner
642: .seealso: PetscDSGetCoordinateDimension(), PetscDSGetNumFields(), PetscDSCreate()
643: @*/
644: PetscErrorCode PetscDSGetSpatialDimension(PetscDS prob, PetscInt *dim)
645: {
651: *dim = 0;
652: if (prob->Nf) {
653: PetscObject obj;
654: PetscClassId id;
656: PetscDSGetDiscretization(prob, 0, &obj);
657: if (obj) {
658: PetscObjectGetClassId(obj, &id);
659: if (id == PETSCFE_CLASSID) {PetscFEGetSpatialDimension((PetscFE) obj, dim);}
660: else if (id == PETSCFV_CLASSID) {PetscFVGetSpatialDimension((PetscFV) obj, dim);}
661: else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", 0);
662: }
663: }
664: return(0);
665: }
667: /*@
668: PetscDSGetCoordinateDimension - Returns the coordinate dimension of the DS, meaning the dimension of the space into which the discretiaztions are embedded
670: Not collective
672: Input Parameter:
673: . prob - The PetscDS object
675: Output Parameter:
676: . dimEmbed - The coordinate dimension
678: Level: beginner
680: .seealso: PetscDSSetCoordinateDimension(), PetscDSGetSpatialDimension(), PetscDSGetNumFields(), PetscDSCreate()
681: @*/
682: PetscErrorCode PetscDSGetCoordinateDimension(PetscDS prob, PetscInt *dimEmbed)
683: {
687: if (prob->dimEmbed < 0) SETERRQ(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONGSTATE, "No coordinate dimension set for this DS");
688: *dimEmbed = prob->dimEmbed;
689: return(0);
690: }
692: /*@
693: PetscDSSetCoordinateDimension - Set the coordinate dimension of the DS, meaning the dimension of the space into which the discretiaztions are embedded
695: Logically collective on prob
697: Input Parameters:
698: + prob - The PetscDS object
699: - dimEmbed - The coordinate dimension
701: Level: beginner
703: .seealso: PetscDSGetCoordinateDimension(), PetscDSGetSpatialDimension(), PetscDSGetNumFields(), PetscDSCreate()
704: @*/
705: PetscErrorCode PetscDSSetCoordinateDimension(PetscDS prob, PetscInt dimEmbed)
706: {
709: if (dimEmbed < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Coordinate dimension must be non-negative, not %D", dimEmbed);
710: prob->dimEmbed = dimEmbed;
711: return(0);
712: }
714: /*@
715: PetscDSGetHybrid - Returns the flag for a hybrid (cohesive) cell
717: Not collective
719: Input Parameter:
720: . prob - The PetscDS object
722: Output Parameter:
723: . isHybrid - The flag
725: Level: developer
727: .seealso: PetscDSSetHybrid(), PetscDSCreate()
728: @*/
729: PetscErrorCode PetscDSGetHybrid(PetscDS prob, PetscBool *isHybrid)
730: {
734: *isHybrid = prob->isHybrid;
735: return(0);
736: }
738: /*@
739: PetscDSSetHybrid - Set the flag for a hybrid (cohesive) cell
741: Not collective
743: Input Parameters:
744: + prob - The PetscDS object
745: - isHybrid - The flag
747: Level: developer
749: .seealso: PetscDSGetHybrid(), PetscDSCreate()
750: @*/
751: PetscErrorCode PetscDSSetHybrid(PetscDS prob, PetscBool isHybrid)
752: {
755: prob->isHybrid = isHybrid;
756: return(0);
757: }
759: /*@
760: PetscDSGetTotalDimension - Returns the total size of the approximation space for this system
762: Not collective
764: Input Parameter:
765: . prob - The PetscDS object
767: Output Parameter:
768: . dim - The total problem dimension
770: Level: beginner
772: .seealso: PetscDSGetNumFields(), PetscDSCreate()
773: @*/
774: PetscErrorCode PetscDSGetTotalDimension(PetscDS prob, PetscInt *dim)
775: {
780: PetscDSSetUp(prob);
782: *dim = prob->totDim;
783: return(0);
784: }
786: /*@
787: PetscDSGetTotalComponents - Returns the total number of components in this system
789: Not collective
791: Input Parameter:
792: . prob - The PetscDS object
794: Output Parameter:
795: . dim - The total number of components
797: Level: beginner
799: .seealso: PetscDSGetNumFields(), PetscDSCreate()
800: @*/
801: PetscErrorCode PetscDSGetTotalComponents(PetscDS prob, PetscInt *Nc)
802: {
807: PetscDSSetUp(prob);
809: *Nc = prob->totComp;
810: return(0);
811: }
813: /*@
814: PetscDSGetDiscretization - Returns the discretization object for the given field
816: Not collective
818: Input Parameters:
819: + prob - The PetscDS object
820: - f - The field number
822: Output Parameter:
823: . disc - The discretization object
825: Level: beginner
827: .seealso: PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
828: @*/
829: PetscErrorCode PetscDSGetDiscretization(PetscDS prob, PetscInt f, PetscObject *disc)
830: {
834: if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
835: *disc = prob->disc[f];
836: return(0);
837: }
839: /*@
840: PetscDSSetDiscretization - Sets the discretization object for the given field
842: Not collective
844: Input Parameters:
845: + prob - The PetscDS object
846: . f - The field number
847: - disc - The discretization object
849: Level: beginner
851: .seealso: PetscDSGetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
852: @*/
853: PetscErrorCode PetscDSSetDiscretization(PetscDS prob, PetscInt f, PetscObject disc)
854: {
860: if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
861: PetscDSEnlarge_Static(prob, f+1);
862: PetscObjectDereference(prob->disc[f]);
863: prob->disc[f] = disc;
864: PetscObjectReference(disc);
865: if (disc) {
866: PetscClassId id;
868: PetscObjectGetClassId(disc, &id);
869: if (id == PETSCFE_CLASSID) {
870: PetscDSSetImplicit(prob, f, PETSC_TRUE);
871: } else if (id == PETSCFV_CLASSID) {
872: PetscDSSetImplicit(prob, f, PETSC_FALSE);
873: }
874: PetscDSSetJetDegree(prob, f, 1);
875: }
876: return(0);
877: }
879: /*@
880: PetscDSGetWeakForm - Returns the weak form object
882: Not collective
884: Input Parameter:
885: . ds - The PetscDS object
887: Output Parameter:
888: . wf - The weak form object
890: Level: beginner
892: .seealso: PetscDSSetWeakForm(), PetscDSGetNumFields(), PetscDSCreate()
893: @*/
894: PetscErrorCode PetscDSGetWeakForm(PetscDS ds, PetscWeakForm *wf)
895: {
899: *wf = ds->wf;
900: return(0);
901: }
903: /*@
904: PetscDSSetWeakForm - Sets the weak form object
906: Not collective
908: Input Parameters:
909: + ds - The PetscDS object
910: - wf - The weak form object
912: Level: beginner
914: .seealso: PetscDSGetWeakForm(), PetscDSGetNumFields(), PetscDSCreate()
915: @*/
916: PetscErrorCode PetscDSSetWeakForm(PetscDS ds, PetscWeakForm wf)
917: {
923: PetscObjectDereference((PetscObject) ds->wf);
924: ds->wf = wf;
925: PetscObjectReference((PetscObject) wf);
926: PetscWeakFormSetNumFields(wf, ds->Nf);
927: return(0);
928: }
930: /*@
931: PetscDSAddDiscretization - Adds a discretization object
933: Not collective
935: Input Parameters:
936: + prob - The PetscDS object
937: - disc - The boundary discretization object
939: Level: beginner
941: .seealso: PetscDSGetDiscretization(), PetscDSSetDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
942: @*/
943: PetscErrorCode PetscDSAddDiscretization(PetscDS prob, PetscObject disc)
944: {
948: PetscDSSetDiscretization(prob, prob->Nf, disc);
949: return(0);
950: }
952: /*@
953: PetscDSGetQuadrature - Returns the quadrature, which must agree for all fields in the DS
955: Not collective
957: Input Parameter:
958: . prob - The PetscDS object
960: Output Parameter:
961: . q - The quadrature object
963: Level: intermediate
965: .seealso: PetscDSSetImplicit(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
966: @*/
967: PetscErrorCode PetscDSGetQuadrature(PetscDS prob, PetscQuadrature *q)
968: {
969: PetscObject obj;
970: PetscClassId id;
974: *q = NULL;
975: if (!prob->Nf) return(0);
976: PetscDSGetDiscretization(prob, 0, &obj);
977: PetscObjectGetClassId(obj, &id);
978: if (id == PETSCFE_CLASSID) {PetscFEGetQuadrature((PetscFE) obj, q);}
979: else if (id == PETSCFV_CLASSID) {PetscFVGetQuadrature((PetscFV) obj, q);}
980: else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", 0);
981: return(0);
982: }
984: /*@
985: PetscDSGetImplicit - Returns the flag for implicit solve for this field. This is just a guide for IMEX
987: Not collective
989: Input Parameters:
990: + prob - The PetscDS object
991: - f - The field number
993: Output Parameter:
994: . implicit - The flag indicating what kind of solve to use for this field
996: Level: developer
998: .seealso: PetscDSSetImplicit(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
999: @*/
1000: PetscErrorCode PetscDSGetImplicit(PetscDS prob, PetscInt f, PetscBool *implicit)
1001: {
1005: if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
1006: *implicit = prob->implicit[f];
1007: return(0);
1008: }
1010: /*@
1011: PetscDSSetImplicit - Set the flag for implicit solve for this field. This is just a guide for IMEX
1013: Not collective
1015: Input Parameters:
1016: + prob - The PetscDS object
1017: . f - The field number
1018: - implicit - The flag indicating what kind of solve to use for this field
1020: Level: developer
1022: .seealso: PetscDSGetImplicit(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
1023: @*/
1024: PetscErrorCode PetscDSSetImplicit(PetscDS prob, PetscInt f, PetscBool implicit)
1025: {
1028: if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
1029: prob->implicit[f] = implicit;
1030: return(0);
1031: }
1033: /*@
1034: PetscDSGetJetDegree - Returns the highest derivative for this field equation, or the k-jet that the discretization needs to tabulate.
1036: Not collective
1038: Input Parameters:
1039: + ds - The PetscDS object
1040: - f - The field number
1042: Output Parameter:
1043: . k - The highest derivative we need to tabulate
1045: Level: developer
1047: .seealso: PetscDSSetJetDegree(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
1048: @*/
1049: PetscErrorCode PetscDSGetJetDegree(PetscDS ds, PetscInt f, PetscInt *k)
1050: {
1054: if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
1055: *k = ds->jetDegree[f];
1056: return(0);
1057: }
1059: /*@
1060: PetscDSSetJetDegree - Set the highest derivative for this field equation, or the k-jet that the discretization needs to tabulate.
1062: Not collective
1064: Input Parameters:
1065: + ds - The PetscDS object
1066: . f - The field number
1067: - k - The highest derivative we need to tabulate
1069: Level: developer
1071: .seealso: PetscDSGetJetDegree(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
1072: @*/
1073: PetscErrorCode PetscDSSetJetDegree(PetscDS ds, PetscInt f, PetscInt k)
1074: {
1077: if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
1078: ds->jetDegree[f] = k;
1079: return(0);
1080: }
1082: PetscErrorCode PetscDSGetObjective(PetscDS ds, PetscInt f,
1083: void (**obj)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1084: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1085: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1086: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[]))
1087: {
1088: PetscPointFunc *tmp;
1089: PetscInt n;
1090: PetscErrorCode ierr;
1095: if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
1096: PetscWeakFormGetObjective(ds->wf, NULL, 0, f, 0, &n, &tmp);
1097: *obj = tmp ? tmp[0] : NULL;
1098: return(0);
1099: }
1101: PetscErrorCode PetscDSSetObjective(PetscDS ds, PetscInt f,
1102: void (*obj)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1103: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1104: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1105: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[]))
1106: {
1112: if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1113: PetscWeakFormSetIndexObjective(ds->wf, NULL, 0, f, 0, 0, obj);
1114: return(0);
1115: }
1117: /*@C
1118: PetscDSGetResidual - Get the pointwise residual function for a given test field
1120: Not collective
1122: Input Parameters:
1123: + ds - The PetscDS
1124: - f - The test field number
1126: Output Parameters:
1127: + f0 - integrand for the test function term
1128: - f1 - integrand for the test function gradient term
1130: Note: We are using a first order FEM model for the weak form:
1132: \int_\Omega \phi f_0(u, u_t, \nabla u, x, t) + \nabla\phi \cdot {\vec f}_1(u, u_t, \nabla u, x, t)
1134: The calling sequence for the callbacks f0 and f1 is given by:
1136: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1137: $ const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1138: $ const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1139: $ PetscReal t, const PetscReal x[], PetscScalar f0[])
1141: + dim - the spatial dimension
1142: . Nf - the number of fields
1143: . uOff - the offset into u[] and u_t[] for each field
1144: . uOff_x - the offset into u_x[] for each field
1145: . u - each field evaluated at the current point
1146: . u_t - the time derivative of each field evaluated at the current point
1147: . u_x - the gradient of each field evaluated at the current point
1148: . aOff - the offset into a[] and a_t[] for each auxiliary field
1149: . aOff_x - the offset into a_x[] for each auxiliary field
1150: . a - each auxiliary field evaluated at the current point
1151: . a_t - the time derivative of each auxiliary field evaluated at the current point
1152: . a_x - the gradient of auxiliary each field evaluated at the current point
1153: . t - current time
1154: . x - coordinates of the current point
1155: . numConstants - number of constant parameters
1156: . constants - constant parameters
1157: - f0 - output values at the current point
1159: Level: intermediate
1161: .seealso: PetscDSSetResidual()
1162: @*/
1163: PetscErrorCode PetscDSGetResidual(PetscDS ds, PetscInt f,
1164: void (**f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1165: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1166: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1167: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]),
1168: void (**f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1169: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1170: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1171: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
1172: {
1173: PetscPointFunc *tmp0, *tmp1;
1174: PetscInt n0, n1;
1175: PetscErrorCode ierr;
1179: if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
1180: PetscWeakFormGetResidual(ds->wf, NULL, 0, f, 0, &n0, &tmp0, &n1, &tmp1);
1181: *f0 = tmp0 ? tmp0[0] : NULL;
1182: *f1 = tmp1 ? tmp1[0] : NULL;
1183: return(0);
1184: }
1186: /*@C
1187: PetscDSSetResidual - Set the pointwise residual function for a given test field
1189: Not collective
1191: Input Parameters:
1192: + ds - The PetscDS
1193: . f - The test field number
1194: . f0 - integrand for the test function term
1195: - f1 - integrand for the test function gradient term
1197: Note: We are using a first order FEM model for the weak form:
1199: \int_\Omega \phi f_0(u, u_t, \nabla u, x, t) + \nabla\phi \cdot {\vec f}_1(u, u_t, \nabla u, x, t)
1201: The calling sequence for the callbacks f0 and f1 is given by:
1203: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1204: $ const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1205: $ const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1206: $ PetscReal t, const PetscReal x[], PetscScalar f0[])
1208: + dim - the spatial dimension
1209: . Nf - the number of fields
1210: . uOff - the offset into u[] and u_t[] for each field
1211: . uOff_x - the offset into u_x[] for each field
1212: . u - each field evaluated at the current point
1213: . u_t - the time derivative of each field evaluated at the current point
1214: . u_x - the gradient of each field evaluated at the current point
1215: . aOff - the offset into a[] and a_t[] for each auxiliary field
1216: . aOff_x - the offset into a_x[] for each auxiliary field
1217: . a - each auxiliary field evaluated at the current point
1218: . a_t - the time derivative of each auxiliary field evaluated at the current point
1219: . a_x - the gradient of auxiliary each field evaluated at the current point
1220: . t - current time
1221: . x - coordinates of the current point
1222: . numConstants - number of constant parameters
1223: . constants - constant parameters
1224: - f0 - output values at the current point
1226: Level: intermediate
1228: .seealso: PetscDSGetResidual()
1229: @*/
1230: PetscErrorCode PetscDSSetResidual(PetscDS ds, PetscInt f,
1231: void (*f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1232: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1233: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1234: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]),
1235: void (*f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1236: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1237: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1238: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
1239: {
1246: if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1247: PetscWeakFormSetIndexResidual(ds->wf, NULL, 0, f, 0, 0, f0, 0, f1);
1248: return(0);
1249: }
1251: /*@C
1252: PetscDSHasJacobian - Signals that Jacobian functions have been set
1254: Not collective
1256: Input Parameter:
1257: . prob - The PetscDS
1259: Output Parameter:
1260: . hasJac - flag that pointwise function for the Jacobian has been set
1262: Level: intermediate
1264: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1265: @*/
1266: PetscErrorCode PetscDSHasJacobian(PetscDS ds, PetscBool *hasJac)
1267: {
1272: PetscWeakFormHasJacobian(ds->wf, hasJac);
1273: return(0);
1274: }
1276: /*@C
1277: PetscDSGetJacobian - Get the pointwise Jacobian function for given test and basis field
1279: Not collective
1281: Input Parameters:
1282: + ds - The PetscDS
1283: . f - The test field number
1284: - g - The field number
1286: Output Parameters:
1287: + g0 - integrand for the test and basis function term
1288: . g1 - integrand for the test function and basis function gradient term
1289: . g2 - integrand for the test function gradient and basis function term
1290: - g3 - integrand for the test function gradient and basis function gradient term
1292: Note: We are using a first order FEM model for the weak form:
1294: \int_\Omega \phi g_0(u, u_t, \nabla u, x, t) \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi
1296: The calling sequence for the callbacks g0, g1, g2 and g3 is given by:
1298: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1299: $ const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1300: $ const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1301: $ PetscReal t, const PetscReal u_tShift, const PetscReal x[], PetscScalar g0[])
1303: + dim - the spatial dimension
1304: . Nf - the number of fields
1305: . uOff - the offset into u[] and u_t[] for each field
1306: . uOff_x - the offset into u_x[] for each field
1307: . u - each field evaluated at the current point
1308: . u_t - the time derivative of each field evaluated at the current point
1309: . u_x - the gradient of each field evaluated at the current point
1310: . aOff - the offset into a[] and a_t[] for each auxiliary field
1311: . aOff_x - the offset into a_x[] for each auxiliary field
1312: . a - each auxiliary field evaluated at the current point
1313: . a_t - the time derivative of each auxiliary field evaluated at the current point
1314: . a_x - the gradient of auxiliary each field evaluated at the current point
1315: . t - current time
1316: . u_tShift - the multiplier a for dF/dU_t
1317: . x - coordinates of the current point
1318: . numConstants - number of constant parameters
1319: . constants - constant parameters
1320: - g0 - output values at the current point
1322: Level: intermediate
1324: .seealso: PetscDSSetJacobian()
1325: @*/
1326: PetscErrorCode PetscDSGetJacobian(PetscDS ds, PetscInt f, PetscInt g,
1327: void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1328: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1329: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1330: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1331: void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1332: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1333: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1334: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1335: void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1336: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1337: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1338: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1339: void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1340: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1341: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1342: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1343: {
1344: PetscPointJac *tmp0, *tmp1, *tmp2, *tmp3;
1345: PetscInt n0, n1, n2, n3;
1350: if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
1351: if ((g < 0) || (g >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, ds->Nf);
1352: PetscWeakFormGetJacobian(ds->wf, NULL, 0, f, g, 0, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3);
1353: *g0 = tmp0 ? tmp0[0] : NULL;
1354: *g1 = tmp1 ? tmp1[0] : NULL;
1355: *g2 = tmp2 ? tmp2[0] : NULL;
1356: *g3 = tmp3 ? tmp3[0] : NULL;
1357: return(0);
1358: }
1360: /*@C
1361: PetscDSSetJacobian - Set the pointwise Jacobian function for given test and basis fields
1363: Not collective
1365: Input Parameters:
1366: + ds - The PetscDS
1367: . f - The test field number
1368: . g - The field number
1369: . g0 - integrand for the test and basis function term
1370: . g1 - integrand for the test function and basis function gradient term
1371: . g2 - integrand for the test function gradient and basis function term
1372: - g3 - integrand for the test function gradient and basis function gradient term
1374: Note: We are using a first order FEM model for the weak form:
1376: \int_\Omega \phi g_0(u, u_t, \nabla u, x, t) \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi
1378: The calling sequence for the callbacks g0, g1, g2 and g3 is given by:
1380: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1381: $ const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1382: $ const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1383: $ PetscReal t, const PetscReal x[], PetscScalar g0[])
1385: + dim - the spatial dimension
1386: . Nf - the number of fields
1387: . uOff - the offset into u[] and u_t[] for each field
1388: . uOff_x - the offset into u_x[] for each field
1389: . u - each field evaluated at the current point
1390: . u_t - the time derivative of each field evaluated at the current point
1391: . u_x - the gradient of each field evaluated at the current point
1392: . aOff - the offset into a[] and a_t[] for each auxiliary field
1393: . aOff_x - the offset into a_x[] for each auxiliary field
1394: . a - each auxiliary field evaluated at the current point
1395: . a_t - the time derivative of each auxiliary field evaluated at the current point
1396: . a_x - the gradient of auxiliary each field evaluated at the current point
1397: . t - current time
1398: . u_tShift - the multiplier a for dF/dU_t
1399: . x - coordinates of the current point
1400: . numConstants - number of constant parameters
1401: . constants - constant parameters
1402: - g0 - output values at the current point
1404: Level: intermediate
1406: .seealso: PetscDSGetJacobian()
1407: @*/
1408: PetscErrorCode PetscDSSetJacobian(PetscDS ds, PetscInt f, PetscInt g,
1409: void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1410: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1411: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1412: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1413: void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1414: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1415: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1416: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1417: void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1418: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1419: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1420: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1421: void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1422: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1423: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1424: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1425: {
1434: if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1435: if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
1436: PetscWeakFormSetIndexJacobian(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3);
1437: return(0);
1438: }
1440: /*@C
1441: PetscDSUseJacobianPreconditioner - Whether to construct a Jacobian preconditioner
1443: Not collective
1445: Input Parameters:
1446: + prob - The PetscDS
1447: - useJacPre - flag that enables construction of a Jacobian preconditioner
1449: Level: intermediate
1451: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1452: @*/
1453: PetscErrorCode PetscDSUseJacobianPreconditioner(PetscDS prob, PetscBool useJacPre)
1454: {
1457: prob->useJacPre = useJacPre;
1458: return(0);
1459: }
1461: /*@C
1462: PetscDSHasJacobianPreconditioner - Signals that a Jacobian preconditioner matrix has been set
1464: Not collective
1466: Input Parameter:
1467: . prob - The PetscDS
1469: Output Parameter:
1470: . hasJacPre - flag that pointwise function for Jacobian preconditioner matrix has been set
1472: Level: intermediate
1474: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1475: @*/
1476: PetscErrorCode PetscDSHasJacobianPreconditioner(PetscDS ds, PetscBool *hasJacPre)
1477: {
1482: *hasJacPre = PETSC_FALSE;
1483: if (!ds->useJacPre) return(0);
1484: PetscWeakFormHasJacobianPreconditioner(ds->wf, hasJacPre);
1485: return(0);
1486: }
1488: /*@C
1489: PetscDSGetJacobianPreconditioner - Get the pointwise Jacobian preconditioner function for given test and basis field. If this is missing, the system matrix is used to build the preconditioner.
1491: Not collective
1493: Input Parameters:
1494: + ds - The PetscDS
1495: . f - The test field number
1496: - g - The field number
1498: Output Parameters:
1499: + g0 - integrand for the test and basis function term
1500: . g1 - integrand for the test function and basis function gradient term
1501: . g2 - integrand for the test function gradient and basis function term
1502: - g3 - integrand for the test function gradient and basis function gradient term
1504: Note: We are using a first order FEM model for the weak form:
1506: \int_\Omega \phi g_0(u, u_t, \nabla u, x, t) \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi
1508: The calling sequence for the callbacks g0, g1, g2 and g3 is given by:
1510: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1511: $ const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1512: $ const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1513: $ PetscReal t, const PetscReal u_tShift, const PetscReal x[], PetscScalar g0[])
1515: + dim - the spatial dimension
1516: . Nf - the number of fields
1517: . uOff - the offset into u[] and u_t[] for each field
1518: . uOff_x - the offset into u_x[] for each field
1519: . u - each field evaluated at the current point
1520: . u_t - the time derivative of each field evaluated at the current point
1521: . u_x - the gradient of each field evaluated at the current point
1522: . aOff - the offset into a[] and a_t[] for each auxiliary field
1523: . aOff_x - the offset into a_x[] for each auxiliary field
1524: . a - each auxiliary field evaluated at the current point
1525: . a_t - the time derivative of each auxiliary field evaluated at the current point
1526: . a_x - the gradient of auxiliary each field evaluated at the current point
1527: . t - current time
1528: . u_tShift - the multiplier a for dF/dU_t
1529: . x - coordinates of the current point
1530: . numConstants - number of constant parameters
1531: . constants - constant parameters
1532: - g0 - output values at the current point
1534: Level: intermediate
1536: .seealso: PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
1537: @*/
1538: PetscErrorCode PetscDSGetJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g,
1539: void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1540: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1541: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1542: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1543: void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1544: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1545: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1546: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1547: void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1548: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1549: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1550: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1551: void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1552: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1553: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1554: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1555: {
1556: PetscPointJac *tmp0, *tmp1, *tmp2, *tmp3;
1557: PetscInt n0, n1, n2, n3;
1562: if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
1563: if ((g < 0) || (g >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, ds->Nf);
1564: PetscWeakFormGetJacobianPreconditioner(ds->wf, NULL, 0, f, g, 0, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3);
1565: *g0 = tmp0 ? tmp0[0] : NULL;
1566: *g1 = tmp1 ? tmp1[0] : NULL;
1567: *g2 = tmp2 ? tmp2[0] : NULL;
1568: *g3 = tmp3 ? tmp3[0] : NULL;
1569: return(0);
1570: }
1572: /*@C
1573: PetscDSSetJacobianPreconditioner - Set the pointwise Jacobian preconditioner function for given test and basis fields. If this is missing, the system matrix is used to build the preconditioner.
1575: Not collective
1577: Input Parameters:
1578: + ds - The PetscDS
1579: . f - The test field number
1580: . g - The field number
1581: . g0 - integrand for the test and basis function term
1582: . g1 - integrand for the test function and basis function gradient term
1583: . g2 - integrand for the test function gradient and basis function term
1584: - g3 - integrand for the test function gradient and basis function gradient term
1586: Note: We are using a first order FEM model for the weak form:
1588: \int_\Omega \phi g_0(u, u_t, \nabla u, x, t) \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi
1590: The calling sequence for the callbacks g0, g1, g2 and g3 is given by:
1592: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1593: $ const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1594: $ const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1595: $ PetscReal t, const PetscReal x[], PetscScalar g0[])
1597: + dim - the spatial dimension
1598: . Nf - the number of fields
1599: . uOff - the offset into u[] and u_t[] for each field
1600: . uOff_x - the offset into u_x[] for each field
1601: . u - each field evaluated at the current point
1602: . u_t - the time derivative of each field evaluated at the current point
1603: . u_x - the gradient of each field evaluated at the current point
1604: . aOff - the offset into a[] and a_t[] for each auxiliary field
1605: . aOff_x - the offset into a_x[] for each auxiliary field
1606: . a - each auxiliary field evaluated at the current point
1607: . a_t - the time derivative of each auxiliary field evaluated at the current point
1608: . a_x - the gradient of auxiliary each field evaluated at the current point
1609: . t - current time
1610: . u_tShift - the multiplier a for dF/dU_t
1611: . x - coordinates of the current point
1612: . numConstants - number of constant parameters
1613: . constants - constant parameters
1614: - g0 - output values at the current point
1616: Level: intermediate
1618: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobian()
1619: @*/
1620: PetscErrorCode PetscDSSetJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g,
1621: void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1622: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1623: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1624: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1625: void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1626: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1627: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1628: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1629: void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1630: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1631: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1632: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1633: void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1634: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1635: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1636: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1637: {
1646: if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1647: if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
1648: PetscWeakFormSetIndexJacobianPreconditioner(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3);
1649: return(0);
1650: }
1652: /*@C
1653: PetscDSHasDynamicJacobian - Signals that a dynamic Jacobian, dF/du_t, has been set
1655: Not collective
1657: Input Parameter:
1658: . ds - The PetscDS
1660: Output Parameter:
1661: . hasDynJac - flag that pointwise function for dynamic Jacobian has been set
1663: Level: intermediate
1665: .seealso: PetscDSGetDynamicJacobian(), PetscDSSetDynamicJacobian(), PetscDSGetJacobian()
1666: @*/
1667: PetscErrorCode PetscDSHasDynamicJacobian(PetscDS ds, PetscBool *hasDynJac)
1668: {
1673: PetscWeakFormHasDynamicJacobian(ds->wf, hasDynJac);
1674: return(0);
1675: }
1677: /*@C
1678: PetscDSGetDynamicJacobian - Get the pointwise dynamic Jacobian, dF/du_t, function for given test and basis field
1680: Not collective
1682: Input Parameters:
1683: + ds - The PetscDS
1684: . f - The test field number
1685: - g - The field number
1687: Output Parameters:
1688: + g0 - integrand for the test and basis function term
1689: . g1 - integrand for the test function and basis function gradient term
1690: . g2 - integrand for the test function gradient and basis function term
1691: - g3 - integrand for the test function gradient and basis function gradient term
1693: Note: We are using a first order FEM model for the weak form:
1695: \int_\Omega \phi g_0(u, u_t, \nabla u, x, t) \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi
1697: The calling sequence for the callbacks g0, g1, g2 and g3 is given by:
1699: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1700: $ const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1701: $ const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1702: $ PetscReal t, const PetscReal u_tShift, const PetscReal x[], PetscScalar g0[])
1704: + dim - the spatial dimension
1705: . Nf - the number of fields
1706: . uOff - the offset into u[] and u_t[] for each field
1707: . uOff_x - the offset into u_x[] for each field
1708: . u - each field evaluated at the current point
1709: . u_t - the time derivative of each field evaluated at the current point
1710: . u_x - the gradient of each field evaluated at the current point
1711: . aOff - the offset into a[] and a_t[] for each auxiliary field
1712: . aOff_x - the offset into a_x[] for each auxiliary field
1713: . a - each auxiliary field evaluated at the current point
1714: . a_t - the time derivative of each auxiliary field evaluated at the current point
1715: . a_x - the gradient of auxiliary each field evaluated at the current point
1716: . t - current time
1717: . u_tShift - the multiplier a for dF/dU_t
1718: . x - coordinates of the current point
1719: . numConstants - number of constant parameters
1720: . constants - constant parameters
1721: - g0 - output values at the current point
1723: Level: intermediate
1725: .seealso: PetscDSSetJacobian()
1726: @*/
1727: PetscErrorCode PetscDSGetDynamicJacobian(PetscDS ds, PetscInt f, PetscInt g,
1728: void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1729: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1730: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1731: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1732: void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1733: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1734: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1735: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1736: void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1737: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1738: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1739: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1740: void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1741: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1742: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1743: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1744: {
1745: PetscPointJac *tmp0, *tmp1, *tmp2, *tmp3;
1746: PetscInt n0, n1, n2, n3;
1751: if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
1752: if ((g < 0) || (g >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, ds->Nf);
1753: PetscWeakFormGetDynamicJacobian(ds->wf, NULL, 0, f, g, 0, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3);
1754: *g0 = tmp0 ? tmp0[0] : NULL;
1755: *g1 = tmp1 ? tmp1[0] : NULL;
1756: *g2 = tmp2 ? tmp2[0] : NULL;
1757: *g3 = tmp3 ? tmp3[0] : NULL;
1758: return(0);
1759: }
1761: /*@C
1762: PetscDSSetDynamicJacobian - Set the pointwise dynamic Jacobian, dF/du_t, function for given test and basis fields
1764: Not collective
1766: Input Parameters:
1767: + ds - The PetscDS
1768: . f - The test field number
1769: . g - The field number
1770: . g0 - integrand for the test and basis function term
1771: . g1 - integrand for the test function and basis function gradient term
1772: . g2 - integrand for the test function gradient and basis function term
1773: - g3 - integrand for the test function gradient and basis function gradient term
1775: Note: We are using a first order FEM model for the weak form:
1777: \int_\Omega \phi g_0(u, u_t, \nabla u, x, t) \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi
1779: The calling sequence for the callbacks g0, g1, g2 and g3 is given by:
1781: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1782: $ const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1783: $ const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1784: $ PetscReal t, const PetscReal x[], PetscScalar g0[])
1786: + dim - the spatial dimension
1787: . Nf - the number of fields
1788: . uOff - the offset into u[] and u_t[] for each field
1789: . uOff_x - the offset into u_x[] for each field
1790: . u - each field evaluated at the current point
1791: . u_t - the time derivative of each field evaluated at the current point
1792: . u_x - the gradient of each field evaluated at the current point
1793: . aOff - the offset into a[] and a_t[] for each auxiliary field
1794: . aOff_x - the offset into a_x[] for each auxiliary field
1795: . a - each auxiliary field evaluated at the current point
1796: . a_t - the time derivative of each auxiliary field evaluated at the current point
1797: . a_x - the gradient of auxiliary each field evaluated at the current point
1798: . t - current time
1799: . u_tShift - the multiplier a for dF/dU_t
1800: . x - coordinates of the current point
1801: . numConstants - number of constant parameters
1802: . constants - constant parameters
1803: - g0 - output values at the current point
1805: Level: intermediate
1807: .seealso: PetscDSGetJacobian()
1808: @*/
1809: PetscErrorCode PetscDSSetDynamicJacobian(PetscDS ds, PetscInt f, PetscInt g,
1810: void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1811: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1812: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1813: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
1814: void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1815: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1816: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1817: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
1818: void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1819: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1820: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1821: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
1822: void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1823: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1824: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1825: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1826: {
1835: if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1836: if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
1837: PetscWeakFormSetIndexDynamicJacobian(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3);
1838: return(0);
1839: }
1841: /*@C
1842: PetscDSGetRiemannSolver - Returns the Riemann solver for the given field
1844: Not collective
1846: Input Parameters:
1847: + ds - The PetscDS object
1848: - f - The field number
1850: Output Parameter:
1851: . r - Riemann solver
1853: Calling sequence for r:
1855: $ r(PetscInt dim, PetscInt Nf, const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscScalar flux[], void *ctx)
1857: + dim - The spatial dimension
1858: . Nf - The number of fields
1859: . x - The coordinates at a point on the interface
1860: . n - The normal vector to the interface
1861: . uL - The state vector to the left of the interface
1862: . uR - The state vector to the right of the interface
1863: . flux - output array of flux through the interface
1864: . numConstants - number of constant parameters
1865: . constants - constant parameters
1866: - ctx - optional user context
1868: Level: intermediate
1870: .seealso: PetscDSSetRiemannSolver()
1871: @*/
1872: PetscErrorCode PetscDSGetRiemannSolver(PetscDS ds, PetscInt f,
1873: void (**r)(PetscInt dim, PetscInt Nf, const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscInt numConstants, const PetscScalar constants[], PetscScalar flux[], void *ctx))
1874: {
1875: PetscRiemannFunc *tmp;
1876: PetscInt n;
1877: PetscErrorCode ierr;
1882: if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
1883: PetscWeakFormGetRiemannSolver(ds->wf, NULL, 0, f, 0, &n, &tmp);
1884: *r = tmp ? tmp[0] : NULL;
1885: return(0);
1886: }
1888: /*@C
1889: PetscDSSetRiemannSolver - Sets the Riemann solver for the given field
1891: Not collective
1893: Input Parameters:
1894: + ds - The PetscDS object
1895: . f - The field number
1896: - r - Riemann solver
1898: Calling sequence for r:
1900: $ r(PetscInt dim, PetscInt Nf, const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscScalar flux[], void *ctx)
1902: + dim - The spatial dimension
1903: . Nf - The number of fields
1904: . x - The coordinates at a point on the interface
1905: . n - The normal vector to the interface
1906: . uL - The state vector to the left of the interface
1907: . uR - The state vector to the right of the interface
1908: . flux - output array of flux through the interface
1909: . numConstants - number of constant parameters
1910: . constants - constant parameters
1911: - ctx - optional user context
1913: Level: intermediate
1915: .seealso: PetscDSGetRiemannSolver()
1916: @*/
1917: PetscErrorCode PetscDSSetRiemannSolver(PetscDS ds, PetscInt f,
1918: void (*r)(PetscInt dim, PetscInt Nf, const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscInt numConstants, const PetscScalar constants[], PetscScalar flux[], void *ctx))
1919: {
1925: if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
1926: PetscWeakFormSetIndexRiemannSolver(ds->wf, NULL, 0, f, 0, 0, r);
1927: return(0);
1928: }
1930: /*@C
1931: PetscDSGetUpdate - Get the pointwise update function for a given field
1933: Not collective
1935: Input Parameters:
1936: + ds - The PetscDS
1937: - f - The field number
1939: Output Parameter:
1940: . update - update function
1942: Note: The calling sequence for the callback update is given by:
1944: $ update(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1945: $ const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1946: $ const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1947: $ PetscReal t, const PetscReal x[], PetscScalar uNew[])
1949: + dim - the spatial dimension
1950: . Nf - the number of fields
1951: . uOff - the offset into u[] and u_t[] for each field
1952: . uOff_x - the offset into u_x[] for each field
1953: . u - each field evaluated at the current point
1954: . u_t - the time derivative of each field evaluated at the current point
1955: . u_x - the gradient of each field evaluated at the current point
1956: . aOff - the offset into a[] and a_t[] for each auxiliary field
1957: . aOff_x - the offset into a_x[] for each auxiliary field
1958: . a - each auxiliary field evaluated at the current point
1959: . a_t - the time derivative of each auxiliary field evaluated at the current point
1960: . a_x - the gradient of auxiliary each field evaluated at the current point
1961: . t - current time
1962: . x - coordinates of the current point
1963: - uNew - new value for field at the current point
1965: Level: intermediate
1967: .seealso: PetscDSSetUpdate(), PetscDSSetResidual()
1968: @*/
1969: PetscErrorCode PetscDSGetUpdate(PetscDS ds, PetscInt f,
1970: void (**update)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1971: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1972: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1973: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar uNew[]))
1974: {
1977: if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
1979: return(0);
1980: }
1982: /*@C
1983: PetscDSSetUpdate - Set the pointwise update function for a given field
1985: Not collective
1987: Input Parameters:
1988: + ds - The PetscDS
1989: . f - The field number
1990: - update - update function
1992: Note: The calling sequence for the callback update is given by:
1994: $ update(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1995: $ const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1996: $ const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1997: $ PetscReal t, const PetscReal x[], PetscScalar uNew[])
1999: + dim - the spatial dimension
2000: . Nf - the number of fields
2001: . uOff - the offset into u[] and u_t[] for each field
2002: . uOff_x - the offset into u_x[] for each field
2003: . u - each field evaluated at the current point
2004: . u_t - the time derivative of each field evaluated at the current point
2005: . u_x - the gradient of each field evaluated at the current point
2006: . aOff - the offset into a[] and a_t[] for each auxiliary field
2007: . aOff_x - the offset into a_x[] for each auxiliary field
2008: . a - each auxiliary field evaluated at the current point
2009: . a_t - the time derivative of each auxiliary field evaluated at the current point
2010: . a_x - the gradient of auxiliary each field evaluated at the current point
2011: . t - current time
2012: . x - coordinates of the current point
2013: - uNew - new field values at the current point
2015: Level: intermediate
2017: .seealso: PetscDSGetResidual()
2018: @*/
2019: PetscErrorCode PetscDSSetUpdate(PetscDS ds, PetscInt f,
2020: void (*update)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2021: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2022: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2023: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar uNew[]))
2024: {
2030: if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2031: PetscDSEnlarge_Static(ds, f+1);
2032: ds->update[f] = update;
2033: return(0);
2034: }
2036: PetscErrorCode PetscDSGetContext(PetscDS ds, PetscInt f, void *ctx)
2037: {
2040: if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
2042: *(void**)ctx = ds->ctx[f];
2043: return(0);
2044: }
2046: PetscErrorCode PetscDSSetContext(PetscDS ds, PetscInt f, void *ctx)
2047: {
2052: if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2053: PetscDSEnlarge_Static(ds, f+1);
2054: ds->ctx[f] = ctx;
2055: return(0);
2056: }
2058: /*@C
2059: PetscDSGetBdResidual - Get the pointwise boundary residual function for a given test field
2061: Not collective
2063: Input Parameters:
2064: + ds - The PetscDS
2065: - f - The test field number
2067: Output Parameters:
2068: + f0 - boundary integrand for the test function term
2069: - f1 - boundary integrand for the test function gradient term
2071: Note: We are using a first order FEM model for the weak form:
2073: \int_\Gamma \phi {\vec f}_0(u, u_t, \nabla u, x, t) \cdot \hat n + \nabla\phi \cdot {\overleftrightarrow f}_1(u, u_t, \nabla u, x, t) \cdot \hat n
2075: The calling sequence for the callbacks f0 and f1 is given by:
2077: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2078: $ const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2079: $ const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2080: $ PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[])
2082: + dim - the spatial dimension
2083: . Nf - the number of fields
2084: . uOff - the offset into u[] and u_t[] for each field
2085: . uOff_x - the offset into u_x[] for each field
2086: . u - each field evaluated at the current point
2087: . u_t - the time derivative of each field evaluated at the current point
2088: . u_x - the gradient of each field evaluated at the current point
2089: . aOff - the offset into a[] and a_t[] for each auxiliary field
2090: . aOff_x - the offset into a_x[] for each auxiliary field
2091: . a - each auxiliary field evaluated at the current point
2092: . a_t - the time derivative of each auxiliary field evaluated at the current point
2093: . a_x - the gradient of auxiliary each field evaluated at the current point
2094: . t - current time
2095: . x - coordinates of the current point
2096: . n - unit normal at the current point
2097: . numConstants - number of constant parameters
2098: . constants - constant parameters
2099: - f0 - output values at the current point
2101: Level: intermediate
2103: .seealso: PetscDSSetBdResidual()
2104: @*/
2105: PetscErrorCode PetscDSGetBdResidual(PetscDS ds, PetscInt f,
2106: void (**f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2107: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2108: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2109: PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]),
2110: void (**f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2111: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2112: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2113: PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
2114: {
2115: PetscBdPointFunc *tmp0, *tmp1;
2116: PetscInt n0, n1;
2117: PetscErrorCode ierr;
2121: if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
2122: PetscWeakFormGetBdResidual(ds->wf, NULL, 0, f, 0, &n0, &tmp0, &n1, &tmp1);
2123: *f0 = tmp0 ? tmp0[0] : NULL;
2124: *f1 = tmp1 ? tmp1[0] : NULL;
2125: return(0);
2126: }
2128: /*@C
2129: PetscDSSetBdResidual - Get the pointwise boundary residual function for a given test field
2131: Not collective
2133: Input Parameters:
2134: + ds - The PetscDS
2135: . f - The test field number
2136: . f0 - boundary integrand for the test function term
2137: - f1 - boundary integrand for the test function gradient term
2139: Note: We are using a first order FEM model for the weak form:
2141: \int_\Gamma \phi {\vec f}_0(u, u_t, \nabla u, x, t) \cdot \hat n + \nabla\phi \cdot {\overleftrightarrow f}_1(u, u_t, \nabla u, x, t) \cdot \hat n
2143: The calling sequence for the callbacks f0 and f1 is given by:
2145: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2146: $ const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2147: $ const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2148: $ PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[])
2150: + dim - the spatial dimension
2151: . Nf - the number of fields
2152: . uOff - the offset into u[] and u_t[] for each field
2153: . uOff_x - the offset into u_x[] for each field
2154: . u - each field evaluated at the current point
2155: . u_t - the time derivative of each field evaluated at the current point
2156: . u_x - the gradient of each field evaluated at the current point
2157: . aOff - the offset into a[] and a_t[] for each auxiliary field
2158: . aOff_x - the offset into a_x[] for each auxiliary field
2159: . a - each auxiliary field evaluated at the current point
2160: . a_t - the time derivative of each auxiliary field evaluated at the current point
2161: . a_x - the gradient of auxiliary each field evaluated at the current point
2162: . t - current time
2163: . x - coordinates of the current point
2164: . n - unit normal at the current point
2165: . numConstants - number of constant parameters
2166: . constants - constant parameters
2167: - f0 - output values at the current point
2169: Level: intermediate
2171: .seealso: PetscDSGetBdResidual()
2172: @*/
2173: PetscErrorCode PetscDSSetBdResidual(PetscDS ds, PetscInt f,
2174: void (*f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2175: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2176: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2177: PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]),
2178: void (*f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2179: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2180: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2181: PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
2182: {
2187: if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2188: PetscWeakFormSetIndexBdResidual(ds->wf, NULL, 0, f, 0, 0, f0, 0, f1);
2189: return(0);
2190: }
2192: /*@
2193: PetscDSHasBdJacobian - Signals that boundary Jacobian functions have been set
2195: Not collective
2197: Input Parameter:
2198: . ds - The PetscDS
2200: Output Parameter:
2201: . hasBdJac - flag that pointwise function for the boundary Jacobian has been set
2203: Level: intermediate
2205: .seealso: PetscDSHasJacobian(), PetscDSSetBdJacobian(), PetscDSGetBdJacobian()
2206: @*/
2207: PetscErrorCode PetscDSHasBdJacobian(PetscDS ds, PetscBool *hasBdJac)
2208: {
2214: PetscWeakFormHasBdJacobian(ds->wf, hasBdJac);
2215: return(0);
2216: }
2218: /*@C
2219: PetscDSGetBdJacobian - Get the pointwise boundary Jacobian function for given test and basis field
2221: Not collective
2223: Input Parameters:
2224: + ds - The PetscDS
2225: . f - The test field number
2226: - g - The field number
2228: Output Parameters:
2229: + g0 - integrand for the test and basis function term
2230: . g1 - integrand for the test function and basis function gradient term
2231: . g2 - integrand for the test function gradient and basis function term
2232: - g3 - integrand for the test function gradient and basis function gradient term
2234: Note: We are using a first order FEM model for the weak form:
2236: \int_\Gamma \phi {\vec g}_0(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \cdot \hat n \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \hat n \cdot \nabla \psi
2238: The calling sequence for the callbacks g0, g1, g2 and g3 is given by:
2240: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2241: $ const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2242: $ const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2243: $ PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar g0[])
2245: + dim - the spatial dimension
2246: . Nf - the number of fields
2247: . uOff - the offset into u[] and u_t[] for each field
2248: . uOff_x - the offset into u_x[] for each field
2249: . u - each field evaluated at the current point
2250: . u_t - the time derivative of each field evaluated at the current point
2251: . u_x - the gradient of each field evaluated at the current point
2252: . aOff - the offset into a[] and a_t[] for each auxiliary field
2253: . aOff_x - the offset into a_x[] for each auxiliary field
2254: . a - each auxiliary field evaluated at the current point
2255: . a_t - the time derivative of each auxiliary field evaluated at the current point
2256: . a_x - the gradient of auxiliary each field evaluated at the current point
2257: . t - current time
2258: . u_tShift - the multiplier a for dF/dU_t
2259: . x - coordinates of the current point
2260: . n - normal at the current point
2261: . numConstants - number of constant parameters
2262: . constants - constant parameters
2263: - g0 - output values at the current point
2265: Level: intermediate
2267: .seealso: PetscDSSetBdJacobian()
2268: @*/
2269: PetscErrorCode PetscDSGetBdJacobian(PetscDS ds, PetscInt f, PetscInt g,
2270: void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2271: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2272: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2273: PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
2274: void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2275: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2276: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2277: PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
2278: void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2279: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2280: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2281: PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
2282: void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2283: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2284: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2285: PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
2286: {
2287: PetscBdPointJac *tmp0, *tmp1, *tmp2, *tmp3;
2288: PetscInt n0, n1, n2, n3;
2289: PetscErrorCode ierr;
2293: if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
2294: if ((g < 0) || (g >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, ds->Nf);
2295: PetscWeakFormGetBdJacobian(ds->wf, NULL, 0, f, g, 0, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3);
2296: *g0 = tmp0 ? tmp0[0] : NULL;
2297: *g1 = tmp1 ? tmp1[0] : NULL;
2298: *g2 = tmp2 ? tmp2[0] : NULL;
2299: *g3 = tmp3 ? tmp3[0] : NULL;
2300: return(0);
2301: }
2303: /*@C
2304: PetscDSSetBdJacobian - Set the pointwise boundary Jacobian function for given test and basis field
2306: Not collective
2308: Input Parameters:
2309: + ds - The PetscDS
2310: . f - The test field number
2311: . g - The field number
2312: . g0 - integrand for the test and basis function term
2313: . g1 - integrand for the test function and basis function gradient term
2314: . g2 - integrand for the test function gradient and basis function term
2315: - g3 - integrand for the test function gradient and basis function gradient term
2317: Note: We are using a first order FEM model for the weak form:
2319: \int_\Gamma \phi {\vec g}_0(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \cdot \hat n \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \hat n \cdot \nabla \psi
2321: The calling sequence for the callbacks g0, g1, g2 and g3 is given by:
2323: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2324: $ const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2325: $ const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2326: $ PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar g0[])
2328: + dim - the spatial dimension
2329: . Nf - the number of fields
2330: . uOff - the offset into u[] and u_t[] for each field
2331: . uOff_x - the offset into u_x[] for each field
2332: . u - each field evaluated at the current point
2333: . u_t - the time derivative of each field evaluated at the current point
2334: . u_x - the gradient of each field evaluated at the current point
2335: . aOff - the offset into a[] and a_t[] for each auxiliary field
2336: . aOff_x - the offset into a_x[] for each auxiliary field
2337: . a - each auxiliary field evaluated at the current point
2338: . a_t - the time derivative of each auxiliary field evaluated at the current point
2339: . a_x - the gradient of auxiliary each field evaluated at the current point
2340: . t - current time
2341: . u_tShift - the multiplier a for dF/dU_t
2342: . x - coordinates of the current point
2343: . n - normal at the current point
2344: . numConstants - number of constant parameters
2345: . constants - constant parameters
2346: - g0 - output values at the current point
2348: Level: intermediate
2350: .seealso: PetscDSGetBdJacobian()
2351: @*/
2352: PetscErrorCode PetscDSSetBdJacobian(PetscDS ds, PetscInt f, PetscInt g,
2353: void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2354: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2355: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2356: PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
2357: void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2358: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2359: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2360: PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
2361: void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2362: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2363: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2364: PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
2365: void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2366: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2367: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2368: PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
2369: {
2378: if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2379: if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
2380: PetscWeakFormSetIndexBdJacobian(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3);
2381: return(0);
2382: }
2384: /*@
2385: PetscDSHasBdJacobianPreconditioner - Signals that boundary Jacobian preconditioner functions have been set
2387: Not collective
2389: Input Parameter:
2390: . ds - The PetscDS
2392: Output Parameter:
2393: . hasBdJac - flag that pointwise function for the boundary Jacobian preconditioner has been set
2395: Level: intermediate
2397: .seealso: PetscDSHasJacobian(), PetscDSSetBdJacobian(), PetscDSGetBdJacobian()
2398: @*/
2399: PetscErrorCode PetscDSHasBdJacobianPreconditioner(PetscDS ds, PetscBool *hasBdJacPre)
2400: {
2406: PetscWeakFormHasBdJacobianPreconditioner(ds->wf, hasBdJacPre);
2407: return(0);
2408: }
2410: /*@C
2411: PetscDSGetBdJacobianPreconditioner - Get the pointwise boundary Jacobian preconditioner function for given test and basis field
2413: Not collective
2415: Input Parameters:
2416: + ds - The PetscDS
2417: . f - The test field number
2418: - g - The field number
2420: Output Parameters:
2421: + g0 - integrand for the test and basis function term
2422: . g1 - integrand for the test function and basis function gradient term
2423: . g2 - integrand for the test function gradient and basis function term
2424: - g3 - integrand for the test function gradient and basis function gradient term
2426: Note: We are using a first order FEM model for the weak form:
2428: \int_\Gamma \phi {\vec g}_0(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \cdot \hat n \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \hat n \cdot \nabla \psi
2430: The calling sequence for the callbacks g0, g1, g2 and g3 is given by:
2432: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2433: $ const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2434: $ const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2435: $ PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
2437: + dim - the spatial dimension
2438: . Nf - the number of fields
2439: . NfAux - the number of auxiliary fields
2440: . uOff - the offset into u[] and u_t[] for each field
2441: . uOff_x - the offset into u_x[] for each field
2442: . u - each field evaluated at the current point
2443: . u_t - the time derivative of each field evaluated at the current point
2444: . u_x - the gradient of each field evaluated at the current point
2445: . aOff - the offset into a[] and a_t[] for each auxiliary field
2446: . aOff_x - the offset into a_x[] for each auxiliary field
2447: . a - each auxiliary field evaluated at the current point
2448: . a_t - the time derivative of each auxiliary field evaluated at the current point
2449: . a_x - the gradient of auxiliary each field evaluated at the current point
2450: . t - current time
2451: . u_tShift - the multiplier a for dF/dU_t
2452: . x - coordinates of the current point
2453: . n - normal at the current point
2454: . numConstants - number of constant parameters
2455: . constants - constant parameters
2456: - g0 - output values at the current point
2458: This is not yet available in Fortran.
2460: Level: intermediate
2462: .seealso: PetscDSSetBdJacobianPreconditioner()
2463: @*/
2464: PetscErrorCode PetscDSGetBdJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g,
2465: void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2466: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2467: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2468: PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
2469: void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2470: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2471: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2472: PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
2473: void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2474: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2475: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2476: PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
2477: void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2478: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2479: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2480: PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
2481: {
2482: PetscBdPointJac *tmp0, *tmp1, *tmp2, *tmp3;
2483: PetscInt n0, n1, n2, n3;
2484: PetscErrorCode ierr;
2488: if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
2489: if ((g < 0) || (g >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, ds->Nf);
2490: PetscWeakFormGetBdJacobianPreconditioner(ds->wf, NULL, 0, f, g, 0, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3);
2491: *g0 = tmp0 ? tmp0[0] : NULL;
2492: *g1 = tmp1 ? tmp1[0] : NULL;
2493: *g2 = tmp2 ? tmp2[0] : NULL;
2494: *g3 = tmp3 ? tmp3[0] : NULL;
2495: return(0);
2496: }
2498: /*@C
2499: PetscDSSetBdJacobianPreconditioner - Set the pointwise boundary Jacobian preconditioner function for given test and basis field
2501: Not collective
2503: Input Parameters:
2504: + ds - The PetscDS
2505: . f - The test field number
2506: . g - The field number
2507: . g0 - integrand for the test and basis function term
2508: . g1 - integrand for the test function and basis function gradient term
2509: . g2 - integrand for the test function gradient and basis function term
2510: - g3 - integrand for the test function gradient and basis function gradient term
2512: Note: We are using a first order FEM model for the weak form:
2514: \int_\Gamma \phi {\vec g}_0(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \cdot \hat n \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \hat n \cdot \nabla \psi
2516: The calling sequence for the callbacks g0, g1, g2 and g3 is given by:
2518: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2519: $ const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2520: $ const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2521: $ PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
2523: + dim - the spatial dimension
2524: . Nf - the number of fields
2525: . NfAux - the number of auxiliary fields
2526: . uOff - the offset into u[] and u_t[] for each field
2527: . uOff_x - the offset into u_x[] for each field
2528: . u - each field evaluated at the current point
2529: . u_t - the time derivative of each field evaluated at the current point
2530: . u_x - the gradient of each field evaluated at the current point
2531: . aOff - the offset into a[] and a_t[] for each auxiliary field
2532: . aOff_x - the offset into a_x[] for each auxiliary field
2533: . a - each auxiliary field evaluated at the current point
2534: . a_t - the time derivative of each auxiliary field evaluated at the current point
2535: . a_x - the gradient of auxiliary each field evaluated at the current point
2536: . t - current time
2537: . u_tShift - the multiplier a for dF/dU_t
2538: . x - coordinates of the current point
2539: . n - normal at the current point
2540: . numConstants - number of constant parameters
2541: . constants - constant parameters
2542: - g0 - output values at the current point
2544: This is not yet available in Fortran.
2546: Level: intermediate
2548: .seealso: PetscDSGetBdJacobianPreconditioner()
2549: @*/
2550: PetscErrorCode PetscDSSetBdJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g,
2551: void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2552: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2553: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2554: PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
2555: void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2556: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2557: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2558: PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
2559: void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2560: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2561: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2562: PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
2563: void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2564: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2565: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2566: PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
2567: {
2576: if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2577: if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
2578: PetscWeakFormSetIndexBdJacobianPreconditioner(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3);
2579: return(0);
2580: }
2582: /*@C
2583: PetscDSGetExactSolution - Get the pointwise exact solution function for a given test field
2585: Not collective
2587: Input Parameters:
2588: + prob - The PetscDS
2589: - f - The test field number
2591: Output Parameters:
2592: + exactSol - exact solution for the test field
2593: - exactCtx - exact solution context
2595: Note: The calling sequence for the solution functions is given by:
2597: $ sol(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx)
2599: + dim - the spatial dimension
2600: . t - current time
2601: . x - coordinates of the current point
2602: . Nc - the number of field components
2603: . u - the solution field evaluated at the current point
2604: - ctx - a user context
2606: Level: intermediate
2608: .seealso: PetscDSSetExactSolution(), PetscDSGetExactSolutionTimeDerivative()
2609: @*/
2610: PetscErrorCode PetscDSGetExactSolution(PetscDS prob, PetscInt f, PetscErrorCode (**sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void **ctx)
2611: {
2614: if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
2617: return(0);
2618: }
2620: /*@C
2621: PetscDSSetExactSolution - Set the pointwise exact solution function for a given test field
2623: Not collective
2625: Input Parameters:
2626: + prob - The PetscDS
2627: . f - The test field number
2628: . sol - solution function for the test fields
2629: - ctx - solution context or NULL
2631: Note: The calling sequence for solution functions is given by:
2633: $ sol(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx)
2635: + dim - the spatial dimension
2636: . t - current time
2637: . x - coordinates of the current point
2638: . Nc - the number of field components
2639: . u - the solution field evaluated at the current point
2640: - ctx - a user context
2642: Level: intermediate
2644: .seealso: PetscDSGetExactSolution()
2645: @*/
2646: PetscErrorCode PetscDSSetExactSolution(PetscDS prob, PetscInt f, PetscErrorCode (*sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void *ctx)
2647: {
2652: if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2653: PetscDSEnlarge_Static(prob, f+1);
2656: return(0);
2657: }
2659: /*@C
2660: PetscDSGetExactSolutionTimeDerivative - Get the pointwise time derivative of the exact solution function for a given test field
2662: Not collective
2664: Input Parameters:
2665: + prob - The PetscDS
2666: - f - The test field number
2668: Output Parameters:
2669: + exactSol - time derivative of the exact solution for the test field
2670: - exactCtx - time derivative of the exact solution context
2672: Note: The calling sequence for the solution functions is given by:
2674: $ sol(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx)
2676: + dim - the spatial dimension
2677: . t - current time
2678: . x - coordinates of the current point
2679: . Nc - the number of field components
2680: . u - the solution field evaluated at the current point
2681: - ctx - a user context
2683: Level: intermediate
2685: .seealso: PetscDSSetExactSolutionTimeDerivative(), PetscDSGetExactSolution()
2686: @*/
2687: PetscErrorCode PetscDSGetExactSolutionTimeDerivative(PetscDS prob, PetscInt f, PetscErrorCode (**sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void **ctx)
2688: {
2691: if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
2694: return(0);
2695: }
2697: /*@C
2698: PetscDSSetExactSolutionTimeDerivative - Set the pointwise time derivative of the exact solution function for a given test field
2700: Not collective
2702: Input Parameters:
2703: + prob - The PetscDS
2704: . f - The test field number
2705: . sol - time derivative of the solution function for the test fields
2706: - ctx - time derivative of the solution context or NULL
2708: Note: The calling sequence for solution functions is given by:
2710: $ sol(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx)
2712: + dim - the spatial dimension
2713: . t - current time
2714: . x - coordinates of the current point
2715: . Nc - the number of field components
2716: . u - the solution field evaluated at the current point
2717: - ctx - a user context
2719: Level: intermediate
2721: .seealso: PetscDSGetExactSolutionTimeDerivative(), PetscDSSetExactSolution()
2722: @*/
2723: PetscErrorCode PetscDSSetExactSolutionTimeDerivative(PetscDS prob, PetscInt f, PetscErrorCode (*sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void *ctx)
2724: {
2729: if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
2730: PetscDSEnlarge_Static(prob, f+1);
2733: return(0);
2734: }
2736: /*@C
2737: PetscDSGetConstants - Returns the array of constants passed to point functions
2739: Not collective
2741: Input Parameter:
2742: . prob - The PetscDS object
2744: Output Parameters:
2745: + numConstants - The number of constants
2746: - constants - The array of constants, NULL if there are none
2748: Level: intermediate
2750: .seealso: PetscDSSetConstants(), PetscDSCreate()
2751: @*/
2752: PetscErrorCode PetscDSGetConstants(PetscDS prob, PetscInt *numConstants, const PetscScalar *constants[])
2753: {
2758: return(0);
2759: }
2761: /*@C
2762: PetscDSSetConstants - Set the array of constants passed to point functions
2764: Not collective
2766: Input Parameters:
2767: + prob - The PetscDS object
2768: . numConstants - The number of constants
2769: - constants - The array of constants, NULL if there are none
2771: Level: intermediate
2773: .seealso: PetscDSGetConstants(), PetscDSCreate()
2774: @*/
2775: PetscErrorCode PetscDSSetConstants(PetscDS prob, PetscInt numConstants, PetscScalar constants[])
2776: {
2781: if (numConstants != prob->numConstants) {
2782: PetscFree(prob->constants);
2783: prob->numConstants = numConstants;
2784: if (prob->numConstants) {
2785: PetscMalloc1(prob->numConstants, &prob->constants);
2786: } else {
2787: prob->constants = NULL;
2788: }
2789: }
2790: if (prob->numConstants) {
2792: PetscArraycpy(prob->constants, constants, prob->numConstants);
2793: }
2794: return(0);
2795: }
2797: /*@
2798: PetscDSGetFieldIndex - Returns the index of the given field
2800: Not collective
2802: Input Parameters:
2803: + prob - The PetscDS object
2804: - disc - The discretization object
2806: Output Parameter:
2807: . f - The field number
2809: Level: beginner
2811: .seealso: PetscGetDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
2812: @*/
2813: PetscErrorCode PetscDSGetFieldIndex(PetscDS prob, PetscObject disc, PetscInt *f)
2814: {
2815: PetscInt g;
2820: *f = -1;
2821: for (g = 0; g < prob->Nf; ++g) {if (disc == prob->disc[g]) break;}
2822: if (g == prob->Nf) SETERRQ(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Field not found in PetscDS.");
2823: *f = g;
2824: return(0);
2825: }
2827: /*@
2828: PetscDSGetFieldSize - Returns the size of the given field in the full space basis
2830: Not collective
2832: Input Parameters:
2833: + prob - The PetscDS object
2834: - f - The field number
2836: Output Parameter:
2837: . size - The size
2839: Level: beginner
2841: .seealso: PetscDSGetFieldOffset(), PetscDSGetNumFields(), PetscDSCreate()
2842: @*/
2843: PetscErrorCode PetscDSGetFieldSize(PetscDS prob, PetscInt f, PetscInt *size)
2844: {
2850: if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
2851: PetscDSSetUp(prob);
2852: *size = prob->Nb[f];
2853: return(0);
2854: }
2856: /*@
2857: PetscDSGetFieldOffset - Returns the offset of the given field in the full space basis
2859: Not collective
2861: Input Parameters:
2862: + prob - The PetscDS object
2863: - f - The field number
2865: Output Parameter:
2866: . off - The offset
2868: Level: beginner
2870: .seealso: PetscDSGetFieldSize(), PetscDSGetNumFields(), PetscDSCreate()
2871: @*/
2872: PetscErrorCode PetscDSGetFieldOffset(PetscDS prob, PetscInt f, PetscInt *off)
2873: {
2874: PetscInt size, g;
2880: if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
2881: *off = 0;
2882: for (g = 0; g < f; ++g) {
2883: PetscDSGetFieldSize(prob, g, &size);
2884: *off += size;
2885: }
2886: return(0);
2887: }
2889: /*@
2890: PetscDSGetDimensions - Returns the size of the approximation space for each field on an evaluation point
2892: Not collective
2894: Input Parameter:
2895: . prob - The PetscDS object
2897: Output Parameter:
2898: . dimensions - The number of dimensions
2900: Level: beginner
2902: .seealso: PetscDSGetComponentOffsets(), PetscDSGetNumFields(), PetscDSCreate()
2903: @*/
2904: PetscErrorCode PetscDSGetDimensions(PetscDS prob, PetscInt *dimensions[])
2905: {
2910: PetscDSSetUp(prob);
2912: *dimensions = prob->Nb;
2913: return(0);
2914: }
2916: /*@
2917: PetscDSGetComponents - Returns the number of components for each field on an evaluation point
2919: Not collective
2921: Input Parameter:
2922: . prob - The PetscDS object
2924: Output Parameter:
2925: . components - The number of components
2927: Level: beginner
2929: .seealso: PetscDSGetComponentOffsets(), PetscDSGetNumFields(), PetscDSCreate()
2930: @*/
2931: PetscErrorCode PetscDSGetComponents(PetscDS prob, PetscInt *components[])
2932: {
2937: PetscDSSetUp(prob);
2939: *components = prob->Nc;
2940: return(0);
2941: }
2943: /*@
2944: PetscDSGetComponentOffset - Returns the offset of the given field on an evaluation point
2946: Not collective
2948: Input Parameters:
2949: + prob - The PetscDS object
2950: - f - The field number
2952: Output Parameter:
2953: . off - The offset
2955: Level: beginner
2957: .seealso: PetscDSGetNumFields(), PetscDSCreate()
2958: @*/
2959: PetscErrorCode PetscDSGetComponentOffset(PetscDS prob, PetscInt f, PetscInt *off)
2960: {
2966: if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
2967: PetscDSSetUp(prob);
2968: *off = prob->off[f];
2969: return(0);
2970: }
2972: /*@
2973: PetscDSGetComponentOffsets - Returns the offset of each field on an evaluation point
2975: Not collective
2977: Input Parameter:
2978: . prob - The PetscDS object
2980: Output Parameter:
2981: . offsets - The offsets
2983: Level: beginner
2985: .seealso: PetscDSGetNumFields(), PetscDSCreate()
2986: @*/
2987: PetscErrorCode PetscDSGetComponentOffsets(PetscDS prob, PetscInt *offsets[])
2988: {
2994: PetscDSSetUp(prob);
2995: *offsets = prob->off;
2996: return(0);
2997: }
2999: /*@
3000: PetscDSGetComponentDerivativeOffsets - Returns the offset of each field derivative on an evaluation point
3002: Not collective
3004: Input Parameter:
3005: . prob - The PetscDS object
3007: Output Parameter:
3008: . offsets - The offsets
3010: Level: beginner
3012: .seealso: PetscDSGetNumFields(), PetscDSCreate()
3013: @*/
3014: PetscErrorCode PetscDSGetComponentDerivativeOffsets(PetscDS prob, PetscInt *offsets[])
3015: {
3021: PetscDSSetUp(prob);
3022: *offsets = prob->offDer;
3023: return(0);
3024: }
3026: /*@C
3027: PetscDSGetTabulation - Return the basis tabulation at quadrature points for the volume discretization
3029: Not collective
3031: Input Parameter:
3032: . prob - The PetscDS object
3034: Output Parameter:
3035: . T - The basis function and derivatives tabulation at quadrature points for each field
3037: Level: intermediate
3039: .seealso: PetscDSCreate()
3040: @*/
3041: PetscErrorCode PetscDSGetTabulation(PetscDS prob, PetscTabulation *T[])
3042: {
3048: PetscDSSetUp(prob);
3049: *T = prob->T;
3050: return(0);
3051: }
3053: /*@C
3054: PetscDSGetFaceTabulation - Return the basis tabulation at quadrature points on the faces
3056: Not collective
3058: Input Parameter:
3059: . prob - The PetscDS object
3061: Output Parameter:
3062: . Tf - The basis function and derivative tabulation on each local face at quadrature points for each and field
3064: Level: intermediate
3066: .seealso: PetscDSGetTabulation(), PetscDSCreate()
3067: @*/
3068: PetscErrorCode PetscDSGetFaceTabulation(PetscDS prob, PetscTabulation *Tf[])
3069: {
3075: PetscDSSetUp(prob);
3076: *Tf = prob->Tf;
3077: return(0);
3078: }
3080: PetscErrorCode PetscDSGetEvaluationArrays(PetscDS prob, PetscScalar **u, PetscScalar **u_t, PetscScalar **u_x)
3081: {
3086: PetscDSSetUp(prob);
3090: return(0);
3091: }
3093: PetscErrorCode PetscDSGetWeakFormArrays(PetscDS prob, PetscScalar **f0, PetscScalar **f1, PetscScalar **g0, PetscScalar **g1, PetscScalar **g2, PetscScalar **g3)
3094: {
3099: PetscDSSetUp(prob);
3106: return(0);
3107: }
3109: PetscErrorCode PetscDSGetWorkspace(PetscDS prob, PetscReal **x, PetscScalar **basisReal, PetscScalar **basisDerReal, PetscScalar **testReal, PetscScalar **testDerReal)
3110: {
3115: PetscDSSetUp(prob);
3121: return(0);
3122: }
3124: /*@C
3125: PetscDSAddBoundary - Add a boundary condition to the model. The pointwise functions are used to provide boundary values for essential boundary conditions. In FEM, they are acting upon by dual basis functionals to generate FEM coefficients which are fixed. Natural boundary conditions signal to PETSc that boundary integrals should be performaed, using the kernels from PetscDSSetBdResidual().
3127: Collective on ds
3129: Input Parameters:
3130: + ds - The PetscDS object
3131: . type - The type of condition, e.g. DM_BC_ESSENTIAL/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
3132: . name - The BC name
3133: . label - The label defining constrained points
3134: . Nv - The number of DMLabel values for constrained points
3135: . values - An array of label values for constrained points
3136: . field - The field to constrain
3137: . Nc - The number of constrained field components (0 will constrain all fields)
3138: . comps - An array of constrained component numbers
3139: . bcFunc - A pointwise function giving boundary values
3140: . bcFunc_t - A pointwise function giving the time derivative of the boundary values, or NULL
3141: - ctx - An optional user context for bcFunc
3143: Output Parameters:
3144: - bd - The boundary number
3146: Options Database Keys:
3147: + -bc_<boundary name> <num> - Overrides the boundary ids
3148: - -bc_<boundary name>_comp <num> - Overrides the boundary components
3150: Note:
3151: Both bcFunc abd bcFunc_t will depend on the boundary condition type. If the type if DM_BC_ESSENTIAL, Then the calling sequence is:
3153: $ bcFunc(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar bcval[])
3155: If the type is DM_BC_ESSENTIAL_FIELD or other _FIELD value, then the calling sequence is:
3157: $ bcFunc(PetscInt dim, PetscInt Nf, PetscInt NfAux,
3158: $ const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
3159: $ const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
3160: $ PetscReal time, const PetscReal x[], PetscScalar bcval[])
3162: + dim - the spatial dimension
3163: . Nf - the number of fields
3164: . uOff - the offset into u[] and u_t[] for each field
3165: . uOff_x - the offset into u_x[] for each field
3166: . u - each field evaluated at the current point
3167: . u_t - the time derivative of each field evaluated at the current point
3168: . u_x - the gradient of each field evaluated at the current point
3169: . aOff - the offset into a[] and a_t[] for each auxiliary field
3170: . aOff_x - the offset into a_x[] for each auxiliary field
3171: . a - each auxiliary field evaluated at the current point
3172: . a_t - the time derivative of each auxiliary field evaluated at the current point
3173: . a_x - the gradient of auxiliary each field evaluated at the current point
3174: . t - current time
3175: . x - coordinates of the current point
3176: . numConstants - number of constant parameters
3177: . constants - constant parameters
3178: - bcval - output values at the current point
3180: Level: developer
3182: .seealso: PetscDSAddBoundaryByName(), PetscDSGetBoundary(), PetscDSSetResidual(), PetscDSSetBdResidual()
3183: @*/
3184: PetscErrorCode PetscDSAddBoundary(PetscDS ds, DMBoundaryConditionType type, const char name[], DMLabel label, PetscInt Nv, const PetscInt values[], PetscInt field, PetscInt Nc, const PetscInt comps[], void (*bcFunc)(void), void (*bcFunc_t)(void), void *ctx, PetscInt *bd)
3185: {
3186: DSBoundary head = ds->boundary, b;
3187: PetscInt n = 0;
3188: const char *lname;
3199: if (Nc > 0) {
3200: PetscInt *fcomps;
3201: PetscInt c;
3203: PetscDSGetComponents(ds, &fcomps);
3204: if (Nc > fcomps[field]) SETERRQ3(PetscObjectComm((PetscObject) ds), PETSC_ERR_ARG_OUTOFRANGE, "Number of constrained components %D > %D components for field %D", Nc, fcomps[field], field);
3205: for (c = 0; c < Nc; ++c) {
3206: if (comps[c] < 0 || comps[c] >= fcomps[field]) SETERRQ4(PetscObjectComm((PetscObject) ds), PETSC_ERR_ARG_OUTOFRANGE, "Constrained component[%D] %D not in [0, %D) components for field %D", c, comps[c], fcomps[field], field);
3207: }
3208: }
3209: PetscNew(&b);
3210: PetscStrallocpy(name, (char **) &b->name);
3211: PetscWeakFormCreate(PETSC_COMM_SELF, &b->wf);
3212: PetscWeakFormSetNumFields(b->wf, ds->Nf);
3213: PetscMalloc1(Nv, &b->values);
3214: if (Nv) {PetscArraycpy(b->values, values, Nv);}
3215: PetscMalloc1(Nc, &b->comps);
3216: if (Nc) {PetscArraycpy(b->comps, comps, Nc);}
3217: PetscObjectGetName((PetscObject) label, &lname);
3218: PetscStrallocpy(lname, (char **) &b->lname);
3219: b->type = type;
3220: b->label = label;
3221: b->Nv = Nv;
3222: b->field = field;
3223: b->Nc = Nc;
3224: b->func = bcFunc;
3225: b->func_t = bcFunc_t;
3226: b->ctx = ctx;
3227: b->next = NULL;
3228: /* Append to linked list so that we can preserve the order */
3229: if (!head) ds->boundary = b;
3230: while (head) {
3231: if (!head->next) {
3232: head->next = b;
3233: head = b;
3234: }
3235: head = head->next;
3236: ++n;
3237: }
3239: return(0);
3240: }
3242: /*@C
3243: PetscDSAddBoundaryByName - Add a boundary condition to the model. The pointwise functions are used to provide boundary values for essential boundary conditions. In FEM, they are acting upon by dual basis functionals to generate FEM coefficients which are fixed. Natural boundary conditions signal to PETSc that boundary integrals should be performaed, using the kernels from PetscDSSetBdResidual().
3245: Collective on ds
3247: Input Parameters:
3248: + ds - The PetscDS object
3249: . type - The type of condition, e.g. DM_BC_ESSENTIAL/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
3250: . name - The BC name
3251: . lname - The naem of the label defining constrained points
3252: . Nv - The number of DMLabel values for constrained points
3253: . values - An array of label values for constrained points
3254: . field - The field to constrain
3255: . Nc - The number of constrained field components (0 will constrain all fields)
3256: . comps - An array of constrained component numbers
3257: . bcFunc - A pointwise function giving boundary values
3258: . bcFunc_t - A pointwise function giving the time derivative of the boundary values, or NULL
3259: - ctx - An optional user context for bcFunc
3261: Output Parameters:
3262: - bd - The boundary number
3264: Options Database Keys:
3265: + -bc_<boundary name> <num> - Overrides the boundary ids
3266: - -bc_<boundary name>_comp <num> - Overrides the boundary components
3268: Note:
3269: This function should only be used with DMForest currently, since labels cannot be defined before the underlygin Plex is built.
3271: Both bcFunc abd bcFunc_t will depend on the boundary condition type. If the type if DM_BC_ESSENTIAL, Then the calling sequence is:
3273: $ bcFunc(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar bcval[])
3275: If the type is DM_BC_ESSENTIAL_FIELD or other _FIELD value, then the calling sequence is:
3277: $ bcFunc(PetscInt dim, PetscInt Nf, PetscInt NfAux,
3278: $ const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
3279: $ const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
3280: $ PetscReal time, const PetscReal x[], PetscScalar bcval[])
3282: + dim - the spatial dimension
3283: . Nf - the number of fields
3284: . uOff - the offset into u[] and u_t[] for each field
3285: . uOff_x - the offset into u_x[] for each field
3286: . u - each field evaluated at the current point
3287: . u_t - the time derivative of each field evaluated at the current point
3288: . u_x - the gradient of each field evaluated at the current point
3289: . aOff - the offset into a[] and a_t[] for each auxiliary field
3290: . aOff_x - the offset into a_x[] for each auxiliary field
3291: . a - each auxiliary field evaluated at the current point
3292: . a_t - the time derivative of each auxiliary field evaluated at the current point
3293: . a_x - the gradient of auxiliary each field evaluated at the current point
3294: . t - current time
3295: . x - coordinates of the current point
3296: . numConstants - number of constant parameters
3297: . constants - constant parameters
3298: - bcval - output values at the current point
3300: Level: developer
3302: .seealso: PetscDSAddBoundary(), PetscDSGetBoundary(), PetscDSSetResidual(), PetscDSSetBdResidual()
3303: @*/
3304: PetscErrorCode PetscDSAddBoundaryByName(PetscDS ds, DMBoundaryConditionType type, const char name[], const char lname[], PetscInt Nv, const PetscInt values[], PetscInt field, PetscInt Nc, const PetscInt comps[], void (*bcFunc)(void), void (*bcFunc_t)(void), void *ctx, PetscInt *bd)
3305: {
3306: DSBoundary head = ds->boundary, b;
3307: PetscInt n = 0;
3318: PetscNew(&b);
3319: PetscStrallocpy(name, (char **) &b->name);
3320: PetscWeakFormCreate(PETSC_COMM_SELF, &b->wf);
3321: PetscWeakFormSetNumFields(b->wf, ds->Nf);
3322: PetscMalloc1(Nv, &b->values);
3323: if (Nv) {PetscArraycpy(b->values, values, Nv);}
3324: PetscMalloc1(Nc, &b->comps);
3325: if (Nc) {PetscArraycpy(b->comps, comps, Nc);}
3326: PetscStrallocpy(lname, (char **) &b->lname);
3327: b->type = type;
3328: b->label = NULL;
3329: b->Nv = Nv;
3330: b->field = field;
3331: b->Nc = Nc;
3332: b->func = bcFunc;
3333: b->func_t = bcFunc_t;
3334: b->ctx = ctx;
3335: b->next = NULL;
3336: /* Append to linked list so that we can preserve the order */
3337: if (!head) ds->boundary = b;
3338: while (head) {
3339: if (!head->next) {
3340: head->next = b;
3341: head = b;
3342: }
3343: head = head->next;
3344: ++n;
3345: }
3347: return(0);
3348: }
3350: /*@C
3351: PetscDSUpdateBoundary - Change a boundary condition for the model. The pointwise functions are used to provide boundary values for essential boundary conditions. In FEM, they are acting upon by dual basis functionals to generate FEM coefficients which are fixed. Natural boundary conditions signal to PETSc that boundary integrals should be performaed, using the kernels from PetscDSSetBdResidual().
3353: Input Parameters:
3354: + ds - The PetscDS object
3355: . bd - The boundary condition number
3356: . type - The type of condition, e.g. DM_BC_ESSENTIAL/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
3357: . name - The BC name
3358: . label - The label defining constrained points
3359: . Nv - The number of DMLabel ids for constrained points
3360: . values - An array of ids for constrained points
3361: . field - The field to constrain
3362: . Nc - The number of constrained field components
3363: . comps - An array of constrained component numbers
3364: . bcFunc - A pointwise function giving boundary values
3365: . bcFunc_t - A pointwise function giving the time derivative of the boundary values, or NULL
3366: - ctx - An optional user context for bcFunc
3368: Note:
3369: The boundary condition number is the order in which it was registered. The user can get the number of boundary conditions from PetscDSGetNumBoundary(). See PetscDSAddBoundary() for a description of the calling sequences for the callbacks.
3371: Level: developer
3373: .seealso: PetscDSAddBoundary(), PetscDSGetBoundary(), PetscDSGetNumBoundary()
3374: @*/
3375: PetscErrorCode PetscDSUpdateBoundary(PetscDS ds, PetscInt bd, DMBoundaryConditionType type, const char name[], DMLabel label, PetscInt Nv, const PetscInt values[], PetscInt field, PetscInt Nc, const PetscInt comps[], void (*bcFunc)(void), void (*bcFunc_t)(void), void *ctx)
3376: {
3377: DSBoundary b = ds->boundary;
3378: PetscInt n = 0;
3383: while (b) {
3384: if (n == bd) break;
3385: b = b->next;
3386: ++n;
3387: }
3388: if (!b) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %d is not in [0, %d)", bd, n);
3389: if (name) {
3390: PetscFree(b->name);
3391: PetscStrallocpy(name, (char **) &b->name);
3392: }
3393: b->type = type;
3394: if (label) {
3395: const char *name;
3397: b->label = label;
3398: PetscFree(b->lname);
3399: PetscObjectGetName((PetscObject) label, &name);
3400: PetscStrallocpy(name, (char **) &b->lname);
3401: }
3402: if (Nv >= 0) {
3403: b->Nv = Nv;
3404: PetscFree(b->values);
3405: PetscMalloc1(Nv, &b->values);
3406: if (Nv) {PetscArraycpy(b->values, values, Nv);}
3407: }
3408: if (field >= 0) b->field = field;
3409: if (Nc >= 0) {
3410: b->Nc = Nc;
3411: PetscFree(b->comps);
3412: PetscMalloc1(Nc, &b->comps);
3413: if (Nc) {PetscArraycpy(b->comps, comps, Nc);}
3414: }
3415: if (bcFunc) b->func = bcFunc;
3416: if (bcFunc_t) b->func_t = bcFunc_t;
3417: if (ctx) b->ctx = ctx;
3418: return(0);
3419: }
3421: /*@
3422: PetscDSGetNumBoundary - Get the number of registered BC
3424: Input Parameters:
3425: . ds - The PetscDS object
3427: Output Parameters:
3428: . numBd - The number of BC
3430: Level: intermediate
3432: .seealso: PetscDSAddBoundary(), PetscDSGetBoundary()
3433: @*/
3434: PetscErrorCode PetscDSGetNumBoundary(PetscDS ds, PetscInt *numBd)
3435: {
3436: DSBoundary b = ds->boundary;
3441: *numBd = 0;
3442: while (b) {++(*numBd); b = b->next;}
3443: return(0);
3444: }
3446: /*@C
3447: PetscDSGetBoundary - Gets a boundary condition to the model
3449: Input Parameters:
3450: + ds - The PetscDS object
3451: - bd - The BC number
3453: Output Parameters:
3454: + wf - The PetscWeakForm holding the pointwise functions
3455: . type - The type of condition, e.g. DM_BC_ESSENTIAL/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
3456: . name - The BC name
3457: . label - The label defining constrained points
3458: . Nv - The number of DMLabel ids for constrained points
3459: . values - An array of ids for constrained points
3460: . field - The field to constrain
3461: . Nc - The number of constrained field components
3462: . comps - An array of constrained component numbers
3463: . bcFunc - A pointwise function giving boundary values
3464: . bcFunc_t - A pointwise function giving the time derivative of the boundary values
3465: - ctx - An optional user context for bcFunc
3467: Options Database Keys:
3468: + -bc_<boundary name> <num> - Overrides the boundary ids
3469: - -bc_<boundary name>_comp <num> - Overrides the boundary components
3471: Level: developer
3473: .seealso: PetscDSAddBoundary()
3474: @*/
3475: PetscErrorCode PetscDSGetBoundary(PetscDS ds, PetscInt bd, PetscWeakForm *wf, DMBoundaryConditionType *type, const char *name[], DMLabel *label, PetscInt *Nv, const PetscInt *values[], PetscInt *field, PetscInt *Nc, const PetscInt *comps[], void (**func)(void), void (**func_t)(void), void **ctx)
3476: {
3477: DSBoundary b = ds->boundary;
3478: PetscInt n = 0;
3482: while (b) {
3483: if (n == bd) break;
3484: b = b->next;
3485: ++n;
3486: }
3487: if (!b) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %d is not in [0, %d)", bd, n);
3488: if (wf) {
3490: *wf = b->wf;
3491: }
3492: if (type) {
3494: *type = b->type;
3495: }
3496: if (name) {
3498: *name = b->name;
3499: }
3500: if (label) {
3502: *label = b->label;
3503: }
3504: if (Nv) {
3506: *Nv = b->Nv;
3507: }
3508: if (values) {
3510: *values = b->values;
3511: }
3512: if (field) {
3514: *field = b->field;
3515: }
3516: if (Nc) {
3518: *Nc = b->Nc;
3519: }
3520: if (comps) {
3522: *comps = b->comps;
3523: }
3524: if (func) {
3526: *func = b->func;
3527: }
3528: if (func_t) {
3530: *func_t = b->func_t;
3531: }
3532: if (ctx) {
3534: *ctx = b->ctx;
3535: }
3536: return(0);
3537: }
3539: static PetscErrorCode DSBoundaryDuplicate_Internal(DSBoundary b, DSBoundary *bNew)
3540: {
3544: PetscNew(bNew);
3545: PetscWeakFormCreate(PETSC_COMM_SELF, &(*bNew)->wf);
3546: PetscWeakFormCopy(b->wf, (*bNew)->wf);
3547: PetscStrallocpy(b->name,(char **) &((*bNew)->name));
3548: PetscStrallocpy(b->lname,(char **) &((*bNew)->lname));
3549: (*bNew)->type = b->type;
3550: (*bNew)->label = b->label;
3551: (*bNew)->Nv = b->Nv;
3552: PetscMalloc1(b->Nv, &(*bNew)->values);
3553: PetscArraycpy((*bNew)->values, b->values, b->Nv);
3554: (*bNew)->field = b->field;
3555: (*bNew)->Nc = b->Nc;
3556: PetscMalloc1(b->Nc, &(*bNew)->comps);
3557: PetscArraycpy((*bNew)->comps, b->comps, b->Nc);
3558: (*bNew)->func = b->func;
3559: (*bNew)->func_t = b->func_t;
3560: (*bNew)->ctx = b->ctx;
3561: return(0);
3562: }
3564: /*@
3565: PetscDSCopyBoundary - Copy all boundary condition objects to the new problem
3567: Not collective
3569: Input Parameters:
3570: + ds - The source PetscDS object
3571: . numFields - The number of selected fields, or PETSC_DEFAULT for all fields
3572: - fields - The selected fields, or NULL for all fields
3574: Output Parameter:
3575: . newds - The target PetscDS, now with a copy of the boundary conditions
3577: Level: intermediate
3579: .seealso: PetscDSCopyEquations(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3580: @*/
3581: PetscErrorCode PetscDSCopyBoundary(PetscDS ds, PetscInt numFields, const PetscInt fields[], PetscDS newds)
3582: {
3583: DSBoundary b, *lastnext;
3589: if (ds == newds) return(0);
3590: PetscDSDestroyBoundary(newds);
3591: lastnext = &(newds->boundary);
3592: for (b = ds->boundary; b; b = b->next) {
3593: DSBoundary bNew;
3594: PetscInt fieldNew = -1;
3596: if (numFields > 0 && fields) {
3597: PetscInt f;
3599: for (f = 0; f < numFields; ++f) if (b->field == fields[f]) break;
3600: if (f == numFields) continue;
3601: fieldNew = f;
3602: }
3603: DSBoundaryDuplicate_Internal(b, &bNew);
3604: bNew->field = fieldNew < 0 ? b->field : fieldNew;
3605: *lastnext = bNew;
3606: lastnext = &(bNew->next);
3607: }
3608: return(0);
3609: }
3611: /*@
3612: PetscDSDestroyBoundary - Remove all DMBoundary objects from the PetscDS
3614: Not collective
3616: Input Parameter:
3617: . ds - The PetscDS object
3619: Level: intermediate
3621: .seealso: PetscDSCopyBoundary(), PetscDSCopyEquations()
3622: @*/
3623: PetscErrorCode PetscDSDestroyBoundary(PetscDS ds)
3624: {
3625: DSBoundary next = ds->boundary;
3629: while (next) {
3630: DSBoundary b = next;
3632: next = b->next;
3633: PetscWeakFormDestroy(&b->wf);
3634: PetscFree(b->name);
3635: PetscFree(b->lname);
3636: PetscFree(b->values);
3637: PetscFree(b->comps);
3638: PetscFree(b);
3639: }
3640: return(0);
3641: }
3643: /*@
3644: PetscDSSelectDiscretizations - Copy discretizations to the new problem with different field layout
3646: Not collective
3648: Input Parameters:
3649: + prob - The PetscDS object
3650: . numFields - Number of new fields
3651: - fields - Old field number for each new field
3653: Output Parameter:
3654: . newprob - The PetscDS copy
3656: Level: intermediate
3658: .seealso: PetscDSSelectEquations(), PetscDSCopyBoundary(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3659: @*/
3660: PetscErrorCode PetscDSSelectDiscretizations(PetscDS prob, PetscInt numFields, const PetscInt fields[], PetscDS newprob)
3661: {
3662: PetscInt Nf, Nfn, fn;
3669: PetscDSGetNumFields(prob, &Nf);
3670: PetscDSGetNumFields(newprob, &Nfn);
3671: numFields = numFields < 0 ? Nf : numFields;
3672: for (fn = 0; fn < numFields; ++fn) {
3673: const PetscInt f = fields ? fields[fn] : fn;
3674: PetscObject disc;
3676: if (f >= Nf) continue;
3677: PetscDSGetDiscretization(prob, f, &disc);
3678: PetscDSSetDiscretization(newprob, fn, disc);
3679: }
3680: return(0);
3681: }
3683: /*@
3684: PetscDSSelectEquations - Copy pointwise function pointers to the new problem with different field layout
3686: Not collective
3688: Input Parameters:
3689: + prob - The PetscDS object
3690: . numFields - Number of new fields
3691: - fields - Old field number for each new field
3693: Output Parameter:
3694: . newprob - The PetscDS copy
3696: Level: intermediate
3698: .seealso: PetscDSSelectDiscretizations(), PetscDSCopyBoundary(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3699: @*/
3700: PetscErrorCode PetscDSSelectEquations(PetscDS prob, PetscInt numFields, const PetscInt fields[], PetscDS newprob)
3701: {
3702: PetscInt Nf, Nfn, fn, gn;
3709: PetscDSGetNumFields(prob, &Nf);
3710: PetscDSGetNumFields(newprob, &Nfn);
3711: if (numFields > Nfn) SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_SIZ, "Number of fields %D to transfer must not be greater then the total number of fields %D", numFields, Nfn);
3712: for (fn = 0; fn < numFields; ++fn) {
3713: const PetscInt f = fields ? fields[fn] : fn;
3714: PetscPointFunc obj;
3715: PetscPointFunc f0, f1;
3716: PetscBdPointFunc f0Bd, f1Bd;
3717: PetscRiemannFunc r;
3719: if (f >= Nf) continue;
3720: PetscDSGetObjective(prob, f, &obj);
3721: PetscDSGetResidual(prob, f, &f0, &f1);
3722: PetscDSGetBdResidual(prob, f, &f0Bd, &f1Bd);
3723: PetscDSGetRiemannSolver(prob, f, &r);
3724: PetscDSSetObjective(newprob, fn, obj);
3725: PetscDSSetResidual(newprob, fn, f0, f1);
3726: PetscDSSetBdResidual(newprob, fn, f0Bd, f1Bd);
3727: PetscDSSetRiemannSolver(newprob, fn, r);
3728: for (gn = 0; gn < numFields; ++gn) {
3729: const PetscInt g = fields ? fields[gn] : gn;
3730: PetscPointJac g0, g1, g2, g3;
3731: PetscPointJac g0p, g1p, g2p, g3p;
3732: PetscBdPointJac g0Bd, g1Bd, g2Bd, g3Bd;
3734: if (g >= Nf) continue;
3735: PetscDSGetJacobian(prob, f, g, &g0, &g1, &g2, &g3);
3736: PetscDSGetJacobianPreconditioner(prob, f, g, &g0p, &g1p, &g2p, &g3p);
3737: PetscDSGetBdJacobian(prob, f, g, &g0Bd, &g1Bd, &g2Bd, &g3Bd);
3738: PetscDSSetJacobian(newprob, fn, gn, g0, g1, g2, g3);
3739: PetscDSSetJacobianPreconditioner(newprob, fn, gn, g0p, g1p, g2p, g3p);
3740: PetscDSSetBdJacobian(newprob, fn, gn, g0Bd, g1Bd, g2Bd, g3Bd);
3741: }
3742: }
3743: return(0);
3744: }
3746: /*@
3747: PetscDSCopyEquations - Copy all pointwise function pointers to the new problem
3749: Not collective
3751: Input Parameter:
3752: . prob - The PetscDS object
3754: Output Parameter:
3755: . newprob - The PetscDS copy
3757: Level: intermediate
3759: .seealso: PetscDSCopyBoundary(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3760: @*/
3761: PetscErrorCode PetscDSCopyEquations(PetscDS prob, PetscDS newprob)
3762: {
3763: PetscWeakForm wf, newwf;
3764: PetscInt Nf, Ng;
3770: PetscDSGetNumFields(prob, &Nf);
3771: PetscDSGetNumFields(newprob, &Ng);
3772: if (Nf != Ng) SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_SIZ, "Number of fields must match %D != %D", Nf, Ng);
3773: PetscDSGetWeakForm(prob, &wf);
3774: PetscDSGetWeakForm(newprob, &newwf);
3775: PetscWeakFormCopy(wf, newwf);
3776: return(0);
3777: }
3779: /*@
3780: PetscDSCopyConstants - Copy all constants to the new problem
3782: Not collective
3784: Input Parameter:
3785: . prob - The PetscDS object
3787: Output Parameter:
3788: . newprob - The PetscDS copy
3790: Level: intermediate
3792: .seealso: PetscDSCopyBoundary(), PetscDSCopyEquations(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3793: @*/
3794: PetscErrorCode PetscDSCopyConstants(PetscDS prob, PetscDS newprob)
3795: {
3796: PetscInt Nc;
3797: const PetscScalar *constants;
3798: PetscErrorCode ierr;
3803: PetscDSGetConstants(prob, &Nc, &constants);
3804: PetscDSSetConstants(newprob, Nc, (PetscScalar *) constants);
3805: return(0);
3806: }
3808: /*@
3809: PetscDSCopyExactSolutions - Copy all exact solutions to the new problem
3811: Not collective
3813: Input Parameter:
3814: . ds - The PetscDS object
3816: Output Parameter:
3817: . newds - The PetscDS copy
3819: Level: intermediate
3821: .seealso: PetscDSCopyBoundary(), PetscDSCopyEquations(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
3822: @*/
3823: PetscErrorCode PetscDSCopyExactSolutions(PetscDS ds, PetscDS newds)
3824: {
3825: PetscSimplePointFunc sol;
3826: void *ctx;
3827: PetscInt Nf, f;
3828: PetscErrorCode ierr;
3833: PetscDSGetNumFields(ds, &Nf);
3834: for (f = 0; f < Nf; ++f) {
3835: PetscDSGetExactSolution(ds, f, &sol, &ctx);
3836: PetscDSSetExactSolution(newds, f, sol, ctx);
3837: PetscDSGetExactSolutionTimeDerivative(ds, f, &sol, &ctx);
3838: PetscDSSetExactSolutionTimeDerivative(newds, f, sol, ctx);
3839: }
3840: return(0);
3841: }
3843: PetscErrorCode PetscDSGetHeightSubspace(PetscDS prob, PetscInt height, PetscDS *subprob)
3844: {
3845: PetscInt dim, Nf, f;
3851: if (height == 0) {*subprob = prob; return(0);}
3852: PetscDSGetNumFields(prob, &Nf);
3853: PetscDSGetSpatialDimension(prob, &dim);
3854: if (height > dim) SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_OUTOFRANGE, "DS can only handle height in [0, %D], not %D", dim, height);
3855: if (!prob->subprobs) {PetscCalloc1(dim, &prob->subprobs);}
3856: if (!prob->subprobs[height-1]) {
3857: PetscInt cdim;
3859: PetscDSCreate(PetscObjectComm((PetscObject) prob), &prob->subprobs[height-1]);
3860: PetscDSGetCoordinateDimension(prob, &cdim);
3861: PetscDSSetCoordinateDimension(prob->subprobs[height-1], cdim);
3862: for (f = 0; f < Nf; ++f) {
3863: PetscFE subfe;
3864: PetscObject obj;
3865: PetscClassId id;
3867: PetscDSGetDiscretization(prob, f, &obj);
3868: PetscObjectGetClassId(obj, &id);
3869: if (id == PETSCFE_CLASSID) {PetscFEGetHeightSubspace((PetscFE) obj, height, &subfe);}
3870: else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unsupported discretization type for field %d", f);
3871: PetscDSSetDiscretization(prob->subprobs[height-1], f, (PetscObject) subfe);
3872: }
3873: }
3874: *subprob = prob->subprobs[height-1];
3875: return(0);
3876: }
3878: PetscErrorCode PetscDSGetDiscType_Internal(PetscDS ds, PetscInt f, PetscDiscType *disctype)
3879: {
3880: PetscObject obj;
3881: PetscClassId id;
3882: PetscInt Nf;
3888: *disctype = PETSC_DISC_NONE;
3889: PetscDSGetNumFields(ds, &Nf);
3890: if (f >= Nf) SETERRQ2(PetscObjectComm((PetscObject) ds), PETSC_ERR_ARG_SIZ, "Field %D must be in [0, %D)", f, Nf);
3891: PetscDSGetDiscretization(ds, f, &obj);
3892: if (obj) {
3893: PetscObjectGetClassId(obj, &id);
3894: if (id == PETSCFE_CLASSID) *disctype = PETSC_DISC_FE;
3895: else *disctype = PETSC_DISC_FV;
3896: }
3897: return(0);
3898: }
3900: static PetscErrorCode PetscDSDestroy_Basic(PetscDS ds)
3901: {
3905: PetscFree(ds->data);
3906: return(0);
3907: }
3909: static PetscErrorCode PetscDSInitialize_Basic(PetscDS ds)
3910: {
3912: ds->ops->setfromoptions = NULL;
3913: ds->ops->setup = NULL;
3914: ds->ops->view = NULL;
3915: ds->ops->destroy = PetscDSDestroy_Basic;
3916: return(0);
3917: }
3919: /*MC
3920: PETSCDSBASIC = "basic" - A discrete system with pointwise residual and boundary residual functions
3922: Level: intermediate
3924: .seealso: PetscDSType, PetscDSCreate(), PetscDSSetType()
3925: M*/
3927: PETSC_EXTERN PetscErrorCode PetscDSCreate_Basic(PetscDS ds)
3928: {
3929: PetscDS_Basic *b;
3934: PetscNewLog(ds, &b);
3935: ds->data = b;
3937: PetscDSInitialize_Basic(ds);
3938: return(0);
3939: }