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;

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

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

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

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

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

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

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

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

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

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

349: /*TEST

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

426: TEST*/