Actual source code: dmfieldshell.c
1: #include <petsc/private/dmfieldimpl.h>
3: typedef struct _n_DMField_Shell {
4: PetscCtx ctx;
5: PetscErrorCode (*destroy)(DMField);
6: } DMField_Shell;
8: /*@
9: DMFieldShellGetContext - Retrieve the user-supplied context associated with a `DMFIELDSHELL`.
11: Not Collective
13: Input Parameter:
14: . field - the `DMField` of type `DMFIELDSHELL`
16: Output Parameter:
17: . ctx - the context pointer that was passed to `DMFieldCreateShell()`
19: Level: intermediate
21: .seealso: `DMField`, `DMFIELDSHELL`, `DMFieldCreateShell()`
22: @*/
23: PetscErrorCode DMFieldShellGetContext(DMField field, PetscCtxRt ctx)
24: {
25: PetscBool flg;
27: PetscFunctionBegin;
29: PetscAssertPointer(ctx, 2);
30: PetscCall(PetscObjectTypeCompare((PetscObject)field, DMFIELDSHELL, &flg));
31: PetscCheck(flg, PetscObjectComm((PetscObject)field), PETSC_ERR_SUP, "Cannot get context from non-shell shield");
32: *(void **)ctx = ((DMField_Shell *)field->data)->ctx;
33: PetscFunctionReturn(PETSC_SUCCESS);
34: }
36: static PetscErrorCode DMFieldDestroy_Shell(DMField field)
37: {
38: DMField_Shell *shell = (DMField_Shell *)field->data;
40: PetscFunctionBegin;
41: if (shell->destroy) PetscCall((*shell->destroy)(field));
42: PetscCall(PetscFree(field->data));
43: PetscFunctionReturn(PETSC_SUCCESS);
44: }
46: /*@C
47: DMFieldShellEvaluateFEDefault - Default finite-element evaluation for a `DMFIELDSHELL` that maps the quadrature points to real space using the coordinate `DMField` and then calls `DMFieldEvaluate()`.
49: Not Collective
51: Input Parameters:
52: + field - the `DMField` of type `DMFIELDSHELL`
53: . pointIS - the `IS` of mesh points at which to evaluate
54: . quad - the reference-element quadrature
55: - type - `PETSC_SCALAR` or `PETSC_REAL`
57: Output Parameters:
58: + B - values at quadrature points, or `NULL`
59: . D - derivatives at quadrature points, or `NULL`
60: - H - Hessians at quadrature points, or `NULL`
62: Level: developer
64: Note:
65: Intended to be registered as the FE evaluation callback via `DMFieldShellSetEvaluateFE()` when the shell only supplies a bulk `DMFieldEvaluate()` implementation.
67: .seealso: `DMField`, `DMFIELDSHELL`, `DMFieldShellSetEvaluateFE()`, `DMFieldShellEvaluateFVDefault()`, `DMFieldEvaluate()`
68: @*/
69: PetscErrorCode DMFieldShellEvaluateFEDefault(DMField field, IS pointIS, PetscQuadrature quad, PetscDataType type, void *B, void *D, void *H)
70: {
71: DM dm = field->dm;
72: DMField coordField;
73: PetscFEGeom *geom;
74: Vec pushforward;
75: PetscInt dimC, dim, numPoints, Nq, p, Nc;
76: PetscScalar *pfArray;
78: PetscFunctionBegin;
79: Nc = field->numComponents;
80: PetscCall(DMGetCoordinateField(dm, &coordField));
81: PetscCall(DMFieldCreateFEGeom(coordField, pointIS, quad, PETSC_FEGEOM_BASIC, &geom));
82: PetscCall(DMGetCoordinateDim(dm, &dimC));
83: PetscCall(PetscQuadratureGetData(quad, &dim, NULL, &Nq, NULL, NULL));
84: PetscCall(ISGetLocalSize(pointIS, &numPoints));
85: PetscCall(PetscMalloc1(dimC * Nq * numPoints, &pfArray));
86: for (p = 0; p < numPoints * dimC * Nq; p++) pfArray[p] = geom->v[p];
87: PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)pointIS), dimC, dimC * Nq * numPoints, PETSC_DETERMINE, pfArray, &pushforward));
88: PetscCall(DMFieldEvaluate(field, pushforward, type, B, D, H));
89: /* TODO: handle covariant/contravariant pullbacks */
90: if (D) {
91: if (type == PETSC_SCALAR) {
92: PetscScalar *sD = (PetscScalar *)D;
93: PetscInt q;
95: for (p = 0; p < numPoints * Nq; p++) {
96: for (q = 0; q < Nc; q++) {
97: PetscScalar d[3];
99: PetscInt i, j;
101: for (i = 0; i < dimC; i++) d[i] = sD[(p * Nc + q) * dimC + i];
102: for (i = 0; i < dimC; i++) sD[(p * Nc + q) * dimC + i] = 0.;
103: for (i = 0; i < dimC; i++) {
104: for (j = 0; j < dimC; j++) sD[(p * Nc + q) * dimC + i] += geom->J[(p * dimC + j) * dimC + i] * d[j];
105: }
106: }
107: }
108: } else {
109: PetscReal *rD = (PetscReal *)D;
110: PetscInt q;
112: for (p = 0; p < numPoints * Nq; p++) {
113: for (q = 0; q < Nc; q++) {
114: PetscReal d[3];
116: PetscInt i, j;
118: for (i = 0; i < dimC; i++) d[i] = rD[(p * Nc + q) * dimC + i];
119: for (i = 0; i < dimC; i++) rD[(p * Nc + q) * dimC + i] = 0.;
120: for (i = 0; i < dimC; i++) {
121: for (j = 0; j < dimC; j++) rD[(p * Nc + q) * dimC + i] += geom->J[(p * dimC + j) * dimC + i] * d[j];
122: }
123: }
124: }
125: }
126: }
127: if (H) {
128: if (type == PETSC_SCALAR) {
129: PetscScalar *sH = (PetscScalar *)H;
130: PetscInt q;
132: for (p = 0; p < numPoints * Nq; p++) {
133: for (q = 0; q < Nc; q++) {
134: PetscScalar d[3][3];
136: PetscInt i, j, k, l;
138: for (i = 0; i < dimC; i++)
139: for (j = 0; j < dimC; j++) d[i][j] = sH[((p * Nc + q) * dimC + i) * dimC + j];
140: for (i = 0; i < dimC; i++)
141: for (j = 0; j < dimC; j++) sH[((p * Nc + q) * dimC + i) * dimC + j] = 0.;
142: for (i = 0; i < dimC; i++) {
143: for (j = 0; j < dimC; j++) {
144: for (k = 0; k < dimC; k++) {
145: for (l = 0; l < dimC; l++) sH[((p * Nc + q) * dimC + i) * dimC + j] += geom->J[(p * dimC + k) * dimC + i] * geom->J[(p * dimC + l) * dimC + j] * d[k][l];
146: }
147: }
148: }
149: }
150: }
151: } else {
152: PetscReal *rH = (PetscReal *)H;
153: PetscInt q;
155: for (p = 0; p < numPoints * Nq; p++) {
156: for (q = 0; q < Nc; q++) {
157: PetscReal d[3][3];
159: PetscInt i, j, k, l;
161: for (i = 0; i < dimC; i++)
162: for (j = 0; j < dimC; j++) d[i][j] = rH[((p * Nc + q) * dimC + i) * dimC + j];
163: for (i = 0; i < dimC; i++)
164: for (j = 0; j < dimC; j++) rH[((p * Nc + q) * dimC + i) * dimC + j] = 0.;
165: for (i = 0; i < dimC; i++) {
166: for (j = 0; j < dimC; j++) {
167: for (k = 0; k < dimC; k++) {
168: for (l = 0; l < dimC; l++) rH[((p * Nc + q) * dimC + i) * dimC + j] += geom->J[(p * dimC + k) * dimC + i] * geom->J[(p * dimC + l) * dimC + j] * d[k][l];
169: }
170: }
171: }
172: }
173: }
174: }
175: }
176: PetscCall(VecDestroy(&pushforward));
177: PetscCall(PetscFree(pfArray));
178: PetscCall(PetscFEGeomDestroy(&geom));
179: PetscFunctionReturn(PETSC_SUCCESS);
180: }
182: /*@C
183: DMFieldShellEvaluateFVDefault - Default finite-volume evaluation for a `DMFIELDSHELL` that samples at cell centroids using the coordinate `DMField`'s default quadrature and calls `DMFieldEvaluate()`.
185: Not Collective
187: Input Parameters:
188: + field - the `DMField` of type `DMFIELDSHELL`
189: . pointIS - the `IS` of mesh cells at which to evaluate
190: - type - `PETSC_SCALAR` or `PETSC_REAL`
192: Output Parameters:
193: + B - cell-averaged values, or `NULL`
194: . D - cell-averaged derivatives, or `NULL`
195: - H - cell-averaged Hessians, or `NULL`
197: Level: developer
199: Note:
200: Intended to be registered as the FV evaluation callback via `DMFieldShellSetEvaluateFV()` when the shell only supplies a bulk `DMFieldEvaluate()` implementation.
202: .seealso: `DMField`, `DMFIELDSHELL`, `DMFieldShellSetEvaluateFV()`, `DMFieldShellEvaluateFEDefault()`, `DMFieldEvaluate()`
203: @*/
204: PetscErrorCode DMFieldShellEvaluateFVDefault(DMField field, IS pointIS, PetscDataType type, void *B, void *D, void *H)
205: {
206: DM dm = field->dm;
207: DMField coordField;
208: PetscFEGeom *geom;
209: Vec pushforward;
210: PetscInt dimC, dim, numPoints, Nq, p;
211: PetscScalar *pfArray;
212: PetscQuadrature quad;
213: MPI_Comm comm;
215: PetscFunctionBegin;
216: PetscCall(PetscObjectGetComm((PetscObject)field, &comm));
217: PetscCall(DMGetDimension(dm, &dim));
218: PetscCall(DMGetCoordinateDim(dm, &dimC));
219: PetscCall(DMGetCoordinateField(dm, &coordField));
220: PetscCall(DMFieldGetFVQuadrature_Internal(coordField, pointIS, &quad));
221: PetscCheck(quad, comm, PETSC_ERR_ARG_WRONGSTATE, "coordinate field must have default quadrature for FV computation");
222: PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL));
223: PetscCheck(Nq == 1, comm, PETSC_ERR_ARG_WRONGSTATE, "quadrature must have only one point");
224: PetscCall(DMFieldCreateFEGeom(coordField, pointIS, quad, PETSC_FEGEOM_BASIC, &geom));
225: PetscCall(ISGetLocalSize(pointIS, &numPoints));
226: PetscCall(PetscMalloc1(dimC * numPoints, &pfArray));
227: for (p = 0; p < numPoints * dimC; p++) pfArray[p] = geom->v[p];
228: PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)pointIS), dimC, dimC * numPoints, PETSC_DETERMINE, pfArray, &pushforward));
229: PetscCall(DMFieldEvaluate(field, pushforward, type, B, D, H));
230: PetscCall(PetscQuadratureDestroy(&quad));
231: PetscCall(VecDestroy(&pushforward));
232: PetscCall(PetscFree(pfArray));
233: PetscCall(PetscFEGeomDestroy(&geom));
234: PetscFunctionReturn(PETSC_SUCCESS);
235: }
237: /*@C
238: DMFieldShellSetDestroy - Register a destroy callback that will be invoked when a `DMFIELDSHELL` is destroyed.
240: Logically Collective
242: Input Parameters:
243: + field - the `DMField` of type `DMFIELDSHELL`
244: - destroy - the destroy routine, called before the shell's own data is freed
246: Calling sequence of `destroy`:
247: . field - the `DMField` of type `DMFIELDSHELL` being destroyed
249: Level: intermediate
251: .seealso: `DMField`, `DMFIELDSHELL`, `DMFieldCreateShell()`, `DMFieldDestroy()`
252: @*/
253: PetscErrorCode DMFieldShellSetDestroy(DMField field, PetscErrorCode (*destroy)(DMField field))
254: {
255: DMField_Shell *shell = (DMField_Shell *)field->data;
257: PetscFunctionBegin;
259: shell->destroy = destroy;
260: PetscFunctionReturn(PETSC_SUCCESS);
261: }
263: /*@C
264: DMFieldShellSetEvaluate - Register the routine that evaluates a `DMFIELDSHELL` at an arbitrary set of real-space points supplied as a `Vec` of coordinates.
266: Logically Collective
268: Input Parameters:
269: + field - the `DMField` of type `DMFIELDSHELL`
270: - evaluate - the evaluation callback
272: Calling sequence of `evaluate`:
273: + field - the `DMField` of type `DMFIELDSHELL`
274: . u - the points at which to evaluate the field, as a `Vec` of coordinates of size d x n
275: . dtype - `PETSC_SCALAR` or `PETSC_REAL`
276: . B - array of field values at each point, or `NULL`
277: . D - array of field spatial derivatives at each point, or `NULL`
278: - H - array of field spatial Hessians at each point, or `NULL`
280: Level: intermediate
282: .seealso: `DMField`, `DMFIELDSHELL`, `DMFieldCreateShell()`, `DMFieldEvaluate()`, `DMFieldShellSetEvaluateFE()`, `DMFieldShellSetEvaluateFV()`
283: @*/
284: PetscErrorCode DMFieldShellSetEvaluate(DMField field, PetscErrorCode (*evaluate)(DMField field, Vec u, PetscDataType dtype, void *B, void *D, void *H))
285: {
286: PetscFunctionBegin;
288: field->ops->evaluate = evaluate;
289: PetscFunctionReturn(PETSC_SUCCESS);
290: }
292: /*@C
293: DMFieldShellSetEvaluateFE - Register the routine that evaluates a `DMFIELDSHELL` at finite-element quadrature points over a set of mesh points.
295: Logically Collective
297: Input Parameters:
298: + field - the `DMField` of type `DMFIELDSHELL`
299: - evaluateFE - the FE evaluation callback
301: Calling sequence of `evaluateFE`:
302: + field - the `DMField` of type `DMFIELDSHELL`
303: . is - the `IS` of mesh cells on which to evaluate the field
304: . quad - the reference-cell `PetscQuadrature` supplying the evaluation points
305: . dtype - `PETSC_SCALAR` or `PETSC_REAL`
306: . B - array of field values at each quadrature point, or `NULL`
307: . D - array of field reference derivatives at each quadrature point, or `NULL`
308: - H - array of field reference Hessians at each quadrature point, or `NULL`
310: Level: intermediate
312: Note:
313: If the shell only supplies a generic `DMFieldEvaluate()` via `DMFieldShellSetEvaluate()`, pass `DMFieldShellEvaluateFEDefault()` here.
315: .seealso: `DMField`, `DMFIELDSHELL`, `DMFieldCreateShell()`, `DMFieldEvaluateFE()`, `DMFieldShellEvaluateFEDefault()`, `DMFieldShellSetEvaluateFV()`
316: @*/
317: PetscErrorCode DMFieldShellSetEvaluateFE(DMField field, PetscErrorCode (*evaluateFE)(DMField field, IS is, PetscQuadrature quad, PetscDataType dtype, void *B, void *D, void *H))
318: {
319: PetscFunctionBegin;
321: field->ops->evaluateFE = evaluateFE;
322: PetscFunctionReturn(PETSC_SUCCESS);
323: }
325: /*@C
326: DMFieldShellSetEvaluateFV - Register the routine that evaluates a `DMFIELDSHELL` as cell averages over a set of mesh cells.
328: Logically Collective
330: Input Parameters:
331: + field - the `DMField` of type `DMFIELDSHELL`
332: - evaluateFV - the FV evaluation callback
334: Calling sequence of `evaluateFV`:
335: + field - the `DMField` of type `DMFIELDSHELL`
336: . is - the `IS` of mesh cells on which to evaluate the field
337: . dtype - `PETSC_SCALAR` or `PETSC_REAL`
338: . B - array of cell-averaged field values, or `NULL`
339: . D - array of cell-averaged field derivatives, or `NULL`
340: - H - array of cell-averaged field Hessians, or `NULL`
342: Level: intermediate
344: Note:
345: If the shell only supplies a generic `DMFieldEvaluate()` via `DMFieldShellSetEvaluate()`, pass `DMFieldShellEvaluateFVDefault()` here.
347: .seealso: `DMField`, `DMFIELDSHELL`, `DMFieldCreateShell()`, `DMFieldEvaluateFV()`, `DMFieldShellEvaluateFVDefault()`, `DMFieldShellSetEvaluateFE()`
348: @*/
349: PetscErrorCode DMFieldShellSetEvaluateFV(DMField field, PetscErrorCode (*evaluateFV)(DMField field, IS is, PetscDataType dtype, void *B, void *D, void *H))
350: {
351: PetscFunctionBegin;
353: field->ops->evaluateFV = evaluateFV;
354: PetscFunctionReturn(PETSC_SUCCESS);
355: }
357: /*@C
358: DMFieldShellSetGetDegree - Register the routine that reports the polynomial degree bounds of a `DMFIELDSHELL` over a set of mesh points.
360: Logically Collective
362: Input Parameters:
363: + field - the `DMField` of type `DMFIELDSHELL`
364: - getDegree - callback that returns the minimum and maximum polynomial degrees of the field over the given point `IS`
366: Calling sequence of `getDegree`:
367: + field - the `DMField` of type `DMFIELDSHELL`
368: . is - the `IS` of mesh points over which the degree bounds are requested
369: . minDegree - the degree of the largest polynomial space contained in the field on each element
370: - maxDegree - the largest degree of the smallest polynomial space containing the field on any element
372: Level: intermediate
374: .seealso: `DMField`, `DMFIELDSHELL`, `DMFieldCreateShell()`, `DMFieldGetDegree()`
375: @*/
376: PetscErrorCode DMFieldShellSetGetDegree(DMField field, PetscErrorCode (*getDegree)(DMField field, IS is, PetscInt *minDegree, PetscInt *maxDegree))
377: {
378: PetscFunctionBegin;
380: field->ops->getDegree = getDegree;
381: PetscFunctionReturn(PETSC_SUCCESS);
382: }
384: /*@C
385: DMFieldShellSetCreateDefaultQuadrature - Register the routine that supplies a default `PetscQuadrature` sufficient to integrate a `DMFIELDSHELL` exactly over a set of mesh points.
387: Logically Collective
389: Input Parameters:
390: + field - the `DMField` of type `DMFIELDSHELL`
391: - create - callback that returns a newly created `PetscQuadrature` for the given point `IS`
393: Calling sequence of `create`:
394: + f - the `DMField` of type `DMFIELDSHELL`
395: . is - the `IS` of mesh points over which the field will be integrated
396: - quad - the newly created `PetscQuadrature`
398: Level: intermediate
400: .seealso: `DMField`, `DMFIELDSHELL`, `DMFieldCreateShell()`, `DMFieldCreateDefaultQuadrature()`
401: @*/
402: PetscErrorCode DMFieldShellSetCreateDefaultQuadrature(DMField field, PetscErrorCode (*create)(DMField f, IS is, PetscQuadrature *quad))
403: {
404: PetscFunctionBegin;
406: field->ops->createDefaultQuadrature = create;
407: PetscFunctionReturn(PETSC_SUCCESS);
408: }
410: static PetscErrorCode DMFieldInitialize_Shell(DMField field)
411: {
412: PetscFunctionBegin;
413: field->ops->destroy = DMFieldDestroy_Shell;
414: field->ops->evaluate = NULL;
415: field->ops->evaluateFE = DMFieldShellEvaluateFEDefault;
416: field->ops->evaluateFV = DMFieldShellEvaluateFVDefault;
417: field->ops->getDegree = NULL;
418: field->ops->createDefaultQuadrature = NULL;
419: field->ops->view = NULL;
420: PetscFunctionReturn(PETSC_SUCCESS);
421: }
423: PETSC_INTERN PetscErrorCode DMFieldCreate_Shell(DMField field)
424: {
425: DMField_Shell *shell;
427: PetscFunctionBegin;
428: PetscCall(PetscNew(&shell));
429: field->data = shell;
430: PetscCall(DMFieldInitialize_Shell(field));
431: PetscFunctionReturn(PETSC_SUCCESS);
432: }
434: /*@
435: DMFieldCreateShell - Create a `DMFIELDSHELL`, a `DMField` whose evaluation is implemented entirely by user-supplied callbacks.
437: Collective
439: Input Parameters:
440: + dm - the `DM` on which the field lives
441: . numComponents - the number of components of the field
442: . continuity - the continuity of the field (e.g. `DMFIELD_VERTEX`)
443: - ctx - optional application context returned by `DMFieldShellGetContext()`
445: Output Parameter:
446: . field - the newly created `DMField` of type `DMFIELDSHELL`
448: Level: intermediate
450: Note:
451: After creation the user must register the desired evaluation callbacks with `DMFieldShellSetEvaluate()`, `DMFieldShellSetEvaluateFE()`, `DMFieldShellSetEvaluateFV()`, and optionally `DMFieldShellSetDestroy()`, `DMFieldShellSetGetDegree()`, and `DMFieldShellSetCreateDefaultQuadrature()`.
453: .seealso: `DMField`, `DMFIELDSHELL`, `DMFieldShellGetContext()`, `DMFieldShellSetEvaluate()`, `DMFieldShellSetEvaluateFE()`, `DMFieldShellSetEvaluateFV()`, `DMFieldShellSetDestroy()`
454: @*/
455: PetscErrorCode DMFieldCreateShell(DM dm, PetscInt numComponents, DMFieldContinuity continuity, PetscCtx ctx, DMField *field)
456: {
457: DMField b;
458: DMField_Shell *shell;
460: PetscFunctionBegin;
462: if (ctx) PetscAssertPointer(ctx, 4);
463: PetscAssertPointer(field, 5);
464: PetscCall(DMFieldCreate(dm, numComponents, continuity, &b));
465: PetscCall(DMFieldSetType(b, DMFIELDSHELL));
466: shell = (DMField_Shell *)b->data;
467: shell->ctx = ctx;
468: *field = b;
469: PetscFunctionReturn(PETSC_SUCCESS);
470: }