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; No Fortran Support
26: Input Parameters:
27: + sname - The name of a new user-defined creation routine
28: - function - The creation routine itself
30: Sample usage:
31: .vb
32: PetscDSRegister("my_ds", MyPetscDSCreate);
33: .ve
35: Then, your PetscDS type can be chosen with the procedural interface via
36: .vb
37: PetscDSCreate(MPI_Comm, PetscDS *);
38: PetscDSSetType(PetscDS, "my_ds");
39: .ve
40: or at runtime via the option
41: .vb
42: -petscds_type my_ds
43: .ve
45: Level: advanced
47: Note:
48: `PetscDSRegister()` may be called multiple times to add several user-defined `PetscDSs`
50: .seealso: `PetscDSType`, `PetscDS`, `PetscDSRegisterAll()`, `PetscDSRegisterDestroy()`
51: @*/
52: PetscErrorCode PetscDSRegister(const char sname[], PetscErrorCode (*function)(PetscDS))
53: {
54: PetscFunctionBegin;
55: PetscCall(PetscFunctionListAdd(&PetscDSList, sname, function));
56: PetscFunctionReturn(PETSC_SUCCESS);
57: }
59: /*@C
60: PetscDSSetType - Builds a particular `PetscDS`
62: Collective; No Fortran Support
64: Input Parameters:
65: + prob - The `PetscDS` object
66: - name - The `PetscDSType`
68: Options Database Key:
69: . -petscds_type <type> - Sets the PetscDS type; use -help for a list of available types
71: Level: intermediate
73: .seealso: `PetscDSType`, `PetscDS`, `PetscDSGetType()`, `PetscDSCreate()`
74: @*/
75: PetscErrorCode PetscDSSetType(PetscDS prob, PetscDSType name)
76: {
77: PetscErrorCode (*r)(PetscDS);
78: PetscBool match;
80: PetscFunctionBegin;
82: PetscCall(PetscObjectTypeCompare((PetscObject)prob, name, &match));
83: if (match) PetscFunctionReturn(PETSC_SUCCESS);
85: PetscCall(PetscDSRegisterAll());
86: PetscCall(PetscFunctionListFind(PetscDSList, name, &r));
87: PetscCheck(r, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscDS type: %s", name);
89: PetscTryTypeMethod(prob, destroy);
90: prob->ops->destroy = NULL;
92: PetscCall((*r)(prob));
93: PetscCall(PetscObjectChangeTypeName((PetscObject)prob, name));
94: PetscFunctionReturn(PETSC_SUCCESS);
95: }
97: /*@C
98: PetscDSGetType - Gets the `PetscDSType` name (as a string) from the `PetscDS`
100: Not Collective; No Fortran Support
102: Input Parameter:
103: . prob - The `PetscDS`
105: Output Parameter:
106: . name - The `PetscDSType` name
108: Level: intermediate
110: .seealso: `PetscDSType`, `PetscDS`, `PetscDSSetType()`, `PetscDSCreate()`
111: @*/
112: PetscErrorCode PetscDSGetType(PetscDS prob, PetscDSType *name)
113: {
114: PetscFunctionBegin;
117: PetscCall(PetscDSRegisterAll());
118: *name = ((PetscObject)prob)->type_name;
119: PetscFunctionReturn(PETSC_SUCCESS);
120: }
122: static PetscErrorCode PetscDSView_Ascii(PetscDS ds, PetscViewer viewer)
123: {
124: PetscViewerFormat format;
125: const PetscScalar *constants;
126: PetscInt Nf, numConstants, f;
128: PetscFunctionBegin;
129: PetscCall(PetscDSGetNumFields(ds, &Nf));
130: PetscCall(PetscViewerGetFormat(viewer, &format));
131: PetscCall(PetscViewerASCIIPrintf(viewer, "Discrete System with %" PetscInt_FMT " fields\n", Nf));
132: PetscCall(PetscViewerASCIIPushTab(viewer));
133: PetscCall(PetscViewerASCIIPrintf(viewer, " cell total dim %" PetscInt_FMT " total comp %" PetscInt_FMT "\n", ds->totDim, ds->totComp));
134: if (ds->isCohesive) PetscCall(PetscViewerASCIIPrintf(viewer, " cohesive cell\n"));
135: for (f = 0; f < Nf; ++f) {
136: DSBoundary b;
137: PetscObject obj;
138: PetscClassId id;
139: PetscQuadrature q;
140: const char *name;
141: PetscInt Nc, Nq, Nqc;
143: PetscCall(PetscDSGetDiscretization(ds, f, &obj));
144: PetscCall(PetscObjectGetClassId(obj, &id));
145: PetscCall(PetscObjectGetName(obj, &name));
146: PetscCall(PetscViewerASCIIPrintf(viewer, "Field %s", name ? name : "<unknown>"));
147: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
148: if (id == PETSCFE_CLASSID) {
149: PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
150: PetscCall(PetscFEGetQuadrature((PetscFE)obj, &q));
151: PetscCall(PetscViewerASCIIPrintf(viewer, " FEM"));
152: } else if (id == PETSCFV_CLASSID) {
153: PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc));
154: PetscCall(PetscFVGetQuadrature((PetscFV)obj, &q));
155: PetscCall(PetscViewerASCIIPrintf(viewer, " FVM"));
156: } else SETERRQ(PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %" PetscInt_FMT, f);
157: if (Nc > 1) PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT " components", Nc));
158: else PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT " component ", Nc));
159: if (ds->implicit[f]) PetscCall(PetscViewerASCIIPrintf(viewer, " (implicit)"));
160: else PetscCall(PetscViewerASCIIPrintf(viewer, " (explicit)"));
161: if (q) {
162: PetscCall(PetscQuadratureGetData(q, NULL, &Nqc, &Nq, NULL, NULL));
163: PetscCall(PetscViewerASCIIPrintf(viewer, " (Nq %" PetscInt_FMT " Nqc %" PetscInt_FMT ")", Nq, Nqc));
164: }
165: PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT "-jet", ds->jetDegree[f]));
166: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
167: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
168: PetscCall(PetscViewerASCIIPushTab(viewer));
169: if (id == PETSCFE_CLASSID) PetscCall(PetscFEView((PetscFE)obj, viewer));
170: else if (id == PETSCFV_CLASSID) PetscCall(PetscFVView((PetscFV)obj, viewer));
171: PetscCall(PetscViewerASCIIPopTab(viewer));
173: for (b = ds->boundary; b; b = b->next) {
174: char *name;
175: PetscInt c, i;
177: if (b->field != f) continue;
178: PetscCall(PetscViewerASCIIPushTab(viewer));
179: PetscCall(PetscViewerASCIIPrintf(viewer, "Boundary %s (%s) %s\n", b->name, b->lname, DMBoundaryConditionTypes[b->type]));
180: if (!b->Nc) {
181: PetscCall(PetscViewerASCIIPrintf(viewer, " all components\n"));
182: } else {
183: PetscCall(PetscViewerASCIIPrintf(viewer, " components: "));
184: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
185: for (c = 0; c < b->Nc; ++c) {
186: if (c > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ", "));
187: PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT, b->comps[c]));
188: }
189: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
190: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
191: }
192: PetscCall(PetscViewerASCIIPrintf(viewer, " values: "));
193: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
194: for (i = 0; i < b->Nv; ++i) {
195: if (i > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ", "));
196: PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT, b->values[i]));
197: }
198: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
199: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
200: #if defined(__clang__)
201: PETSC_PRAGMA_DIAGNOSTIC_IGNORED_BEGIN("-Wformat-pedantic");
202: #elif defined(__GNUC__) || defined(__GNUG__)
203: PETSC_PRAGMA_DIAGNOSTIC_IGNORED_BEGIN("-Wformat");
204: #endif
205: if (b->func) {
206: PetscCall(PetscDLAddr(b->func, &name));
207: if (name) PetscCall(PetscViewerASCIIPrintf(viewer, " func: %s\n", name));
208: else PetscCall(PetscViewerASCIIPrintf(viewer, " func: %p\n", b->func));
209: PetscCall(PetscFree(name));
210: }
211: if (b->func_t) {
212: PetscCall(PetscDLAddr(b->func_t, &name));
213: if (name) PetscCall(PetscViewerASCIIPrintf(viewer, " func_t: %s\n", name));
214: else PetscCall(PetscViewerASCIIPrintf(viewer, " func_t: %p\n", b->func_t));
215: PetscCall(PetscFree(name));
216: }
217: PETSC_PRAGMA_DIAGNOSTIC_IGNORED_END();
218: PetscCall(PetscWeakFormView(b->wf, viewer));
219: PetscCall(PetscViewerASCIIPopTab(viewer));
220: }
221: }
222: PetscCall(PetscDSGetConstants(ds, &numConstants, &constants));
223: if (numConstants) {
224: PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " constants\n", numConstants));
225: PetscCall(PetscViewerASCIIPushTab(viewer));
226: for (f = 0; f < numConstants; ++f) PetscCall(PetscViewerASCIIPrintf(viewer, "%g\n", (double)PetscRealPart(constants[f])));
227: PetscCall(PetscViewerASCIIPopTab(viewer));
228: }
229: PetscCall(PetscWeakFormView(ds->wf, viewer));
230: PetscCall(PetscViewerASCIIPopTab(viewer));
231: PetscFunctionReturn(PETSC_SUCCESS);
232: }
234: /*@C
235: PetscDSViewFromOptions - View a `PetscDS` based on values in the options database
237: Collective
239: Input Parameters:
240: + A - the `PetscDS` object
241: . obj - Optional object that provides the options prefix used in the search
242: - name - command line option
244: Level: intermediate
246: .seealso: `PetscDSType`, `PetscDS`, `PetscDSView()`, `PetscObjectViewFromOptions()`, `PetscDSCreate()`
247: @*/
248: PetscErrorCode PetscDSViewFromOptions(PetscDS A, PetscObject obj, const char name[])
249: {
250: PetscFunctionBegin;
252: PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
253: PetscFunctionReturn(PETSC_SUCCESS);
254: }
256: /*@C
257: PetscDSView - Views a `PetscDS`
259: Collective
261: Input Parameters:
262: + prob - the `PetscDS` object to view
263: - v - the viewer
265: Level: developer
267: .seealso: `PetscDSType`, `PetscDS`, `PetscViewer`, `PetscDSDestroy()`, `PetscDSViewFromOptions()`
268: @*/
269: PetscErrorCode PetscDSView(PetscDS prob, PetscViewer v)
270: {
271: PetscBool iascii;
273: PetscFunctionBegin;
275: if (!v) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)prob), &v));
277: PetscCall(PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &iascii));
278: if (iascii) PetscCall(PetscDSView_Ascii(prob, v));
279: PetscTryTypeMethod(prob, view, v);
280: PetscFunctionReturn(PETSC_SUCCESS);
281: }
283: /*@
284: PetscDSSetFromOptions - sets parameters in a `PetscDS` from the options database
286: Collective
288: Input Parameter:
289: . prob - the `PetscDS` object to set options for
291: Options Database Keys:
292: + -petscds_type <type> - Set the `PetscDS` type
293: . -petscds_view <view opt> - View the `PetscDS`
294: . -petscds_jac_pre - Turn formation of a separate Jacobian preconditioner on or off
295: . -bc_<name> <ids> - Specify a list of label ids for a boundary condition
296: - -bc_<name>_comp <comps> - Specify a list of field components to constrain for a boundary condition
298: Level: intermediate
300: .seealso: `PetscDS`, `PetscDSView()`
301: @*/
302: PetscErrorCode PetscDSSetFromOptions(PetscDS prob)
303: {
304: DSBoundary b;
305: const char *defaultType;
306: char name[256];
307: PetscBool flg;
309: PetscFunctionBegin;
311: if (!((PetscObject)prob)->type_name) {
312: defaultType = PETSCDSBASIC;
313: } else {
314: defaultType = ((PetscObject)prob)->type_name;
315: }
316: PetscCall(PetscDSRegisterAll());
318: PetscObjectOptionsBegin((PetscObject)prob);
319: for (b = prob->boundary; b; b = b->next) {
320: char optname[1024];
321: PetscInt ids[1024], len = 1024;
322: PetscBool flg;
324: PetscCall(PetscSNPrintf(optname, sizeof(optname), "-bc_%s", b->name));
325: PetscCall(PetscMemzero(ids, sizeof(ids)));
326: PetscCall(PetscOptionsIntArray(optname, "List of boundary IDs", "", ids, &len, &flg));
327: if (flg) {
328: b->Nv = len;
329: PetscCall(PetscFree(b->values));
330: PetscCall(PetscMalloc1(len, &b->values));
331: PetscCall(PetscArraycpy(b->values, ids, len));
332: PetscCall(PetscWeakFormRewriteKeys(b->wf, b->label, len, b->values));
333: }
334: len = 1024;
335: PetscCall(PetscSNPrintf(optname, sizeof(optname), "-bc_%s_comp", b->name));
336: PetscCall(PetscMemzero(ids, sizeof(ids)));
337: PetscCall(PetscOptionsIntArray(optname, "List of boundary field components", "", ids, &len, &flg));
338: if (flg) {
339: b->Nc = len;
340: PetscCall(PetscFree(b->comps));
341: PetscCall(PetscMalloc1(len, &b->comps));
342: PetscCall(PetscArraycpy(b->comps, ids, len));
343: }
344: }
345: PetscCall(PetscOptionsFList("-petscds_type", "Discrete System", "PetscDSSetType", PetscDSList, defaultType, name, 256, &flg));
346: if (flg) {
347: PetscCall(PetscDSSetType(prob, name));
348: } else if (!((PetscObject)prob)->type_name) {
349: PetscCall(PetscDSSetType(prob, defaultType));
350: }
351: PetscCall(PetscOptionsBool("-petscds_jac_pre", "Discrete System", "PetscDSUseJacobianPreconditioner", prob->useJacPre, &prob->useJacPre, &flg));
352: PetscCall(PetscOptionsBool("-petscds_force_quad", "Discrete System", "PetscDSSetForceQuad", prob->forceQuad, &prob->forceQuad, &flg));
353: PetscTryTypeMethod(prob, setfromoptions);
354: /* process any options handlers added with PetscObjectAddOptionsHandler() */
355: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)prob, PetscOptionsObject));
356: PetscOptionsEnd();
357: if (prob->Nf) PetscCall(PetscDSViewFromOptions(prob, NULL, "-petscds_view"));
358: PetscFunctionReturn(PETSC_SUCCESS);
359: }
361: /*@C
362: PetscDSSetUp - Construct data structures for the `PetscDS`
364: Collective
366: Input Parameter:
367: . prob - the `PetscDS` object to setup
369: Level: developer
371: .seealso: `PetscDS`, `PetscDSView()`, `PetscDSDestroy()`
372: @*/
373: PetscErrorCode PetscDSSetUp(PetscDS prob)
374: {
375: const PetscInt Nf = prob->Nf;
376: PetscBool hasH = PETSC_FALSE;
377: PetscInt maxOrder[4] = {-1, -1, -1, -1};
378: PetscInt dim, dimEmbed, NbMax = 0, NcMax = 0, NqMax = 0, NsMax = 1, f;
380: PetscFunctionBegin;
382: if (prob->setup) PetscFunctionReturn(PETSC_SUCCESS);
383: /* Calculate sizes */
384: PetscCall(PetscDSGetSpatialDimension(prob, &dim));
385: PetscCall(PetscDSGetCoordinateDimension(prob, &dimEmbed));
386: prob->totDim = prob->totComp = 0;
387: PetscCall(PetscMalloc2(Nf, &prob->Nc, Nf, &prob->Nb));
388: PetscCall(PetscCalloc2(Nf + 1, &prob->off, Nf + 1, &prob->offDer));
389: PetscCall(PetscCalloc6(Nf + 1, &prob->offCohesive[0], Nf + 1, &prob->offCohesive[1], Nf + 1, &prob->offCohesive[2], Nf + 1, &prob->offDerCohesive[0], Nf + 1, &prob->offDerCohesive[1], Nf + 1, &prob->offDerCohesive[2]));
390: PetscCall(PetscMalloc2(Nf, &prob->T, Nf, &prob->Tf));
391: if (prob->forceQuad) {
392: // Note: This assumes we have one kind of cell at each dimension.
393: // We can fix this by having quadrature hold the celltype
394: PetscQuadrature maxQuad[4] = {NULL, NULL, NULL, NULL};
396: for (f = 0; f < Nf; ++f) {
397: PetscObject obj;
398: PetscClassId id;
399: PetscQuadrature q = NULL, fq = NULL;
400: PetscInt dim = -1, order = -1, forder = -1;
402: PetscCall(PetscDSGetDiscretization(prob, f, &obj));
403: if (!obj) continue;
404: PetscCall(PetscObjectGetClassId(obj, &id));
405: if (id == PETSCFE_CLASSID) {
406: PetscFE fe = (PetscFE)obj;
408: PetscCall(PetscFEGetQuadrature(fe, &q));
409: PetscCall(PetscFEGetFaceQuadrature(fe, &fq));
410: } else if (id == PETSCFV_CLASSID) {
411: PetscFV fv = (PetscFV)obj;
413: PetscCall(PetscFVGetQuadrature(fv, &q));
414: }
415: if (q) {
416: PetscCall(PetscQuadratureGetData(q, &dim, NULL, NULL, NULL, NULL));
417: PetscCall(PetscQuadratureGetOrder(q, &order));
418: if (order > maxOrder[dim]) {
419: maxOrder[dim] = order;
420: maxQuad[dim] = q;
421: }
422: }
423: if (fq) {
424: PetscCall(PetscQuadratureGetData(fq, &dim, NULL, NULL, NULL, NULL));
425: PetscCall(PetscQuadratureGetOrder(fq, &forder));
426: if (forder > maxOrder[dim]) {
427: maxOrder[dim] = forder;
428: maxQuad[dim] = fq;
429: }
430: }
431: }
432: for (f = 0; f < Nf; ++f) {
433: PetscObject obj;
434: PetscClassId id;
435: PetscQuadrature q;
436: PetscInt dim;
438: PetscCall(PetscDSGetDiscretization(prob, f, &obj));
439: if (!obj) continue;
440: PetscCall(PetscObjectGetClassId(obj, &id));
441: if (id == PETSCFE_CLASSID) {
442: PetscFE fe = (PetscFE)obj;
444: PetscCall(PetscFEGetQuadrature(fe, &q));
445: PetscCall(PetscQuadratureGetData(q, &dim, NULL, NULL, NULL, NULL));
446: PetscCall(PetscFESetQuadrature(fe, maxQuad[dim]));
447: PetscCall(PetscFESetFaceQuadrature(fe, maxQuad[dim - 1]));
448: } else if (id == PETSCFV_CLASSID) {
449: PetscFV fv = (PetscFV)obj;
451: PetscCall(PetscFVGetQuadrature(fv, &q));
452: PetscCall(PetscQuadratureGetData(q, &dim, NULL, NULL, NULL, NULL));
453: PetscCall(PetscFVSetQuadrature(fv, maxQuad[dim]));
454: }
455: }
456: }
457: for (f = 0; f < Nf; ++f) {
458: PetscObject obj;
459: PetscClassId id;
460: PetscQuadrature q = NULL;
461: PetscInt Nq = 0, Nb, Nc;
463: PetscCall(PetscDSGetDiscretization(prob, f, &obj));
464: if (prob->jetDegree[f] > 1) hasH = PETSC_TRUE;
465: if (!obj) {
466: /* Empty mesh */
467: Nb = Nc = 0;
468: prob->T[f] = prob->Tf[f] = NULL;
469: } else {
470: PetscCall(PetscObjectGetClassId(obj, &id));
471: if (id == PETSCFE_CLASSID) {
472: PetscFE fe = (PetscFE)obj;
474: PetscCall(PetscFEGetQuadrature(fe, &q));
475: {
476: PetscQuadrature fq;
477: PetscInt dim, order;
479: PetscCall(PetscQuadratureGetData(q, &dim, NULL, NULL, NULL, NULL));
480: PetscCall(PetscQuadratureGetOrder(q, &order));
481: if (maxOrder[dim] < 0) maxOrder[dim] = order;
482: PetscCheck(order == maxOrder[dim], PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Field %" PetscInt_FMT " cell quadrature order %" PetscInt_FMT " != %" PetscInt_FMT " DS cell quadrature order", f, order, maxOrder[dim]);
483: PetscCall(PetscFEGetFaceQuadrature(fe, &fq));
484: if (fq) {
485: PetscCall(PetscQuadratureGetData(fq, &dim, NULL, NULL, NULL, NULL));
486: PetscCall(PetscQuadratureGetOrder(fq, &order));
487: if (maxOrder[dim] < 0) maxOrder[dim] = order;
488: PetscCheck(order == maxOrder[dim], PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Field %" PetscInt_FMT " face quadrature order %" PetscInt_FMT " != %" PetscInt_FMT " DS face quadrature order", f, order, maxOrder[dim]);
489: }
490: }
491: PetscCall(PetscFEGetDimension(fe, &Nb));
492: PetscCall(PetscFEGetNumComponents(fe, &Nc));
493: PetscCall(PetscFEGetCellTabulation(fe, prob->jetDegree[f], &prob->T[f]));
494: PetscCall(PetscFEGetFaceTabulation(fe, prob->jetDegree[f], &prob->Tf[f]));
495: } else if (id == PETSCFV_CLASSID) {
496: PetscFV fv = (PetscFV)obj;
498: PetscCall(PetscFVGetQuadrature(fv, &q));
499: PetscCall(PetscFVGetNumComponents(fv, &Nc));
500: Nb = Nc;
501: PetscCall(PetscFVGetCellTabulation(fv, &prob->T[f]));
502: /* TODO: should PetscFV also have face tabulation? Otherwise there will be a null pointer in prob->basisFace */
503: } else SETERRQ(PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %" PetscInt_FMT, f);
504: }
505: prob->Nc[f] = Nc;
506: prob->Nb[f] = Nb;
507: prob->off[f + 1] = Nc + prob->off[f];
508: prob->offDer[f + 1] = Nc * dim + prob->offDer[f];
509: prob->offCohesive[0][f + 1] = (prob->cohesive[f] ? Nc : Nc * 2) + prob->offCohesive[0][f];
510: prob->offDerCohesive[0][f + 1] = (prob->cohesive[f] ? Nc : Nc * 2) * dimEmbed + prob->offDerCohesive[0][f];
511: prob->offCohesive[1][f] = (prob->cohesive[f] ? 0 : Nc) + prob->offCohesive[0][f];
512: prob->offDerCohesive[1][f] = (prob->cohesive[f] ? 0 : Nc) * dimEmbed + prob->offDerCohesive[0][f];
513: prob->offCohesive[2][f + 1] = (prob->cohesive[f] ? Nc : Nc * 2) + prob->offCohesive[2][f];
514: prob->offDerCohesive[2][f + 1] = (prob->cohesive[f] ? Nc : Nc * 2) * dimEmbed + prob->offDerCohesive[2][f];
515: if (q) PetscCall(PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL));
516: NqMax = PetscMax(NqMax, Nq);
517: NbMax = PetscMax(NbMax, Nb);
518: NcMax = PetscMax(NcMax, Nc);
519: prob->totDim += Nb;
520: prob->totComp += Nc;
521: /* There are two faces for all fields on a cohesive cell, except for cohesive fields */
522: if (prob->isCohesive && !prob->cohesive[f]) prob->totDim += Nb;
523: }
524: prob->offCohesive[1][Nf] = prob->offCohesive[0][Nf];
525: prob->offDerCohesive[1][Nf] = prob->offDerCohesive[0][Nf];
526: /* Allocate works space */
527: NsMax = 2; /* A non-cohesive discretizations can be used on a cohesive cell, so we need this extra workspace for all DS */
528: PetscCall(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));
529: PetscCall(PetscMalloc5(dimEmbed, &prob->x, NbMax * NcMax, &prob->basisReal, NbMax * NcMax * dimEmbed, &prob->basisDerReal, NbMax * NcMax, &prob->testReal, NbMax * NcMax * dimEmbed, &prob->testDerReal));
530: PetscCall(PetscMalloc6(NsMax * NqMax * NcMax, &prob->f0, NsMax * NqMax * NcMax * dimEmbed, &prob->f1, NsMax * NsMax * NqMax * NcMax * NcMax, &prob->g0, NsMax * NsMax * NqMax * NcMax * NcMax * dimEmbed, &prob->g1, NsMax * NsMax * NqMax * NcMax * NcMax * dimEmbed,
531: &prob->g2, NsMax * NsMax * NqMax * NcMax * NcMax * dimEmbed * dimEmbed, &prob->g3));
532: PetscTryTypeMethod(prob, setup);
533: prob->setup = PETSC_TRUE;
534: PetscFunctionReturn(PETSC_SUCCESS);
535: }
537: static PetscErrorCode PetscDSDestroyStructs_Static(PetscDS prob)
538: {
539: PetscFunctionBegin;
540: PetscCall(PetscFree2(prob->Nc, prob->Nb));
541: PetscCall(PetscFree2(prob->off, prob->offDer));
542: PetscCall(PetscFree6(prob->offCohesive[0], prob->offCohesive[1], prob->offCohesive[2], prob->offDerCohesive[0], prob->offDerCohesive[1], prob->offDerCohesive[2]));
543: PetscCall(PetscFree2(prob->T, prob->Tf));
544: PetscCall(PetscFree3(prob->u, prob->u_t, prob->u_x));
545: PetscCall(PetscFree5(prob->x, prob->basisReal, prob->basisDerReal, prob->testReal, prob->testDerReal));
546: PetscCall(PetscFree6(prob->f0, prob->f1, prob->g0, prob->g1, prob->g2, prob->g3));
547: PetscFunctionReturn(PETSC_SUCCESS);
548: }
550: static PetscErrorCode PetscDSEnlarge_Static(PetscDS prob, PetscInt NfNew)
551: {
552: PetscObject *tmpd;
553: PetscBool *tmpi;
554: PetscInt *tmpk;
555: PetscBool *tmpc;
556: PetscPointFunc *tmpup;
557: PetscSimplePointFunc *tmpexactSol, *tmpexactSol_t;
558: void **tmpexactCtx, **tmpexactCtx_t;
559: void **tmpctx;
560: PetscInt Nf = prob->Nf, f;
562: PetscFunctionBegin;
563: if (Nf >= NfNew) PetscFunctionReturn(PETSC_SUCCESS);
564: prob->setup = PETSC_FALSE;
565: PetscCall(PetscDSDestroyStructs_Static(prob));
566: PetscCall(PetscMalloc4(NfNew, &tmpd, NfNew, &tmpi, NfNew, &tmpc, NfNew, &tmpk));
567: for (f = 0; f < Nf; ++f) {
568: tmpd[f] = prob->disc[f];
569: tmpi[f] = prob->implicit[f];
570: tmpc[f] = prob->cohesive[f];
571: tmpk[f] = prob->jetDegree[f];
572: }
573: for (f = Nf; f < NfNew; ++f) {
574: tmpd[f] = NULL;
575: tmpi[f] = PETSC_TRUE, tmpc[f] = PETSC_FALSE;
576: tmpk[f] = 1;
577: }
578: PetscCall(PetscFree4(prob->disc, prob->implicit, prob->cohesive, prob->jetDegree));
579: PetscCall(PetscWeakFormSetNumFields(prob->wf, NfNew));
580: prob->Nf = NfNew;
581: prob->disc = tmpd;
582: prob->implicit = tmpi;
583: prob->cohesive = tmpc;
584: prob->jetDegree = tmpk;
585: PetscCall(PetscCalloc2(NfNew, &tmpup, NfNew, &tmpctx));
586: for (f = 0; f < Nf; ++f) tmpup[f] = prob->update[f];
587: for (f = 0; f < Nf; ++f) tmpctx[f] = prob->ctx[f];
588: for (f = Nf; f < NfNew; ++f) tmpup[f] = NULL;
589: for (f = Nf; f < NfNew; ++f) tmpctx[f] = NULL;
590: PetscCall(PetscFree2(prob->update, prob->ctx));
591: prob->update = tmpup;
592: prob->ctx = tmpctx;
593: PetscCall(PetscCalloc4(NfNew, &tmpexactSol, NfNew, &tmpexactCtx, NfNew, &tmpexactSol_t, NfNew, &tmpexactCtx_t));
594: for (f = 0; f < Nf; ++f) tmpexactSol[f] = prob->exactSol[f];
595: for (f = 0; f < Nf; ++f) tmpexactCtx[f] = prob->exactCtx[f];
596: for (f = 0; f < Nf; ++f) tmpexactSol_t[f] = prob->exactSol_t[f];
597: for (f = 0; f < Nf; ++f) tmpexactCtx_t[f] = prob->exactCtx_t[f];
598: for (f = Nf; f < NfNew; ++f) tmpexactSol[f] = NULL;
599: for (f = Nf; f < NfNew; ++f) tmpexactCtx[f] = NULL;
600: for (f = Nf; f < NfNew; ++f) tmpexactSol_t[f] = NULL;
601: for (f = Nf; f < NfNew; ++f) tmpexactCtx_t[f] = NULL;
602: PetscCall(PetscFree4(prob->exactSol, prob->exactCtx, prob->exactSol_t, prob->exactCtx_t));
603: prob->exactSol = tmpexactSol;
604: prob->exactCtx = tmpexactCtx;
605: prob->exactSol_t = tmpexactSol_t;
606: prob->exactCtx_t = tmpexactCtx_t;
607: PetscFunctionReturn(PETSC_SUCCESS);
608: }
610: /*@
611: PetscDSDestroy - Destroys a `PetscDS` object
613: Collective
615: Input Parameter:
616: . prob - the `PetscDS` object to destroy
618: Level: developer
620: .seealso: `PetscDSView()`
621: @*/
622: PetscErrorCode PetscDSDestroy(PetscDS *ds)
623: {
624: PetscInt f;
626: PetscFunctionBegin;
627: if (!*ds) PetscFunctionReturn(PETSC_SUCCESS);
630: if (--((PetscObject)(*ds))->refct > 0) {
631: *ds = NULL;
632: PetscFunctionReturn(PETSC_SUCCESS);
633: }
634: ((PetscObject)(*ds))->refct = 0;
635: if ((*ds)->subprobs) {
636: PetscInt dim, d;
638: PetscCall(PetscDSGetSpatialDimension(*ds, &dim));
639: for (d = 0; d < dim; ++d) PetscCall(PetscDSDestroy(&(*ds)->subprobs[d]));
640: }
641: PetscCall(PetscFree((*ds)->subprobs));
642: PetscCall(PetscDSDestroyStructs_Static(*ds));
643: for (f = 0; f < (*ds)->Nf; ++f) PetscCall(PetscObjectDereference((*ds)->disc[f]));
644: PetscCall(PetscFree4((*ds)->disc, (*ds)->implicit, (*ds)->cohesive, (*ds)->jetDegree));
645: PetscCall(PetscWeakFormDestroy(&(*ds)->wf));
646: PetscCall(PetscFree2((*ds)->update, (*ds)->ctx));
647: PetscCall(PetscFree4((*ds)->exactSol, (*ds)->exactCtx, (*ds)->exactSol_t, (*ds)->exactCtx_t));
648: PetscTryTypeMethod((*ds), destroy);
649: PetscCall(PetscDSDestroyBoundary(*ds));
650: PetscCall(PetscFree((*ds)->constants));
651: for (PetscInt c = 0; c < DM_NUM_POLYTOPES; ++c) {
652: const PetscInt Na = DMPolytopeTypeGetNumArrangments((DMPolytopeType)c);
653: if ((*ds)->quadPerm[c])
654: for (PetscInt o = 0; o < Na; ++o) PetscCall(ISDestroy(&(*ds)->quadPerm[c][o]));
655: PetscCall(PetscFree((*ds)->quadPerm[c]));
656: (*ds)->quadPerm[c] = NULL;
657: }
658: PetscCall(PetscHeaderDestroy(ds));
659: PetscFunctionReturn(PETSC_SUCCESS);
660: }
662: /*@
663: PetscDSCreate - Creates an empty `PetscDS` object. The type can then be set with `PetscDSSetType()`.
665: Collective
667: Input Parameter:
668: . comm - The communicator for the `PetscDS` object
670: Output Parameter:
671: . ds - The `PetscDS` object
673: Level: beginner
675: .seealso: `PetscDS`, `PetscDSSetType()`, `PETSCDSBASIC`, `PetscDSType`
676: @*/
677: PetscErrorCode PetscDSCreate(MPI_Comm comm, PetscDS *ds)
678: {
679: PetscDS p;
681: PetscFunctionBegin;
683: *ds = NULL;
684: PetscCall(PetscDSInitializePackage());
686: PetscCall(PetscHeaderCreate(p, PETSCDS_CLASSID, "PetscDS", "Discrete System", "PetscDS", comm, PetscDSDestroy, PetscDSView));
688: p->Nf = 0;
689: p->setup = PETSC_FALSE;
690: p->numConstants = 0;
691: p->constants = NULL;
692: p->dimEmbed = -1;
693: p->useJacPre = PETSC_TRUE;
694: p->forceQuad = PETSC_TRUE;
695: PetscCall(PetscWeakFormCreate(comm, &p->wf));
696: PetscCall(PetscArrayzero(p->quadPerm, DM_NUM_POLYTOPES));
698: *ds = p;
699: PetscFunctionReturn(PETSC_SUCCESS);
700: }
702: /*@
703: PetscDSGetNumFields - Returns the number of fields in the `PetscDS`
705: Not Collective
707: Input Parameter:
708: . prob - The `PetscDS` object
710: Output Parameter:
711: . Nf - The number of fields
713: Level: beginner
715: .seealso: `PetscDS`, `PetscDSGetSpatialDimension()`, `PetscDSCreate()`
716: @*/
717: PetscErrorCode PetscDSGetNumFields(PetscDS prob, PetscInt *Nf)
718: {
719: PetscFunctionBegin;
722: *Nf = prob->Nf;
723: PetscFunctionReturn(PETSC_SUCCESS);
724: }
726: /*@
727: PetscDSGetSpatialDimension - Returns the spatial dimension of the `PetscDS`, meaning the topological dimension of the discretizations
729: Not Collective
731: Input Parameter:
732: . prob - The `PetscDS` object
734: Output Parameter:
735: . dim - The spatial dimension
737: Level: beginner
739: .seealso: `PetscDS`, `PetscDSGetCoordinateDimension()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
740: @*/
741: PetscErrorCode PetscDSGetSpatialDimension(PetscDS prob, PetscInt *dim)
742: {
743: PetscFunctionBegin;
746: *dim = 0;
747: if (prob->Nf) {
748: PetscObject obj;
749: PetscClassId id;
751: PetscCall(PetscDSGetDiscretization(prob, 0, &obj));
752: if (obj) {
753: PetscCall(PetscObjectGetClassId(obj, &id));
754: if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetSpatialDimension((PetscFE)obj, dim));
755: else if (id == PETSCFV_CLASSID) PetscCall(PetscFVGetSpatialDimension((PetscFV)obj, dim));
756: else SETERRQ(PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", 0);
757: }
758: }
759: PetscFunctionReturn(PETSC_SUCCESS);
760: }
762: /*@
763: PetscDSGetCoordinateDimension - Returns the coordinate dimension of the `PetscDS`, meaning the dimension of the space into which the discretiaztions are embedded
765: Not Collective
767: Input Parameter:
768: . prob - The `PetscDS` object
770: Output Parameter:
771: . dimEmbed - The coordinate dimension
773: Level: beginner
775: .seealso: `PetscDS`, `PetscDSSetCoordinateDimension()`, `PetscDSGetSpatialDimension()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
776: @*/
777: PetscErrorCode PetscDSGetCoordinateDimension(PetscDS prob, PetscInt *dimEmbed)
778: {
779: PetscFunctionBegin;
782: PetscCheck(prob->dimEmbed >= 0, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONGSTATE, "No coordinate dimension set for this DS");
783: *dimEmbed = prob->dimEmbed;
784: PetscFunctionReturn(PETSC_SUCCESS);
785: }
787: /*@
788: PetscDSSetCoordinateDimension - Set the coordinate dimension of the `PetscDS`, meaning the dimension of the space into which the discretiaztions are embedded
790: Logically Collective
792: Input Parameters:
793: + prob - The `PetscDS` object
794: - dimEmbed - The coordinate dimension
796: Level: beginner
798: .seealso: `PetscDS`, `PetscDSGetCoordinateDimension()`, `PetscDSGetSpatialDimension()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
799: @*/
800: PetscErrorCode PetscDSSetCoordinateDimension(PetscDS prob, PetscInt dimEmbed)
801: {
802: PetscFunctionBegin;
804: PetscCheck(dimEmbed >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Coordinate dimension must be non-negative, not %" PetscInt_FMT, dimEmbed);
805: prob->dimEmbed = dimEmbed;
806: PetscFunctionReturn(PETSC_SUCCESS);
807: }
809: /*@
810: PetscDSGetForceQuad - Returns the flag to force matching quadratures among the field discretizations
812: Not collective
814: Input Parameter:
815: . prob - The `PetscDS` object
817: Output Parameter:
818: . forceQuad - The flag
820: Level: intermediate
822: .seealso: `PetscDS`, `PetscDSSetForceQuad()`, `PetscDSGetDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
823: @*/
824: PetscErrorCode PetscDSGetForceQuad(PetscDS ds, PetscBool *forceQuad)
825: {
826: PetscFunctionBegin;
829: *forceQuad = ds->forceQuad;
830: PetscFunctionReturn(PETSC_SUCCESS);
831: }
833: /*@
834: PetscDSSetForceQuad - Set the flag to force matching quadratures among the field discretizations
836: Logically collective on ds
838: Input Parameters:
839: + ds - The `PetscDS` object
840: - forceQuad - The flag
842: Level: intermediate
844: .seealso: `PetscDS`, `PetscDSGetForceQuad()`, `PetscDSGetDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
845: @*/
846: PetscErrorCode PetscDSSetForceQuad(PetscDS ds, PetscBool forceQuad)
847: {
848: PetscFunctionBegin;
850: ds->forceQuad = forceQuad;
851: PetscFunctionReturn(PETSC_SUCCESS);
852: }
854: /*@
855: PetscDSIsCohesive - Returns the flag indicating that this `PetscDS` is for a cohesive cell
857: Not Collective
859: Input Parameter:
860: . ds - The `PetscDS` object
862: Output Parameter:
863: . isCohesive - The flag
865: Level: developer
867: .seealso: `PetscDS`, `PetscDSGetNumCohesive()`, `PetscDSGetCohesive()`, `PetscDSSetCohesive()`, `PetscDSCreate()`
868: @*/
869: PetscErrorCode PetscDSIsCohesive(PetscDS ds, PetscBool *isCohesive)
870: {
871: PetscFunctionBegin;
874: *isCohesive = ds->isCohesive;
875: PetscFunctionReturn(PETSC_SUCCESS);
876: }
878: /*@
879: PetscDSGetNumCohesive - Returns the number of cohesive fields, meaning those defined on the interior of a cohesive cell
881: Not Collective
883: Input Parameter:
884: . ds - The `PetscDS` object
886: Output Parameter:
887: . numCohesive - The number of cohesive fields
889: Level: developer
891: .seealso: `PetscDS`, `PetscDSSetCohesive()`, `PetscDSCreate()`
892: @*/
893: PetscErrorCode PetscDSGetNumCohesive(PetscDS ds, PetscInt *numCohesive)
894: {
895: PetscInt f;
897: PetscFunctionBegin;
900: *numCohesive = 0;
901: for (f = 0; f < ds->Nf; ++f) *numCohesive += ds->cohesive[f] ? 1 : 0;
902: PetscFunctionReturn(PETSC_SUCCESS);
903: }
905: /*@
906: PetscDSGetCohesive - Returns the flag indicating that a field is cohesive, meaning it is defined on the interior of a cohesive cell
908: Not Collective
910: Input Parameters:
911: + ds - The `PetscDS` object
912: - f - The field index
914: Output Parameter:
915: . isCohesive - The flag
917: Level: developer
919: .seealso: `PetscDS`, `PetscDSSetCohesive()`, `PetscDSIsCohesive()`, `PetscDSCreate()`
920: @*/
921: PetscErrorCode PetscDSGetCohesive(PetscDS ds, PetscInt f, PetscBool *isCohesive)
922: {
923: PetscFunctionBegin;
926: PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
927: *isCohesive = ds->cohesive[f];
928: PetscFunctionReturn(PETSC_SUCCESS);
929: }
931: /*@
932: PetscDSSetCohesive - Set the flag indicating that a field is cohesive, meaning it is defined on the interior of a cohesive cell
934: Not Collective
936: Input Parameters:
937: + ds - The `PetscDS` object
938: . f - The field index
939: - isCohesive - The flag for a cohesive field
941: Level: developer
943: .seealso: `PetscDS`, `PetscDSGetCohesive()`, `PetscDSIsCohesive()`, `PetscDSCreate()`
944: @*/
945: PetscErrorCode PetscDSSetCohesive(PetscDS ds, PetscInt f, PetscBool isCohesive)
946: {
947: PetscInt i;
949: PetscFunctionBegin;
951: PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
952: ds->cohesive[f] = isCohesive;
953: ds->isCohesive = PETSC_FALSE;
954: for (i = 0; i < ds->Nf; ++i) ds->isCohesive = ds->isCohesive || ds->cohesive[f] ? PETSC_TRUE : PETSC_FALSE;
955: PetscFunctionReturn(PETSC_SUCCESS);
956: }
958: /*@
959: PetscDSGetTotalDimension - Returns the total size of the approximation space for this system
961: Not Collective
963: Input Parameter:
964: . prob - The `PetscDS` object
966: Output Parameter:
967: . dim - The total problem dimension
969: Level: beginner
971: .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
972: @*/
973: PetscErrorCode PetscDSGetTotalDimension(PetscDS prob, PetscInt *dim)
974: {
975: PetscFunctionBegin;
977: PetscCall(PetscDSSetUp(prob));
979: *dim = prob->totDim;
980: PetscFunctionReturn(PETSC_SUCCESS);
981: }
983: /*@
984: PetscDSGetTotalComponents - Returns the total number of components in this system
986: Not Collective
988: Input Parameter:
989: . prob - The `PetscDS` object
991: Output Parameter:
992: . dim - The total number of components
994: Level: beginner
996: .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
997: @*/
998: PetscErrorCode PetscDSGetTotalComponents(PetscDS prob, PetscInt *Nc)
999: {
1000: PetscFunctionBegin;
1002: PetscCall(PetscDSSetUp(prob));
1004: *Nc = prob->totComp;
1005: PetscFunctionReturn(PETSC_SUCCESS);
1006: }
1008: /*@
1009: PetscDSGetDiscretization - Returns the discretization object for the given field
1011: Not Collective
1013: Input Parameters:
1014: + prob - The `PetscDS` object
1015: - f - The field number
1017: Output Parameter:
1018: . disc - The discretization object
1020: Level: beginner
1022: .seealso: `PetscDS`, `PetscFE`, `PetscFV`, `PetscDSSetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1023: @*/
1024: PetscErrorCode PetscDSGetDiscretization(PetscDS prob, PetscInt f, PetscObject *disc)
1025: {
1026: PetscFunctionBeginHot;
1029: PetscCheck(!(f < 0) && !(f >= prob->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, prob->Nf);
1030: *disc = prob->disc[f];
1031: PetscFunctionReturn(PETSC_SUCCESS);
1032: }
1034: /*@
1035: PetscDSSetDiscretization - Sets the discretization object for the given field
1037: Not Collective
1039: Input Parameters:
1040: + prob - The `PetscDS` object
1041: . f - The field number
1042: - disc - The discretization object
1044: Level: beginner
1046: .seealso: `PetscDS`, `PetscFE`, `PetscFV`, `PetscDSGetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1047: @*/
1048: PetscErrorCode PetscDSSetDiscretization(PetscDS prob, PetscInt f, PetscObject disc)
1049: {
1050: PetscFunctionBegin;
1053: PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
1054: PetscCall(PetscDSEnlarge_Static(prob, f + 1));
1055: PetscCall(PetscObjectDereference(prob->disc[f]));
1056: prob->disc[f] = disc;
1057: PetscCall(PetscObjectReference(disc));
1058: if (disc) {
1059: PetscClassId id;
1061: PetscCall(PetscObjectGetClassId(disc, &id));
1062: if (id == PETSCFE_CLASSID) {
1063: PetscCall(PetscDSSetImplicit(prob, f, PETSC_TRUE));
1064: } else if (id == PETSCFV_CLASSID) {
1065: PetscCall(PetscDSSetImplicit(prob, f, PETSC_FALSE));
1066: }
1067: PetscCall(PetscDSSetJetDegree(prob, f, 1));
1068: }
1069: PetscFunctionReturn(PETSC_SUCCESS);
1070: }
1072: /*@
1073: PetscDSGetWeakForm - Returns the weak form object
1075: Not Collective
1077: Input Parameter:
1078: . ds - The `PetscDS` object
1080: Output Parameter:
1081: . wf - The weak form object
1083: Level: beginner
1085: .seealso: `PetscWeakForm`, `PetscDSSetWeakForm()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1086: @*/
1087: PetscErrorCode PetscDSGetWeakForm(PetscDS ds, PetscWeakForm *wf)
1088: {
1089: PetscFunctionBegin;
1092: *wf = ds->wf;
1093: PetscFunctionReturn(PETSC_SUCCESS);
1094: }
1096: /*@
1097: PetscDSSetWeakForm - Sets the weak form object
1099: Not Collective
1101: Input Parameters:
1102: + ds - The `PetscDS` object
1103: - wf - The weak form object
1105: Level: beginner
1107: .seealso: `PetscWeakForm`, `PetscDSGetWeakForm()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1108: @*/
1109: PetscErrorCode PetscDSSetWeakForm(PetscDS ds, PetscWeakForm wf)
1110: {
1111: PetscFunctionBegin;
1114: PetscCall(PetscObjectDereference((PetscObject)ds->wf));
1115: ds->wf = wf;
1116: PetscCall(PetscObjectReference((PetscObject)wf));
1117: PetscCall(PetscWeakFormSetNumFields(wf, ds->Nf));
1118: PetscFunctionReturn(PETSC_SUCCESS);
1119: }
1121: /*@
1122: PetscDSAddDiscretization - Adds a discretization object
1124: Not Collective
1126: Input Parameters:
1127: + prob - The `PetscDS` object
1128: - disc - The boundary discretization object
1130: Level: beginner
1132: .seealso: `PetscWeakForm`, `PetscDSGetDiscretization()`, `PetscDSSetDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1133: @*/
1134: PetscErrorCode PetscDSAddDiscretization(PetscDS prob, PetscObject disc)
1135: {
1136: PetscFunctionBegin;
1137: PetscCall(PetscDSSetDiscretization(prob, prob->Nf, disc));
1138: PetscFunctionReturn(PETSC_SUCCESS);
1139: }
1141: /*@
1142: PetscDSGetQuadrature - Returns the quadrature, which must agree for all fields in the `PetscDS`
1144: Not Collective
1146: Input Parameter:
1147: . prob - The `PetscDS` object
1149: Output Parameter:
1150: . q - The quadrature object
1152: Level: intermediate
1154: .seealso: `PetscDS`, `PetscQuadrature`, `PetscDSSetImplicit()`, `PetscDSSetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1155: @*/
1156: PetscErrorCode PetscDSGetQuadrature(PetscDS prob, PetscQuadrature *q)
1157: {
1158: PetscObject obj;
1159: PetscClassId id;
1161: PetscFunctionBegin;
1162: *q = NULL;
1163: if (!prob->Nf) PetscFunctionReturn(PETSC_SUCCESS);
1164: PetscCall(PetscDSGetDiscretization(prob, 0, &obj));
1165: PetscCall(PetscObjectGetClassId(obj, &id));
1166: if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetQuadrature((PetscFE)obj, q));
1167: else if (id == PETSCFV_CLASSID) PetscCall(PetscFVGetQuadrature((PetscFV)obj, q));
1168: else SETERRQ(PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", 0);
1169: PetscFunctionReturn(PETSC_SUCCESS);
1170: }
1172: /*@
1173: PetscDSGetImplicit - Returns the flag for implicit solve for this field. This is just a guide for `TSIMEX`
1175: Not Collective
1177: Input Parameters:
1178: + prob - The `PetscDS` object
1179: - f - The field number
1181: Output Parameter:
1182: . implicit - The flag indicating what kind of solve to use for this field
1184: Level: developer
1186: .seealso: `TSIMEX`, `PetscDS`, `PetscDSSetImplicit()`, `PetscDSSetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1187: @*/
1188: PetscErrorCode PetscDSGetImplicit(PetscDS prob, PetscInt f, PetscBool *implicit)
1189: {
1190: PetscFunctionBegin;
1193: PetscCheck(!(f < 0) && !(f >= prob->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, prob->Nf);
1194: *implicit = prob->implicit[f];
1195: PetscFunctionReturn(PETSC_SUCCESS);
1196: }
1198: /*@
1199: PetscDSSetImplicit - Set the flag for implicit solve for this field. This is just a guide for `TSIMEX`
1201: Not Collective
1203: Input Parameters:
1204: + prob - The `PetscDS` object
1205: . f - The field number
1206: - implicit - The flag indicating what kind of solve to use for this field
1208: Level: developer
1210: .seealso: `TSIMEX`, `PetscDSGetImplicit()`, `PetscDSSetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1211: @*/
1212: PetscErrorCode PetscDSSetImplicit(PetscDS prob, PetscInt f, PetscBool implicit)
1213: {
1214: PetscFunctionBegin;
1216: PetscCheck(!(f < 0) && !(f >= prob->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, prob->Nf);
1217: prob->implicit[f] = implicit;
1218: PetscFunctionReturn(PETSC_SUCCESS);
1219: }
1221: /*@
1222: PetscDSGetJetDegree - Returns the highest derivative for this field equation, or the k-jet that the discretization needs to tabulate.
1224: Not Collective
1226: Input Parameters:
1227: + ds - The `PetscDS` object
1228: - f - The field number
1230: Output Parameter:
1231: . k - The highest derivative we need to tabulate
1233: Level: developer
1235: .seealso: `PetscDS`, `PetscDSSetJetDegree()`, `PetscDSSetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1236: @*/
1237: PetscErrorCode PetscDSGetJetDegree(PetscDS ds, PetscInt f, PetscInt *k)
1238: {
1239: PetscFunctionBegin;
1242: PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
1243: *k = ds->jetDegree[f];
1244: PetscFunctionReturn(PETSC_SUCCESS);
1245: }
1247: /*@
1248: PetscDSSetJetDegree - Set the highest derivative for this field equation, or the k-jet that the discretization needs to tabulate.
1250: Not Collective
1252: Input Parameters:
1253: + ds - The `PetscDS` object
1254: . f - The field number
1255: - k - The highest derivative we need to tabulate
1257: Level: developer
1259: .seealso: ``PetscDS`, PetscDSGetJetDegree()`, `PetscDSSetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
1260: @*/
1261: PetscErrorCode PetscDSSetJetDegree(PetscDS ds, PetscInt f, PetscInt k)
1262: {
1263: PetscFunctionBegin;
1265: PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
1266: ds->jetDegree[f] = k;
1267: PetscFunctionReturn(PETSC_SUCCESS);
1268: }
1270: /*@C
1271: PetscDSGetObjective - Get the pointwise objective function for a given test field
1273: Not Collective
1275: Input Parameters:
1276: + ds - The `PetscDS`
1277: - f - The test field number
1279: Output Parameters:
1280: . obj - integrand for the test function term
1282: Calling sequence of `obj`:
1283: .vb
1284: void f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1285: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1286: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1287: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[])
1288: .ve
1289: + dim - the spatial dimension
1290: . Nf - the number of fields
1291: . uOff - the offset into u[] and u_t[] for each field
1292: . uOff_x - the offset into u_x[] for each field
1293: . u - each field evaluated at the current point
1294: . u_t - the time derivative of each field evaluated at the current point
1295: . u_x - the gradient of each field evaluated at the current point
1296: . aOff - the offset into a[] and a_t[] for each auxiliary field
1297: . aOff_x - the offset into a_x[] for each auxiliary field
1298: . a - each auxiliary field evaluated at the current point
1299: . a_t - the time derivative of each auxiliary field evaluated at the current point
1300: . a_x - the gradient of auxiliary each field evaluated at the current point
1301: . t - current time
1302: . x - coordinates of the current point
1303: . numConstants - number of constant parameters
1304: . constants - constant parameters
1305: - obj - output values at the current point
1307: Level: intermediate
1309: Note:
1310: We are using a first order FEM model for the weak form: \int_\Omega \phi obj(u, u_t, \nabla u, x, t)
1312: .seealso: `PetscDS`, `PetscDSSetObjective()`, `PetscDSGetResidual()`
1313: @*/
1314: PetscErrorCode PetscDSGetObjective(PetscDS ds, PetscInt f, void (**obj)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[]))
1315: {
1316: PetscPointFunc *tmp;
1317: PetscInt n;
1319: PetscFunctionBegin;
1322: PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
1323: PetscCall(PetscWeakFormGetObjective(ds->wf, NULL, 0, f, 0, &n, &tmp));
1324: *obj = tmp ? tmp[0] : NULL;
1325: PetscFunctionReturn(PETSC_SUCCESS);
1326: }
1328: /*@C
1329: PetscDSSetObjective - Set the pointwise objective function for a given test field
1331: Not Collective
1333: Input Parameters:
1334: + ds - The `PetscDS`
1335: . f - The test field number
1336: - obj - integrand for the test function term
1338: Calling sequence of `obj`:
1339: .vb
1340: void f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1341: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1342: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1343: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[])
1344: .ve
1345: + dim - the spatial dimension
1346: . Nf - the number of fields
1347: . uOff - the offset into u[] and u_t[] for each field
1348: . uOff_x - the offset into u_x[] for each field
1349: . u - each field evaluated at the current point
1350: . u_t - the time derivative of each field evaluated at the current point
1351: . u_x - the gradient of each field evaluated at the current point
1352: . aOff - the offset into a[] and a_t[] for each auxiliary field
1353: . aOff_x - the offset into a_x[] for each auxiliary field
1354: . a - each auxiliary field evaluated at the current point
1355: . a_t - the time derivative of each auxiliary field evaluated at the current point
1356: . a_x - the gradient of auxiliary each field evaluated at the current point
1357: . t - current time
1358: . x - coordinates of the current point
1359: . numConstants - number of constant parameters
1360: . constants - constant parameters
1361: - obj - output values at the current point
1363: Level: intermediate
1365: Note:
1366: We are using a first order FEM model for the weak form: \int_\Omega \phi obj(u, u_t, \nabla u, x, t)
1368: .seealso: `PetscDS`, `PetscDSGetObjective()`, `PetscDSSetResidual()`
1369: @*/
1370: PetscErrorCode PetscDSSetObjective(PetscDS ds, PetscInt f, void (*obj)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[]))
1371: {
1372: PetscFunctionBegin;
1375: PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
1376: PetscCall(PetscWeakFormSetIndexObjective(ds->wf, NULL, 0, f, 0, 0, obj));
1377: PetscFunctionReturn(PETSC_SUCCESS);
1378: }
1380: /*@C
1381: PetscDSGetResidual - Get the pointwise residual function for a given test field
1383: Not Collective
1385: Input Parameters:
1386: + ds - The `PetscDS`
1387: - f - The test field number
1389: Output Parameters:
1390: + f0 - integrand for the test function term
1391: - f1 - integrand for the test function gradient term
1393: Calling sequence of `f0` and `f1`:
1394: .vb
1395: void f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1396: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1397: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1398: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
1399: .ve
1400: + dim - the spatial dimension
1401: . Nf - the number of fields
1402: . uOff - the offset into u[] and u_t[] for each field
1403: . uOff_x - the offset into u_x[] for each field
1404: . u - each field evaluated at the current point
1405: . u_t - the time derivative of each field evaluated at the current point
1406: . u_x - the gradient of each field evaluated at the current point
1407: . aOff - the offset into a[] and a_t[] for each auxiliary field
1408: . aOff_x - the offset into a_x[] for each auxiliary field
1409: . a - each auxiliary field evaluated at the current point
1410: . a_t - the time derivative of each auxiliary field evaluated at the current point
1411: . a_x - the gradient of auxiliary each field evaluated at the current point
1412: . t - current time
1413: . x - coordinates of the current point
1414: . numConstants - number of constant parameters
1415: . constants - constant parameters
1416: - f0 - output values at the current point
1418: Level: intermediate
1420: Note:
1421: We are using a first order FEM model for the weak form: \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)
1423: .seealso: `PetscDS`, `PetscDSSetResidual()`
1424: @*/
1425: PetscErrorCode PetscDSGetResidual(PetscDS ds, PetscInt f, void (**f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]), void (**f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
1426: {
1427: PetscPointFunc *tmp0, *tmp1;
1428: PetscInt n0, n1;
1430: PetscFunctionBegin;
1432: PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
1433: PetscCall(PetscWeakFormGetResidual(ds->wf, NULL, 0, f, 0, &n0, &tmp0, &n1, &tmp1));
1434: *f0 = tmp0 ? tmp0[0] : NULL;
1435: *f1 = tmp1 ? tmp1[0] : NULL;
1436: PetscFunctionReturn(PETSC_SUCCESS);
1437: }
1439: /*@C
1440: PetscDSSetResidual - Set the pointwise residual function for a given test field
1442: Not Collective
1444: Input Parameters:
1445: + ds - The `PetscDS`
1446: . f - The test field number
1447: . f0 - integrand for the test function term
1448: - f1 - integrand for the test function gradient term
1450: Calling sequence of `f0` and `f1`:
1451: .vb
1452: void f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1453: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1454: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1455: PetscReal t, const PetscReal x[], PetscScalar f0[])
1456: .ve
1457: + dim - the spatial dimension
1458: . Nf - the number of fields
1459: . uOff - the offset into u[] and u_t[] for each field
1460: . uOff_x - the offset into u_x[] for each field
1461: . u - each field evaluated at the current point
1462: . u_t - the time derivative of each field evaluated at the current point
1463: . u_x - the gradient of each field evaluated at the current point
1464: . aOff - the offset into a[] and a_t[] for each auxiliary field
1465: . aOff_x - the offset into a_x[] for each auxiliary field
1466: . a - each auxiliary field evaluated at the current point
1467: . a_t - the time derivative of each auxiliary field evaluated at the current point
1468: . a_x - the gradient of auxiliary each field evaluated at the current point
1469: . t - current time
1470: . x - coordinates of the current point
1471: . numConstants - number of constant parameters
1472: . constants - constant parameters
1473: - f0 - output values at the current point
1475: Level: intermediate
1477: Note:
1478: We are using a first order FEM model for the weak form: \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)
1480: .seealso: `PetscDS`, `PetscDSGetResidual()`
1481: @*/
1482: PetscErrorCode PetscDSSetResidual(PetscDS ds, PetscInt f, void (*f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]), void (*f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
1483: {
1484: PetscFunctionBegin;
1488: PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
1489: PetscCall(PetscWeakFormSetIndexResidual(ds->wf, NULL, 0, f, 0, 0, f0, 0, f1));
1490: PetscFunctionReturn(PETSC_SUCCESS);
1491: }
1493: /*@C
1494: PetscDSGetRHSResidual - Get the pointwise RHS residual function for explicit timestepping for a given test field
1496: Not Collective
1498: Input Parameters:
1499: + ds - The `PetscDS`
1500: - f - The test field number
1502: Output Parameters:
1503: + f0 - integrand for the test function term
1504: - f1 - integrand for the test function gradient term
1506: Calling sequence of `f0` and `f1`:
1507: .vb
1508: void f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1509: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1510: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1511: PetscReal t, const PetscReal x[], PetscScalar f0[])
1512: .ve
1513: + dim - the spatial dimension
1514: . Nf - the number of fields
1515: . uOff - the offset into u[] and u_t[] for each field
1516: . uOff_x - the offset into u_x[] for each field
1517: . u - each field evaluated at the current point
1518: . u_t - the time derivative of each field evaluated at the current point
1519: . u_x - the gradient of each field evaluated at the current point
1520: . aOff - the offset into a[] and a_t[] for each auxiliary field
1521: . aOff_x - the offset into a_x[] for each auxiliary field
1522: . a - each auxiliary field evaluated at the current point
1523: . a_t - the time derivative of each auxiliary field evaluated at the current point
1524: . a_x - the gradient of auxiliary each field evaluated at the current point
1525: . t - current time
1526: . x - coordinates of the current point
1527: . numConstants - number of constant parameters
1528: . constants - constant parameters
1529: - f0 - output values at the current point
1531: Level: intermediate
1533: Note:
1534: We are using a first order FEM model for the weak form: \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)
1536: .seealso: `PetscDS`, `PetscDSSetRHSResidual()`
1537: @*/
1538: PetscErrorCode PetscDSGetRHSResidual(PetscDS ds, PetscInt f, void (**f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]), void (**f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
1539: {
1540: PetscPointFunc *tmp0, *tmp1;
1541: PetscInt n0, n1;
1543: PetscFunctionBegin;
1545: PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
1546: PetscCall(PetscWeakFormGetResidual(ds->wf, NULL, 0, f, 100, &n0, &tmp0, &n1, &tmp1));
1547: *f0 = tmp0 ? tmp0[0] : NULL;
1548: *f1 = tmp1 ? tmp1[0] : NULL;
1549: PetscFunctionReturn(PETSC_SUCCESS);
1550: }
1552: /*@C
1553: PetscDSSetRHSResidual - Set the pointwise residual function for explicit timestepping for a given test field
1555: Not Collective
1557: Input Parameters:
1558: + ds - The `PetscDS`
1559: . f - The test field number
1560: . f0 - integrand for the test function term
1561: - f1 - integrand for the test function gradient term
1563: Clling sequence for the callbacks f0 and f1:
1564: .vb
1565: f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1566: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1567: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1568: PetscReal t, const PetscReal x[], PetscScalar f0[])
1569: .ve
1570: + dim - the spatial dimension
1571: . Nf - the number of fields
1572: . uOff - the offset into u[] and u_t[] for each field
1573: . uOff_x - the offset into u_x[] for each field
1574: . u - each field evaluated at the current point
1575: . u_t - the time derivative of each field evaluated at the current point
1576: . u_x - the gradient of each field evaluated at the current point
1577: . aOff - the offset into a[] and a_t[] for each auxiliary field
1578: . aOff_x - the offset into a_x[] for each auxiliary field
1579: . a - each auxiliary field evaluated at the current point
1580: . a_t - the time derivative of each auxiliary field evaluated at the current point
1581: . a_x - the gradient of auxiliary each field evaluated at the current point
1582: . t - current time
1583: . x - coordinates of the current point
1584: . numConstants - number of constant parameters
1585: . constants - constant parameters
1586: - f0 - output values at the current point
1588: Level: intermediate
1590: Note:
1591: We are using a first order FEM model for the weak form: \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)
1593: .seealso: `PetscDS`, `PetscDSGetResidual()`
1594: @*/
1595: PetscErrorCode PetscDSSetRHSResidual(PetscDS ds, PetscInt f, void (*f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]), void (*f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
1596: {
1597: PetscFunctionBegin;
1601: PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
1602: PetscCall(PetscWeakFormSetIndexResidual(ds->wf, NULL, 0, f, 100, 0, f0, 0, f1));
1603: PetscFunctionReturn(PETSC_SUCCESS);
1604: }
1606: /*@C
1607: PetscDSHasJacobian - Checks that the Jacobian functions have been set
1609: Not Collective
1611: Input Parameter:
1612: . prob - The `PetscDS`
1614: Output Parameter:
1615: . hasJac - flag that pointwise function for the Jacobian has been set
1617: Level: intermediate
1619: .seealso: `PetscDS`, `PetscDSGetJacobianPreconditioner()`, `PetscDSSetJacobianPreconditioner()`, `PetscDSGetJacobian()`
1620: @*/
1621: PetscErrorCode PetscDSHasJacobian(PetscDS ds, PetscBool *hasJac)
1622: {
1623: PetscFunctionBegin;
1625: PetscCall(PetscWeakFormHasJacobian(ds->wf, hasJac));
1626: PetscFunctionReturn(PETSC_SUCCESS);
1627: }
1629: /*@C
1630: PetscDSGetJacobian - Get the pointwise Jacobian function for given test and basis field
1632: Not Collective
1634: Input Parameters:
1635: + ds - The `PetscDS`
1636: . f - The test field number
1637: - g - The field number
1639: Output Parameters:
1640: + g0 - integrand for the test and basis function term
1641: . g1 - integrand for the test function and basis function gradient term
1642: . g2 - integrand for the test function gradient and basis function term
1643: - g3 - integrand for the test function gradient and basis function gradient term
1645: Calling sequence of `g0`, `g1`, `g2` and `g3`:
1646: .vb
1647: void g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1648: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1649: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1650: PetscReal t, const PetscReal u_tShift, const PetscReal x[], PetscScalar g0[])
1651: .ve
1652: + dim - the spatial dimension
1653: . Nf - the number of fields
1654: . uOff - the offset into u[] and u_t[] for each field
1655: . uOff_x - the offset into u_x[] for each field
1656: . u - each field evaluated at the current point
1657: . u_t - the time derivative of each field evaluated at the current point
1658: . u_x - the gradient of each field evaluated at the current point
1659: . aOff - the offset into a[] and a_t[] for each auxiliary field
1660: . aOff_x - the offset into a_x[] for each auxiliary field
1661: . a - each auxiliary field evaluated at the current point
1662: . a_t - the time derivative of each auxiliary field evaluated at the current point
1663: . a_x - the gradient of auxiliary each field evaluated at the current point
1664: . t - current time
1665: . u_tShift - the multiplier a for dF/dU_t
1666: . x - coordinates of the current point
1667: . numConstants - number of constant parameters
1668: . constants - constant parameters
1669: - g0 - output values at the current point
1671: Level: intermediate
1673: Note:
1674: We are using a first order FEM model for the weak form:
1675: \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
1677: .seealso: `PetscDS`, `PetscDSSetJacobian()`
1678: @*/
1679: PetscErrorCode PetscDSGetJacobian(PetscDS ds, PetscInt f, PetscInt g, void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]), void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]), void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]), void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1680: {
1681: PetscPointJac *tmp0, *tmp1, *tmp2, *tmp3;
1682: PetscInt n0, n1, n2, n3;
1684: PetscFunctionBegin;
1686: PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
1687: PetscCheck(!(g < 0) && !(g >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", g, ds->Nf);
1688: PetscCall(PetscWeakFormGetJacobian(ds->wf, NULL, 0, f, g, 0, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3));
1689: *g0 = tmp0 ? tmp0[0] : NULL;
1690: *g1 = tmp1 ? tmp1[0] : NULL;
1691: *g2 = tmp2 ? tmp2[0] : NULL;
1692: *g3 = tmp3 ? tmp3[0] : NULL;
1693: PetscFunctionReturn(PETSC_SUCCESS);
1694: }
1696: /*@C
1697: PetscDSSetJacobian - Set the pointwise Jacobian function for given test and basis fields
1699: Not Collective
1701: Input Parameters:
1702: + ds - The `PetscDS`
1703: . f - The test field number
1704: . g - The field number
1705: . g0 - integrand for the test and basis function term
1706: . g1 - integrand for the test function and basis function gradient term
1707: . g2 - integrand for the test function gradient and basis function term
1708: - g3 - integrand for the test function gradient and basis function gradient term
1710: Calling sequence of `g0`, `g1`, `g2` and `g3`:
1711: .vb
1712: void g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1713: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1714: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1715: PetscReal t, const PetscReal x[], PetscScalar g0[])
1716: .ve
1717: + dim - the spatial dimension
1718: . Nf - the number of fields
1719: . uOff - the offset into u[] and u_t[] for each field
1720: . uOff_x - the offset into u_x[] for each field
1721: . u - each field evaluated at the current point
1722: . u_t - the time derivative of each field evaluated at the current point
1723: . u_x - the gradient of each field evaluated at the current point
1724: . aOff - the offset into a[] and a_t[] for each auxiliary field
1725: . aOff_x - the offset into a_x[] for each auxiliary field
1726: . a - each auxiliary field evaluated at the current point
1727: . a_t - the time derivative of each auxiliary field evaluated at the current point
1728: . a_x - the gradient of auxiliary each field evaluated at the current point
1729: . t - current time
1730: . u_tShift - the multiplier a for dF/dU_t
1731: . x - coordinates of the current point
1732: . numConstants - number of constant parameters
1733: . constants - constant parameters
1734: - g0 - output values at the current point
1736: Level: intermediate
1738: Note:
1739: We are using a first order FEM model for the weak form:
1740: \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
1742: .seealso: `PetscDS`, `PetscDSGetJacobian()`
1743: @*/
1744: PetscErrorCode PetscDSSetJacobian(PetscDS ds, PetscInt f, PetscInt g, void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]), void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]), void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]), void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1745: {
1746: PetscFunctionBegin;
1752: PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
1753: PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", g);
1754: PetscCall(PetscWeakFormSetIndexJacobian(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3));
1755: PetscFunctionReturn(PETSC_SUCCESS);
1756: }
1758: /*@C
1759: PetscDSUseJacobianPreconditioner - Set whether to construct a Jacobian preconditioner
1761: Not Collective
1763: Input Parameters:
1764: + prob - The `PetscDS`
1765: - useJacPre - flag that enables construction of a Jacobian preconditioner
1767: Level: intermediate
1769: .seealso: `PetscDS`, `PetscDSGetJacobianPreconditioner()`, `PetscDSSetJacobianPreconditioner()`, `PetscDSGetJacobian()`
1770: @*/
1771: PetscErrorCode PetscDSUseJacobianPreconditioner(PetscDS prob, PetscBool useJacPre)
1772: {
1773: PetscFunctionBegin;
1775: prob->useJacPre = useJacPre;
1776: PetscFunctionReturn(PETSC_SUCCESS);
1777: }
1779: /*@C
1780: PetscDSHasJacobianPreconditioner - Checks if a Jacobian preconditioner matrix has been set
1782: Not Collective
1784: Input Parameter:
1785: . prob - The `PetscDS`
1787: Output Parameter:
1788: . hasJacPre - flag that pointwise function for Jacobian preconditioner matrix has been set
1790: Level: intermediate
1792: .seealso: `PetscDS`, `PetscDSGetJacobianPreconditioner()`, `PetscDSSetJacobianPreconditioner()`, `PetscDSGetJacobian()`
1793: @*/
1794: PetscErrorCode PetscDSHasJacobianPreconditioner(PetscDS ds, PetscBool *hasJacPre)
1795: {
1796: PetscFunctionBegin;
1798: *hasJacPre = PETSC_FALSE;
1799: if (!ds->useJacPre) PetscFunctionReturn(PETSC_SUCCESS);
1800: PetscCall(PetscWeakFormHasJacobianPreconditioner(ds->wf, hasJacPre));
1801: PetscFunctionReturn(PETSC_SUCCESS);
1802: }
1804: /*@C
1805: PetscDSGetJacobianPreconditioner - Get the pointwise Jacobian preconditioner function for given test and basis field. If this is missing,
1806: the system matrix is used to build the preconditioner.
1808: Not Collective
1810: Input Parameters:
1811: + ds - The `PetscDS`
1812: . f - The test field number
1813: - g - The field number
1815: Output Parameters:
1816: + g0 - integrand for the test and basis function term
1817: . g1 - integrand for the test function and basis function gradient term
1818: . g2 - integrand for the test function gradient and basis function term
1819: - g3 - integrand for the test function gradient and basis function gradient term
1821: Calling sequence of `g0`, `g1`, `g2` and `g3`:
1822: .vb
1823: void g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1824: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1825: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1826: PetscReal t, const PetscReal u_tShift, const PetscReal x[], PetscScalar g0[])
1827: .ve
1828: + dim - the spatial dimension
1829: . Nf - the number of fields
1830: . uOff - the offset into u[] and u_t[] for each field
1831: . uOff_x - the offset into u_x[] for each field
1832: . u - each field evaluated at the current point
1833: . u_t - the time derivative of each field evaluated at the current point
1834: . u_x - the gradient of each field evaluated at the current point
1835: . aOff - the offset into a[] and a_t[] for each auxiliary field
1836: . aOff_x - the offset into a_x[] for each auxiliary field
1837: . a - each auxiliary field evaluated at the current point
1838: . a_t - the time derivative of each auxiliary field evaluated at the current point
1839: . a_x - the gradient of auxiliary each field evaluated at the current point
1840: . t - current time
1841: . u_tShift - the multiplier a for dF/dU_t
1842: . x - coordinates of the current point
1843: . numConstants - number of constant parameters
1844: . constants - constant parameters
1845: - g0 - output values at the current point
1847: Level: intermediate
1849: Note:
1850: We are using a first order FEM model for the weak form:
1851: \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
1853: .seealso: `PetscDS`, `PetscDSSetJacobianPreconditioner()`, `PetscDSGetJacobian()`
1854: @*/
1855: PetscErrorCode PetscDSGetJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g, void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]), void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]), void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]), void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1856: {
1857: PetscPointJac *tmp0, *tmp1, *tmp2, *tmp3;
1858: PetscInt n0, n1, n2, n3;
1860: PetscFunctionBegin;
1862: PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
1863: PetscCheck(!(g < 0) && !(g >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", g, ds->Nf);
1864: PetscCall(PetscWeakFormGetJacobianPreconditioner(ds->wf, NULL, 0, f, g, 0, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3));
1865: *g0 = tmp0 ? tmp0[0] : NULL;
1866: *g1 = tmp1 ? tmp1[0] : NULL;
1867: *g2 = tmp2 ? tmp2[0] : NULL;
1868: *g3 = tmp3 ? tmp3[0] : NULL;
1869: PetscFunctionReturn(PETSC_SUCCESS);
1870: }
1872: /*@C
1873: PetscDSSetJacobianPreconditioner - Set the pointwise Jacobian preconditioner function for given test and basis fields.
1874: If this is missing, the system matrix is used to build the preconditioner.
1876: Not Collective
1878: Input Parameters:
1879: + ds - The `PetscDS`
1880: . f - The test field number
1881: . g - The field number
1882: . g0 - integrand for the test and basis function term
1883: . g1 - integrand for the test function and basis function gradient term
1884: . g2 - integrand for the test function gradient and basis function term
1885: - g3 - integrand for the test function gradient and basis function gradient term
1887: Calling sequence of `g0`, `g1`, `g2` and `g3`:
1888: .vb
1889: void g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1890: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1891: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1892: PetscReal t, const PetscReal x[], PetscScalar g0[])
1893: .ve
1894: + dim - the spatial dimension
1895: . Nf - the number of fields
1896: . uOff - the offset into u[] and u_t[] for each field
1897: . uOff_x - the offset into u_x[] for each field
1898: . u - each field evaluated at the current point
1899: . u_t - the time derivative of each field evaluated at the current point
1900: . u_x - the gradient of each field evaluated at the current point
1901: . aOff - the offset into a[] and a_t[] for each auxiliary field
1902: . aOff_x - the offset into a_x[] for each auxiliary field
1903: . a - each auxiliary field evaluated at the current point
1904: . a_t - the time derivative of each auxiliary field evaluated at the current point
1905: . a_x - the gradient of auxiliary each field evaluated at the current point
1906: . t - current time
1907: . u_tShift - the multiplier a for dF/dU_t
1908: . x - coordinates of the current point
1909: . numConstants - number of constant parameters
1910: . constants - constant parameters
1911: - g0 - output values at the current point
1913: Level: intermediate
1915: Note:
1916: We are using a first order FEM model for the weak form:
1917: \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
1919: .seealso: `PetscDS`, `PetscDSGetJacobianPreconditioner()`, `PetscDSSetJacobian()`
1920: @*/
1921: PetscErrorCode PetscDSSetJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g, void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]), void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]), void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]), void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
1922: {
1923: PetscFunctionBegin;
1929: PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
1930: PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", g);
1931: PetscCall(PetscWeakFormSetIndexJacobianPreconditioner(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3));
1932: PetscFunctionReturn(PETSC_SUCCESS);
1933: }
1935: /*@C
1936: PetscDSHasDynamicJacobian - Signals that a dynamic Jacobian, dF/du_t, has been set
1938: Not Collective
1940: Input Parameter:
1941: . ds - The `PetscDS`
1943: Output Parameter:
1944: . hasDynJac - flag that pointwise function for dynamic Jacobian has been set
1946: Level: intermediate
1948: .seealso: `PetscDS`, `PetscDSGetDynamicJacobian()`, `PetscDSSetDynamicJacobian()`, `PetscDSGetJacobian()`
1949: @*/
1950: PetscErrorCode PetscDSHasDynamicJacobian(PetscDS ds, PetscBool *hasDynJac)
1951: {
1952: PetscFunctionBegin;
1954: PetscCall(PetscWeakFormHasDynamicJacobian(ds->wf, hasDynJac));
1955: PetscFunctionReturn(PETSC_SUCCESS);
1956: }
1958: /*@C
1959: PetscDSGetDynamicJacobian - Get the pointwise dynamic Jacobian, dF/du_t, function for given test and basis field
1961: Not Collective
1963: Input Parameters:
1964: + ds - The `PetscDS`
1965: . f - The test field number
1966: - g - The field number
1968: Output Parameters:
1969: + g0 - integrand for the test and basis function term
1970: . g1 - integrand for the test function and basis function gradient term
1971: . g2 - integrand for the test function gradient and basis function term
1972: - g3 - integrand for the test function gradient and basis function gradient term
1974: Calling sequence of `g0`, `g1`, `g2` and `g3`:
1975: .vb
1976: void g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1977: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1978: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1979: PetscReal t, const PetscReal u_tShift, const PetscReal x[], PetscScalar g0[])
1980: .ve
1981: + dim - the spatial dimension
1982: . Nf - the number of fields
1983: . uOff - the offset into u[] and u_t[] for each field
1984: . uOff_x - the offset into u_x[] for each field
1985: . u - each field evaluated at the current point
1986: . u_t - the time derivative of each field evaluated at the current point
1987: . u_x - the gradient of each field evaluated at the current point
1988: . aOff - the offset into a[] and a_t[] for each auxiliary field
1989: . aOff_x - the offset into a_x[] for each auxiliary field
1990: . a - each auxiliary field evaluated at the current point
1991: . a_t - the time derivative of each auxiliary field evaluated at the current point
1992: . a_x - the gradient of auxiliary each field evaluated at the current point
1993: . t - current time
1994: . u_tShift - the multiplier a for dF/dU_t
1995: . x - coordinates of the current point
1996: . numConstants - number of constant parameters
1997: . constants - constant parameters
1998: - g0 - output values at the current point
2000: Level: intermediate
2002: Note:
2003: We are using a first order FEM model for the weak form:
2004: \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
2006: .seealso: `PetscDS`, `PetscDSSetJacobian()`
2007: @*/
2008: PetscErrorCode PetscDSGetDynamicJacobian(PetscDS ds, PetscInt f, PetscInt g, void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]), void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]), void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]), void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
2009: {
2010: PetscPointJac *tmp0, *tmp1, *tmp2, *tmp3;
2011: PetscInt n0, n1, n2, n3;
2013: PetscFunctionBegin;
2015: PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
2016: PetscCheck(!(g < 0) && !(g >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", g, ds->Nf);
2017: PetscCall(PetscWeakFormGetDynamicJacobian(ds->wf, NULL, 0, f, g, 0, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3));
2018: *g0 = tmp0 ? tmp0[0] : NULL;
2019: *g1 = tmp1 ? tmp1[0] : NULL;
2020: *g2 = tmp2 ? tmp2[0] : NULL;
2021: *g3 = tmp3 ? tmp3[0] : NULL;
2022: PetscFunctionReturn(PETSC_SUCCESS);
2023: }
2025: /*@C
2026: PetscDSSetDynamicJacobian - Set the pointwise dynamic Jacobian, dF/du_t, function for given test and basis fields
2028: Not Collective
2030: Input Parameters:
2031: + ds - The `PetscDS`
2032: . f - The test field number
2033: . g - The field number
2034: . g0 - integrand for the test and basis function term
2035: . g1 - integrand for the test function and basis function gradient term
2036: . g2 - integrand for the test function gradient and basis function term
2037: - g3 - integrand for the test function gradient and basis function gradient term
2039: Calling sequence of `g0`, `g1`, `g2` and `g3`:
2040: .vb
2041: void g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2042: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2043: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2044: PetscReal t, const PetscReal x[], PetscScalar g0[])
2045: .ve
2046: + dim - the spatial dimension
2047: . Nf - the number of fields
2048: . uOff - the offset into u[] and u_t[] for each field
2049: . uOff_x - the offset into u_x[] for each field
2050: . u - each field evaluated at the current point
2051: . u_t - the time derivative of each field evaluated at the current point
2052: . u_x - the gradient of each field evaluated at the current point
2053: . aOff - the offset into a[] and a_t[] for each auxiliary field
2054: . aOff_x - the offset into a_x[] for each auxiliary field
2055: . a - each auxiliary field evaluated at the current point
2056: . a_t - the time derivative of each auxiliary field evaluated at the current point
2057: . a_x - the gradient of auxiliary each field evaluated at the current point
2058: . t - current time
2059: . u_tShift - the multiplier a for dF/dU_t
2060: . x - coordinates of the current point
2061: . numConstants - number of constant parameters
2062: . constants - constant parameters
2063: - g0 - output values at the current point
2065: Level: intermediate
2067: Note:
2068: We are using a first order FEM model for the weak form:
2069: \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
2071: .seealso: `PetscDS`, `PetscDSGetJacobian()`
2072: @*/
2073: PetscErrorCode PetscDSSetDynamicJacobian(PetscDS ds, PetscInt f, PetscInt g, void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]), void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]), void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]), void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
2074: {
2075: PetscFunctionBegin;
2081: PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2082: PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", g);
2083: PetscCall(PetscWeakFormSetIndexDynamicJacobian(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3));
2084: PetscFunctionReturn(PETSC_SUCCESS);
2085: }
2087: /*@C
2088: PetscDSGetRiemannSolver - Returns the Riemann solver for the given field
2090: Not Collective
2092: Input Parameters:
2093: + ds - The `PetscDS` object
2094: - f - The field number
2096: Output Parameter:
2097: . r - Riemann solver
2099: Calling sequence of `r`:
2100: .vb
2101: void r(PetscInt dim, PetscInt Nf, const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscScalar flux[], void *ctx)
2102: .ve
2103: + dim - The spatial dimension
2104: . Nf - The number of fields
2105: . x - The coordinates at a point on the interface
2106: . n - The normal vector to the interface
2107: . uL - The state vector to the left of the interface
2108: . uR - The state vector to the right of the interface
2109: . flux - output array of flux through the interface
2110: . numConstants - number of constant parameters
2111: . constants - constant parameters
2112: - ctx - optional user context
2114: Level: intermediate
2116: .seealso: `PetscDS`, `PetscDSSetRiemannSolver()`
2117: @*/
2118: PetscErrorCode PetscDSGetRiemannSolver(PetscDS ds, PetscInt f, 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))
2119: {
2120: PetscRiemannFunc *tmp;
2121: PetscInt n;
2123: PetscFunctionBegin;
2126: PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
2127: PetscCall(PetscWeakFormGetRiemannSolver(ds->wf, NULL, 0, f, 0, &n, &tmp));
2128: *r = tmp ? tmp[0] : NULL;
2129: PetscFunctionReturn(PETSC_SUCCESS);
2130: }
2132: /*@C
2133: PetscDSSetRiemannSolver - Sets the Riemann solver for the given field
2135: Not Collective
2137: Input Parameters:
2138: + ds - The `PetscDS` object
2139: . f - The field number
2140: - r - Riemann solver
2142: Calling sequence of `r`:
2143: .vb
2144: void r(PetscInt dim, PetscInt Nf, const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscScalar flux[], void *ctx)
2145: .ve
2146: + dim - The spatial dimension
2147: . Nf - The number of fields
2148: . x - The coordinates at a point on the interface
2149: . n - The normal vector to the interface
2150: . uL - The state vector to the left of the interface
2151: . uR - The state vector to the right of the interface
2152: . flux - output array of flux through the interface
2153: . numConstants - number of constant parameters
2154: . constants - constant parameters
2155: - ctx - optional user context
2157: Level: intermediate
2159: .seealso: `PetscDS`, `PetscDSGetRiemannSolver()`
2160: @*/
2161: PetscErrorCode PetscDSSetRiemannSolver(PetscDS ds, PetscInt f, 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))
2162: {
2163: PetscFunctionBegin;
2166: PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2167: PetscCall(PetscWeakFormSetIndexRiemannSolver(ds->wf, NULL, 0, f, 0, 0, r));
2168: PetscFunctionReturn(PETSC_SUCCESS);
2169: }
2171: /*@C
2172: PetscDSGetUpdate - Get the pointwise update function for a given field
2174: Not Collective
2176: Input Parameters:
2177: + ds - The `PetscDS`
2178: - f - The field number
2180: Output Parameter:
2181: . update - update function
2183: Calling sequence of `update`:
2184: .vb
2185: void update(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2186: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2187: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2188: PetscReal t, const PetscReal x[], PetscScalar uNew[])
2189: .ve
2190: + dim - the spatial dimension
2191: . Nf - the number of fields
2192: . uOff - the offset into u[] and u_t[] for each field
2193: . uOff_x - the offset into u_x[] for each field
2194: . u - each field evaluated at the current point
2195: . u_t - the time derivative of each field evaluated at the current point
2196: . u_x - the gradient of each field evaluated at the current point
2197: . aOff - the offset into a[] and a_t[] for each auxiliary field
2198: . aOff_x - the offset into a_x[] for each auxiliary field
2199: . a - each auxiliary field evaluated at the current point
2200: . a_t - the time derivative of each auxiliary field evaluated at the current point
2201: . a_x - the gradient of auxiliary each field evaluated at the current point
2202: . t - current time
2203: . x - coordinates of the current point
2204: - uNew - new value for field at the current point
2206: Level: intermediate
2208: .seealso: `PetscDS`, `PetscDSSetUpdate()`, `PetscDSSetResidual()`
2209: @*/
2210: PetscErrorCode PetscDSGetUpdate(PetscDS ds, PetscInt f, void (**update)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar uNew[]))
2211: {
2212: PetscFunctionBegin;
2214: PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
2215: if (update) {
2217: *update = ds->update[f];
2218: }
2219: PetscFunctionReturn(PETSC_SUCCESS);
2220: }
2222: /*@C
2223: PetscDSSetUpdate - Set the pointwise update function for a given field
2225: Not Collective
2227: Input Parameters:
2228: + ds - The `PetscDS`
2229: . f - The field number
2230: - update - update function
2232: Calling sequence of `update`:
2233: .vb
2234: void update(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2235: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2236: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2237: PetscReal t, const PetscReal x[], PetscScalar uNew[])
2238: .ve
2239: + dim - the spatial dimension
2240: . Nf - the number of fields
2241: . uOff - the offset into u[] and u_t[] for each field
2242: . uOff_x - the offset into u_x[] for each field
2243: . u - each field evaluated at the current point
2244: . u_t - the time derivative of each field evaluated at the current point
2245: . u_x - the gradient of each field evaluated at the current point
2246: . aOff - the offset into a[] and a_t[] for each auxiliary field
2247: . aOff_x - the offset into a_x[] for each auxiliary field
2248: . a - each auxiliary field evaluated at the current point
2249: . a_t - the time derivative of each auxiliary field evaluated at the current point
2250: . a_x - the gradient of auxiliary each field evaluated at the current point
2251: . t - current time
2252: . x - coordinates of the current point
2253: - uNew - new field values at the current point
2255: Level: intermediate
2257: .seealso: `PetscDS`, `PetscDSGetResidual()`
2258: @*/
2259: PetscErrorCode PetscDSSetUpdate(PetscDS ds, PetscInt f, void (*update)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar uNew[]))
2260: {
2261: PetscFunctionBegin;
2264: PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2265: PetscCall(PetscDSEnlarge_Static(ds, f + 1));
2266: ds->update[f] = update;
2267: PetscFunctionReturn(PETSC_SUCCESS);
2268: }
2270: PetscErrorCode PetscDSGetContext(PetscDS ds, PetscInt f, void *ctx)
2271: {
2272: PetscFunctionBegin;
2274: PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
2276: *(void **)ctx = ds->ctx[f];
2277: PetscFunctionReturn(PETSC_SUCCESS);
2278: }
2280: PetscErrorCode PetscDSSetContext(PetscDS ds, PetscInt f, void *ctx)
2281: {
2282: PetscFunctionBegin;
2284: PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2285: PetscCall(PetscDSEnlarge_Static(ds, f + 1));
2286: ds->ctx[f] = ctx;
2287: PetscFunctionReturn(PETSC_SUCCESS);
2288: }
2290: /*@C
2291: PetscDSGetBdResidual - Get the pointwise boundary residual function for a given test field
2293: Not Collective
2295: Input Parameters:
2296: + ds - The PetscDS
2297: - f - The test field number
2299: Output Parameters:
2300: + f0 - boundary integrand for the test function term
2301: - f1 - boundary integrand for the test function gradient term
2303: Calling sequence of `f0` and `f1`:
2304: .vb
2305: void f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2306: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2307: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2308: PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[])
2309: .ve
2310: + dim - the spatial dimension
2311: . Nf - the number of fields
2312: . uOff - the offset into u[] and u_t[] for each field
2313: . uOff_x - the offset into u_x[] for each field
2314: . u - each field evaluated at the current point
2315: . u_t - the time derivative of each field evaluated at the current point
2316: . u_x - the gradient of each field evaluated at the current point
2317: . aOff - the offset into a[] and a_t[] for each auxiliary field
2318: . aOff_x - the offset into a_x[] for each auxiliary field
2319: . a - each auxiliary field evaluated at the current point
2320: . a_t - the time derivative of each auxiliary field evaluated at the current point
2321: . a_x - the gradient of auxiliary each field evaluated at the current point
2322: . t - current time
2323: . x - coordinates of the current point
2324: . n - unit normal at the current point
2325: . numConstants - number of constant parameters
2326: . constants - constant parameters
2327: - f0 - output values at the current point
2329: Level: intermediate
2331: Note:
2332: We are using a first order FEM model for the weak form:
2333: \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
2335: .seealso: `PetscDS`, `PetscDSSetBdResidual()`
2336: @*/
2337: PetscErrorCode PetscDSGetBdResidual(PetscDS ds, PetscInt f, void (**f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]), void (**f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
2338: {
2339: PetscBdPointFunc *tmp0, *tmp1;
2340: PetscInt n0, n1;
2342: PetscFunctionBegin;
2344: PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
2345: PetscCall(PetscWeakFormGetBdResidual(ds->wf, NULL, 0, f, 0, &n0, &tmp0, &n1, &tmp1));
2346: *f0 = tmp0 ? tmp0[0] : NULL;
2347: *f1 = tmp1 ? tmp1[0] : NULL;
2348: PetscFunctionReturn(PETSC_SUCCESS);
2349: }
2351: /*@C
2352: PetscDSSetBdResidual - Get the pointwise boundary residual function for a given test field
2354: Not Collective
2356: Input Parameters:
2357: + ds - The `PetscDS`
2358: . f - The test field number
2359: . f0 - boundary integrand for the test function term
2360: - f1 - boundary integrand for the test function gradient term
2362: Calling sequence of `f0` and `f1`:
2363: .vb
2364: void f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2365: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2366: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2367: PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[])
2368: .ve
2369: + dim - the spatial dimension
2370: . Nf - the number of fields
2371: . uOff - the offset into u[] and u_t[] for each field
2372: . uOff_x - the offset into u_x[] for each field
2373: . u - each field evaluated at the current point
2374: . u_t - the time derivative of each field evaluated at the current point
2375: . u_x - the gradient of each field evaluated at the current point
2376: . aOff - the offset into a[] and a_t[] for each auxiliary field
2377: . aOff_x - the offset into a_x[] for each auxiliary field
2378: . a - each auxiliary field evaluated at the current point
2379: . a_t - the time derivative of each auxiliary field evaluated at the current point
2380: . a_x - the gradient of auxiliary each field evaluated at the current point
2381: . t - current time
2382: . x - coordinates of the current point
2383: . n - unit normal at the current point
2384: . numConstants - number of constant parameters
2385: . constants - constant parameters
2386: - f0 - output values at the current point
2388: Level: intermediate
2390: Note:
2391: We are using a first order FEM model for the weak form:
2392: \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
2394: .seealso: `PetscDS`, `PetscDSGetBdResidual()`
2395: @*/
2396: PetscErrorCode PetscDSSetBdResidual(PetscDS ds, PetscInt f, void (*f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]), void (*f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
2397: {
2398: PetscFunctionBegin;
2400: PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2401: PetscCall(PetscWeakFormSetIndexBdResidual(ds->wf, NULL, 0, f, 0, 0, f0, 0, f1));
2402: PetscFunctionReturn(PETSC_SUCCESS);
2403: }
2405: /*@
2406: PetscDSHasBdJacobian - Indicates that boundary Jacobian functions have been set
2408: Not Collective
2410: Input Parameter:
2411: . ds - The `PetscDS`
2413: Output Parameter:
2414: . hasBdJac - flag that pointwise function for the boundary Jacobian has been set
2416: Level: intermediate
2418: .seealso: `PetscDS`, `PetscDSHasJacobian()`, `PetscDSSetBdJacobian()`, `PetscDSGetBdJacobian()`
2419: @*/
2420: PetscErrorCode PetscDSHasBdJacobian(PetscDS ds, PetscBool *hasBdJac)
2421: {
2422: PetscFunctionBegin;
2425: PetscCall(PetscWeakFormHasBdJacobian(ds->wf, hasBdJac));
2426: PetscFunctionReturn(PETSC_SUCCESS);
2427: }
2429: /*@C
2430: PetscDSGetBdJacobian - Get the pointwise boundary Jacobian function for given test and basis field
2432: Not Collective
2434: Input Parameters:
2435: + ds - The `PetscDS`
2436: . f - The test field number
2437: - g - The field number
2439: Output Parameters:
2440: + g0 - integrand for the test and basis function term
2441: . g1 - integrand for the test function and basis function gradient term
2442: . g2 - integrand for the test function gradient and basis function term
2443: - g3 - integrand for the test function gradient and basis function gradient term
2445: Calling sequence of `g0`, `g1`, `g2` and `g3`:
2446: .vb
2447: void g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2448: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2449: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2450: PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar g0[])
2451: .ve
2452: + dim - the spatial dimension
2453: . Nf - the number of fields
2454: . uOff - the offset into u[] and u_t[] for each field
2455: . uOff_x - the offset into u_x[] for each field
2456: . u - each field evaluated at the current point
2457: . u_t - the time derivative of each field evaluated at the current point
2458: . u_x - the gradient of each field evaluated at the current point
2459: . aOff - the offset into a[] and a_t[] for each auxiliary field
2460: . aOff_x - the offset into a_x[] for each auxiliary field
2461: . a - each auxiliary field evaluated at the current point
2462: . a_t - the time derivative of each auxiliary field evaluated at the current point
2463: . a_x - the gradient of auxiliary each field evaluated at the current point
2464: . t - current time
2465: . u_tShift - the multiplier a for dF/dU_t
2466: . x - coordinates of the current point
2467: . n - normal at the current point
2468: . numConstants - number of constant parameters
2469: . constants - constant parameters
2470: - g0 - output values at the current point
2472: Level: intermediate
2474: Note:
2475: We are using a first order FEM model for the weak form:
2476: \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
2478: .seealso: `PetscDS`, `PetscDSSetBdJacobian()`
2479: @*/
2480: PetscErrorCode PetscDSGetBdJacobian(PetscDS ds, PetscInt f, PetscInt g, void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]), void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]), void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]), void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 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;
2485: PetscFunctionBegin;
2487: PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
2488: PetscCheck(!(g < 0) && !(g >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", g, ds->Nf);
2489: PetscCall(PetscWeakFormGetBdJacobian(ds->wf, NULL, 0, f, g, 0, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3));
2490: *g0 = tmp0 ? tmp0[0] : NULL;
2491: *g1 = tmp1 ? tmp1[0] : NULL;
2492: *g2 = tmp2 ? tmp2[0] : NULL;
2493: *g3 = tmp3 ? tmp3[0] : NULL;
2494: PetscFunctionReturn(PETSC_SUCCESS);
2495: }
2497: /*@C
2498: PetscDSSetBdJacobian - Set the pointwise boundary Jacobian function for given test and basis field
2500: Not Collective
2502: Input Parameters:
2503: + ds - The PetscDS
2504: . f - The test field number
2505: . g - The field number
2506: . g0 - integrand for the test and basis function term
2507: . g1 - integrand for the test function and basis function gradient term
2508: . g2 - integrand for the test function gradient and basis function term
2509: - g3 - integrand for the test function gradient and basis function gradient term
2511: Calling sequence of `g0`, `g1`, `g2` and `g3`:
2512: .vb
2513: void g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2514: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2515: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2516: PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar g0[])
2517: .ve
2518: + dim - the spatial dimension
2519: . Nf - the number of fields
2520: . uOff - the offset into u[] and u_t[] for each field
2521: . uOff_x - the offset into u_x[] for each field
2522: . u - each field evaluated at the current point
2523: . u_t - the time derivative of each field evaluated at the current point
2524: . u_x - the gradient of each field evaluated at the current point
2525: . aOff - the offset into a[] and a_t[] for each auxiliary field
2526: . aOff_x - the offset into a_x[] for each auxiliary field
2527: . a - each auxiliary field evaluated at the current point
2528: . a_t - the time derivative of each auxiliary field evaluated at the current point
2529: . a_x - the gradient of auxiliary each field evaluated at the current point
2530: . t - current time
2531: . u_tShift - the multiplier a for dF/dU_t
2532: . x - coordinates of the current point
2533: . n - normal at the current point
2534: . numConstants - number of constant parameters
2535: . constants - constant parameters
2536: - g0 - output values at the current point
2538: Level: intermediate
2540: Note:
2541: We are using a first order FEM model for the weak form:
2542: \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
2544: .seealso: `PetscDS`, `PetscDSGetBdJacobian()`
2545: @*/
2546: PetscErrorCode PetscDSSetBdJacobian(PetscDS ds, PetscInt f, PetscInt g, void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]), void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]), void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]), void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
2547: {
2548: PetscFunctionBegin;
2554: PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2555: PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", g);
2556: PetscCall(PetscWeakFormSetIndexBdJacobian(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3));
2557: PetscFunctionReturn(PETSC_SUCCESS);
2558: }
2560: /*@
2561: PetscDSHasBdJacobianPreconditioner - Signals that boundary Jacobian preconditioner functions have been set
2563: Not Collective
2565: Input Parameter:
2566: . ds - The `PetscDS`
2568: Output Parameter:
2569: . hasBdJac - flag that pointwise function for the boundary Jacobian preconditioner has been set
2571: Level: intermediate
2573: .seealso: `PetscDS`, `PetscDSHasJacobian()`, `PetscDSSetBdJacobian()`, `PetscDSGetBdJacobian()`
2574: @*/
2575: PetscErrorCode PetscDSHasBdJacobianPreconditioner(PetscDS ds, PetscBool *hasBdJacPre)
2576: {
2577: PetscFunctionBegin;
2580: PetscCall(PetscWeakFormHasBdJacobianPreconditioner(ds->wf, hasBdJacPre));
2581: PetscFunctionReturn(PETSC_SUCCESS);
2582: }
2584: /*@C
2585: PetscDSGetBdJacobianPreconditioner - Get the pointwise boundary Jacobian preconditioner function for given test and basis field
2587: Not Collective; No Fortran Support
2589: Input Parameters:
2590: + ds - The `PetscDS`
2591: . f - The test field number
2592: - g - The field number
2594: Output Parameters:
2595: + g0 - integrand for the test and basis function term
2596: . g1 - integrand for the test function and basis function gradient term
2597: . g2 - integrand for the test function gradient and basis function term
2598: - g3 - integrand for the test function gradient and basis function gradient term
2600: Calling sequence of `g0`, `g1`, `g2` and `g3`:
2601: .vb
2602: void g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2603: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2604: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2605: PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
2606: .ve
2607: + dim - the spatial dimension
2608: . Nf - the number of fields
2609: . NfAux - the number of auxiliary fields
2610: . uOff - the offset into u[] and u_t[] for each field
2611: . uOff_x - the offset into u_x[] for each field
2612: . u - each field evaluated at the current point
2613: . u_t - the time derivative of each field evaluated at the current point
2614: . u_x - the gradient of each field evaluated at the current point
2615: . aOff - the offset into a[] and a_t[] for each auxiliary field
2616: . aOff_x - the offset into a_x[] for each auxiliary field
2617: . a - each auxiliary field evaluated at the current point
2618: . a_t - the time derivative of each auxiliary field evaluated at the current point
2619: . a_x - the gradient of auxiliary each field evaluated at the current point
2620: . t - current time
2621: . u_tShift - the multiplier a for dF/dU_t
2622: . x - coordinates of the current point
2623: . n - normal at the current point
2624: . numConstants - number of constant parameters
2625: . constants - constant parameters
2626: - g0 - output values at the current point
2628: Level: intermediate
2630: Note:
2631: We are using a first order FEM model for the weak form:
2632: \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
2634: .seealso: `PetscDS`, `PetscDSSetBdJacobianPreconditioner()`
2635: @*/
2636: PetscErrorCode PetscDSGetBdJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g, void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]), void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]), void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]), void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
2637: {
2638: PetscBdPointJac *tmp0, *tmp1, *tmp2, *tmp3;
2639: PetscInt n0, n1, n2, n3;
2641: PetscFunctionBegin;
2643: PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
2644: PetscCheck(!(g < 0) && !(g >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", g, ds->Nf);
2645: PetscCall(PetscWeakFormGetBdJacobianPreconditioner(ds->wf, NULL, 0, f, g, 0, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3));
2646: *g0 = tmp0 ? tmp0[0] : NULL;
2647: *g1 = tmp1 ? tmp1[0] : NULL;
2648: *g2 = tmp2 ? tmp2[0] : NULL;
2649: *g3 = tmp3 ? tmp3[0] : NULL;
2650: PetscFunctionReturn(PETSC_SUCCESS);
2651: }
2653: /*@C
2654: PetscDSSetBdJacobianPreconditioner - Set the pointwise boundary Jacobian preconditioner function for given test and basis field
2656: Not Collective; No Fortran Support
2658: Input Parameters:
2659: + ds - The `PetscDS`
2660: . f - The test field number
2661: . g - The field number
2662: . g0 - integrand for the test and basis function term
2663: . g1 - integrand for the test function and basis function gradient term
2664: . g2 - integrand for the test function gradient and basis function term
2665: - g3 - integrand for the test function gradient and basis function gradient term
2667: Calling sequence of `g0`, `g1`, `g2` and `g3`:
2668: .vb
2669: void g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
2670: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
2671: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
2672: PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
2673: .ve
2674: + dim - the spatial dimension
2675: . Nf - the number of fields
2676: . NfAux - the number of auxiliary fields
2677: . uOff - the offset into u[] and u_t[] for each field
2678: . uOff_x - the offset into u_x[] for each field
2679: . u - each field evaluated at the current point
2680: . u_t - the time derivative of each field evaluated at the current point
2681: . u_x - the gradient of each field evaluated at the current point
2682: . aOff - the offset into a[] and a_t[] for each auxiliary field
2683: . aOff_x - the offset into a_x[] for each auxiliary field
2684: . a - each auxiliary field evaluated at the current point
2685: . a_t - the time derivative of each auxiliary field evaluated at the current point
2686: . a_x - the gradient of auxiliary each field evaluated at the current point
2687: . t - current time
2688: . u_tShift - the multiplier a for dF/dU_t
2689: . x - coordinates of the current point
2690: . n - normal at the current point
2691: . numConstants - number of constant parameters
2692: . constants - constant parameters
2693: - g0 - output values at the current point
2695: Level: intermediate
2697: Note:
2698: We are using a first order FEM model for the weak form:
2699: \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
2701: .seealso: `PetscDS`, `PetscDSGetBdJacobianPreconditioner()`
2702: @*/
2703: PetscErrorCode PetscDSSetBdJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g, void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]), void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]), void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]), void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
2704: {
2705: PetscFunctionBegin;
2711: PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2712: PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", g);
2713: PetscCall(PetscWeakFormSetIndexBdJacobianPreconditioner(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3));
2714: PetscFunctionReturn(PETSC_SUCCESS);
2715: }
2717: /*@C
2718: PetscDSGetExactSolution - Get the pointwise exact solution function for a given test field
2720: Not Collective
2722: Input Parameters:
2723: + prob - The PetscDS
2724: - f - The test field number
2726: Output Parameters:
2727: + exactSol - exact solution for the test field
2728: - exactCtx - exact solution context
2730: Calling sequence of `exactSol`:
2731: .vb
2732: PetscErrorCode sol(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx)
2733: .ve
2734: + dim - the spatial dimension
2735: . t - current time
2736: . x - coordinates of the current point
2737: . Nc - the number of field components
2738: . u - the solution field evaluated at the current point
2739: - ctx - a user context
2741: Level: intermediate
2743: .seealso: `PetscDS`, `PetscDSSetExactSolution()`, `PetscDSGetExactSolutionTimeDerivative()`
2744: @*/
2745: PetscErrorCode PetscDSGetExactSolution(PetscDS prob, PetscInt f, PetscErrorCode (**sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void **ctx)
2746: {
2747: PetscFunctionBegin;
2749: PetscCheck(!(f < 0) && !(f >= prob->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, prob->Nf);
2750: if (sol) {
2752: *sol = prob->exactSol[f];
2753: }
2754: if (ctx) {
2756: *ctx = prob->exactCtx[f];
2757: }
2758: PetscFunctionReturn(PETSC_SUCCESS);
2759: }
2761: /*@C
2762: PetscDSSetExactSolution - Set the pointwise exact solution function for a given test field
2764: Not Collective
2766: Input Parameters:
2767: + prob - The `PetscDS`
2768: . f - The test field number
2769: . sol - solution function for the test fields
2770: - ctx - solution context or `NULL`
2772: Calling sequence of `sol`:
2773: .vb
2774: PetscErrorCode sol(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx)
2775: .ve
2776: + dim - the spatial dimension
2777: . t - current time
2778: . x - coordinates of the current point
2779: . Nc - the number of field components
2780: . u - the solution field evaluated at the current point
2781: - ctx - a user context
2783: Level: intermediate
2785: .seealso: `PetscDS`, `PetscDSGetExactSolution()`
2786: @*/
2787: PetscErrorCode PetscDSSetExactSolution(PetscDS prob, PetscInt f, PetscErrorCode (*sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void *ctx)
2788: {
2789: PetscFunctionBegin;
2791: PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2792: PetscCall(PetscDSEnlarge_Static(prob, f + 1));
2793: if (sol) {
2795: prob->exactSol[f] = sol;
2796: }
2797: if (ctx) {
2799: prob->exactCtx[f] = ctx;
2800: }
2801: PetscFunctionReturn(PETSC_SUCCESS);
2802: }
2804: /*@C
2805: PetscDSGetExactSolutionTimeDerivative - Get the pointwise time derivative of the exact solution function for a given test field
2807: Not Collective
2809: Input Parameters:
2810: + prob - The `PetscDS`
2811: - f - The test field number
2813: Output Parameters:
2814: + exactSol - time derivative of the exact solution for the test field
2815: - exactCtx - time derivative of the exact solution context
2817: Calling sequence of `exactSol`:
2818: .vb
2819: PetscErrorCode sol(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx)
2820: .ve
2821: + dim - the spatial dimension
2822: . t - current time
2823: . x - coordinates of the current point
2824: . Nc - the number of field components
2825: . u - the solution field evaluated at the current point
2826: - ctx - a user context
2828: Level: intermediate
2830: .seealso: `PetscDS`, `PetscDSSetExactSolutionTimeDerivative()`, `PetscDSGetExactSolution()`
2831: @*/
2832: PetscErrorCode PetscDSGetExactSolutionTimeDerivative(PetscDS prob, PetscInt f, PetscErrorCode (**sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void **ctx)
2833: {
2834: PetscFunctionBegin;
2836: PetscCheck(!(f < 0) && !(f >= prob->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, prob->Nf);
2837: if (sol) {
2839: *sol = prob->exactSol_t[f];
2840: }
2841: if (ctx) {
2843: *ctx = prob->exactCtx_t[f];
2844: }
2845: PetscFunctionReturn(PETSC_SUCCESS);
2846: }
2848: /*@C
2849: PetscDSSetExactSolutionTimeDerivative - Set the pointwise time derivative of the exact solution function for a given test field
2851: Not Collective
2853: Input Parameters:
2854: + prob - The `PetscDS`
2855: . f - The test field number
2856: . sol - time derivative of the solution function for the test fields
2857: - ctx - time derivative of the solution context or `NULL`
2859: Calling sequence of `sol`:
2860: .vb
2861: PetscErrorCode sol(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx)
2862: .ve
2863: + dim - the spatial dimension
2864: . t - current time
2865: . x - coordinates of the current point
2866: . Nc - the number of field components
2867: . u - the solution field evaluated at the current point
2868: - ctx - a user context
2870: Level: intermediate
2872: .seealso: `PetscDS`, `PetscDSGetExactSolutionTimeDerivative()`, `PetscDSSetExactSolution()`
2873: @*/
2874: PetscErrorCode PetscDSSetExactSolutionTimeDerivative(PetscDS prob, PetscInt f, PetscErrorCode (*sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void *ctx)
2875: {
2876: PetscFunctionBegin;
2878: PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
2879: PetscCall(PetscDSEnlarge_Static(prob, f + 1));
2880: if (sol) {
2882: prob->exactSol_t[f] = sol;
2883: }
2884: if (ctx) {
2886: prob->exactCtx_t[f] = ctx;
2887: }
2888: PetscFunctionReturn(PETSC_SUCCESS);
2889: }
2891: /*@C
2892: PetscDSGetConstants - Returns the array of constants passed to point functions
2894: Not Collective
2896: Input Parameter:
2897: . prob - The `PetscDS` object
2899: Output Parameters:
2900: + numConstants - The number of constants
2901: - constants - The array of constants, NULL if there are none
2903: Level: intermediate
2905: .seealso: `PetscDS`, `PetscDSSetConstants()`, `PetscDSCreate()`
2906: @*/
2907: PetscErrorCode PetscDSGetConstants(PetscDS prob, PetscInt *numConstants, const PetscScalar *constants[])
2908: {
2909: PetscFunctionBegin;
2911: if (numConstants) {
2913: *numConstants = prob->numConstants;
2914: }
2915: if (constants) {
2917: *constants = prob->constants;
2918: }
2919: PetscFunctionReturn(PETSC_SUCCESS);
2920: }
2922: /*@C
2923: PetscDSSetConstants - Set the array of constants passed to point functions
2925: Not Collective
2927: Input Parameters:
2928: + prob - The `PetscDS` object
2929: . numConstants - The number of constants
2930: - constants - The array of constants, NULL if there are none
2932: Level: intermediate
2934: .seealso: `PetscDS`, `PetscDSGetConstants()`, `PetscDSCreate()`
2935: @*/
2936: PetscErrorCode PetscDSSetConstants(PetscDS prob, PetscInt numConstants, PetscScalar constants[])
2937: {
2938: PetscFunctionBegin;
2940: if (numConstants != prob->numConstants) {
2941: PetscCall(PetscFree(prob->constants));
2942: prob->numConstants = numConstants;
2943: if (prob->numConstants) {
2944: PetscCall(PetscMalloc1(prob->numConstants, &prob->constants));
2945: } else {
2946: prob->constants = NULL;
2947: }
2948: }
2949: if (prob->numConstants) {
2951: PetscCall(PetscArraycpy(prob->constants, constants, prob->numConstants));
2952: }
2953: PetscFunctionReturn(PETSC_SUCCESS);
2954: }
2956: /*@
2957: PetscDSGetFieldIndex - Returns the index of the given field
2959: Not Collective
2961: Input Parameters:
2962: + prob - The `PetscDS` object
2963: - disc - The discretization object
2965: Output Parameter:
2966: . f - The field number
2968: Level: beginner
2970: .seealso: `PetscDS`, `PetscGetDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
2971: @*/
2972: PetscErrorCode PetscDSGetFieldIndex(PetscDS prob, PetscObject disc, PetscInt *f)
2973: {
2974: PetscInt g;
2976: PetscFunctionBegin;
2979: *f = -1;
2980: for (g = 0; g < prob->Nf; ++g) {
2981: if (disc == prob->disc[g]) break;
2982: }
2983: PetscCheck(g != prob->Nf, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONG, "Field not found in PetscDS.");
2984: *f = g;
2985: PetscFunctionReturn(PETSC_SUCCESS);
2986: }
2988: /*@
2989: PetscDSGetFieldSize - Returns the size of the given field in the full space basis
2991: Not Collective
2993: Input Parameters:
2994: + prob - The `PetscDS` object
2995: - f - The field number
2997: Output Parameter:
2998: . size - The size
3000: Level: beginner
3002: .seealso: `PetscDS`, `PetscDSGetFieldOffset()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3003: @*/
3004: PetscErrorCode PetscDSGetFieldSize(PetscDS prob, PetscInt f, PetscInt *size)
3005: {
3006: PetscFunctionBegin;
3009: PetscCheck(!(f < 0) && !(f >= prob->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, prob->Nf);
3010: PetscCall(PetscDSSetUp(prob));
3011: *size = prob->Nb[f];
3012: PetscFunctionReturn(PETSC_SUCCESS);
3013: }
3015: /*@
3016: PetscDSGetFieldOffset - Returns the offset of the given field in the full space basis
3018: Not Collective
3020: Input Parameters:
3021: + prob - The `PetscDS` object
3022: - f - The field number
3024: Output Parameter:
3025: . off - The offset
3027: Level: beginner
3029: .seealso: `PetscDS`, `PetscDSGetFieldSize()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3030: @*/
3031: PetscErrorCode PetscDSGetFieldOffset(PetscDS prob, PetscInt f, PetscInt *off)
3032: {
3033: PetscInt size, g;
3035: PetscFunctionBegin;
3038: PetscCheck(!(f < 0) && !(f >= prob->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, prob->Nf);
3039: *off = 0;
3040: for (g = 0; g < f; ++g) {
3041: PetscCall(PetscDSGetFieldSize(prob, g, &size));
3042: *off += size;
3043: }
3044: PetscFunctionReturn(PETSC_SUCCESS);
3045: }
3047: /*@
3048: PetscDSGetFieldOffsetCohesive - Returns the offset of the given field in the full space basis on a cohesive cell
3050: Not Collective
3052: Input Parameters:
3053: + prob - The `PetscDS` object
3054: - f - The field number
3056: Output Parameter:
3057: . off - The offset
3059: Level: beginner
3061: .seealso: `PetscDS`, `PetscDSGetFieldSize()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3062: @*/
3063: PetscErrorCode PetscDSGetFieldOffsetCohesive(PetscDS ds, PetscInt f, PetscInt *off)
3064: {
3065: PetscInt size, g;
3067: PetscFunctionBegin;
3070: PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
3071: *off = 0;
3072: for (g = 0; g < f; ++g) {
3073: PetscBool cohesive;
3075: PetscCall(PetscDSGetCohesive(ds, g, &cohesive));
3076: PetscCall(PetscDSGetFieldSize(ds, g, &size));
3077: *off += cohesive ? size : size * 2;
3078: }
3079: PetscFunctionReturn(PETSC_SUCCESS);
3080: }
3082: /*@
3083: PetscDSGetDimensions - Returns the size of the approximation space for each field on an evaluation point
3085: Not Collective
3087: Input Parameter:
3088: . prob - The `PetscDS` object
3090: Output Parameter:
3091: . dimensions - The number of dimensions
3093: Level: beginner
3095: .seealso: `PetscDS`, `PetscDSGetComponentOffsets()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3096: @*/
3097: PetscErrorCode PetscDSGetDimensions(PetscDS prob, PetscInt *dimensions[])
3098: {
3099: PetscFunctionBegin;
3101: PetscCall(PetscDSSetUp(prob));
3103: *dimensions = prob->Nb;
3104: PetscFunctionReturn(PETSC_SUCCESS);
3105: }
3107: /*@
3108: PetscDSGetComponents - Returns the number of components for each field on an evaluation point
3110: Not Collective
3112: Input Parameter:
3113: . prob - The `PetscDS` object
3115: Output Parameter:
3116: . components - The number of components
3118: Level: beginner
3120: .seealso: `PetscDS`, `PetscDSGetComponentOffsets()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3121: @*/
3122: PetscErrorCode PetscDSGetComponents(PetscDS prob, PetscInt *components[])
3123: {
3124: PetscFunctionBegin;
3126: PetscCall(PetscDSSetUp(prob));
3128: *components = prob->Nc;
3129: PetscFunctionReturn(PETSC_SUCCESS);
3130: }
3132: /*@
3133: PetscDSGetComponentOffset - Returns the offset of the given field on an evaluation point
3135: Not Collective
3137: Input Parameters:
3138: + prob - The `PetscDS` object
3139: - f - The field number
3141: Output Parameter:
3142: . off - The offset
3144: Level: beginner
3146: .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3147: @*/
3148: PetscErrorCode PetscDSGetComponentOffset(PetscDS prob, PetscInt f, PetscInt *off)
3149: {
3150: PetscFunctionBegin;
3153: PetscCheck(!(f < 0) && !(f >= prob->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, prob->Nf);
3154: PetscCall(PetscDSSetUp(prob));
3155: *off = prob->off[f];
3156: PetscFunctionReturn(PETSC_SUCCESS);
3157: }
3159: /*@
3160: PetscDSGetComponentOffsets - Returns the offset of each field on an evaluation point
3162: Not Collective
3164: Input Parameter:
3165: . prob - The `PetscDS` object
3167: Output Parameter:
3168: . offsets - The offsets
3170: Level: beginner
3172: .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3173: @*/
3174: PetscErrorCode PetscDSGetComponentOffsets(PetscDS prob, PetscInt *offsets[])
3175: {
3176: PetscFunctionBegin;
3179: PetscCall(PetscDSSetUp(prob));
3180: *offsets = prob->off;
3181: PetscFunctionReturn(PETSC_SUCCESS);
3182: }
3184: /*@
3185: PetscDSGetComponentDerivativeOffsets - Returns the offset of each field derivative on an evaluation point
3187: Not Collective
3189: Input Parameter:
3190: . prob - The `PetscDS` object
3192: Output Parameter:
3193: . offsets - The offsets
3195: Level: beginner
3197: .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3198: @*/
3199: PetscErrorCode PetscDSGetComponentDerivativeOffsets(PetscDS prob, PetscInt *offsets[])
3200: {
3201: PetscFunctionBegin;
3204: PetscCall(PetscDSSetUp(prob));
3205: *offsets = prob->offDer;
3206: PetscFunctionReturn(PETSC_SUCCESS);
3207: }
3209: /*@
3210: PetscDSGetComponentOffsetsCohesive - Returns the offset of each field on an evaluation point
3212: Not Collective
3214: Input Parameters:
3215: + ds - The `PetscDS` object
3216: - s - The cohesive side, 0 for negative, 1 for positive, 2 for cohesive
3218: Output Parameter:
3219: . offsets - The offsets
3221: Level: beginner
3223: .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3224: @*/
3225: PetscErrorCode PetscDSGetComponentOffsetsCohesive(PetscDS ds, PetscInt s, PetscInt *offsets[])
3226: {
3227: PetscFunctionBegin;
3230: PetscCheck(ds->isCohesive, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cohesive offsets are only valid for a cohesive DS");
3231: PetscCheck(!(s < 0) && !(s > 2), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cohesive side %" PetscInt_FMT " is not in [0, 2]", s);
3232: PetscCall(PetscDSSetUp(ds));
3233: *offsets = ds->offCohesive[s];
3234: PetscFunctionReturn(PETSC_SUCCESS);
3235: }
3237: /*@
3238: PetscDSGetComponentDerivativeOffsetsCohesive - Returns the offset of each field derivative on an evaluation point
3240: Not Collective
3242: Input Parameters:
3243: + ds - The `PetscDS` object
3244: - s - The cohesive side, 0 for negative, 1 for positive, 2 for cohesive
3246: Output Parameter:
3247: . offsets - The offsets
3249: Level: beginner
3251: .seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
3252: @*/
3253: PetscErrorCode PetscDSGetComponentDerivativeOffsetsCohesive(PetscDS ds, PetscInt s, PetscInt *offsets[])
3254: {
3255: PetscFunctionBegin;
3258: PetscCheck(ds->isCohesive, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cohesive offsets are only valid for a cohesive DS");
3259: PetscCheck(!(s < 0) && !(s > 2), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cohesive side %" PetscInt_FMT " is not in [0, 2]", s);
3260: PetscCall(PetscDSSetUp(ds));
3261: *offsets = ds->offDerCohesive[s];
3262: PetscFunctionReturn(PETSC_SUCCESS);
3263: }
3265: /*@C
3266: PetscDSGetTabulation - Return the basis tabulation at quadrature points for the volume discretization
3268: Not Collective
3270: Input Parameter:
3271: . prob - The `PetscDS` object
3273: Output Parameter:
3274: . T - The basis function and derivatives tabulation at quadrature points for each field
3276: Level: intermediate
3278: .seealso: `PetscDS`, `PetscTabulation`, `PetscDSCreate()`
3279: @*/
3280: PetscErrorCode PetscDSGetTabulation(PetscDS prob, PetscTabulation *T[])
3281: {
3282: PetscFunctionBegin;
3285: PetscCall(PetscDSSetUp(prob));
3286: *T = prob->T;
3287: PetscFunctionReturn(PETSC_SUCCESS);
3288: }
3290: /*@C
3291: PetscDSGetFaceTabulation - Return the basis tabulation at quadrature points on the faces
3293: Not Collective
3295: Input Parameter:
3296: . prob - The `PetscDS` object
3298: Output Parameter:
3299: . Tf - The basis function and derivative tabulation on each local face at quadrature points for each and field
3301: Level: intermediate
3303: .seealso: `PetscTabulation`, `PetscDS`, `PetscDSGetTabulation()`, `PetscDSCreate()`
3304: @*/
3305: PetscErrorCode PetscDSGetFaceTabulation(PetscDS prob, PetscTabulation *Tf[])
3306: {
3307: PetscFunctionBegin;
3310: PetscCall(PetscDSSetUp(prob));
3311: *Tf = prob->Tf;
3312: PetscFunctionReturn(PETSC_SUCCESS);
3313: }
3315: PetscErrorCode PetscDSGetEvaluationArrays(PetscDS prob, PetscScalar **u, PetscScalar **u_t, PetscScalar **u_x)
3316: {
3317: PetscFunctionBegin;
3319: PetscCall(PetscDSSetUp(prob));
3320: if (u) {
3322: *u = prob->u;
3323: }
3324: if (u_t) {
3326: *u_t = prob->u_t;
3327: }
3328: if (u_x) {
3330: *u_x = prob->u_x;
3331: }
3332: PetscFunctionReturn(PETSC_SUCCESS);
3333: }
3335: PetscErrorCode PetscDSGetWeakFormArrays(PetscDS prob, PetscScalar **f0, PetscScalar **f1, PetscScalar **g0, PetscScalar **g1, PetscScalar **g2, PetscScalar **g3)
3336: {
3337: PetscFunctionBegin;
3339: PetscCall(PetscDSSetUp(prob));
3340: if (f0) {
3342: *f0 = prob->f0;
3343: }
3344: if (f1) {
3346: *f1 = prob->f1;
3347: }
3348: if (g0) {
3350: *g0 = prob->g0;
3351: }
3352: if (g1) {
3354: *g1 = prob->g1;
3355: }
3356: if (g2) {
3358: *g2 = prob->g2;
3359: }
3360: if (g3) {
3362: *g3 = prob->g3;
3363: }
3364: PetscFunctionReturn(PETSC_SUCCESS);
3365: }
3367: PetscErrorCode PetscDSGetWorkspace(PetscDS prob, PetscReal **x, PetscScalar **basisReal, PetscScalar **basisDerReal, PetscScalar **testReal, PetscScalar **testDerReal)
3368: {
3369: PetscFunctionBegin;
3371: PetscCall(PetscDSSetUp(prob));
3372: if (x) {
3374: *x = prob->x;
3375: }
3376: if (basisReal) {
3378: *basisReal = prob->basisReal;
3379: }
3380: if (basisDerReal) {
3382: *basisDerReal = prob->basisDerReal;
3383: }
3384: if (testReal) {
3386: *testReal = prob->testReal;
3387: }
3388: if (testDerReal) {
3390: *testDerReal = prob->testDerReal;
3391: }
3392: PetscFunctionReturn(PETSC_SUCCESS);
3393: }
3395: /*@C
3396: PetscDSAddBoundary - Add a boundary condition to the model. The pointwise functions are used to provide boundary values for essential boundary conditions.
3397: 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
3398: integrals should be performed, using the kernels from `PetscDSSetBdResidual()`.
3400: Collective
3402: Input Parameters:
3403: + ds - The PetscDS object
3404: . type - The type of condition, e.g. `DM_BC_ESSENTIAL`/`DM_BC_ESSENTIAL_FIELD` (Dirichlet), or `DM_BC_NATURAL` (Neumann)
3405: . name - The BC name
3406: . label - The label defining constrained points
3407: . Nv - The number of `DMLabel` values for constrained points
3408: . values - An array of label values for constrained points
3409: . field - The field to constrain
3410: . Nc - The number of constrained field components (0 will constrain all fields)
3411: . comps - An array of constrained component numbers
3412: . bcFunc - A pointwise function giving boundary values
3413: . bcFunc_t - A pointwise function giving the time derivative of the boundary values, or NULL
3414: - ctx - An optional user context for bcFunc
3416: Output Parameter:
3417: - bd - The boundary number
3419: Options Database Keys:
3420: + -bc_<boundary name> <num> - Overrides the boundary ids
3421: - -bc_<boundary name>_comp <num> - Overrides the boundary components
3423: Level: developer
3425: Note:
3426: Both `bcFunc` and `bcFunc_t` will depend on the boundary condition type. If the type if `DM_BC_ESSENTIAL`, Then the calling sequence is:
3428: $ void bcFunc(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar bcval[])
3430: If the type is `DM_BC_ESSENTIAL_FIELD` or other _FIELD value, then the calling sequence is:
3431: .vb
3432: void bcFunc(PetscInt dim, PetscInt Nf, PetscInt NfAux,
3433: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
3434: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
3435: PetscReal time, const PetscReal x[], PetscScalar bcval[])
3436: .ve
3437: + dim - the spatial dimension
3438: . Nf - the number of fields
3439: . uOff - the offset into u[] and u_t[] for each field
3440: . uOff_x - the offset into u_x[] for each field
3441: . u - each field evaluated at the current point
3442: . u_t - the time derivative of each field evaluated at the current point
3443: . u_x - the gradient of each field evaluated at the current point
3444: . aOff - the offset into a[] and a_t[] for each auxiliary field
3445: . aOff_x - the offset into a_x[] for each auxiliary field
3446: . a - each auxiliary field evaluated at the current point
3447: . a_t - the time derivative of each auxiliary field evaluated at the current point
3448: . a_x - the gradient of auxiliary each field evaluated at the current point
3449: . t - current time
3450: . x - coordinates of the current point
3451: . numConstants - number of constant parameters
3452: . constants - constant parameters
3453: - bcval - output values at the current point
3455: .seealso: `PetscDS`, `PetscWeakForm`, `DMLabel`, `DMBoundaryConditionType`, `PetscDSAddBoundaryByName()`, `PetscDSGetBoundary()`, `PetscDSSetResidual()`, `PetscDSSetBdResidual()`
3456: @*/
3457: 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)
3458: {
3459: DSBoundary head = ds->boundary, b;
3460: PetscInt n = 0;
3461: const char *lname;
3463: PetscFunctionBegin;
3471: PetscCheck(field >= 0 && field < ds->Nf, PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_OUTOFRANGE, "Field %" PetscInt_FMT " is not in [0, %" PetscInt_FMT ")", field, ds->Nf);
3472: if (Nc > 0) {
3473: PetscInt *fcomps;
3474: PetscInt c;
3476: PetscCall(PetscDSGetComponents(ds, &fcomps));
3477: PetscCheck(Nc <= fcomps[field], PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_OUTOFRANGE, "Number of constrained components %" PetscInt_FMT " > %" PetscInt_FMT " components for field %" PetscInt_FMT, Nc, fcomps[field], field);
3478: for (c = 0; c < Nc; ++c) {
3479: PetscCheck(comps[c] >= 0 && comps[c] < fcomps[field], PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_OUTOFRANGE, "Constrained component[%" PetscInt_FMT "] %" PetscInt_FMT " not in [0, %" PetscInt_FMT ") components for field %" PetscInt_FMT, c, comps[c], fcomps[field], field);
3480: }
3481: }
3482: PetscCall(PetscNew(&b));
3483: PetscCall(PetscStrallocpy(name, (char **)&b->name));
3484: PetscCall(PetscWeakFormCreate(PETSC_COMM_SELF, &b->wf));
3485: PetscCall(PetscWeakFormSetNumFields(b->wf, ds->Nf));
3486: PetscCall(PetscMalloc1(Nv, &b->values));
3487: if (Nv) PetscCall(PetscArraycpy(b->values, values, Nv));
3488: PetscCall(PetscMalloc1(Nc, &b->comps));
3489: if (Nc) PetscCall(PetscArraycpy(b->comps, comps, Nc));
3490: PetscCall(PetscObjectGetName((PetscObject)label, &lname));
3491: PetscCall(PetscStrallocpy(lname, (char **)&b->lname));
3492: b->type = type;
3493: b->label = label;
3494: b->Nv = Nv;
3495: b->field = field;
3496: b->Nc = Nc;
3497: b->func = bcFunc;
3498: b->func_t = bcFunc_t;
3499: b->ctx = ctx;
3500: b->next = NULL;
3501: /* Append to linked list so that we can preserve the order */
3502: if (!head) ds->boundary = b;
3503: while (head) {
3504: if (!head->next) {
3505: head->next = b;
3506: head = b;
3507: }
3508: head = head->next;
3509: ++n;
3510: }
3511: if (bd) {
3513: *bd = n;
3514: }
3515: PetscFunctionReturn(PETSC_SUCCESS);
3516: }
3518: /*@C
3519: PetscDSAddBoundaryByName - Add a boundary condition to the model. The pointwise functions are used to provide boundary values for essential boundary conditions.
3520: In FEM, they are acting upon by dual basis functionals to generate FEM coefficients which are fixed. Natural boundary conditions signal to PETSc that
3521: boundary integrals should be performed, using the kernels from `PetscDSSetBdResidual()`.
3523: Collective
3525: Input Parameters:
3526: + ds - The `PetscDS` object
3527: . type - The type of condition, e.g. `DM_BC_ESSENTIAL`/`DM_BC_ESSENTIAL_FIELD` (Dirichlet), or `DM_BC_NATURAL` (Neumann)
3528: . name - The BC name
3529: . lname - The naem of the label defining constrained points
3530: . Nv - The number of `DMLabel` values for constrained points
3531: . values - An array of label values for constrained points
3532: . field - The field to constrain
3533: . Nc - The number of constrained field components (0 will constrain all fields)
3534: . comps - An array of constrained component numbers
3535: . bcFunc - A pointwise function giving boundary values
3536: . bcFunc_t - A pointwise function giving the time derivative of the boundary values, or NULL
3537: - ctx - An optional user context for bcFunc
3539: Output Parameter:
3540: - bd - The boundary number
3542: Options Database Keys:
3543: + -bc_<boundary name> <num> - Overrides the boundary ids
3544: - -bc_<boundary name>_comp <num> - Overrides the boundary components
3546: Calling Sequence of `bcFunc` and `bcFunc_t`:
3547: If the type is `DM_BC_ESSENTIAL`
3548: .vb
3549: void bcFunc(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar bcval[])
3550: .ve
3551: If the type is `DM_BC_ESSENTIAL_FIELD` or other _FIELD value,
3552: .vb
3553: void bcFunc(PetscInt dim, PetscInt Nf, PetscInt NfAux,
3554: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
3555: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
3556: PetscReal time, const PetscReal x[], PetscScalar bcval[])
3557: .ve
3558: + dim - the spatial dimension
3559: . Nf - the number of fields
3560: . uOff - the offset into u[] and u_t[] for each field
3561: . uOff_x - the offset into u_x[] for each field
3562: . u - each field evaluated at the current point
3563: . u_t - the time derivative of each field evaluated at the current point
3564: . u_x - the gradient of each field evaluated at the current point
3565: . aOff - the offset into a[] and a_t[] for each auxiliary field
3566: . aOff_x - the offset into a_x[] for each auxiliary field
3567: . a - each auxiliary field evaluated at the current point
3568: . a_t - the time derivative of each auxiliary field evaluated at the current point
3569: . a_x - the gradient of auxiliary each field evaluated at the current point
3570: . t - current time
3571: . x - coordinates of the current point
3572: . numConstants - number of constant parameters
3573: . constants - constant parameters
3574: - bcval - output values at the current point
3576: Level: developer
3578: Note:
3579: This function should only be used with `DMFOREST` currently, since labels cannot be defined before the underlying `DMPLEX` is built.
3581: .seealso: `PetscDS`, `PetscWeakForm`, `DMLabel`, `DMBoundaryConditionType`, `PetscDSAddBoundary()`, `PetscDSGetBoundary()`, `PetscDSSetResidual()`, `PetscDSSetBdResidual()`
3582: @*/
3583: 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)
3584: {
3585: DSBoundary head = ds->boundary, b;
3586: PetscInt n = 0;
3588: PetscFunctionBegin;
3596: PetscCall(PetscNew(&b));
3597: PetscCall(PetscStrallocpy(name, (char **)&b->name));
3598: PetscCall(PetscWeakFormCreate(PETSC_COMM_SELF, &b->wf));
3599: PetscCall(PetscWeakFormSetNumFields(b->wf, ds->Nf));
3600: PetscCall(PetscMalloc1(Nv, &b->values));
3601: if (Nv) PetscCall(PetscArraycpy(b->values, values, Nv));
3602: PetscCall(PetscMalloc1(Nc, &b->comps));
3603: if (Nc) PetscCall(PetscArraycpy(b->comps, comps, Nc));
3604: PetscCall(PetscStrallocpy(lname, (char **)&b->lname));
3605: b->type = type;
3606: b->label = NULL;
3607: b->Nv = Nv;
3608: b->field = field;
3609: b->Nc = Nc;
3610: b->func = bcFunc;
3611: b->func_t = bcFunc_t;
3612: b->ctx = ctx;
3613: b->next = NULL;
3614: /* Append to linked list so that we can preserve the order */
3615: if (!head) ds->boundary = b;
3616: while (head) {
3617: if (!head->next) {
3618: head->next = b;
3619: head = b;
3620: }
3621: head = head->next;
3622: ++n;
3623: }
3624: if (bd) {
3626: *bd = n;
3627: }
3628: PetscFunctionReturn(PETSC_SUCCESS);
3629: }
3631: /*@C
3632: PetscDSUpdateBoundary - Change a boundary condition for the model. The pointwise functions are used to provide boundary values for essential boundary conditions.
3633: 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
3634: should be performed, using the kernels from `PetscDSSetBdResidual()`.
3636: Input Parameters:
3637: + ds - The `PetscDS` object
3638: . bd - The boundary condition number
3639: . type - The type of condition, e.g. `DM_BC_ESSENTIAL`/`DM_BC_ESSENTIAL_FIELD` (Dirichlet), or `DM_BC_NATURAL` (Neumann)
3640: . name - The BC name
3641: . label - The label defining constrained points
3642: . Nv - The number of `DMLabel` ids for constrained points
3643: . values - An array of ids for constrained points
3644: . field - The field to constrain
3645: . Nc - The number of constrained field components
3646: . comps - An array of constrained component numbers
3647: . bcFunc - A pointwise function giving boundary values
3648: . bcFunc_t - A pointwise function giving the time derivative of the boundary values, or NULL
3649: - ctx - An optional user context for bcFunc
3651: Level: developer
3653: Note:
3654: The boundary condition number is the order in which it was registered. The user can get the number of boundary conditions from `PetscDSGetNumBoundary()`.
3655: See `PetscDSAddBoundary()` for a description of the calling sequences for the callbacks.
3657: .seealso: `PetscDS`, `PetscWeakForm`, `DMBoundaryConditionType`, `PetscDSAddBoundary()`, `PetscDSGetBoundary()`, `PetscDSGetNumBoundary()`, `DMLabel`
3658: @*/
3659: 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)
3660: {
3661: DSBoundary b = ds->boundary;
3662: PetscInt n = 0;
3664: PetscFunctionBegin;
3666: while (b) {
3667: if (n == bd) break;
3668: b = b->next;
3669: ++n;
3670: }
3671: PetscCheck(b, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %" PetscInt_FMT " is not in [0, %" PetscInt_FMT ")", bd, n);
3672: if (name) {
3673: PetscCall(PetscFree(b->name));
3674: PetscCall(PetscStrallocpy(name, (char **)&b->name));
3675: }
3676: b->type = type;
3677: if (label) {
3678: const char *name;
3680: b->label = label;
3681: PetscCall(PetscFree(b->lname));
3682: PetscCall(PetscObjectGetName((PetscObject)label, &name));
3683: PetscCall(PetscStrallocpy(name, (char **)&b->lname));
3684: }
3685: if (Nv >= 0) {
3686: b->Nv = Nv;
3687: PetscCall(PetscFree(b->values));
3688: PetscCall(PetscMalloc1(Nv, &b->values));
3689: if (Nv) PetscCall(PetscArraycpy(b->values, values, Nv));
3690: }
3691: if (field >= 0) b->field = field;
3692: if (Nc >= 0) {
3693: b->Nc = Nc;
3694: PetscCall(PetscFree(b->comps));
3695: PetscCall(PetscMalloc1(Nc, &b->comps));
3696: if (Nc) PetscCall(PetscArraycpy(b->comps, comps, Nc));
3697: }
3698: if (bcFunc) b->func = bcFunc;
3699: if (bcFunc_t) b->func_t = bcFunc_t;
3700: if (ctx) b->ctx = ctx;
3701: PetscFunctionReturn(PETSC_SUCCESS);
3702: }
3704: /*@
3705: PetscDSGetNumBoundary - Get the number of registered BC
3707: Input Parameter:
3708: . ds - The `PetscDS` object
3710: Output Parameter:
3711: . numBd - The number of BC
3713: Level: intermediate
3715: .seealso: `PetscDS`, `PetscDSAddBoundary()`, `PetscDSGetBoundary()`
3716: @*/
3717: PetscErrorCode PetscDSGetNumBoundary(PetscDS ds, PetscInt *numBd)
3718: {
3719: DSBoundary b = ds->boundary;
3721: PetscFunctionBegin;
3724: *numBd = 0;
3725: while (b) {
3726: ++(*numBd);
3727: b = b->next;
3728: }
3729: PetscFunctionReturn(PETSC_SUCCESS);
3730: }
3732: /*@C
3733: PetscDSGetBoundary - Gets a boundary condition to the model
3735: Input Parameters:
3736: + ds - The `PetscDS` object
3737: - bd - The BC number
3739: Output Parameters:
3740: + wf - The `PetscWeakForm` holding the pointwise functions
3741: . type - The type of condition, e.g. `DM_BC_ESSENTIAL`/`DM_BC_ESSENTIAL_FIELD` (Dirichlet), or `DM_BC_NATURAL` (Neumann)
3742: . name - The BC name
3743: . label - The label defining constrained points
3744: . Nv - The number of `DMLabel` ids for constrained points
3745: . values - An array of ids for constrained points
3746: . field - The field to constrain
3747: . Nc - The number of constrained field components
3748: . comps - An array of constrained component numbers
3749: . bcFunc - A pointwise function giving boundary values
3750: . bcFunc_t - A pointwise function giving the time derivative of the boundary values
3751: - ctx - An optional user context for bcFunc
3753: Options Database Keys:
3754: + -bc_<boundary name> <num> - Overrides the boundary ids
3755: - -bc_<boundary name>_comp <num> - Overrides the boundary components
3757: Level: developer
3759: .seealso: `PetscDS`, `PetscWeakForm`, `DMBoundaryConditionType`, `PetscDSAddBoundary()`, `DMLabel`
3760: @*/
3761: 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)
3762: {
3763: DSBoundary b = ds->boundary;
3764: PetscInt n = 0;
3766: PetscFunctionBegin;
3768: while (b) {
3769: if (n == bd) break;
3770: b = b->next;
3771: ++n;
3772: }
3773: PetscCheck(b, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %" PetscInt_FMT " is not in [0, %" PetscInt_FMT ")", bd, n);
3774: if (wf) {
3776: *wf = b->wf;
3777: }
3778: if (type) {
3780: *type = b->type;
3781: }
3782: if (name) {
3784: *name = b->name;
3785: }
3786: if (label) {
3788: *label = b->label;
3789: }
3790: if (Nv) {
3792: *Nv = b->Nv;
3793: }
3794: if (values) {
3796: *values = b->values;
3797: }
3798: if (field) {
3800: *field = b->field;
3801: }
3802: if (Nc) {
3804: *Nc = b->Nc;
3805: }
3806: if (comps) {
3808: *comps = b->comps;
3809: }
3810: if (func) {
3812: *func = b->func;
3813: }
3814: if (func_t) {
3816: *func_t = b->func_t;
3817: }
3818: if (ctx) {
3820: *ctx = b->ctx;
3821: }
3822: PetscFunctionReturn(PETSC_SUCCESS);
3823: }
3825: static PetscErrorCode DSBoundaryDuplicate_Internal(DSBoundary b, DSBoundary *bNew)
3826: {
3827: PetscFunctionBegin;
3828: PetscCall(PetscNew(bNew));
3829: PetscCall(PetscWeakFormCreate(PETSC_COMM_SELF, &(*bNew)->wf));
3830: PetscCall(PetscWeakFormCopy(b->wf, (*bNew)->wf));
3831: PetscCall(PetscStrallocpy(b->name, (char **)&((*bNew)->name)));
3832: PetscCall(PetscStrallocpy(b->lname, (char **)&((*bNew)->lname)));
3833: (*bNew)->type = b->type;
3834: (*bNew)->label = b->label;
3835: (*bNew)->Nv = b->Nv;
3836: PetscCall(PetscMalloc1(b->Nv, &(*bNew)->values));
3837: PetscCall(PetscArraycpy((*bNew)->values, b->values, b->Nv));
3838: (*bNew)->field = b->field;
3839: (*bNew)->Nc = b->Nc;
3840: PetscCall(PetscMalloc1(b->Nc, &(*bNew)->comps));
3841: PetscCall(PetscArraycpy((*bNew)->comps, b->comps, b->Nc));
3842: (*bNew)->func = b->func;
3843: (*bNew)->func_t = b->func_t;
3844: (*bNew)->ctx = b->ctx;
3845: PetscFunctionReturn(PETSC_SUCCESS);
3846: }
3848: /*@
3849: PetscDSCopyBoundary - Copy all boundary condition objects to the new problem
3851: Not Collective
3853: Input Parameters:
3854: + ds - The source `PetscDS` object
3855: . numFields - The number of selected fields, or `PETSC_DEFAULT` for all fields
3856: - fields - The selected fields, or NULL for all fields
3858: Output Parameter:
3859: . newds - The target `PetscDS`, now with a copy of the boundary conditions
3861: Level: intermediate
3863: .seealso: `PetscDS`, `DMBoundary`, `PetscDSCopyEquations()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
3864: @*/
3865: PetscErrorCode PetscDSCopyBoundary(PetscDS ds, PetscInt numFields, const PetscInt fields[], PetscDS newds)
3866: {
3867: DSBoundary b, *lastnext;
3869: PetscFunctionBegin;
3872: if (ds == newds) PetscFunctionReturn(PETSC_SUCCESS);
3873: PetscCall(PetscDSDestroyBoundary(newds));
3874: lastnext = &(newds->boundary);
3875: for (b = ds->boundary; b; b = b->next) {
3876: DSBoundary bNew;
3877: PetscInt fieldNew = -1;
3879: if (numFields > 0 && fields) {
3880: PetscInt f;
3882: for (f = 0; f < numFields; ++f)
3883: if (b->field == fields[f]) break;
3884: if (f == numFields) continue;
3885: fieldNew = f;
3886: }
3887: PetscCall(DSBoundaryDuplicate_Internal(b, &bNew));
3888: bNew->field = fieldNew < 0 ? b->field : fieldNew;
3889: *lastnext = bNew;
3890: lastnext = &(bNew->next);
3891: }
3892: PetscFunctionReturn(PETSC_SUCCESS);
3893: }
3895: /*@
3896: PetscDSDestroyBoundary - Remove all `DMBoundary` objects from the `PetscDS`
3898: Not Collective
3900: Input Parameter:
3901: . ds - The `PetscDS` object
3903: Level: intermediate
3905: .seealso: `PetscDS`, `DMBoundary`, `PetscDSCopyBoundary()`, `PetscDSCopyEquations()`
3906: @*/
3907: PetscErrorCode PetscDSDestroyBoundary(PetscDS ds)
3908: {
3909: DSBoundary next = ds->boundary;
3911: PetscFunctionBegin;
3912: while (next) {
3913: DSBoundary b = next;
3915: next = b->next;
3916: PetscCall(PetscWeakFormDestroy(&b->wf));
3917: PetscCall(PetscFree(b->name));
3918: PetscCall(PetscFree(b->lname));
3919: PetscCall(PetscFree(b->values));
3920: PetscCall(PetscFree(b->comps));
3921: PetscCall(PetscFree(b));
3922: }
3923: PetscFunctionReturn(PETSC_SUCCESS);
3924: }
3926: /*@
3927: PetscDSSelectDiscretizations - Copy discretizations to the new problem with different field layout
3929: Not Collective
3931: Input Parameters:
3932: + prob - The `PetscDS` object
3933: . numFields - Number of new fields
3934: - fields - Old field number for each new field
3936: Output Parameter:
3937: . newprob - The `PetscDS` copy
3939: Level: intermediate
3941: .seealso: `PetscDS`, `PetscDSSelectEquations()`, `PetscDSCopyBoundary()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
3942: @*/
3943: PetscErrorCode PetscDSSelectDiscretizations(PetscDS prob, PetscInt numFields, const PetscInt fields[], PetscDS newprob)
3944: {
3945: PetscInt Nf, Nfn, fn;
3947: PetscFunctionBegin;
3951: PetscCall(PetscDSGetNumFields(prob, &Nf));
3952: PetscCall(PetscDSGetNumFields(newprob, &Nfn));
3953: numFields = numFields < 0 ? Nf : numFields;
3954: for (fn = 0; fn < numFields; ++fn) {
3955: const PetscInt f = fields ? fields[fn] : fn;
3956: PetscObject disc;
3958: if (f >= Nf) continue;
3959: PetscCall(PetscDSGetDiscretization(prob, f, &disc));
3960: PetscCall(PetscDSSetDiscretization(newprob, fn, disc));
3961: }
3962: PetscFunctionReturn(PETSC_SUCCESS);
3963: }
3965: /*@
3966: PetscDSSelectEquations - Copy pointwise function pointers to the new problem with different field layout
3968: Not Collective
3970: Input Parameters:
3971: + prob - The `PetscDS` object
3972: . numFields - Number of new fields
3973: - fields - Old field number for each new field
3975: Output Parameter:
3976: . newprob - The `PetscDS` copy
3978: Level: intermediate
3980: .seealso: `PetscDS`, `PetscDSSelectDiscretizations()`, `PetscDSCopyBoundary()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
3981: @*/
3982: PetscErrorCode PetscDSSelectEquations(PetscDS prob, PetscInt numFields, const PetscInt fields[], PetscDS newprob)
3983: {
3984: PetscInt Nf, Nfn, fn, gn;
3986: PetscFunctionBegin;
3990: PetscCall(PetscDSGetNumFields(prob, &Nf));
3991: PetscCall(PetscDSGetNumFields(newprob, &Nfn));
3992: PetscCheck(numFields <= Nfn, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_SIZ, "Number of fields %" PetscInt_FMT " to transfer must not be greater then the total number of fields %" PetscInt_FMT, numFields, Nfn);
3993: for (fn = 0; fn < numFields; ++fn) {
3994: const PetscInt f = fields ? fields[fn] : fn;
3995: PetscPointFunc obj;
3996: PetscPointFunc f0, f1;
3997: PetscBdPointFunc f0Bd, f1Bd;
3998: PetscRiemannFunc r;
4000: if (f >= Nf) continue;
4001: PetscCall(PetscDSGetObjective(prob, f, &obj));
4002: PetscCall(PetscDSGetResidual(prob, f, &f0, &f1));
4003: PetscCall(PetscDSGetBdResidual(prob, f, &f0Bd, &f1Bd));
4004: PetscCall(PetscDSGetRiemannSolver(prob, f, &r));
4005: PetscCall(PetscDSSetObjective(newprob, fn, obj));
4006: PetscCall(PetscDSSetResidual(newprob, fn, f0, f1));
4007: PetscCall(PetscDSSetBdResidual(newprob, fn, f0Bd, f1Bd));
4008: PetscCall(PetscDSSetRiemannSolver(newprob, fn, r));
4009: for (gn = 0; gn < numFields; ++gn) {
4010: const PetscInt g = fields ? fields[gn] : gn;
4011: PetscPointJac g0, g1, g2, g3;
4012: PetscPointJac g0p, g1p, g2p, g3p;
4013: PetscBdPointJac g0Bd, g1Bd, g2Bd, g3Bd;
4015: if (g >= Nf) continue;
4016: PetscCall(PetscDSGetJacobian(prob, f, g, &g0, &g1, &g2, &g3));
4017: PetscCall(PetscDSGetJacobianPreconditioner(prob, f, g, &g0p, &g1p, &g2p, &g3p));
4018: PetscCall(PetscDSGetBdJacobian(prob, f, g, &g0Bd, &g1Bd, &g2Bd, &g3Bd));
4019: PetscCall(PetscDSSetJacobian(newprob, fn, gn, g0, g1, g2, g3));
4020: PetscCall(PetscDSSetJacobianPreconditioner(newprob, fn, gn, g0p, g1p, g2p, g3p));
4021: PetscCall(PetscDSSetBdJacobian(newprob, fn, gn, g0Bd, g1Bd, g2Bd, g3Bd));
4022: }
4023: }
4024: PetscFunctionReturn(PETSC_SUCCESS);
4025: }
4027: /*@
4028: PetscDSCopyEquations - Copy all pointwise function pointers to another `PetscDS`
4030: Not Collective
4032: Input Parameter:
4033: . prob - The `PetscDS` object
4035: Output Parameter:
4036: . newprob - The `PetscDS` copy
4038: Level: intermediate
4040: .seealso: `PetscDS`, `PetscDSCopyBoundary()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
4041: @*/
4042: PetscErrorCode PetscDSCopyEquations(PetscDS prob, PetscDS newprob)
4043: {
4044: PetscWeakForm wf, newwf;
4045: PetscInt Nf, Ng;
4047: PetscFunctionBegin;
4050: PetscCall(PetscDSGetNumFields(prob, &Nf));
4051: PetscCall(PetscDSGetNumFields(newprob, &Ng));
4052: PetscCheck(Nf == Ng, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_SIZ, "Number of fields must match %" PetscInt_FMT " != %" PetscInt_FMT, Nf, Ng);
4053: PetscCall(PetscDSGetWeakForm(prob, &wf));
4054: PetscCall(PetscDSGetWeakForm(newprob, &newwf));
4055: PetscCall(PetscWeakFormCopy(wf, newwf));
4056: PetscFunctionReturn(PETSC_SUCCESS);
4057: }
4059: /*@
4060: PetscDSCopyConstants - Copy all constants to another `PetscDS`
4062: Not Collective
4064: Input Parameter:
4065: . prob - The `PetscDS` object
4067: Output Parameter:
4068: . newprob - The `PetscDS` copy
4070: Level: intermediate
4072: .seealso: `PetscDS`, `PetscDSCopyBoundary()`, `PetscDSCopyEquations()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
4073: @*/
4074: PetscErrorCode PetscDSCopyConstants(PetscDS prob, PetscDS newprob)
4075: {
4076: PetscInt Nc;
4077: const PetscScalar *constants;
4079: PetscFunctionBegin;
4082: PetscCall(PetscDSGetConstants(prob, &Nc, &constants));
4083: PetscCall(PetscDSSetConstants(newprob, Nc, (PetscScalar *)constants));
4084: PetscFunctionReturn(PETSC_SUCCESS);
4085: }
4087: /*@
4088: PetscDSCopyExactSolutions - Copy all exact solutions to another `PetscDS`
4090: Not Collective
4092: Input Parameter:
4093: . ds - The `PetscDS` object
4095: Output Parameter:
4096: . newds - The `PetscDS` copy
4098: Level: intermediate
4100: .seealso: `PetscDS`, `PetscDSCopyBoundary()`, `PetscDSCopyEquations()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
4101: @*/
4102: PetscErrorCode PetscDSCopyExactSolutions(PetscDS ds, PetscDS newds)
4103: {
4104: PetscSimplePointFunc sol;
4105: void *ctx;
4106: PetscInt Nf, f;
4108: PetscFunctionBegin;
4111: PetscCall(PetscDSGetNumFields(ds, &Nf));
4112: for (f = 0; f < Nf; ++f) {
4113: PetscCall(PetscDSGetExactSolution(ds, f, &sol, &ctx));
4114: PetscCall(PetscDSSetExactSolution(newds, f, sol, ctx));
4115: PetscCall(PetscDSGetExactSolutionTimeDerivative(ds, f, &sol, &ctx));
4116: PetscCall(PetscDSSetExactSolutionTimeDerivative(newds, f, sol, ctx));
4117: }
4118: PetscFunctionReturn(PETSC_SUCCESS);
4119: }
4121: PetscErrorCode PetscDSCopy(PetscDS ds, DM dmNew, PetscDS dsNew)
4122: {
4123: DSBoundary b;
4124: PetscInt cdim, Nf, f, d;
4125: PetscBool isCohesive;
4126: void *ctx;
4128: PetscFunctionBegin;
4129: PetscCall(PetscDSCopyConstants(ds, dsNew));
4130: PetscCall(PetscDSCopyExactSolutions(ds, dsNew));
4131: PetscCall(PetscDSSelectDiscretizations(ds, PETSC_DETERMINE, NULL, dsNew));
4132: PetscCall(PetscDSCopyEquations(ds, dsNew));
4133: PetscCall(PetscDSGetNumFields(ds, &Nf));
4134: for (f = 0; f < Nf; ++f) {
4135: PetscCall(PetscDSGetContext(ds, f, &ctx));
4136: PetscCall(PetscDSSetContext(dsNew, f, ctx));
4137: PetscCall(PetscDSGetCohesive(ds, f, &isCohesive));
4138: PetscCall(PetscDSSetCohesive(dsNew, f, isCohesive));
4139: PetscCall(PetscDSGetJetDegree(ds, f, &d));
4140: PetscCall(PetscDSSetJetDegree(dsNew, f, d));
4141: }
4142: if (Nf) {
4143: PetscCall(PetscDSGetCoordinateDimension(ds, &cdim));
4144: PetscCall(PetscDSSetCoordinateDimension(dsNew, cdim));
4145: }
4146: PetscCall(PetscDSCopyBoundary(ds, PETSC_DETERMINE, NULL, dsNew));
4147: for (b = dsNew->boundary; b; b = b->next) {
4148: PetscCall(DMGetLabel(dmNew, b->lname, &b->label));
4149: /* Do not check if label exists here, since p4est calls this for the reference tree which does not have the labels */
4150: //PetscCheck(b->label,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Label %s missing in new DM", name);
4151: }
4152: PetscFunctionReturn(PETSC_SUCCESS);
4153: }
4155: PetscErrorCode PetscDSGetHeightSubspace(PetscDS prob, PetscInt height, PetscDS *subprob)
4156: {
4157: PetscInt dim, Nf, f;
4159: PetscFunctionBegin;
4162: if (height == 0) {
4163: *subprob = prob;
4164: PetscFunctionReturn(PETSC_SUCCESS);
4165: }
4166: PetscCall(PetscDSGetNumFields(prob, &Nf));
4167: PetscCall(PetscDSGetSpatialDimension(prob, &dim));
4168: PetscCheck(height <= dim, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_OUTOFRANGE, "DS can only handle height in [0, %" PetscInt_FMT "], not %" PetscInt_FMT, dim, height);
4169: if (!prob->subprobs) PetscCall(PetscCalloc1(dim, &prob->subprobs));
4170: if (!prob->subprobs[height - 1]) {
4171: PetscInt cdim;
4173: PetscCall(PetscDSCreate(PetscObjectComm((PetscObject)prob), &prob->subprobs[height - 1]));
4174: PetscCall(PetscDSGetCoordinateDimension(prob, &cdim));
4175: PetscCall(PetscDSSetCoordinateDimension(prob->subprobs[height - 1], cdim));
4176: for (f = 0; f < Nf; ++f) {
4177: PetscFE subfe;
4178: PetscObject obj;
4179: PetscClassId id;
4181: PetscCall(PetscDSGetDiscretization(prob, f, &obj));
4182: PetscCall(PetscObjectGetClassId(obj, &id));
4183: if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetHeightSubspace((PetscFE)obj, height, &subfe));
4184: else SETERRQ(PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONG, "Unsupported discretization type for field %" PetscInt_FMT, f);
4185: PetscCall(PetscDSSetDiscretization(prob->subprobs[height - 1], f, (PetscObject)subfe));
4186: }
4187: }
4188: *subprob = prob->subprobs[height - 1];
4189: PetscFunctionReturn(PETSC_SUCCESS);
4190: }
4192: PetscErrorCode PetscDSPermuteQuadPoint(PetscDS ds, PetscInt ornt, PetscInt field, PetscInt q, PetscInt *qperm)
4193: {
4194: IS permIS;
4195: PetscQuadrature quad;
4196: DMPolytopeType ct;
4197: const PetscInt *perm;
4198: PetscInt Na, Nq;
4200: PetscFunctionBeginHot;
4201: PetscCall(PetscFEGetQuadrature((PetscFE)ds->disc[field], &quad));
4202: PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL));
4203: PetscCall(PetscQuadratureGetCellType(quad, &ct));
4204: PetscCheck(q >= 0 && q < Nq, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Quadrature point %" PetscInt_FMT " is not in [0, %" PetscInt_FMT ")", q, Nq);
4205: Na = DMPolytopeTypeGetNumArrangments(ct) / 2;
4206: PetscCheck(ornt >= -Na && ornt < Na, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Orientation %" PetscInt_FMT " of %s is not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", ornt, DMPolytopeTypes[ct], -Na, Na);
4207: if (!ds->quadPerm[(PetscInt)ct]) PetscCall(PetscQuadratureComputePermutations(quad, NULL, &ds->quadPerm[(PetscInt)ct]));
4208: permIS = ds->quadPerm[(PetscInt)ct][ornt + Na];
4209: PetscCall(ISGetIndices(permIS, &perm));
4210: *qperm = perm[q];
4211: PetscCall(ISRestoreIndices(permIS, &perm));
4212: PetscFunctionReturn(PETSC_SUCCESS);
4213: }
4215: PetscErrorCode PetscDSGetDiscType_Internal(PetscDS ds, PetscInt f, PetscDiscType *disctype)
4216: {
4217: PetscObject obj;
4218: PetscClassId id;
4219: PetscInt Nf;
4221: PetscFunctionBegin;
4224: *disctype = PETSC_DISC_NONE;
4225: PetscCall(PetscDSGetNumFields(ds, &Nf));
4226: PetscCheck(f < Nf, PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_SIZ, "Field %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, Nf);
4227: PetscCall(PetscDSGetDiscretization(ds, f, &obj));
4228: if (obj) {
4229: PetscCall(PetscObjectGetClassId(obj, &id));
4230: if (id == PETSCFE_CLASSID) *disctype = PETSC_DISC_FE;
4231: else *disctype = PETSC_DISC_FV;
4232: }
4233: PetscFunctionReturn(PETSC_SUCCESS);
4234: }
4236: static PetscErrorCode PetscDSDestroy_Basic(PetscDS ds)
4237: {
4238: PetscFunctionBegin;
4239: PetscCall(PetscFree(ds->data));
4240: PetscFunctionReturn(PETSC_SUCCESS);
4241: }
4243: static PetscErrorCode PetscDSInitialize_Basic(PetscDS ds)
4244: {
4245: PetscFunctionBegin;
4246: ds->ops->setfromoptions = NULL;
4247: ds->ops->setup = NULL;
4248: ds->ops->view = NULL;
4249: ds->ops->destroy = PetscDSDestroy_Basic;
4250: PetscFunctionReturn(PETSC_SUCCESS);
4251: }
4253: /*MC
4254: PETSCDSBASIC = "basic" - A discrete system with pointwise residual and boundary residual functions
4256: Level: intermediate
4258: .seealso: `PetscDSType`, `PetscDSCreate()`, `PetscDSSetType()`
4259: M*/
4261: PETSC_EXTERN PetscErrorCode PetscDSCreate_Basic(PetscDS ds)
4262: {
4263: PetscDS_Basic *b;
4265: PetscFunctionBegin;
4267: PetscCall(PetscNew(&b));
4268: ds->data = b;
4270: PetscCall(PetscDSInitialize_Basic(ds));
4271: PetscFunctionReturn(PETSC_SUCCESS);
4272: }