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