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