Actual source code: ex1.c

  1: static char help[] = "Demonstration of creating and viewing DMFields objects.\n\n";

  3: #include <petscdmfield.h>
  4: #include <petscdmplex.h>
  5: #include <petscdmda.h>

  7: static PetscErrorCode ViewResults(PetscViewer viewer, PetscInt N, PetscInt dim, PetscScalar *B, PetscScalar *D, PetscScalar *H, PetscReal *rB, PetscReal *rD, PetscReal *rH)
  8: {
  9:   PetscViewerASCIIPrintf(viewer, "B:\n");
 10:   PetscScalarView(N, B, viewer);
 11:   PetscViewerASCIIPrintf(viewer, "D:\n");
 12:   PetscScalarView(N * dim, D, viewer);
 13:   PetscViewerASCIIPrintf(viewer, "H:\n");
 14:   PetscScalarView(N * dim * dim, H, viewer);

 16:   PetscViewerASCIIPrintf(viewer, "rB:\n");
 17:   PetscRealView(N, rB, viewer);
 18:   PetscViewerASCIIPrintf(viewer, "rD:\n");
 19:   PetscRealView(N * dim, rD, viewer);
 20:   PetscViewerASCIIPrintf(viewer, "rH:\n");
 21:   PetscRealView(N * dim * dim, rH, viewer);
 22:   return 0;
 23: }

 25: static PetscErrorCode TestEvaluate(DMField field, PetscInt n, PetscRandom rand)
 26: {
 27:   DM           dm;
 28:   PetscInt     dim, i, nc;
 29:   PetscScalar *B, *D, *H;
 30:   PetscReal   *rB, *rD, *rH;
 31:   Vec          points;
 32:   PetscScalar *array;
 33:   PetscViewer  viewer;
 34:   MPI_Comm     comm;

 36:   comm = PetscObjectComm((PetscObject)field);
 37:   DMFieldGetNumComponents(field, &nc);
 38:   DMFieldGetDM(field, &dm);
 39:   DMGetDimension(dm, &dim);
 40:   VecCreateMPI(PetscObjectComm((PetscObject)field), n * dim, PETSC_DETERMINE, &points);
 41:   VecSetBlockSize(points, dim);
 42:   VecGetArray(points, &array);
 43:   for (i = 0; i < n * dim; i++) PetscRandomGetValue(rand, &array[i]);
 44:   VecRestoreArray(points, &array);
 45:   PetscMalloc6(n * nc, &B, n * nc, &rB, n * nc * dim, &D, n * nc * dim, &rD, n * nc * dim * dim, &H, n * nc * dim * dim, &rH);
 46:   DMFieldEvaluate(field, points, PETSC_SCALAR, B, D, H);
 47:   DMFieldEvaluate(field, points, PETSC_REAL, rB, rD, rH);
 48:   viewer = PETSC_VIEWER_STDOUT_(comm);

 50:   PetscObjectSetName((PetscObject)points, "Test Points");
 51:   VecView(points, viewer);
 52:   ViewResults(viewer, n * nc, dim, B, D, H, rB, rD, rH);

 54:   PetscFree6(B, rB, D, rD, H, rH);
 55:   VecDestroy(&points);
 56:   return 0;
 57: }

 59: static PetscErrorCode TestEvaluateFE(DMField field, PetscInt n, PetscInt cStart, PetscInt cEnd, PetscQuadrature quad, PetscRandom rand)
 60: {
 61:   DM           dm;
 62:   PetscInt     dim, i, nc, nq;
 63:   PetscInt     N;
 64:   PetscScalar *B, *D, *H;
 65:   PetscReal   *rB, *rD, *rH;
 66:   PetscInt    *cells;
 67:   IS           cellIS;
 68:   PetscViewer  viewer;
 69:   MPI_Comm     comm;

 71:   comm = PetscObjectComm((PetscObject)field);
 72:   DMFieldGetNumComponents(field, &nc);
 73:   DMFieldGetDM(field, &dm);
 74:   DMGetDimension(dm, &dim);
 75:   PetscRandomSetInterval(rand, (PetscScalar)cStart, (PetscScalar)cEnd);
 76:   PetscMalloc1(n, &cells);
 77:   for (i = 0; i < n; i++) {
 78:     PetscReal rc;

 80:     PetscRandomGetValueReal(rand, &rc);
 81:     cells[i] = PetscFloorReal(rc);
 82:   }
 83:   ISCreateGeneral(comm, n, cells, PETSC_OWN_POINTER, &cellIS);
 84:   PetscObjectSetName((PetscObject)cellIS, "FE Test Cells");
 85:   PetscQuadratureGetData(quad, NULL, NULL, &nq, NULL, NULL);
 86:   N = n * nq * nc;
 87:   PetscMalloc6(N, &B, N, &rB, N * dim, &D, N * dim, &rD, N * dim * dim, &H, N * dim * dim, &rH);
 88:   DMFieldEvaluateFE(field, cellIS, quad, PETSC_SCALAR, B, D, H);
 89:   DMFieldEvaluateFE(field, cellIS, quad, PETSC_REAL, rB, rD, rH);
 90:   viewer = PETSC_VIEWER_STDOUT_(comm);

 92:   PetscObjectSetName((PetscObject)quad, "Test quadrature");
 93:   PetscQuadratureView(quad, viewer);
 94:   ISView(cellIS, viewer);
 95:   ViewResults(viewer, N, dim, B, D, H, rB, rD, rH);

 97:   PetscFree6(B, rB, D, rD, H, rH);
 98:   ISDestroy(&cellIS);
 99:   return 0;
100: }

102: static PetscErrorCode TestEvaluateFV(DMField field, PetscInt n, PetscInt cStart, PetscInt cEnd, PetscRandom rand)
103: {
104:   DM           dm;
105:   PetscInt     dim, i, nc;
106:   PetscInt     N;
107:   PetscScalar *B, *D, *H;
108:   PetscReal   *rB, *rD, *rH;
109:   PetscInt    *cells;
110:   IS           cellIS;
111:   PetscViewer  viewer;
112:   MPI_Comm     comm;

114:   comm = PetscObjectComm((PetscObject)field);
115:   DMFieldGetNumComponents(field, &nc);
116:   DMFieldGetDM(field, &dm);
117:   DMGetDimension(dm, &dim);
118:   PetscRandomSetInterval(rand, (PetscScalar)cStart, (PetscScalar)cEnd);
119:   PetscMalloc1(n, &cells);
120:   for (i = 0; i < n; i++) {
121:     PetscReal rc;

123:     PetscRandomGetValueReal(rand, &rc);
124:     cells[i] = PetscFloorReal(rc);
125:   }
126:   ISCreateGeneral(comm, n, cells, PETSC_OWN_POINTER, &cellIS);
127:   PetscObjectSetName((PetscObject)cellIS, "FV Test Cells");
128:   N = n * nc;
129:   PetscMalloc6(N, &B, N, &rB, N * dim, &D, N * dim, &rD, N * dim * dim, &H, N * dim * dim, &rH);
130:   DMFieldEvaluateFV(field, cellIS, PETSC_SCALAR, B, D, H);
131:   DMFieldEvaluateFV(field, cellIS, PETSC_REAL, rB, rD, rH);
132:   viewer = PETSC_VIEWER_STDOUT_(comm);

134:   ISView(cellIS, viewer);
135:   ViewResults(viewer, N, dim, B, D, H, rB, rD, rH);

137:   PetscFree6(B, rB, D, rD, H, rH);
138:   ISDestroy(&cellIS);
139:   return 0;
140: }

142: static PetscErrorCode radiusSquared(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx)
143: {
144:   PetscInt  i;
145:   PetscReal r2 = 0.;

147:   for (i = 0; i < dim; i++) r2 += PetscSqr(x[i]);
148:   for (i = 0; i < Nf; i++) u[i] = (i + 1) * r2;
149:   return 0;
150: }

152: static PetscErrorCode TestShellEvaluate(DMField field, Vec points, PetscDataType type, void *B, void *D, void *H)
153: {
154:   Vec                ctxVec = NULL;
155:   const PetscScalar *mult;
156:   PetscInt           dim;
157:   const PetscScalar *x;
158:   PetscInt           Nc, n, i, j, k, l;

160:   DMFieldGetNumComponents(field, &Nc);
161:   DMFieldShellGetContext(field, &ctxVec);
162:   VecGetBlockSize(points, &dim);
163:   VecGetLocalSize(points, &n);
164:   n /= Nc;
165:   VecGetArrayRead(ctxVec, &mult);
166:   VecGetArrayRead(points, &x);
167:   for (i = 0; i < n; i++) {
168:     PetscReal r2 = 0.;

170:     for (j = 0; j < dim; j++) r2 += PetscSqr(PetscRealPart(x[i * dim + j]));
171:     for (j = 0; j < Nc; j++) {
172:       PetscReal m = PetscRealPart(mult[j]);
173:       if (B) {
174:         if (type == PETSC_SCALAR) {
175:           ((PetscScalar *)B)[i * Nc + j] = m * r2;
176:         } else {
177:           ((PetscReal *)B)[i * Nc + j] = m * r2;
178:         }
179:       }
180:       if (D) {
181:         if (type == PETSC_SCALAR) {
182:           for (k = 0; k < dim; k++) ((PetscScalar *)D)[(i * Nc + j) * dim + k] = 2. * m * x[i * dim + k];
183:         } else {
184:           for (k = 0; k < dim; k++) ((PetscReal *)D)[(i * Nc + j) * dim + k] = 2. * m * PetscRealPart(x[i * dim + k]);
185:         }
186:       }
187:       if (H) {
188:         if (type == PETSC_SCALAR) {
189:           for (k = 0; k < dim; k++)
190:             for (l = 0; l < dim; l++) ((PetscScalar *)H)[((i * Nc + j) * dim + k) * dim + l] = (k == l) ? 2. * m : 0.;
191:         } else {
192:           for (k = 0; k < dim; k++)
193:             for (l = 0; l < dim; l++) ((PetscReal *)H)[((i * Nc + j) * dim + k) * dim + l] = (k == l) ? 2. * m : 0.;
194:         }
195:       }
196:     }
197:   }
198:   VecRestoreArrayRead(points, &x);
199:   VecRestoreArrayRead(ctxVec, &mult);
200:   return 0;
201: }

203: static PetscErrorCode TestShellDestroy(DMField field)
204: {
205:   Vec ctxVec = NULL;

207:   DMFieldShellGetContext(field, &ctxVec);
208:   VecDestroy(&ctxVec);
209:   return 0;
210: }

212: int main(int argc, char **argv)
213: {
214:   DM              dm = NULL;
215:   MPI_Comm        comm;
216:   char            type[256] = DMPLEX;
217:   PetscBool       isda, isplex;
218:   PetscInt        dim    = 2;
219:   DMField         field  = NULL;
220:   PetscInt        nc     = 1;
221:   PetscInt        cStart = -1, cEnd = -1;
222:   PetscRandom     rand;
223:   PetscQuadrature quad          = NULL;
224:   PetscInt        pointsPerEdge = 2;
225:   PetscInt        numPoint = 0, numFE = 0, numFV = 0;
226:   PetscBool       testShell = PETSC_FALSE;

229:   PetscInitialize(&argc, &argv, NULL, help);
230:   comm = PETSC_COMM_WORLD;
231:   PetscOptionsBegin(comm, "", "DMField Tutorial Options", "DM");
232:   PetscOptionsFList("-dm_type", "DM implementation on which to define field", "ex1.c", DMList, type, type, 256, NULL);
233:   PetscOptionsRangeInt("-dim", "DM intrinsic dimension", "ex1.c", dim, &dim, NULL, 1, 3);
234:   PetscOptionsBoundedInt("-num_components", "Number of components in field", "ex1.c", nc, &nc, NULL, 1);
235:   PetscOptionsBoundedInt("-num_quad_points", "Number of quadrature points per dimension", "ex1.c", pointsPerEdge, &pointsPerEdge, NULL, 1);
236:   PetscOptionsBoundedInt("-num_point_tests", "Number of test points for DMFieldEvaluate()", "ex1.c", numPoint, &numPoint, NULL, 0);
237:   PetscOptionsBoundedInt("-num_fe_tests", "Number of test cells for DMFieldEvaluateFE()", "ex1.c", numFE, &numFE, NULL, 0);
238:   PetscOptionsBoundedInt("-num_fv_tests", "Number of test cells for DMFieldEvaluateFV()", "ex1.c", numFV, &numFV, NULL, 0);
239:   PetscOptionsBool("-test_shell", "Test the DMFIELDSHELL implementation of DMField", "ex1.c", testShell, &testShell, NULL);
240:   PetscOptionsEnd();

243:   PetscStrncmp(type, DMPLEX, 256, &isplex);
244:   PetscStrncmp(type, DMDA, 256, &isda);

246:   PetscRandomCreate(PETSC_COMM_SELF, &rand);
247:   PetscRandomSetFromOptions(rand);
248:   if (isplex) {
249:     PetscInt  overlap = 0;
250:     Vec       fieldvec;
251:     PetscInt  cells[3] = {3, 3, 3};
252:     PetscBool simplex;
253:     PetscFE   fe;

255:     PetscOptionsBegin(comm, "", "DMField DMPlex Options", "DM");
256:     PetscOptionsBoundedInt("-overlap", "DMPlex parallel overlap", "ex1.c", overlap, &overlap, NULL, 0);
257:     PetscOptionsEnd();
258:     if (0) {
259:       DMPlexCreateBoxMesh(comm, 2, PETSC_TRUE, cells, NULL, NULL, NULL, PETSC_TRUE, &dm);
260:     } else {
261:       DMCreate(comm, &dm);
262:       DMSetType(dm, DMPLEX);
263:       DMSetFromOptions(dm);
264:       CHKMEMQ;
265:     }
266:     DMGetDimension(dm, &dim);
267:     DMPlexIsSimplex(dm, &simplex);
268:     if (simplex) PetscDTStroudConicalQuadrature(dim, 1, pointsPerEdge, -1.0, 1.0, &quad);
269:     else PetscDTGaussTensorQuadrature(dim, 1, pointsPerEdge, -1.0, 1.0, &quad);
270:     DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
271:     if (testShell) {
272:       Vec          ctxVec;
273:       PetscInt     i;
274:       PetscScalar *array;

276:       VecCreateSeq(PETSC_COMM_SELF, nc, &ctxVec);
277:       VecSetUp(ctxVec);
278:       VecGetArray(ctxVec, &array);
279:       for (i = 0; i < nc; i++) array[i] = i + 1.;
280:       VecRestoreArray(ctxVec, &array);
281:       DMFieldCreateShell(dm, nc, DMFIELD_VERTEX, (void *)ctxVec, &field);
282:       DMFieldShellSetEvaluate(field, TestShellEvaluate);
283:       DMFieldShellSetDestroy(field, TestShellDestroy);
284:     } else {
285:       PetscFECreateDefault(PETSC_COMM_SELF, dim, nc, simplex, NULL, PETSC_DEFAULT, &fe);
286:       PetscFESetName(fe, "MyPetscFE");
287:       DMSetField(dm, 0, NULL, (PetscObject)fe);
288:       PetscFEDestroy(&fe);
289:       DMCreateDS(dm);
290:       DMCreateLocalVector(dm, &fieldvec);
291:       {
292:         PetscErrorCode (*func[1])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
293:         void *ctxs[1];

295:         func[0] = radiusSquared;
296:         ctxs[0] = NULL;

298:         DMProjectFunctionLocal(dm, 0.0, func, ctxs, INSERT_ALL_VALUES, fieldvec);
299:       }
300:       DMFieldCreateDS(dm, 0, fieldvec, &field);
301:       VecDestroy(&fieldvec);
302:     }
303:   } else if (isda) {
304:     PetscInt     i;
305:     PetscScalar *cv;

307:     switch (dim) {
308:     case 1:
309:       DMDACreate1d(comm, DM_BOUNDARY_NONE, 3, 1, 1, NULL, &dm);
310:       break;
311:     case 2:
312:       DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, 3, 3, PETSC_DETERMINE, PETSC_DETERMINE, 1, 1, NULL, NULL, &dm);
313:       break;
314:     default:
315:       DMDACreate3d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, 3, 3, 3, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE, 1, 1, NULL, NULL, NULL, &dm);
316:       break;
317:     }
318:     DMSetUp(dm);
319:     DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);
320:     PetscMalloc1(nc * (1 << dim), &cv);
321:     for (i = 0; i < nc * (1 << dim); i++) {
322:       PetscReal rv;

324:       PetscRandomGetValueReal(rand, &rv);
325:       cv[i] = rv;
326:     }
327:     DMFieldCreateDA(dm, nc, cv, &field);
328:     PetscFree(cv);
329:     PetscDTGaussTensorQuadrature(dim, 1, pointsPerEdge, -1.0, 1.0, &quad);
330:   } else SETERRQ(comm, PETSC_ERR_SUP, "This test does not run for DM type %s", type);

332:   PetscObjectSetName((PetscObject)dm, "mesh");
333:   DMViewFromOptions(dm, NULL, "-dm_view");
334:   DMDestroy(&dm);
335:   PetscObjectSetName((PetscObject)field, "field");
336:   PetscObjectViewFromOptions((PetscObject)field, NULL, "-dmfield_view");
337:   if (numPoint) TestEvaluate(field, numPoint, rand);
338:   if (numFE) TestEvaluateFE(field, numFE, cStart, cEnd, quad, rand);
339:   if (numFV) TestEvaluateFV(field, numFV, cStart, cEnd, rand);
340:   DMFieldDestroy(&field);
341:   PetscQuadratureDestroy(&quad);
342:   PetscRandomDestroy(&rand);
343:   PetscFinalize();
344:   return 0;
345: }

347: /*TEST

349:   test:
350:     suffix: da
351:     requires: !complex
352:     args: -dm_type da -dim 2 -num_components 2 -num_point_tests 2 -num_fe_tests 2 -num_fv_tests 2 -dmfield_view

354:   test:
355:     suffix: da_1
356:     requires: !complex
357:     args: -dm_type da -dim 1  -num_fe_tests 2

359:   test:
360:     suffix: da_2
361:     requires: !complex
362:     args: -dm_type da -dim 2  -num_fe_tests 2

364:   test:
365:     suffix: da_3
366:     requires: !complex
367:     args: -dm_type da -dim 3  -num_fe_tests 2

369:   test:
370:     suffix: ds
371:     requires: !complex triangle
372:     args: -dm_type plex -dm_plex_dim 2 -dm_plex_box_faces 3,3 -num_components 2 -num_point_tests 2 -num_fe_tests 2 -num_fv_tests 2 -dmfield_view -petscspace_degree 2 -num_quad_points 1

374:   test:
375:     suffix: ds_simplex_0
376:     requires: !complex triangle
377:     args: -dm_type plex -dm_plex_dim 2 -dm_plex_box_faces 3,3  -num_fe_tests 2  -petscspace_degree 0

379:   test:
380:     suffix: ds_simplex_1
381:     requires: !complex triangle
382:     args: -dm_type plex -dm_plex_dim 2 -dm_plex_box_faces 3,3  -num_fe_tests 2  -petscspace_degree 1

384:   test:
385:     suffix: ds_simplex_2
386:     requires: !complex triangle
387:     args: -dm_type plex -dm_plex_dim 2 -dm_plex_box_faces 3,3  -num_fe_tests 2  -petscspace_degree 2

389:   test:
390:     suffix: ds_tensor_2_0
391:     requires: !complex
392:     args: -dm_type plex -dm_plex_dim 2 -dm_plex_box_faces 3,3  -num_fe_tests 2  -petscspace_poly_tensor 1 -petscspace_degree 0 -dm_plex_simplex 0

394:   test:
395:     suffix: ds_tensor_2_1
396:     requires: !complex
397:     args: -dm_type plex -dm_plex_dim 2 -dm_plex_box_faces 3,3  -num_fe_tests 2  -petscspace_poly_tensor 1 -petscspace_degree 1 -dm_plex_simplex 0

399:   test:
400:     suffix: ds_tensor_2_2
401:     requires: !complex
402:     args: -dm_type plex -dm_plex_dim 2 -dm_plex_box_faces 3,3  -num_fe_tests 2  -petscspace_poly_tensor 1 -petscspace_degree 2 -dm_plex_simplex 0

404:   test:
405:     suffix: ds_tensor_3_0
406:     requires: !complex
407:     args: -dm_type plex -dm_plex_dim 3 -dm_plex_box_faces 3,3,3  -num_fe_tests 2  -petscspace_poly_tensor 1 -petscspace_degree 0 -dm_plex_simplex 0

409:   test:
410:     suffix: ds_tensor_3_1
411:     requires: !complex
412:     args: -dm_type plex -dm_plex_dim 3 -dm_plex_box_faces 3,3,3  -num_fe_tests 2  -petscspace_poly_tensor 1 -petscspace_degree 1 -dm_plex_simplex 0

414:   test:
415:     suffix: ds_tensor_3_2
416:     requires: !complex
417:     args: -dm_type plex -dm_plex_dim 3 -dm_plex_box_faces 3,3,3  -num_fe_tests 2  -petscspace_poly_tensor 1 -petscspace_degree 2 -dm_plex_simplex 0

419:   test:
420:     suffix: shell
421:     requires: !complex triangle
422:     args: -dm_coord_space 0 -dm_type plex -dm_plex_dim 2 -dm_plex_box_faces 3,3 -num_components 2 -num_point_tests 2 -num_fe_tests 2 -num_fv_tests 2 -dmfield_view -num_quad_points 1 -test_shell

424: TEST*/