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: PetscErrorCode DMFieldShellGetContext(DMField field, PetscCtxRt ctx)
9: {
10: PetscBool flg;
12: PetscFunctionBegin;
14: PetscAssertPointer(ctx, 2);
15: PetscCall(PetscObjectTypeCompare((PetscObject)field, DMFIELDSHELL, &flg));
16: PetscCheck(flg, PetscObjectComm((PetscObject)field), PETSC_ERR_SUP, "Cannot get context from non-shell shield");
17: *(void **)ctx = ((DMField_Shell *)field->data)->ctx;
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_FEGEOM_BASIC, &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] = 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;
56: for (p = 0; p < numPoints * Nq; p++) {
57: for (PetscInt q = 0; q < Nc; q++) {
58: PetscScalar d[3];
60: for (PetscInt i = 0; i < dimC; i++) d[i] = sD[(p * Nc + q) * dimC + i];
61: for (PetscInt i = 0; i < dimC; i++) sD[(p * Nc + q) * dimC + i] = 0.;
62: for (PetscInt i = 0; i < dimC; i++) {
63: for (PetscInt j = 0; j < dimC; j++) sD[(p * Nc + q) * dimC + i] += geom->J[(p * dimC + j) * dimC + i] * d[j];
64: }
65: }
66: }
67: } else {
68: PetscReal *rD = (PetscReal *)D;
70: for (p = 0; p < numPoints * Nq; p++) {
71: for (PetscInt q = 0; q < Nc; q++) {
72: PetscReal d[3];
74: for (PetscInt i = 0; i < dimC; i++) d[i] = rD[(p * Nc + q) * dimC + i];
75: for (PetscInt i = 0; i < dimC; i++) rD[(p * Nc + q) * dimC + i] = 0.;
76: for (PetscInt i = 0; i < dimC; i++) {
77: for (PetscInt j = 0; j < dimC; j++) rD[(p * Nc + q) * dimC + i] += geom->J[(p * dimC + j) * dimC + i] * d[j];
78: }
79: }
80: }
81: }
82: }
83: if (H) {
84: if (type == PETSC_SCALAR) {
85: PetscScalar *sH = (PetscScalar *)H;
87: for (p = 0; p < numPoints * Nq; p++) {
88: for (PetscInt q = 0; q < Nc; q++) {
89: PetscScalar d[3][3];
91: PetscInt i, j, k, l;
93: for (i = 0; i < dimC; i++)
94: for (j = 0; j < dimC; j++) d[i][j] = sH[((p * Nc + q) * dimC + i) * dimC + j];
95: for (i = 0; i < dimC; i++)
96: for (j = 0; j < dimC; j++) sH[((p * Nc + q) * dimC + i) * dimC + j] = 0.;
97: for (i = 0; i < dimC; i++) {
98: for (j = 0; j < dimC; j++) {
99: for (k = 0; k < dimC; k++) {
100: 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];
101: }
102: }
103: }
104: }
105: }
106: } else {
107: PetscReal *rH = (PetscReal *)H;
109: for (p = 0; p < numPoints * Nq; p++) {
110: for (PetscInt q = 0; q < Nc; q++) {
111: PetscReal d[3][3];
113: PetscInt i, j, k, l;
115: for (i = 0; i < dimC; i++)
116: for (j = 0; j < dimC; j++) d[i][j] = rH[((p * Nc + q) * dimC + i) * dimC + j];
117: for (i = 0; i < dimC; i++)
118: for (j = 0; j < dimC; j++) rH[((p * Nc + q) * dimC + i) * dimC + j] = 0.;
119: for (i = 0; i < dimC; i++) {
120: for (j = 0; j < dimC; j++) {
121: for (k = 0; k < dimC; k++) {
122: 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];
123: }
124: }
125: }
126: }
127: }
128: }
129: }
130: PetscCall(VecDestroy(&pushforward));
131: PetscCall(PetscFree(pfArray));
132: PetscCall(PetscFEGeomDestroy(&geom));
133: PetscFunctionReturn(PETSC_SUCCESS);
134: }
136: PetscErrorCode DMFieldShellEvaluateFVDefault(DMField field, IS pointIS, PetscDataType type, void *B, void *D, void *H)
137: {
138: DM dm = field->dm;
139: DMField coordField;
140: PetscFEGeom *geom;
141: Vec pushforward;
142: PetscInt dimC, dim, numPoints, Nq, p;
143: PetscScalar *pfArray;
144: PetscQuadrature quad;
145: MPI_Comm comm;
147: PetscFunctionBegin;
148: PetscCall(PetscObjectGetComm((PetscObject)field, &comm));
149: PetscCall(DMGetDimension(dm, &dim));
150: PetscCall(DMGetCoordinateDim(dm, &dimC));
151: PetscCall(DMGetCoordinateField(dm, &coordField));
152: PetscCall(DMFieldGetFVQuadrature_Internal(coordField, pointIS, &quad));
153: PetscCheck(quad, comm, PETSC_ERR_ARG_WRONGSTATE, "coordinate field must have default quadrature for FV computation");
154: PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL));
155: PetscCheck(Nq == 1, comm, PETSC_ERR_ARG_WRONGSTATE, "quadrature must have only one point");
156: PetscCall(DMFieldCreateFEGeom(coordField, pointIS, quad, PETSC_FEGEOM_BASIC, &geom));
157: PetscCall(ISGetLocalSize(pointIS, &numPoints));
158: PetscCall(PetscMalloc1(dimC * numPoints, &pfArray));
159: for (p = 0; p < numPoints * dimC; p++) pfArray[p] = geom->v[p];
160: PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)pointIS), dimC, dimC * numPoints, PETSC_DETERMINE, pfArray, &pushforward));
161: PetscCall(DMFieldEvaluate(field, pushforward, type, B, D, H));
162: PetscCall(PetscQuadratureDestroy(&quad));
163: PetscCall(VecDestroy(&pushforward));
164: PetscCall(PetscFree(pfArray));
165: PetscCall(PetscFEGeomDestroy(&geom));
166: PetscFunctionReturn(PETSC_SUCCESS);
167: }
169: PetscErrorCode DMFieldShellSetDestroy(DMField field, PetscErrorCode (*destroy)(DMField field))
170: {
171: DMField_Shell *shell = (DMField_Shell *)field->data;
173: PetscFunctionBegin;
175: shell->destroy = destroy;
176: PetscFunctionReturn(PETSC_SUCCESS);
177: }
179: PetscErrorCode DMFieldShellSetEvaluate(DMField field, PetscErrorCode (*evaluate)(DMField field, Vec u, PetscDataType dtype, void *B, void *D, void *H))
180: {
181: PetscFunctionBegin;
183: field->ops->evaluate = evaluate;
184: PetscFunctionReturn(PETSC_SUCCESS);
185: }
187: PetscErrorCode DMFieldShellSetEvaluateFE(DMField field, PetscErrorCode (*evaluateFE)(DMField field, IS is, PetscQuadrature quad, PetscDataType dtype, void *B, void *D, void *H))
188: {
189: PetscFunctionBegin;
191: field->ops->evaluateFE = evaluateFE;
192: PetscFunctionReturn(PETSC_SUCCESS);
193: }
195: PetscErrorCode DMFieldShellSetEvaluateFV(DMField field, PetscErrorCode (*evaluateFV)(DMField field, IS is, PetscDataType dtype, void *B, void *D, void *H))
196: {
197: PetscFunctionBegin;
199: field->ops->evaluateFV = evaluateFV;
200: PetscFunctionReturn(PETSC_SUCCESS);
201: }
203: PetscErrorCode DMFieldShellSetGetDegree(DMField field, PetscErrorCode (*getDegree)(DMField field, IS is, PetscInt *minDegree, PetscInt *maxDegree))
204: {
205: PetscFunctionBegin;
207: field->ops->getDegree = getDegree;
208: PetscFunctionReturn(PETSC_SUCCESS);
209: }
211: PetscErrorCode DMFieldShellSetCreateDefaultQuadrature(DMField field, PetscErrorCode (*createDefaultQuadrature)(DMField field, IS is, PetscQuadrature *quad))
212: {
213: PetscFunctionBegin;
215: field->ops->createDefaultQuadrature = createDefaultQuadrature;
216: PetscFunctionReturn(PETSC_SUCCESS);
217: }
219: static PetscErrorCode DMFieldInitialize_Shell(DMField field)
220: {
221: PetscFunctionBegin;
222: field->ops->destroy = DMFieldDestroy_Shell;
223: field->ops->evaluate = NULL;
224: field->ops->evaluateFE = DMFieldShellEvaluateFEDefault;
225: field->ops->evaluateFV = DMFieldShellEvaluateFVDefault;
226: field->ops->getDegree = NULL;
227: field->ops->createDefaultQuadrature = NULL;
228: field->ops->view = NULL;
229: PetscFunctionReturn(PETSC_SUCCESS);
230: }
232: PETSC_INTERN PetscErrorCode DMFieldCreate_Shell(DMField field)
233: {
234: DMField_Shell *shell;
236: PetscFunctionBegin;
237: PetscCall(PetscNew(&shell));
238: field->data = shell;
239: PetscCall(DMFieldInitialize_Shell(field));
240: PetscFunctionReturn(PETSC_SUCCESS);
241: }
243: PetscErrorCode DMFieldCreateShell(DM dm, PetscInt numComponents, DMFieldContinuity continuity, PetscCtx ctx, DMField *field)
244: {
245: DMField b;
246: DMField_Shell *shell;
248: PetscFunctionBegin;
250: if (ctx) PetscAssertPointer(ctx, 4);
251: PetscAssertPointer(field, 5);
252: PetscCall(DMFieldCreate(dm, numComponents, continuity, &b));
253: PetscCall(DMFieldSetType(b, DMFIELDSHELL));
254: shell = (DMField_Shell *)b->data;
255: shell->ctx = ctx;
256: *field = b;
257: PetscFunctionReturn(PETSC_SUCCESS);
258: }