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*/