1 #include <petsc/private/dmfieldimpl.h> /*I "petscdmfield.h" I*/
2
3 typedef struct _n_DMField_Shell {
4 PetscCtx ctx;
5 PetscErrorCode (*destroy)(DMField);
6 } DMField_Shell;
7
DMFieldShellGetContext(DMField field,PetscCtxRt ctx)8 PetscErrorCode DMFieldShellGetContext(DMField field, PetscCtxRt ctx)
9 {
10 PetscBool flg;
11
12 PetscFunctionBegin;
13 PetscValidHeaderSpecific(field, DMFIELD_CLASSID, 1);
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 }
20
DMFieldDestroy_Shell(DMField field)21 static PetscErrorCode DMFieldDestroy_Shell(DMField field)
22 {
23 DMField_Shell *shell = (DMField_Shell *)field->data;
24
25 PetscFunctionBegin;
26 if (shell->destroy) PetscCall((*shell->destroy)(field));
27 PetscCall(PetscFree(field->data));
28 PetscFunctionReturn(PETSC_SUCCESS);
29 }
30
DMFieldShellEvaluateFEDefault(DMField field,IS pointIS,PetscQuadrature quad,PetscDataType type,void * B,void * D,void * H)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;
39
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;
55 PetscInt q;
56
57 for (p = 0; p < numPoints * Nq; p++) {
58 for (q = 0; q < Nc; q++) {
59 PetscScalar d[3];
60
61 PetscInt i, j;
62
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;
73
74 for (p = 0; p < numPoints * Nq; p++) {
75 for (q = 0; q < Nc; q++) {
76 PetscReal d[3];
77
78 PetscInt i, j;
79
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;
93
94 for (p = 0; p < numPoints * Nq; p++) {
95 for (q = 0; q < Nc; q++) {
96 PetscScalar d[3][3];
97
98 PetscInt i, j, k, l;
99
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;
116
117 for (p = 0; p < numPoints * Nq; p++) {
118 for (q = 0; q < Nc; q++) {
119 PetscReal d[3][3];
120
121 PetscInt i, j, k, l;
122
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 }
143
DMFieldShellEvaluateFVDefault(DMField field,IS pointIS,PetscDataType type,void * B,void * D,void * H)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;
154
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_FEGEOM_BASIC, &geom));
165 PetscCall(ISGetLocalSize(pointIS, &numPoints));
166 PetscCall(PetscMalloc1(dimC * numPoints, &pfArray));
167 for (p = 0; p < numPoints * dimC; p++) pfArray[p] = 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 }
176
DMFieldShellSetDestroy(DMField field,PetscErrorCode (* destroy)(DMField))177 PetscErrorCode DMFieldShellSetDestroy(DMField field, PetscErrorCode (*destroy)(DMField))
178 {
179 DMField_Shell *shell = (DMField_Shell *)field->data;
180
181 PetscFunctionBegin;
182 PetscValidHeaderSpecific(field, DMFIELD_CLASSID, 1);
183 shell->destroy = destroy;
184 PetscFunctionReturn(PETSC_SUCCESS);
185 }
186
DMFieldShellSetEvaluate(DMField field,PetscErrorCode (* evaluate)(DMField,Vec,PetscDataType,void *,void *,void *))187 PetscErrorCode DMFieldShellSetEvaluate(DMField field, PetscErrorCode (*evaluate)(DMField, Vec, PetscDataType, void *, void *, void *))
188 {
189 PetscFunctionBegin;
190 PetscValidHeaderSpecific(field, DMFIELD_CLASSID, 1);
191 field->ops->evaluate = evaluate;
192 PetscFunctionReturn(PETSC_SUCCESS);
193 }
194
DMFieldShellSetEvaluateFE(DMField field,PetscErrorCode (* evaluateFE)(DMField,IS,PetscQuadrature,PetscDataType,void *,void *,void *))195 PetscErrorCode DMFieldShellSetEvaluateFE(DMField field, PetscErrorCode (*evaluateFE)(DMField, IS, PetscQuadrature, PetscDataType, void *, void *, void *))
196 {
197 PetscFunctionBegin;
198 PetscValidHeaderSpecific(field, DMFIELD_CLASSID, 1);
199 field->ops->evaluateFE = evaluateFE;
200 PetscFunctionReturn(PETSC_SUCCESS);
201 }
202
DMFieldShellSetEvaluateFV(DMField field,PetscErrorCode (* evaluateFV)(DMField,IS,PetscDataType,void *,void *,void *))203 PetscErrorCode DMFieldShellSetEvaluateFV(DMField field, PetscErrorCode (*evaluateFV)(DMField, IS, PetscDataType, void *, void *, void *))
204 {
205 PetscFunctionBegin;
206 PetscValidHeaderSpecific(field, DMFIELD_CLASSID, 1);
207 field->ops->evaluateFV = evaluateFV;
208 PetscFunctionReturn(PETSC_SUCCESS);
209 }
210
DMFieldShellSetGetDegree(DMField field,PetscErrorCode (* getDegree)(DMField,IS,PetscInt *,PetscInt *))211 PetscErrorCode DMFieldShellSetGetDegree(DMField field, PetscErrorCode (*getDegree)(DMField, IS, PetscInt *, PetscInt *))
212 {
213 PetscFunctionBegin;
214 PetscValidHeaderSpecific(field, DMFIELD_CLASSID, 1);
215 field->ops->getDegree = getDegree;
216 PetscFunctionReturn(PETSC_SUCCESS);
217 }
218
DMFieldShellSetCreateDefaultQuadrature(DMField field,PetscErrorCode (* createDefaultQuadrature)(DMField,IS,PetscQuadrature *))219 PetscErrorCode DMFieldShellSetCreateDefaultQuadrature(DMField field, PetscErrorCode (*createDefaultQuadrature)(DMField, IS, PetscQuadrature *))
220 {
221 PetscFunctionBegin;
222 PetscValidHeaderSpecific(field, DMFIELD_CLASSID, 1);
223 field->ops->createDefaultQuadrature = createDefaultQuadrature;
224 PetscFunctionReturn(PETSC_SUCCESS);
225 }
226
DMFieldInitialize_Shell(DMField field)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 }
239
DMFieldCreate_Shell(DMField field)240 PETSC_INTERN PetscErrorCode DMFieldCreate_Shell(DMField field)
241 {
242 DMField_Shell *shell;
243
244 PetscFunctionBegin;
245 PetscCall(PetscNew(&shell));
246 field->data = shell;
247 PetscCall(DMFieldInitialize_Shell(field));
248 PetscFunctionReturn(PETSC_SUCCESS);
249 }
250
DMFieldCreateShell(DM dm,PetscInt numComponents,DMFieldContinuity continuity,PetscCtx ctx,DMField * field)251 PetscErrorCode DMFieldCreateShell(DM dm, PetscInt numComponents, DMFieldContinuity continuity, PetscCtx ctx, DMField *field)
252 {
253 DMField b;
254 DMField_Shell *shell;
255
256 PetscFunctionBegin;
257 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
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 }
267