1 static const char help[] = "Tests for injecting basis functions";
2
3 #include <petscdmplex.h>
4 #include <petscfe.h>
5 #include <petscds.h>
6
7 typedef struct {
8 PetscInt its; /* Number of replications for timing */
9 } AppCtx;
10
ProcessOptions(MPI_Comm comm,AppCtx * options)11 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
12 {
13 PetscFunctionBeginUser;
14 options->its = 1;
15
16 PetscOptionsBegin(comm, "", "FE Injection Options", "PETSCFE");
17 PetscCall(PetscOptionsInt("-its", "The number of replications for timing", "ex1.c", options->its, &options->its, NULL));
18 PetscOptionsEnd();
19 PetscFunctionReturn(PETSC_SUCCESS);
20 }
21
trig_u(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)22 static PetscErrorCode trig_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
23 {
24 PetscInt d;
25 *u = 0.0;
26 for (d = 0; d < dim; ++d) *u += PetscSinReal(2.0 * PETSC_PI * x[d]);
27 return PETSC_SUCCESS;
28 }
29
f0_trig_u(PetscInt dim,PetscInt Nf,PetscInt NfAux,const PetscInt uOff[],const PetscInt uOff_x[],const PetscScalar u[],const PetscScalar u_t[],const PetscScalar u_x[],const PetscInt aOff[],const PetscInt aOff_x[],const PetscScalar a[],const PetscScalar a_t[],const PetscScalar a_x[],PetscReal t,const PetscReal x[],PetscInt numConstants,const PetscScalar constants[],PetscScalar f0[])30 static void f0_trig_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
31 {
32 PetscInt d;
33 for (d = 0; d < dim; ++d) f0[0] += -4.0 * PetscSqr(PETSC_PI) * PetscSinReal(2.0 * PETSC_PI * x[d]);
34 }
35
f1_u(PetscInt dim,PetscInt Nf,PetscInt NfAux,const PetscInt uOff[],const PetscInt uOff_x[],const PetscScalar u[],const PetscScalar u_t[],const PetscScalar u_x[],const PetscInt aOff[],const PetscInt aOff_x[],const PetscScalar a[],const PetscScalar a_t[],const PetscScalar a_x[],PetscReal t,const PetscReal x[],PetscInt numConstants,const PetscScalar constants[],PetscScalar f1[])36 static void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
37 {
38 PetscInt d;
39 for (d = 0; d < dim; ++d) f1[d] = u_x[d];
40 }
41
g3_uu(PetscInt dim,PetscInt Nf,PetscInt NfAux,const PetscInt uOff[],const PetscInt uOff_x[],const PetscScalar u[],const PetscScalar u_t[],const PetscScalar u_x[],const PetscInt aOff[],const PetscInt aOff_x[],const PetscScalar a[],const PetscScalar a_t[],const PetscScalar a_x[],PetscReal t,PetscReal u_tShift,const PetscReal x[],PetscInt numConstants,const PetscScalar constants[],PetscScalar g3[])42 static void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
43 {
44 PetscInt d;
45 for (d = 0; d < dim; ++d) g3[d * dim + d] = 1.0;
46 }
47
SetupPrimalProblem(DM dm,AppCtx * user)48 static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user)
49 {
50 PetscDS ds;
51 DMLabel label;
52 const PetscInt id = 1;
53
54 PetscFunctionBeginUser;
55 PetscCall(DMGetDS(dm, &ds));
56 PetscCall(PetscDSSetResidual(ds, 0, f0_trig_u, f1_u));
57 PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu));
58 PetscCall(PetscDSSetExactSolution(ds, 0, trig_u, user));
59 PetscCall(DMGetLabel(dm, "marker", &label));
60 if (label) PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)trig_u, NULL, user, NULL));
61 PetscFunctionReturn(PETSC_SUCCESS);
62 }
63
SetupDiscretization(DM dm,const char name[],PetscErrorCode (* setup)(DM,AppCtx *),AppCtx * user)64 static PetscErrorCode SetupDiscretization(DM dm, const char name[], PetscErrorCode (*setup)(DM, AppCtx *), AppCtx *user)
65 {
66 DM cdm = dm;
67 PetscFE fe;
68 char prefix[PETSC_MAX_PATH_LEN];
69 PetscInt dim;
70
71 PetscFunctionBeginUser;
72 PetscCall(DMGetDimension(dm, &dim));
73 PetscCall(PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name));
74 PetscCall(DMCreateFEDefault(dm, dim, name ? prefix : NULL, -1, &fe));
75 PetscCall(PetscObjectSetName((PetscObject)fe, name));
76 /* Set discretization and boundary conditions for each mesh */
77 PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
78 PetscCall(DMCreateDS(dm));
79 PetscCall((*setup)(dm, user));
80 while (cdm) {
81 PetscCall(DMCopyDisc(dm, cdm));
82 PetscCall(DMGetCoarseDM(cdm, &cdm));
83 }
84 PetscCall(PetscFEDestroy(&fe));
85 PetscFunctionReturn(PETSC_SUCCESS);
86 }
87
88 /* PetscObjectContainerCompose() compose requires void ** signature on destructor */
PetscFEGeomDestroy_Void(PetscCtxRt ctx)89 static PetscErrorCode PetscFEGeomDestroy_Void(PetscCtxRt ctx)
90 {
91 return PetscFEGeomDestroy((PetscFEGeom **)ctx);
92 }
93
CellRangeGetFEGeom(IS cellIS,DMField coordField,PetscQuadrature quad,PetscFEGeomMode mode,PetscFEGeom ** geom)94 PetscErrorCode CellRangeGetFEGeom(IS cellIS, DMField coordField, PetscQuadrature quad, PetscFEGeomMode mode, PetscFEGeom **geom)
95 {
96 char composeStr[33] = {0};
97 PetscObjectId id;
98 PetscContainer container;
99
100 PetscFunctionBegin;
101 PetscCall(PetscObjectGetId((PetscObject)quad, &id));
102 PetscCall(PetscSNPrintf(composeStr, 32, "CellRangeGetFEGeom_%" PetscInt64_FMT "\n", id));
103 PetscCall(PetscObjectQuery((PetscObject)cellIS, composeStr, (PetscObject *)&container));
104 if (container) {
105 PetscCall(PetscContainerGetPointer(container, geom));
106 } else {
107 PetscCall(DMFieldCreateFEGeom(coordField, cellIS, quad, mode, geom));
108 PetscCall(PetscObjectContainerCompose((PetscObject)cellIS, composeStr, *geom, PetscFEGeomDestroy_Void));
109 }
110 PetscFunctionReturn(PETSC_SUCCESS);
111 }
112
CellRangeRestoreFEGeom(IS cellIS,DMField coordField,PetscQuadrature quad,PetscBool faceData,PetscFEGeom ** geom)113 PetscErrorCode CellRangeRestoreFEGeom(IS cellIS, DMField coordField, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
114 {
115 PetscFunctionBegin;
116 *geom = NULL;
117 PetscFunctionReturn(PETSC_SUCCESS);
118 }
119
CreateFEGeometry(DM dm,PetscDS ds,IS cellIS,PetscQuadrature * affineQuad,PetscFEGeom ** affineGeom,PetscQuadrature ** quads,PetscFEGeom *** geoms)120 static PetscErrorCode CreateFEGeometry(DM dm, PetscDS ds, IS cellIS, PetscQuadrature *affineQuad, PetscFEGeom **affineGeom, PetscQuadrature **quads, PetscFEGeom ***geoms)
121 {
122 DMField coordField;
123 PetscInt Nf, f, maxDegree;
124
125 PetscFunctionBeginUser;
126 *affineQuad = NULL;
127 *affineGeom = NULL;
128 *quads = NULL;
129 *geoms = NULL;
130 PetscCall(PetscDSGetNumFields(ds, &Nf));
131 PetscCall(DMGetCoordinateField(dm, &coordField));
132 PetscCall(DMFieldGetDegree(coordField, cellIS, NULL, &maxDegree));
133 if (maxDegree <= 1) {
134 PetscCall(DMFieldCreateDefaultQuadrature(coordField, cellIS, affineQuad));
135 if (*affineQuad) PetscCall(CellRangeGetFEGeom(cellIS, coordField, *affineQuad, PETSC_FEGEOM_BASIC, affineGeom));
136 } else {
137 PetscCall(PetscCalloc2(Nf, quads, Nf, geoms));
138 for (f = 0; f < Nf; ++f) {
139 PetscFE fe;
140
141 PetscCall(PetscDSGetDiscretization(ds, f, (PetscObject *)&fe));
142 PetscCall(PetscFEGetQuadrature(fe, &(*quads)[f]));
143 PetscCall(PetscObjectReference((PetscObject)(*quads)[f]));
144 PetscCall(CellRangeGetFEGeom(cellIS, coordField, (*quads)[f], PETSC_FEGEOM_BASIC, &(*geoms)[f]));
145 }
146 }
147 PetscFunctionReturn(PETSC_SUCCESS);
148 }
149
DestroyFEGeometry(DM dm,PetscDS ds,IS cellIS,PetscQuadrature * affineQuad,PetscFEGeom ** affineGeom,PetscQuadrature ** quads,PetscFEGeom *** geoms)150 static PetscErrorCode DestroyFEGeometry(DM dm, PetscDS ds, IS cellIS, PetscQuadrature *affineQuad, PetscFEGeom **affineGeom, PetscQuadrature **quads, PetscFEGeom ***geoms)
151 {
152 DMField coordField;
153 PetscInt Nf, f;
154
155 PetscFunctionBeginUser;
156 PetscCall(PetscDSGetNumFields(ds, &Nf));
157 PetscCall(DMGetCoordinateField(dm, &coordField));
158 if (*affineQuad) {
159 PetscCall(CellRangeRestoreFEGeom(cellIS, coordField, *affineQuad, PETSC_FALSE, affineGeom));
160 PetscCall(PetscQuadratureDestroy(affineQuad));
161 } else {
162 for (f = 0; f < Nf; ++f) {
163 PetscCall(CellRangeRestoreFEGeom(cellIS, coordField, (*quads)[f], PETSC_FALSE, &(*geoms)[f]));
164 PetscCall(PetscQuadratureDestroy(&(*quads)[f]));
165 }
166 PetscCall(PetscFree2(*quads, *geoms));
167 }
168 PetscFunctionReturn(PETSC_SUCCESS);
169 }
170
TestEvaluation(DM dm)171 static PetscErrorCode TestEvaluation(DM dm)
172 {
173 PetscFE fe;
174 PetscSpace sp;
175 PetscReal *points;
176 PetscReal *B, *D, *H;
177 PetscInt dim, Nb, b, Nc, c, Np, p;
178
179 PetscFunctionBeginUser;
180 PetscCall(DMGetDimension(dm, &dim));
181 PetscCall(DMGetField(dm, 0, NULL, (PetscObject *)&fe));
182 Np = 6;
183 PetscCall(PetscMalloc1(Np * dim, &points));
184 if (dim == 3) {
185 points[0] = -1.0;
186 points[1] = -1.0;
187 points[2] = -1.0;
188 points[3] = 1.0;
189 points[4] = -1.0;
190 points[5] = -1.0;
191 points[6] = -1.0;
192 points[7] = 1.0;
193 points[8] = -1.0;
194 points[9] = -1.0;
195 points[10] = -1.0;
196 points[11] = 1.0;
197 points[12] = 1.0;
198 points[13] = -1.0;
199 points[14] = 1.0;
200 points[15] = -1.0;
201 points[16] = 1.0;
202 points[17] = 1.0;
203 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only for 3D right now");
204 PetscCall(PetscFEGetBasisSpace(fe, &sp));
205 PetscCall(PetscSpaceGetDimension(sp, &Nb));
206 PetscCall(PetscSpaceGetNumComponents(sp, &Nc));
207 PetscCall(DMGetWorkArray(dm, Np * Nb * Nc, MPIU_REAL, &B));
208 PetscCall(DMGetWorkArray(dm, Np * Nb * Nc * dim, MPIU_REAL, &D));
209 PetscCall(DMGetWorkArray(dm, Np * Nb * Nc * dim * dim, MPIU_REAL, &H));
210 PetscCall(PetscSpaceEvaluate(sp, Np, points, B, NULL, NULL /*D, H*/));
211 for (p = 0; p < Np; ++p) {
212 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Point %" PetscInt_FMT "\n", p));
213 for (b = 0; b < Nb; ++b) {
214 PetscCall(PetscPrintf(PETSC_COMM_SELF, "B[%" PetscInt_FMT "]:", b));
215 for (c = 0; c < Nc; ++c) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %g", (double)B[(p * Nb + b) * Nc + c]));
216 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
217 #if 0
218 for (c = 0; c < Nc; ++c) {
219 PetscCall(PetscPrintf(PETSC_COMM_SELF, " D[%" PetscInt_FMT ",%" PetscInt_FMT "]:", b, c));
220 for (d = 0; d < dim; ++d) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %g", (double) B[((p*Nb+b)*Nc+c)*dim+d)]));
221 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
222 }
223 #endif
224 }
225 }
226 PetscCall(DMRestoreWorkArray(dm, Np * Nb, MPIU_REAL, &B));
227 PetscCall(DMRestoreWorkArray(dm, Np * Nb * dim, MPIU_REAL, &D));
228 PetscCall(DMRestoreWorkArray(dm, Np * Nb * dim * dim, MPIU_REAL, &H));
229 PetscCall(PetscFree(points));
230 PetscFunctionReturn(PETSC_SUCCESS);
231 }
232
TestIntegration(DM dm,PetscInt cbs,PetscInt its)233 static PetscErrorCode TestIntegration(DM dm, PetscInt cbs, PetscInt its)
234 {
235 PetscDS ds;
236 PetscFEGeom *chunkGeom = NULL;
237 PetscQuadrature affineQuad, *quads = NULL;
238 PetscFEGeom *affineGeom, **geoms = NULL;
239 PetscScalar *u, *elemVec;
240 IS cellIS;
241 PetscInt depth, cStart, cEnd, cell, chunkSize = cbs, Nch = 0, Nf, f, totDim, i, k;
242
243 PetscFunctionBeginUser;
244 PetscCall(DMPlexGetDepth(dm, &depth));
245 PetscCall(DMGetStratumIS(dm, "depth", depth, &cellIS));
246 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
247 PetscCall(DMGetCellDS(dm, cStart, &ds, NULL));
248 PetscCall(PetscDSGetNumFields(ds, &Nf));
249 PetscCall(PetscDSGetTotalDimension(ds, &totDim));
250 PetscCall(CreateFEGeometry(dm, ds, cellIS, &affineQuad, &affineGeom, &quads, &geoms));
251 PetscCall(PetscMalloc2(chunkSize * totDim, &u, chunkSize * totDim, &elemVec));
252 /* Assumptions:
253 - Single field
254 - No input data
255 - No auxiliary data
256 - No time-dependence
257 */
258 for (i = 0; i < its; ++i) {
259 for (cell = cStart; cell < cEnd; cell += chunkSize, ++Nch) {
260 const PetscInt cS = cell, cE = PetscMin(cS + chunkSize, cEnd), Ne = cE - cS;
261
262 PetscCall(PetscArrayzero(elemVec, chunkSize * totDim));
263 /* TODO Replace with DMPlexGetCellFields() */
264 for (k = 0; k < chunkSize * totDim; ++k) u[k] = 1.0;
265 for (f = 0; f < Nf; ++f) {
266 PetscFormKey key;
267 PetscFEGeom *geom = affineGeom ? affineGeom : geoms[f];
268 /* PetscQuadrature quad = affineQuad ? affineQuad : quads[f]; */
269
270 key.label = NULL;
271 key.value = 0;
272 key.field = f;
273 key.part = 0;
274 PetscCall(PetscFEGeomGetChunk(geom, cS, cE, &chunkGeom));
275 PetscCall(PetscFEIntegrateResidual(ds, key, Ne, chunkGeom, u, NULL, NULL, NULL, 0.0, elemVec));
276 }
277 }
278 }
279 PetscCall(PetscFEGeomRestoreChunk(affineGeom, cStart, cEnd, &chunkGeom));
280 PetscCall(DestroyFEGeometry(dm, ds, cellIS, &affineQuad, &affineGeom, &quads, &geoms));
281 PetscCall(ISDestroy(&cellIS));
282 PetscCall(PetscFree2(u, elemVec));
283 PetscFunctionReturn(PETSC_SUCCESS);
284 }
285
TestUnisolvence(DM dm)286 static PetscErrorCode TestUnisolvence(DM dm)
287 {
288 Mat M;
289 Vec v;
290
291 PetscFunctionBeginUser;
292 PetscCall(DMGetLocalVector(dm, &v));
293 PetscCall(DMRestoreLocalVector(dm, &v));
294 PetscCall(DMCreateMassMatrix(dm, dm, &M));
295 PetscCall(MatViewFromOptions(M, NULL, "-mass_view"));
296 PetscCall(MatDestroy(&M));
297 PetscFunctionReturn(PETSC_SUCCESS);
298 }
299
main(int argc,char ** argv)300 int main(int argc, char **argv)
301 {
302 DM dm;
303 AppCtx ctx;
304 PetscMPIInt size;
305
306 PetscFunctionBeginUser;
307 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
308 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
309 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only.");
310 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &ctx));
311 PetscCall(DMCreate(PETSC_COMM_WORLD, &dm));
312 PetscCall(DMSetType(dm, DMPLEX));
313 PetscCall(DMSetFromOptions(dm));
314 PetscCall(PetscObjectSetName((PetscObject)dm, "Mesh"));
315 PetscCall(PetscObjectViewFromOptions((PetscObject)dm, NULL, "-dm_view"));
316 PetscCall(SetupDiscretization(dm, "field", SetupPrimalProblem, &ctx));
317 PetscCall(TestEvaluation(dm));
318 PetscCall(TestIntegration(dm, 1, ctx.its));
319 PetscCall(TestUnisolvence(dm));
320 PetscCall(DMDestroy(&dm));
321 PetscCall(PetscFinalize());
322 return 0;
323 }
324
325 /*TEST
326 test:
327 suffix: 0
328 args: -dm_plex_reference_cell_domain -dm_plex_cell triangular_prism -field_petscspace_degree 0
329
330 test:
331 suffix: 2
332 args: -dm_plex_reference_cell_domain -dm_plex_cell triangular_prism \
333 -field_petscspace_type sum \
334 -field_petscspace_variables 3 \
335 -field_petscspace_components 3 \
336 -field_petscspace_sum_spaces 2 \
337 -field_petscspace_sum_concatenate false \
338 -field_sumcomp_0_petscspace_variables 3 \
339 -field_sumcomp_0_petscspace_components 3 \
340 -field_sumcomp_0_petscspace_degree 1 \
341 -field_sumcomp_1_petscspace_variables 3 \
342 -field_sumcomp_1_petscspace_components 3 \
343 -field_sumcomp_1_petscspace_type wxy \
344 -field_petscdualspace_form_degree 0 \
345 -field_petscdualspace_order 1 \
346 -field_petscdualspace_components 3
347
348 TEST*/
349