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