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;
94: for (p = 0; p < numPoints * Nq; p++) {
95: for (PetscInt q = 0; q < Nc; q++) {
96: PetscScalar d[3];
98: for (PetscInt i = 0; i < dimC; i++) d[i] = sD[(p * Nc + q) * dimC + i];
99: for (PetscInt i = 0; i < dimC; i++) sD[(p * Nc + q) * dimC + i] = 0.;
100: for (PetscInt i = 0; i < dimC; i++) {
101: for (PetscInt j = 0; j < dimC; j++) sD[(p * Nc + q) * dimC + i] += geom->J[(p * dimC + j) * dimC + i] * d[j];
102: }
103: }
104: }
105: } else {
106: PetscReal *rD = (PetscReal *)D;
108: for (p = 0; p < numPoints * Nq; p++) {
109: for (PetscInt q = 0; q < Nc; q++) {
110: PetscReal d[3];
112: for (PetscInt i = 0; i < dimC; i++) d[i] = rD[(p * Nc + q) * dimC + i];
113: for (PetscInt i = 0; i < dimC; i++) rD[(p * Nc + q) * dimC + i] = 0.;
114: for (PetscInt i = 0; i < dimC; i++) {
115: for (PetscInt j = 0; j < dimC; j++) rD[(p * Nc + q) * dimC + i] += geom->J[(p * dimC + j) * dimC + i] * d[j];
116: }
117: }
118: }
119: }
120: }
121: if (H) {
122: if (type == PETSC_SCALAR) {
123: PetscScalar *sH = (PetscScalar *)H;
125: for (p = 0; p < numPoints * Nq; p++) {
126: for (PetscInt q = 0; q < Nc; q++) {
127: PetscScalar d[3][3];
129: PetscInt i, j, k, l;
131: for (i = 0; i < dimC; i++)
132: for (j = 0; j < dimC; j++) d[i][j] = sH[((p * Nc + q) * dimC + i) * dimC + j];
133: for (i = 0; i < dimC; i++)
134: for (j = 0; j < dimC; j++) sH[((p * Nc + q) * dimC + i) * dimC + j] = 0.;
135: for (i = 0; i < dimC; i++) {
136: for (j = 0; j < dimC; j++) {
137: for (k = 0; k < dimC; k++) {
138: 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];
139: }
140: }
141: }
142: }
143: }
144: } else {
145: PetscReal *rH = (PetscReal *)H;
147: for (p = 0; p < numPoints * Nq; p++) {
148: for (PetscInt q = 0; q < Nc; q++) {
149: PetscReal d[3][3];
151: PetscInt i, j, k, l;
153: for (i = 0; i < dimC; i++)
154: for (j = 0; j < dimC; j++) d[i][j] = rH[((p * Nc + q) * dimC + i) * dimC + j];
155: for (i = 0; i < dimC; i++)
156: for (j = 0; j < dimC; j++) rH[((p * Nc + q) * dimC + i) * dimC + j] = 0.;
157: for (i = 0; i < dimC; i++) {
158: for (j = 0; j < dimC; j++) {
159: for (k = 0; k < dimC; k++) {
160: 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];
161: }
162: }
163: }
164: }
165: }
166: }
167: }
168: PetscCall(VecDestroy(&pushforward));
169: PetscCall(PetscFree(pfArray));
170: PetscCall(PetscFEGeomDestroy(&geom));
171: PetscFunctionReturn(PETSC_SUCCESS);
172: }
174: /*@C
175: DMFieldShellEvaluateFVDefault - Default finite-volume evaluation for a `DMFIELDSHELL` that samples at cell centroids using the coordinate `DMField`'s default quadrature and calls `DMFieldEvaluate()`.
177: Not Collective
179: Input Parameters:
180: + field - the `DMField` of type `DMFIELDSHELL`
181: . pointIS - the `IS` of mesh cells at which to evaluate
182: - type - `PETSC_SCALAR` or `PETSC_REAL`
184: Output Parameters:
185: + B - cell-averaged values, or `NULL`
186: . D - cell-averaged derivatives, or `NULL`
187: - H - cell-averaged Hessians, or `NULL`
189: Level: developer
191: Note:
192: Intended to be registered as the FV evaluation callback via `DMFieldShellSetEvaluateFV()` when the shell only supplies a bulk `DMFieldEvaluate()` implementation.
194: .seealso: `DMField`, `DMFIELDSHELL`, `DMFieldShellSetEvaluateFV()`, `DMFieldShellEvaluateFEDefault()`, `DMFieldEvaluate()`
195: @*/
196: PetscErrorCode DMFieldShellEvaluateFVDefault(DMField field, IS pointIS, PetscDataType type, void *B, void *D, void *H)
197: {
198: DM dm = field->dm;
199: DMField coordField;
200: PetscFEGeom *geom;
201: Vec pushforward;
202: PetscInt dimC, dim, numPoints, Nq, p;
203: PetscScalar *pfArray;
204: PetscQuadrature quad;
205: MPI_Comm comm;
207: PetscFunctionBegin;
208: PetscCall(PetscObjectGetComm((PetscObject)field, &comm));
209: PetscCall(DMGetDimension(dm, &dim));
210: PetscCall(DMGetCoordinateDim(dm, &dimC));
211: PetscCall(DMGetCoordinateField(dm, &coordField));
212: PetscCall(DMFieldGetFVQuadrature_Internal(coordField, pointIS, &quad));
213: PetscCheck(quad, comm, PETSC_ERR_ARG_WRONGSTATE, "coordinate field must have default quadrature for FV computation");
214: PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL));
215: PetscCheck(Nq == 1, comm, PETSC_ERR_ARG_WRONGSTATE, "quadrature must have only one point");
216: PetscCall(DMFieldCreateFEGeom(coordField, pointIS, quad, PETSC_FEGEOM_BASIC, &geom));
217: PetscCall(ISGetLocalSize(pointIS, &numPoints));
218: PetscCall(PetscMalloc1(dimC * numPoints, &pfArray));
219: for (p = 0; p < numPoints * dimC; p++) pfArray[p] = geom->v[p];
220: PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)pointIS), dimC, dimC * numPoints, PETSC_DETERMINE, pfArray, &pushforward));
221: PetscCall(DMFieldEvaluate(field, pushforward, type, B, D, H));
222: PetscCall(PetscQuadratureDestroy(&quad));
223: PetscCall(VecDestroy(&pushforward));
224: PetscCall(PetscFree(pfArray));
225: PetscCall(PetscFEGeomDestroy(&geom));
226: PetscFunctionReturn(PETSC_SUCCESS);
227: }
229: /*@C
230: DMFieldShellSetDestroy - Register a destroy callback that will be invoked when a `DMFIELDSHELL` is destroyed.
232: Logically Collective
234: Input Parameters:
235: + field - the `DMField` of type `DMFIELDSHELL`
236: - destroy - the destroy routine, called before the shell's own data is freed
238: Calling sequence of `destroy`:
239: . field - the `DMField` of type `DMFIELDSHELL` being destroyed
241: Level: intermediate
243: .seealso: `DMField`, `DMFIELDSHELL`, `DMFieldCreateShell()`, `DMFieldDestroy()`
244: @*/
245: PetscErrorCode DMFieldShellSetDestroy(DMField field, PetscErrorCode (*destroy)(DMField field))
246: {
247: DMField_Shell *shell = (DMField_Shell *)field->data;
249: PetscFunctionBegin;
251: shell->destroy = destroy;
252: PetscFunctionReturn(PETSC_SUCCESS);
253: }
255: /*@C
256: DMFieldShellSetEvaluate - Register the routine that evaluates a `DMFIELDSHELL` at an arbitrary set of real-space points supplied as a `Vec` of coordinates.
258: Logically Collective
260: Input Parameters:
261: + field - the `DMField` of type `DMFIELDSHELL`
262: - evaluate - the evaluation callback
264: Calling sequence of `evaluate`:
265: + field - the `DMField` of type `DMFIELDSHELL`
266: . u - the points at which to evaluate the field, as a `Vec` of coordinates of size d x n
267: . dtype - `PETSC_SCALAR` or `PETSC_REAL`
268: . B - array of field values at each point, or `NULL`
269: . D - array of field spatial derivatives at each point, or `NULL`
270: - H - array of field spatial Hessians at each point, or `NULL`
272: Level: intermediate
274: .seealso: `DMField`, `DMFIELDSHELL`, `DMFieldCreateShell()`, `DMFieldEvaluate()`, `DMFieldShellSetEvaluateFE()`, `DMFieldShellSetEvaluateFV()`
275: @*/
276: PetscErrorCode DMFieldShellSetEvaluate(DMField field, PetscErrorCode (*evaluate)(DMField field, Vec u, PetscDataType dtype, void *B, void *D, void *H))
277: {
278: PetscFunctionBegin;
280: field->ops->evaluate = evaluate;
281: PetscFunctionReturn(PETSC_SUCCESS);
282: }
284: /*@C
285: DMFieldShellSetEvaluateFE - Register the routine that evaluates a `DMFIELDSHELL` at finite-element quadrature points over a set of mesh points.
287: Logically Collective
289: Input Parameters:
290: + field - the `DMField` of type `DMFIELDSHELL`
291: - evaluateFE - the FE evaluation callback
293: Calling sequence of `evaluateFE`:
294: + field - the `DMField` of type `DMFIELDSHELL`
295: . is - the `IS` of mesh cells on which to evaluate the field
296: . quad - the reference-cell `PetscQuadrature` supplying the evaluation points
297: . dtype - `PETSC_SCALAR` or `PETSC_REAL`
298: . B - array of field values at each quadrature point, or `NULL`
299: . D - array of field reference derivatives at each quadrature point, or `NULL`
300: - H - array of field reference Hessians at each quadrature point, or `NULL`
302: Level: intermediate
304: Note:
305: If the shell only supplies a generic `DMFieldEvaluate()` via `DMFieldShellSetEvaluate()`, pass `DMFieldShellEvaluateFEDefault()` here.
307: .seealso: `DMField`, `DMFIELDSHELL`, `DMFieldCreateShell()`, `DMFieldEvaluateFE()`, `DMFieldShellEvaluateFEDefault()`, `DMFieldShellSetEvaluateFV()`
308: @*/
309: PetscErrorCode DMFieldShellSetEvaluateFE(DMField field, PetscErrorCode (*evaluateFE)(DMField field, IS is, PetscQuadrature quad, PetscDataType dtype, void *B, void *D, void *H))
310: {
311: PetscFunctionBegin;
313: field->ops->evaluateFE = evaluateFE;
314: PetscFunctionReturn(PETSC_SUCCESS);
315: }
317: /*@C
318: DMFieldShellSetEvaluateFV - Register the routine that evaluates a `DMFIELDSHELL` as cell averages over a set of mesh cells.
320: Logically Collective
322: Input Parameters:
323: + field - the `DMField` of type `DMFIELDSHELL`
324: - evaluateFV - the FV evaluation callback
326: Calling sequence of `evaluateFV`:
327: + field - the `DMField` of type `DMFIELDSHELL`
328: . is - the `IS` of mesh cells on which to evaluate the field
329: . dtype - `PETSC_SCALAR` or `PETSC_REAL`
330: . B - array of cell-averaged field values, or `NULL`
331: . D - array of cell-averaged field derivatives, or `NULL`
332: - H - array of cell-averaged field Hessians, or `NULL`
334: Level: intermediate
336: Note:
337: If the shell only supplies a generic `DMFieldEvaluate()` via `DMFieldShellSetEvaluate()`, pass `DMFieldShellEvaluateFVDefault()` here.
339: .seealso: `DMField`, `DMFIELDSHELL`, `DMFieldCreateShell()`, `DMFieldEvaluateFV()`, `DMFieldShellEvaluateFVDefault()`, `DMFieldShellSetEvaluateFE()`
340: @*/
341: PetscErrorCode DMFieldShellSetEvaluateFV(DMField field, PetscErrorCode (*evaluateFV)(DMField field, IS is, PetscDataType dtype, void *B, void *D, void *H))
342: {
343: PetscFunctionBegin;
345: field->ops->evaluateFV = evaluateFV;
346: PetscFunctionReturn(PETSC_SUCCESS);
347: }
349: /*@C
350: DMFieldShellSetGetDegree - Register the routine that reports the polynomial degree bounds of a `DMFIELDSHELL` over a set of mesh points.
352: Logically Collective
354: Input Parameters:
355: + field - the `DMField` of type `DMFIELDSHELL`
356: - getDegree - callback that returns the minimum and maximum polynomial degrees of the field over the given point `IS`
358: Calling sequence of `getDegree`:
359: + field - the `DMField` of type `DMFIELDSHELL`
360: . is - the `IS` of mesh points over which the degree bounds are requested
361: . minDegree - the degree of the largest polynomial space contained in the field on each element
362: - maxDegree - the largest degree of the smallest polynomial space containing the field on any element
364: Level: intermediate
366: .seealso: `DMField`, `DMFIELDSHELL`, `DMFieldCreateShell()`, `DMFieldGetDegree()`
367: @*/
368: PetscErrorCode DMFieldShellSetGetDegree(DMField field, PetscErrorCode (*getDegree)(DMField field, IS is, PetscInt *minDegree, PetscInt *maxDegree))
369: {
370: PetscFunctionBegin;
372: field->ops->getDegree = getDegree;
373: PetscFunctionReturn(PETSC_SUCCESS);
374: }
376: /*@C
377: DMFieldShellSetCreateDefaultQuadrature - Register the routine that supplies a default `PetscQuadrature` sufficient to integrate a `DMFIELDSHELL` exactly over a set of mesh points.
379: Logically Collective
381: Input Parameters:
382: + field - the `DMField` of type `DMFIELDSHELL`
383: - create - callback that returns a newly created `PetscQuadrature` for the given point `IS`
385: Calling sequence of `create`:
386: + f - the `DMField` of type `DMFIELDSHELL`
387: . is - the `IS` of mesh points over which the field will be integrated
388: - quad - the newly created `PetscQuadrature`
390: Level: intermediate
392: .seealso: `DMField`, `DMFIELDSHELL`, `DMFieldCreateShell()`, `DMFieldCreateDefaultQuadrature()`
393: @*/
394: PetscErrorCode DMFieldShellSetCreateDefaultQuadrature(DMField field, PetscErrorCode (*create)(DMField f, IS is, PetscQuadrature *quad))
395: {
396: PetscFunctionBegin;
398: field->ops->createDefaultQuadrature = create;
399: PetscFunctionReturn(PETSC_SUCCESS);
400: }
402: static PetscErrorCode DMFieldInitialize_Shell(DMField field)
403: {
404: PetscFunctionBegin;
405: field->ops->destroy = DMFieldDestroy_Shell;
406: field->ops->evaluate = NULL;
407: field->ops->evaluateFE = DMFieldShellEvaluateFEDefault;
408: field->ops->evaluateFV = DMFieldShellEvaluateFVDefault;
409: field->ops->getDegree = NULL;
410: field->ops->createDefaultQuadrature = NULL;
411: field->ops->view = NULL;
412: PetscFunctionReturn(PETSC_SUCCESS);
413: }
415: PETSC_INTERN PetscErrorCode DMFieldCreate_Shell(DMField field)
416: {
417: DMField_Shell *shell;
419: PetscFunctionBegin;
420: PetscCall(PetscNew(&shell));
421: field->data = shell;
422: PetscCall(DMFieldInitialize_Shell(field));
423: PetscFunctionReturn(PETSC_SUCCESS);
424: }
426: /*@
427: DMFieldCreateShell - Create a `DMFIELDSHELL`, a `DMField` whose evaluation is implemented entirely by user-supplied callbacks.
429: Collective
431: Input Parameters:
432: + dm - the `DM` on which the field lives
433: . numComponents - the number of components of the field
434: . continuity - the continuity of the field (e.g. `DMFIELD_VERTEX`)
435: - ctx - optional application context returned by `DMFieldShellGetContext()`
437: Output Parameter:
438: . field - the newly created `DMField` of type `DMFIELDSHELL`
440: Level: intermediate
442: Note:
443: After creation the user must register the desired evaluation callbacks with `DMFieldShellSetEvaluate()`, `DMFieldShellSetEvaluateFE()`, `DMFieldShellSetEvaluateFV()`, and optionally `DMFieldShellSetDestroy()`, `DMFieldShellSetGetDegree()`, and `DMFieldShellSetCreateDefaultQuadrature()`.
445: .seealso: `DMField`, `DMFIELDSHELL`, `DMFieldShellGetContext()`, `DMFieldShellSetEvaluate()`, `DMFieldShellSetEvaluateFE()`, `DMFieldShellSetEvaluateFV()`, `DMFieldShellSetDestroy()`
446: @*/
447: PetscErrorCode DMFieldCreateShell(DM dm, PetscInt numComponents, DMFieldContinuity continuity, PetscCtx ctx, DMField *field)
448: {
449: DMField b;
450: DMField_Shell *shell;
452: PetscFunctionBegin;
454: if (ctx) PetscAssertPointer(ctx, 4);
455: PetscAssertPointer(field, 5);
456: PetscCall(DMFieldCreate(dm, numComponents, continuity, &b));
457: PetscCall(DMFieldSetType(b, DMFIELDSHELL));
458: shell = (DMField_Shell *)b->data;
459: shell->ctx = ctx;
460: *field = b;
461: PetscFunctionReturn(PETSC_SUCCESS);
462: }