Actual source code: dmfieldshell.c
1: #include <petsc/private/dmfieldimpl.h>
3: typedef struct _n_DMField_Shell {
4: void *ctx;
5: PetscErrorCode (*destroy)(DMField);
6: } DMField_Shell;
8: PetscErrorCode DMFieldShellGetContext(DMField field, void *ctx)
9: {
10: PetscBool flg;
12: PetscFunctionBegin;
14: PetscAssertPointer(ctx, 2);
15: PetscCall(PetscObjectTypeCompare((PetscObject)field, DMFIELDSHELL, &flg));
16: if (flg) *(void **)ctx = ((DMField_Shell *)field->data)->ctx;
17: else SETERRQ(PetscObjectComm((PetscObject)field), PETSC_ERR_SUP, "Cannot get context from non-shell shield");
18: PetscFunctionReturn(PETSC_SUCCESS);
19: }
21: static PetscErrorCode DMFieldDestroy_Shell(DMField field)
22: {
23: DMField_Shell *shell = (DMField_Shell *)field->data;
25: PetscFunctionBegin;
26: if (shell->destroy) PetscCall((*shell->destroy)(field));
27: PetscCall(PetscFree(field->data));
28: PetscFunctionReturn(PETSC_SUCCESS);
29: }
31: PetscErrorCode DMFieldShellEvaluateFEDefault(DMField field, IS pointIS, PetscQuadrature quad, PetscDataType type, void *B, void *D, void *H)
32: {
33: DM dm = field->dm;
34: DMField coordField;
35: PetscFEGeom *geom;
36: Vec pushforward;
37: PetscInt dimC, dim, numPoints, Nq, p, Nc;
38: PetscScalar *pfArray;
40: PetscFunctionBegin;
41: Nc = field->numComponents;
42: PetscCall(DMGetCoordinateField(dm, &coordField));
43: PetscCall(DMFieldCreateFEGeom(coordField, pointIS, quad, PETSC_FALSE, &geom));
44: PetscCall(DMGetCoordinateDim(dm, &dimC));
45: PetscCall(PetscQuadratureGetData(quad, &dim, NULL, &Nq, NULL, NULL));
46: PetscCall(ISGetLocalSize(pointIS, &numPoints));
47: PetscCall(PetscMalloc1(dimC * Nq * numPoints, &pfArray));
48: for (p = 0; p < numPoints * dimC * Nq; p++) pfArray[p] = (PetscScalar)geom->v[p];
49: PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)pointIS), dimC, dimC * Nq * numPoints, PETSC_DETERMINE, pfArray, &pushforward));
50: PetscCall(DMFieldEvaluate(field, pushforward, type, B, D, H));
51: /* TODO: handle covariant/contravariant pullbacks */
52: if (D) {
53: if (type == PETSC_SCALAR) {
54: PetscScalar *sD = (PetscScalar *)D;
55: PetscInt q;
57: for (p = 0; p < numPoints * Nq; p++) {
58: for (q = 0; q < Nc; q++) {
59: PetscScalar d[3];
61: PetscInt i, j;
63: for (i = 0; i < dimC; i++) d[i] = sD[(p * Nc + q) * dimC + i];
64: for (i = 0; i < dimC; i++) sD[(p * Nc + q) * dimC + i] = 0.;
65: for (i = 0; i < dimC; i++) {
66: for (j = 0; j < dimC; j++) sD[(p * Nc + q) * dimC + i] += geom->J[(p * dimC + j) * dimC + i] * d[j];
67: }
68: }
69: }
70: } else {
71: PetscReal *rD = (PetscReal *)D;
72: PetscInt q;
74: for (p = 0; p < numPoints * Nq; p++) {
75: for (q = 0; q < Nc; q++) {
76: PetscReal d[3];
78: PetscInt i, j;
80: for (i = 0; i < dimC; i++) d[i] = rD[(p * Nc + q) * dimC + i];
81: for (i = 0; i < dimC; i++) rD[(p * Nc + q) * dimC + i] = 0.;
82: for (i = 0; i < dimC; i++) {
83: for (j = 0; j < dimC; j++) rD[(p * Nc + q) * dimC + i] += geom->J[(p * dimC + j) * dimC + i] * d[j];
84: }
85: }
86: }
87: }
88: }
89: if (H) {
90: if (type == PETSC_SCALAR) {
91: PetscScalar *sH = (PetscScalar *)H;
92: PetscInt q;
94: for (p = 0; p < numPoints * Nq; p++) {
95: for (q = 0; q < Nc; q++) {
96: PetscScalar d[3][3];
98: PetscInt i, j, k, l;
100: for (i = 0; i < dimC; i++)
101: for (j = 0; j < dimC; j++) d[i][j] = sH[((p * Nc + q) * dimC + i) * dimC + j];
102: for (i = 0; i < dimC; i++)
103: for (j = 0; j < dimC; j++) sH[((p * Nc + q) * dimC + i) * dimC + j] = 0.;
104: for (i = 0; i < dimC; i++) {
105: for (j = 0; j < dimC; j++) {
106: for (k = 0; k < dimC; k++) {
107: 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];
108: }
109: }
110: }
111: }
112: }
113: } else {
114: PetscReal *rH = (PetscReal *)H;
115: PetscInt q;
117: for (p = 0; p < numPoints * Nq; p++) {
118: for (q = 0; q < Nc; q++) {
119: PetscReal d[3][3];
121: PetscInt i, j, k, l;
123: for (i = 0; i < dimC; i++)
124: for (j = 0; j < dimC; j++) d[i][j] = rH[((p * Nc + q) * dimC + i) * dimC + j];
125: for (i = 0; i < dimC; i++)
126: for (j = 0; j < dimC; j++) rH[((p * Nc + q) * dimC + i) * dimC + j] = 0.;
127: for (i = 0; i < dimC; i++) {
128: for (j = 0; j < dimC; j++) {
129: for (k = 0; k < dimC; k++) {
130: 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];
131: }
132: }
133: }
134: }
135: }
136: }
137: }
138: PetscCall(VecDestroy(&pushforward));
139: PetscCall(PetscFree(pfArray));
140: PetscCall(PetscFEGeomDestroy(&geom));
141: PetscFunctionReturn(PETSC_SUCCESS);
142: }
144: PetscErrorCode DMFieldShellEvaluateFVDefault(DMField field, IS pointIS, PetscDataType type, void *B, void *D, void *H)
145: {
146: DM dm = field->dm;
147: DMField coordField;
148: PetscFEGeom *geom;
149: Vec pushforward;
150: PetscInt dimC, dim, numPoints, Nq, p;
151: PetscScalar *pfArray;
152: PetscQuadrature quad;
153: MPI_Comm comm;
155: PetscFunctionBegin;
156: PetscCall(PetscObjectGetComm((PetscObject)field, &comm));
157: PetscCall(DMGetDimension(dm, &dim));
158: PetscCall(DMGetCoordinateDim(dm, &dimC));
159: PetscCall(DMGetCoordinateField(dm, &coordField));
160: PetscCall(DMFieldGetFVQuadrature_Internal(coordField, pointIS, &quad));
161: PetscCheck(quad, comm, PETSC_ERR_ARG_WRONGSTATE, "coordinate field must have default quadrature for FV computation");
162: PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL));
163: PetscCheck(Nq == 1, comm, PETSC_ERR_ARG_WRONGSTATE, "quadrature must have only one point");
164: PetscCall(DMFieldCreateFEGeom(coordField, pointIS, quad, PETSC_FALSE, &geom));
165: PetscCall(ISGetLocalSize(pointIS, &numPoints));
166: PetscCall(PetscMalloc1(dimC * numPoints, &pfArray));
167: for (p = 0; p < numPoints * dimC; p++) pfArray[p] = (PetscScalar)geom->v[p];
168: PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)pointIS), dimC, dimC * numPoints, PETSC_DETERMINE, pfArray, &pushforward));
169: PetscCall(DMFieldEvaluate(field, pushforward, type, B, D, H));
170: PetscCall(PetscQuadratureDestroy(&quad));
171: PetscCall(VecDestroy(&pushforward));
172: PetscCall(PetscFree(pfArray));
173: PetscCall(PetscFEGeomDestroy(&geom));
174: PetscFunctionReturn(PETSC_SUCCESS);
175: }
177: PetscErrorCode DMFieldShellSetDestroy(DMField field, PetscErrorCode (*destroy)(DMField))
178: {
179: DMField_Shell *shell = (DMField_Shell *)field->data;
181: PetscFunctionBegin;
183: shell->destroy = destroy;
184: PetscFunctionReturn(PETSC_SUCCESS);
185: }
187: PetscErrorCode DMFieldShellSetEvaluate(DMField field, PetscErrorCode (*evaluate)(DMField, Vec, PetscDataType, void *, void *, void *))
188: {
189: PetscFunctionBegin;
191: field->ops->evaluate = evaluate;
192: PetscFunctionReturn(PETSC_SUCCESS);
193: }
195: PetscErrorCode DMFieldShellSetEvaluateFE(DMField field, PetscErrorCode (*evaluateFE)(DMField, IS, PetscQuadrature, PetscDataType, void *, void *, void *))
196: {
197: PetscFunctionBegin;
199: field->ops->evaluateFE = evaluateFE;
200: PetscFunctionReturn(PETSC_SUCCESS);
201: }
203: PetscErrorCode DMFieldShellSetEvaluateFV(DMField field, PetscErrorCode (*evaluateFV)(DMField, IS, PetscDataType, void *, void *, void *))
204: {
205: PetscFunctionBegin;
207: field->ops->evaluateFV = evaluateFV;
208: PetscFunctionReturn(PETSC_SUCCESS);
209: }
211: PetscErrorCode DMFieldShellSetGetDegree(DMField field, PetscErrorCode (*getDegree)(DMField, IS, PetscInt *, PetscInt *))
212: {
213: PetscFunctionBegin;
215: field->ops->getDegree = getDegree;
216: PetscFunctionReturn(PETSC_SUCCESS);
217: }
219: PetscErrorCode DMFieldShellSetCreateDefaultQuadrature(DMField field, PetscErrorCode (*createDefaultQuadrature)(DMField, IS, PetscQuadrature *))
220: {
221: PetscFunctionBegin;
223: field->ops->createDefaultQuadrature = createDefaultQuadrature;
224: PetscFunctionReturn(PETSC_SUCCESS);
225: }
227: static PetscErrorCode DMFieldInitialize_Shell(DMField field)
228: {
229: PetscFunctionBegin;
230: field->ops->destroy = DMFieldDestroy_Shell;
231: field->ops->evaluate = NULL;
232: field->ops->evaluateFE = DMFieldShellEvaluateFEDefault;
233: field->ops->evaluateFV = DMFieldShellEvaluateFVDefault;
234: field->ops->getDegree = NULL;
235: field->ops->createDefaultQuadrature = NULL;
236: field->ops->view = NULL;
237: PetscFunctionReturn(PETSC_SUCCESS);
238: }
240: PETSC_INTERN PetscErrorCode DMFieldCreate_Shell(DMField field)
241: {
242: DMField_Shell *shell;
244: PetscFunctionBegin;
245: PetscCall(PetscNew(&shell));
246: field->data = shell;
247: PetscCall(DMFieldInitialize_Shell(field));
248: PetscFunctionReturn(PETSC_SUCCESS);
249: }
251: PetscErrorCode DMFieldCreateShell(DM dm, PetscInt numComponents, DMFieldContinuity continuity, void *ctx, DMField *field)
252: {
253: DMField b;
254: DMField_Shell *shell;
256: PetscFunctionBegin;
258: if (ctx) PetscAssertPointer(ctx, 4);
259: PetscAssertPointer(field, 5);
260: PetscCall(DMFieldCreate(dm, numComponents, continuity, &b));
261: PetscCall(DMFieldSetType(b, DMFIELDSHELL));
262: shell = (DMField_Shell *)b->data;
263: shell->ctx = ctx;
264: *field = b;
265: PetscFunctionReturn(PETSC_SUCCESS);
266: }