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++) {
149:     u[i] = (i + 1) * r2;
150:   }
151:   return 0;
152: }

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

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

172:     for (j = 0; j < dim; j++) {r2 += PetscSqr(PetscRealPart(x[i * dim + j]));}
173:     for (j = 0; j < Nc; j++) {
174:       PetscReal m = PetscRealPart(mult[j]);
175:       if (B) {
176:         if (type == PETSC_SCALAR) {
177:           ((PetscScalar *)B)[i * Nc + j] = m * r2;
178:         } else {
179:           ((PetscReal *)B)[i * Nc + j] = m * r2;
180:         }
181:       }
182:       if (D) {
183:         if (type == PETSC_SCALAR) {
184:           for (k = 0; k < dim; k++) ((PetscScalar *)D)[(i * Nc + j) * dim + k] = 2. * m * x[i * dim + k];
185:         } else {
186:           for (k = 0; k < dim; k++) ((PetscReal   *)D)[(i * Nc + j) * dim + k] = 2. * m * PetscRealPart(x[i * dim + k]);
187:         }
188:       }
189:       if (H) {
190:         if (type == PETSC_SCALAR) {
191:           for (k = 0; k < dim; k++) for (l = 0; l < dim; l++) ((PetscScalar *)H)[((i * Nc + j) * dim + k) * dim + l] = (k == l) ? 2. * m : 0.;
192:         } else {
193:           for (k = 0; k < dim; k++) 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;
227:   PetscErrorCode  ierr;

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) {
269:       PetscDTStroudConicalQuadrature(dim, 1, pointsPerEdge, -1.0, 1.0, &quad);
270:     } else {
271:       PetscDTGaussTensorQuadrature(dim, 1, pointsPerEdge, -1.0, 1.0, &quad);
272:     }
273:     DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
274:     if (testShell) {
275:       Vec ctxVec;
276:       PetscInt i;
277:       PetscScalar *array;

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

298:         func[0] = radiusSquared;
299:         ctxs[0] = NULL;

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

310:     switch (dim) {
311:     case 1:
312:       DMDACreate1d(comm, DM_BOUNDARY_NONE, 3, 1, 1, NULL, &dm);
313:       break;
314:     case 2:
315:       DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, 3, 3, PETSC_DETERMINE, PETSC_DETERMINE, 1, 1, NULL, NULL, &dm);
316:       break;
317:     default:
318:       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);
319:       break;
320:     }
321:     DMSetUp(dm);
322:     DMDAGetHeightStratum(dm,0,&cStart,&cEnd);
323:     PetscMalloc1(nc * (1 << dim),&cv);
324:     for (i = 0; i < nc * (1 << dim); i++) {
325:       PetscReal rv;

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

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

350: /*TEST

352:   test:
353:     suffix: da
354:     requires: !complex
355:     args: -dm_type da -dim 2 -num_components 2 -num_point_tests 2 -num_fe_tests 2 -num_fv_tests 2 -dmfield_view

357:   test:
358:     suffix: da_1
359:     requires: !complex
360:     args: -dm_type da -dim 1  -num_fe_tests 2

362:   test:
363:     suffix: da_2
364:     requires: !complex
365:     args: -dm_type da -dim 2  -num_fe_tests 2

367:   test:
368:     suffix: da_3
369:     requires: !complex
370:     args: -dm_type da -dim 3  -num_fe_tests 2

372:   test:
373:     suffix: ds
374:     requires: !complex triangle
375:     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

377:   test:
378:     suffix: ds_simplex_0
379:     requires: !complex triangle
380:     args: -dm_type plex -dm_plex_dim 2 -dm_plex_box_faces 3,3  -num_fe_tests 2  -petscspace_degree 0

382:   test:
383:     suffix: ds_simplex_1
384:     requires: !complex triangle
385:     args: -dm_type plex -dm_plex_dim 2 -dm_plex_box_faces 3,3  -num_fe_tests 2  -petscspace_degree 1

387:   test:
388:     suffix: ds_simplex_2
389:     requires: !complex triangle
390:     args: -dm_type plex -dm_plex_dim 2 -dm_plex_box_faces 3,3  -num_fe_tests 2  -petscspace_degree 2

392:   test:
393:     suffix: ds_tensor_2_0
394:     requires: !complex
395:     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

397:   test:
398:     suffix: ds_tensor_2_1
399:     requires: !complex
400:     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

402:   test:
403:     suffix: ds_tensor_2_2
404:     requires: !complex
405:     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

407:   test:
408:     suffix: ds_tensor_3_0
409:     requires: !complex
410:     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

412:   test:
413:     suffix: ds_tensor_3_1
414:     requires: !complex
415:     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

417:   test:
418:     suffix: ds_tensor_3_2
419:     requires: !complex
420:     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

422:   test:
423:     suffix: shell
424:     requires: !complex triangle
425:     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

427: TEST*/