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, 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 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: 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;
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, (void (*)(void))trig_u, NULL, user, NULL));
61: PetscFunctionReturn(PETSC_SUCCESS);
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;
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: }
88: static PetscErrorCode PetscContainerUserDestroy_PetscFEGeom(void *ctx)
89: {
90: PetscFEGeom *geom = (PetscFEGeom *)ctx;
92: PetscFunctionBegin;
93: PetscCall(PetscFEGeomDestroy(&geom));
94: PetscFunctionReturn(PETSC_SUCCESS);
95: }
97: PetscErrorCode CellRangeGetFEGeom(IS cellIS, DMField coordField, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
98: {
99: char composeStr[33] = {0};
100: PetscObjectId id;
101: PetscContainer container;
103: PetscFunctionBegin;
104: PetscCall(PetscObjectGetId((PetscObject)quad, &id));
105: PetscCall(PetscSNPrintf(composeStr, 32, "CellRangeGetFEGeom_%" PetscInt64_FMT "\n", id));
106: PetscCall(PetscObjectQuery((PetscObject)cellIS, composeStr, (PetscObject *)&container));
107: if (container) {
108: PetscCall(PetscContainerGetPointer(container, (void **)geom));
109: } else {
110: PetscCall(DMFieldCreateFEGeom(coordField, cellIS, quad, faceData, geom));
111: PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &container));
112: PetscCall(PetscContainerSetPointer(container, (void *)*geom));
113: PetscCall(PetscContainerSetUserDestroy(container, PetscContainerUserDestroy_PetscFEGeom));
114: PetscCall(PetscObjectCompose((PetscObject)cellIS, composeStr, (PetscObject)container));
115: PetscCall(PetscContainerDestroy(&container));
116: }
117: PetscFunctionReturn(PETSC_SUCCESS);
118: }
120: PetscErrorCode CellRangeRestoreFEGeom(IS cellIS, DMField coordField, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
121: {
122: PetscFunctionBegin;
123: *geom = NULL;
124: PetscFunctionReturn(PETSC_SUCCESS);
125: }
127: static PetscErrorCode CreateFEGeometry(DM dm, PetscDS ds, IS cellIS, PetscQuadrature *affineQuad, PetscFEGeom **affineGeom, PetscQuadrature **quads, PetscFEGeom ***geoms)
128: {
129: DMField coordField;
130: PetscInt Nf, f, maxDegree;
132: PetscFunctionBeginUser;
133: *affineQuad = NULL;
134: *affineGeom = NULL;
135: *quads = NULL;
136: *geoms = NULL;
137: PetscCall(PetscDSGetNumFields(ds, &Nf));
138: PetscCall(DMGetCoordinateField(dm, &coordField));
139: PetscCall(DMFieldGetDegree(coordField, cellIS, NULL, &maxDegree));
140: if (maxDegree <= 1) {
141: PetscCall(DMFieldCreateDefaultQuadrature(coordField, cellIS, affineQuad));
142: if (*affineQuad) PetscCall(CellRangeGetFEGeom(cellIS, coordField, *affineQuad, PETSC_FALSE, affineGeom));
143: } else {
144: PetscCall(PetscCalloc2(Nf, quads, Nf, geoms));
145: for (f = 0; f < Nf; ++f) {
146: PetscFE fe;
148: PetscCall(PetscDSGetDiscretization(ds, f, (PetscObject *)&fe));
149: PetscCall(PetscFEGetQuadrature(fe, &(*quads)[f]));
150: PetscCall(PetscObjectReference((PetscObject)(*quads)[f]));
151: PetscCall(CellRangeGetFEGeom(cellIS, coordField, (*quads)[f], PETSC_FALSE, &(*geoms)[f]));
152: }
153: }
154: PetscFunctionReturn(PETSC_SUCCESS);
155: }
157: static PetscErrorCode DestroyFEGeometry(DM dm, PetscDS ds, IS cellIS, PetscQuadrature *affineQuad, PetscFEGeom **affineGeom, PetscQuadrature **quads, PetscFEGeom ***geoms)
158: {
159: DMField coordField;
160: PetscInt Nf, f;
162: PetscFunctionBeginUser;
163: PetscCall(PetscDSGetNumFields(ds, &Nf));
164: PetscCall(DMGetCoordinateField(dm, &coordField));
165: if (*affineQuad) {
166: PetscCall(CellRangeRestoreFEGeom(cellIS, coordField, *affineQuad, PETSC_FALSE, affineGeom));
167: PetscCall(PetscQuadratureDestroy(affineQuad));
168: } else {
169: for (f = 0; f < Nf; ++f) {
170: PetscCall(CellRangeRestoreFEGeom(cellIS, coordField, (*quads)[f], PETSC_FALSE, &(*geoms)[f]));
171: PetscCall(PetscQuadratureDestroy(&(*quads)[f]));
172: }
173: PetscCall(PetscFree2(*quads, *geoms));
174: }
175: PetscFunctionReturn(PETSC_SUCCESS);
176: }
178: static PetscErrorCode TestEvaluation(DM dm)
179: {
180: PetscFE fe;
181: PetscSpace sp;
182: PetscReal *points;
183: PetscReal *B, *D, *H;
184: PetscInt dim, Nb, b, Nc, c, Np, p;
186: PetscFunctionBeginUser;
187: PetscCall(DMGetDimension(dm, &dim));
188: PetscCall(DMGetField(dm, 0, NULL, (PetscObject *)&fe));
189: Np = 6;
190: PetscCall(PetscMalloc1(Np * dim, &points));
191: if (dim == 3) {
192: points[0] = -1.0;
193: points[1] = -1.0;
194: points[2] = -1.0;
195: points[3] = 1.0;
196: points[4] = -1.0;
197: points[5] = -1.0;
198: points[6] = -1.0;
199: points[7] = 1.0;
200: points[8] = -1.0;
201: points[9] = -1.0;
202: points[10] = -1.0;
203: points[11] = 1.0;
204: points[12] = 1.0;
205: points[13] = -1.0;
206: points[14] = 1.0;
207: points[15] = -1.0;
208: points[16] = 1.0;
209: points[17] = 1.0;
210: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only for 3D right now");
211: PetscCall(PetscFEGetBasisSpace(fe, &sp));
212: PetscCall(PetscSpaceGetDimension(sp, &Nb));
213: PetscCall(PetscSpaceGetNumComponents(sp, &Nc));
214: PetscCall(DMGetWorkArray(dm, Np * Nb * Nc, MPIU_REAL, &B));
215: PetscCall(DMGetWorkArray(dm, Np * Nb * Nc * dim, MPIU_REAL, &D));
216: PetscCall(DMGetWorkArray(dm, Np * Nb * Nc * dim * dim, MPIU_REAL, &H));
217: PetscCall(PetscSpaceEvaluate(sp, Np, points, B, NULL, NULL /*D, H*/));
218: for (p = 0; p < Np; ++p) {
219: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Point %" PetscInt_FMT "\n", p));
220: for (b = 0; b < Nb; ++b) {
221: PetscCall(PetscPrintf(PETSC_COMM_SELF, "B[%" PetscInt_FMT "]:", b));
222: for (c = 0; c < Nc; ++c) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %g", (double)B[(p * Nb + b) * Nc + c]));
223: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
224: #if 0
225: for (c = 0; c < Nc; ++c) {
226: PetscCall(PetscPrintf(PETSC_COMM_SELF, " D[%" PetscInt_FMT ",%" PetscInt_FMT "]:", b, c));
227: for (d = 0; d < dim; ++d) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %g", (double) B[((p*Nb+b)*Nc+c)*dim+d)]));
228: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
229: }
230: #endif
231: }
232: }
233: PetscCall(DMRestoreWorkArray(dm, Np * Nb, MPIU_REAL, &B));
234: PetscCall(DMRestoreWorkArray(dm, Np * Nb * dim, MPIU_REAL, &D));
235: PetscCall(DMRestoreWorkArray(dm, Np * Nb * dim * dim, MPIU_REAL, &H));
236: PetscCall(PetscFree(points));
237: PetscFunctionReturn(PETSC_SUCCESS);
238: }
240: static PetscErrorCode TestIntegration(DM dm, PetscInt cbs, PetscInt its)
241: {
242: PetscDS ds;
243: PetscFEGeom *chunkGeom = NULL;
244: PetscQuadrature affineQuad, *quads = NULL;
245: PetscFEGeom *affineGeom, **geoms = NULL;
246: PetscScalar *u, *elemVec;
247: IS cellIS;
248: PetscInt depth, cStart, cEnd, cell, chunkSize = cbs, Nch = 0, Nf, f, totDim, i, k;
250: PetscFunctionBeginUser;
251: PetscCall(DMPlexGetDepth(dm, &depth));
252: PetscCall(DMGetStratumIS(dm, "depth", depth, &cellIS));
253: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
254: PetscCall(DMGetCellDS(dm, cStart, &ds, NULL));
255: PetscCall(PetscDSGetNumFields(ds, &Nf));
256: PetscCall(PetscDSGetTotalDimension(ds, &totDim));
257: PetscCall(CreateFEGeometry(dm, ds, cellIS, &affineQuad, &affineGeom, &quads, &geoms));
258: PetscCall(PetscMalloc2(chunkSize * totDim, &u, chunkSize * totDim, &elemVec));
259: /* Assumptions:
260: - Single field
261: - No input data
262: - No auxiliary data
263: - No time-dependence
264: */
265: for (i = 0; i < its; ++i) {
266: for (cell = cStart; cell < cEnd; cell += chunkSize, ++Nch) {
267: const PetscInt cS = cell, cE = PetscMin(cS + chunkSize, cEnd), Ne = cE - cS;
269: PetscCall(PetscArrayzero(elemVec, chunkSize * totDim));
270: /* TODO Replace with DMPlexGetCellFields() */
271: for (k = 0; k < chunkSize * totDim; ++k) u[k] = 1.0;
272: for (f = 0; f < Nf; ++f) {
273: PetscFormKey key;
274: PetscFEGeom *geom = affineGeom ? affineGeom : geoms[f];
275: /* PetscQuadrature quad = affineQuad ? affineQuad : quads[f]; */
277: key.label = NULL;
278: key.value = 0;
279: key.field = f;
280: key.part = 0;
281: PetscCall(PetscFEGeomGetChunk(geom, cS, cE, &chunkGeom));
282: PetscCall(PetscFEIntegrateResidual(ds, key, Ne, chunkGeom, u, NULL, NULL, NULL, 0.0, elemVec));
283: }
284: }
285: }
286: PetscCall(PetscFEGeomRestoreChunk(affineGeom, cStart, cEnd, &chunkGeom));
287: PetscCall(DestroyFEGeometry(dm, ds, cellIS, &affineQuad, &affineGeom, &quads, &geoms));
288: PetscCall(ISDestroy(&cellIS));
289: PetscCall(PetscFree2(u, elemVec));
290: PetscFunctionReturn(PETSC_SUCCESS);
291: }
293: static PetscErrorCode TestUnisolvence(DM dm)
294: {
295: Mat M;
296: Vec v;
298: PetscFunctionBeginUser;
299: PetscCall(DMGetLocalVector(dm, &v));
300: PetscCall(DMRestoreLocalVector(dm, &v));
301: PetscCall(DMCreateMassMatrix(dm, dm, &M));
302: PetscCall(MatViewFromOptions(M, NULL, "-mass_view"));
303: PetscCall(MatDestroy(&M));
304: PetscFunctionReturn(PETSC_SUCCESS);
305: }
307: int main(int argc, char **argv)
308: {
309: DM dm;
310: AppCtx ctx;
311: PetscMPIInt size;
313: PetscFunctionBeginUser;
314: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
315: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
316: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_SUP, "This is a uniprocessor example only.");
317: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &ctx));
318: PetscCall(DMCreate(PETSC_COMM_WORLD, &dm));
319: PetscCall(DMSetType(dm, DMPLEX));
320: PetscCall(DMSetFromOptions(dm));
321: PetscCall(PetscObjectSetName((PetscObject)dm, "Mesh"));
322: PetscCall(PetscObjectViewFromOptions((PetscObject)dm, NULL, "-dm_view"));
323: PetscCall(SetupDiscretization(dm, "field", SetupPrimalProblem, &ctx));
324: PetscCall(TestEvaluation(dm));
325: PetscCall(TestIntegration(dm, 1, ctx.its));
326: PetscCall(TestUnisolvence(dm));
327: PetscCall(DMDestroy(&dm));
328: PetscCall(PetscFinalize());
329: return 0;
330: }
332: /*TEST
333: test:
334: suffix: 0
335: args: -dm_plex_reference_cell_domain -dm_plex_cell triangular_prism -field_petscspace_degree 0
337: test:
338: suffix: 2
339: args: -dm_plex_reference_cell_domain -dm_plex_cell triangular_prism \
340: -field_petscspace_type sum \
341: -field_petscspace_variables 3 \
342: -field_petscspace_components 3 \
343: -field_petscspace_sum_spaces 2 \
344: -field_petscspace_sum_concatenate false \
345: -field_sumcomp_0_petscspace_variables 3 \
346: -field_sumcomp_0_petscspace_components 3 \
347: -field_sumcomp_0_petscspace_degree 1 \
348: -field_sumcomp_1_petscspace_variables 3 \
349: -field_sumcomp_1_petscspace_components 3 \
350: -field_sumcomp_1_petscspace_type wxy \
351: -field_petscdualspace_form_degree 0 \
352: -field_petscdualspace_order 1 \
353: -field_petscdualspace_components 3
355: TEST*/