Actual source code: ex1.c

petsc-3.14.6 2021-03-30
Report Typos and Errors
  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: {

 12:   PetscViewerASCIIPrintf(viewer,"B:\n");
 13:   PetscScalarView(N,B,viewer);
 14:   PetscViewerASCIIPrintf(viewer,"D:\n");
 15:   PetscScalarView(N*dim,D,viewer);
 16:   PetscViewerASCIIPrintf(viewer,"H:\n");
 17:   PetscScalarView(N*dim*dim,H,viewer);

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

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

 41:   comm = PetscObjectComm((PetscObject)field);
 42:   DMFieldGetNumComponents(field,&nc);
 43:   DMFieldGetDM(field,&dm);
 44:   DMGetDimension(dm,&dim);
 45:   VecCreateMPI(PetscObjectComm((PetscObject)field),n * dim,PETSC_DETERMINE,&points);
 46:   VecSetBlockSize(points,dim);
 47:   VecGetArray(points,&array);
 48:   for (i = 0; i < n * dim; i++) {PetscRandomGetValue(rand,&array[i]);}
 49:   VecRestoreArray(points,&array);
 50:   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);
 51:   DMFieldEvaluate(field,points,PETSC_SCALAR,B,D,H);
 52:   DMFieldEvaluate(field,points,PETSC_REAL,rB,rD,rH);
 53:   viewer = PETSC_VIEWER_STDOUT_(comm);

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

 59:   PetscFree6(B,rB,D,rD,H,rH);
 60:   VecDestroy(&points);
 61:   return(0);
 62: }

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

 78:   comm = PetscObjectComm((PetscObject)field);
 79:   DMFieldGetNumComponents(field,&nc);
 80:   DMFieldGetDM(field,&dm);
 81:   DMGetDimension(dm,&dim);
 82:   PetscRandomSetInterval(rand,(PetscScalar) cStart, (PetscScalar) cEnd);
 83:   PetscMalloc1(n,&cells);
 84:   for (i = 0; i < n; i++) {
 85:     PetscReal rc;

 87:     PetscRandomGetValueReal(rand,&rc);
 88:     cells[i] = PetscFloorReal(rc);
 89:   }
 90:   ISCreateGeneral(comm,n,cells,PETSC_OWN_POINTER,&cellIS);
 91:   PetscObjectSetName((PetscObject)cellIS,"FE Test Cells");
 92:   PetscQuadratureGetData(quad,NULL,NULL,&nq,NULL,NULL);
 93:   N    = n * nq * nc;
 94:   PetscMalloc6(N,&B,N,&rB,N*dim,&D,N*dim,&rD,N*dim*dim,&H,N*dim*dim,&rH);
 95:   DMFieldEvaluateFE(field,cellIS,quad,PETSC_SCALAR,B,D,H);
 96:   DMFieldEvaluateFE(field,cellIS,quad,PETSC_REAL,rB,rD,rH);
 97:   viewer = PETSC_VIEWER_STDOUT_(comm);

 99:   PetscObjectSetName((PetscObject)quad,"Test quadrature");
100:   PetscQuadratureView(quad,viewer);
101:   ISView(cellIS,viewer);
102:   ViewResults(viewer,N,dim,B,D,H,rB,rD,rH);

104:   PetscFree6(B,rB,D,rD,H,rH);
105:   ISDestroy(&cellIS);
106:   return(0);
107: }

109: static PetscErrorCode TestEvaluateFV(DMField field, PetscInt n, PetscInt cStart, PetscInt cEnd, PetscRandom rand)
110: {
111:   DM             dm;
112:   PetscInt       dim, i, nc;
113:   PetscInt       N;
114:   PetscScalar    *B, *D, *H;
115:   PetscReal      *rB, *rD, *rH;
116:   PetscInt       *cells;
117:   IS             cellIS;
118:   PetscViewer    viewer;
119:   MPI_Comm       comm;

123:   comm = PetscObjectComm((PetscObject)field);
124:   DMFieldGetNumComponents(field,&nc);
125:   DMFieldGetDM(field,&dm);
126:   DMGetDimension(dm,&dim);
127:   PetscRandomSetInterval(rand,(PetscScalar) cStart, (PetscScalar) cEnd);
128:   PetscMalloc1(n,&cells);
129:   for (i = 0; i < n; i++) {
130:     PetscReal rc;

132:     PetscRandomGetValueReal(rand,&rc);
133:     cells[i] = PetscFloorReal(rc);
134:   }
135:   ISCreateGeneral(comm,n,cells,PETSC_OWN_POINTER,&cellIS);
136:   PetscObjectSetName((PetscObject)cellIS,"FV Test Cells");
137:   N    = n * nc;
138:   PetscMalloc6(N,&B,N,&rB,N*dim,&D,N*dim,&rD,N*dim*dim,&H,N*dim*dim,&rH);
139:   DMFieldEvaluateFV(field,cellIS,PETSC_SCALAR,B,D,H);
140:   DMFieldEvaluateFV(field,cellIS,PETSC_REAL,rB,rD,rH);
141:   viewer = PETSC_VIEWER_STDOUT_(comm);

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

146:   PetscFree6(B,rB,D,rD,H,rH);
147:   ISDestroy(&cellIS);
148:   return(0);
149: }

151: static PetscErrorCode radiusSquared(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx)
152: {
153:   PetscInt       i;
154:   PetscReal      r2 = 0.;

157:   for (i = 0; i < dim; i++) {r2 += PetscSqr(x[i]);}
158:   for (i = 0; i < Nf; i++) {
159:     u[i] = (i + 1) * r2;
160:   }
161:   return(0);
162: }

164: static PetscErrorCode TestShellEvaluate(DMField field, Vec points, PetscDataType type, void *B, void *D, void *H)
165: {
166:   Vec                ctxVec = NULL;
167:   const PetscScalar *mult;
168:   PetscInt           dim;
169:   const PetscScalar *x;
170:   PetscInt           Nc, n, i, j, k, l;
171:   PetscErrorCode     ierr;

174:   DMFieldGetNumComponents(field, &Nc);
175:   DMFieldShellGetContext(field, (void *) &ctxVec);
176:   VecGetBlockSize(points, &dim);
177:   VecGetLocalSize(points, &n);
178:   n /= Nc;
179:   VecGetArrayRead(ctxVec, &mult);
180:   VecGetArrayRead(points, &x);
181:   for (i = 0; i < n; i++) {
182:     PetscReal r2 = 0.;

184:     for (j = 0; j < dim; j++) {r2 += PetscSqr(PetscRealPart(x[i * dim + j]));}
185:     for (j = 0; j < Nc; j++) {
186:       PetscReal m = PetscRealPart(mult[j]);
187:       if (B) {
188:         if (type == PETSC_SCALAR) {
189:           ((PetscScalar *)B)[i * Nc + j] = m * r2;
190:         } else {
191:           ((PetscReal *)B)[i * Nc + j] = m * r2;
192:         }
193:       }
194:       if (D) {
195:         if (type == PETSC_SCALAR) {
196:           for (k = 0; k < dim; k++) ((PetscScalar *)D)[(i * Nc + j) * dim + k] = 2. * m * x[i * dim + k];
197:         } else {
198:           for (k = 0; k < dim; k++) ((PetscReal   *)D)[(i * Nc + j) * dim + k] = 2. * m * PetscRealPart(x[i * dim + k]);
199:         }
200:       }
201:       if (H) {
202:         if (type == PETSC_SCALAR) {
203:           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.;
204:         } else {
205:           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.;
206:         }
207:       }
208:     }
209:   }
210:   VecRestoreArrayRead(points, &x);
211:   VecRestoreArrayRead(ctxVec, &mult);
212:   return(0);
213: }

215: static PetscErrorCode TestShellDestroy(DMField field)
216: {
217:   Vec                ctxVec = NULL;
218:   PetscErrorCode     ierr;

221:   DMFieldShellGetContext(field, (void *) &ctxVec);
222:   VecDestroy(&ctxVec);
223:   return(0);
224: }

226: int main(int argc, char **argv)
227: {
228:   DM              dm = NULL;
229:   MPI_Comm        comm;
230:   char            type[256] = DMPLEX;
231:   PetscBool       isda, isplex;
232:   PetscInt        dim = 2;
233:   DMField         field = NULL;
234:   PetscInt        nc = 1;
235:   PetscInt        cStart = -1, cEnd = -1;
236:   PetscRandom     rand;
237:   PetscQuadrature quad = NULL;
238:   PetscInt        pointsPerEdge = 2;
239:   PetscInt        numPoint = 0, numFE = 0, numFV = 0;
240:   PetscBool       testShell = PETSC_FALSE;
241:   PetscErrorCode  ierr;

243:   PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr;
244:   comm = PETSC_COMM_WORLD;
245:   PetscOptionsBegin(comm, "", "DMField Tutorial Options", "DM");
246:   PetscOptionsFList("-dm_type","DM implementation on which to define field","ex1.c",DMList,type,type,256,NULL);
247:   PetscOptionsRangeInt("-dim","DM intrinsic dimension", "ex1.c", dim, &dim, NULL,1,3);
248:   PetscOptionsBoundedInt("-num_components","Number of components in field", "ex1.c", nc, &nc, NULL,1);
249:   PetscOptionsBoundedInt("-num_quad_points","Number of quadrature points per dimension", "ex1.c", pointsPerEdge, &pointsPerEdge, NULL,1);
250:   PetscOptionsBoundedInt("-num_point_tests", "Number of test points for DMFieldEvaluate()", "ex1.c", numPoint, &numPoint, NULL,0);
251:   PetscOptionsBoundedInt("-num_fe_tests", "Number of test cells for DMFieldEvaluateFE()", "ex1.c", numFE, &numFE, NULL,0);
252:   PetscOptionsBoundedInt("-num_fv_tests", "Number of test cells for DMFieldEvaluateFV()", "ex1.c", numFV, &numFV, NULL,0);
253:   PetscOptionsBool("-test_shell", "Test the DMFIELDSHELL implementation of DMField", "ex1.c", testShell, &testShell, NULL);
254:   PetscOptionsEnd();

256:   if (dim > 3) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"This examples works for dim <= 3, not %D",dim);
257:   PetscStrncmp(type,DMPLEX,256,&isplex);
258:   PetscStrncmp(type,DMDA,256,&isda);

260:   PetscRandomCreate(PETSC_COMM_SELF,&rand);
261:   PetscRandomSetFromOptions(rand);
262:   if (isplex) {
263:     PetscBool simplex = PETSC_TRUE;
264:     PetscInt  overlap = 0;
265:     Vec       fieldvec;
266:     PetscBool interpolate = PETSC_TRUE;
267:     PetscInt       inttypes[3];
268:     DMBoundaryType types[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE};
269:     PetscInt       cells[3] = {3,3,3};
270:     PetscBool      flags[3] = {PETSC_FALSE, PETSC_FALSE, PETSC_FALSE};
271:     PetscInt       n = 3, i;
272:     PetscFE   fe;

274:     PetscOptionsBegin(comm, "", "DMField DMPlex Options", "DM");
275:     PetscOptionsBool("-simplex","Create a simplicial DMPlex","ex1.c",simplex,&simplex,NULL);
276:     PetscOptionsBoundedInt("-overlap","DMPlex parallel overlap","ex1.c",overlap,&overlap,NULL,0);
277:     PetscOptionsBool("-interpolate","Interpolate the DMPlex","ex1.c",interpolate,&interpolate,NULL);
278:     PetscOptionsEList("-boundary_x","type of boundary in x direction","ex1.c",DMBoundaryTypes,5,DMBoundaryTypes[types[0]],&inttypes[0],&flags[0]);
279:     PetscOptionsEList("-boundary_y","type of boundary in y direction","ex1.c",DMBoundaryTypes,5,DMBoundaryTypes[types[1]],&inttypes[1],&flags[1]);
280:     PetscOptionsEList("-boundary_z","type of boundary in z direction","ex1.c",DMBoundaryTypes,5,DMBoundaryTypes[types[2]],&inttypes[2],&flags[2]);
281:     PetscOptionsIntArray("-cells","Cells per dimension","ex1.c",cells,&n,NULL);
282:     PetscOptionsEnd();
283:     for (i = 0; i < 3; i++) {
284:       if (flags[i]) {types[i] = (DMBoundaryType) inttypes[i];}
285:     }
286:     DMPlexCreateBoxMesh(comm,dim,simplex,cells,NULL,NULL,types,interpolate,&dm);
287:     if (simplex) {
288:       PetscDTStroudConicalQuadrature(dim, 1, pointsPerEdge, -1.0, 1.0, &quad);
289:     } else {
290:       PetscDTGaussTensorQuadrature(dim, 1, pointsPerEdge, -1.0, 1.0, &quad);
291:     }
292:     {
293:       PetscPartitioner part;
294:       DM               dmDist;

296:       DMPlexGetPartitioner(dm,&part);
297:       PetscPartitionerSetType(part,PETSCPARTITIONERSIMPLE);
298:       DMPlexDistribute(dm,overlap,NULL,&dmDist);
299:       if (dmDist) {
300:         DMDestroy(&dm);
301:         dm = dmDist;
302:       }
303:     }
304:     DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
305:     if (testShell) {
306:       Vec ctxVec;
307:       PetscInt i;
308:       PetscScalar *array;

310:       VecCreateSeq(PETSC_COMM_SELF, nc, &ctxVec);
311:       VecSetUp(ctxVec);
312:       VecGetArray(ctxVec,&array);
313:       for (i = 0; i < nc; i++) array[i] = i + 1.;
314:       VecRestoreArray(ctxVec,&array);
315:       DMFieldCreateShell(dm, nc, DMFIELD_VERTEX, (void *) ctxVec, &field);
316:       DMFieldShellSetEvaluate(field, TestShellEvaluate);
317:       DMFieldShellSetDestroy(field, TestShellDestroy);
318:     } else {
319:       PetscFECreateDefault(comm,dim,nc,simplex,NULL,PETSC_DEFAULT,&fe);
320:       PetscFESetName(fe,"MyPetscFE");
321:       DMSetField(dm,0,NULL,(PetscObject)fe);
322:       DMCreateDS(dm);
323:       PetscFEDestroy(&fe);
324:       DMCreateLocalVector(dm,&fieldvec);
325:       {
326:         PetscErrorCode (*func[1]) (PetscInt,PetscReal,const PetscReal [],PetscInt, PetscScalar *,void *);
327:         void            *ctxs[1];

329:         func[0] = radiusSquared;
330:         ctxs[0] = NULL;

332:         DMProjectFunctionLocal(dm,0.0,func,ctxs,INSERT_ALL_VALUES,fieldvec);
333:       }
334:       DMFieldCreateDS(dm,0,fieldvec,&field);
335:       VecDestroy(&fieldvec);
336:     }
337:   } else if (isda) {
338:     PetscInt       i;
339:     PetscScalar    *cv;

341:     switch (dim) {
342:     case 1:
343:       DMDACreate1d(comm, DM_BOUNDARY_NONE, 3, 1, 1, NULL, &dm);
344:       break;
345:     case 2:
346:       DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, 3, 3, PETSC_DETERMINE, PETSC_DETERMINE, 1, 1, NULL, NULL, &dm);
347:       break;
348:     default:
349:       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);
350:       break;
351:     }
352:     DMSetUp(dm);
353:     DMDAGetHeightStratum(dm,0,&cStart,&cEnd);
354:     PetscMalloc1(nc * (1 << dim),&cv);
355:     for (i = 0; i < nc * (1 << dim); i++) {
356:       PetscReal rv;

358:       PetscRandomGetValueReal(rand,&rv);
359:       cv[i] = rv;
360:     }
361:     DMFieldCreateDA(dm,nc,cv,&field);
362:     PetscFree(cv);
363:     PetscDTGaussTensorQuadrature(dim, 1, pointsPerEdge, -1.0, 1.0, &quad);
364:   } else SETERRQ1(comm,PETSC_ERR_SUP,"This test does not run for DM type %s",type);

366:   PetscObjectSetName((PetscObject)dm,"mesh");
367:   DMViewFromOptions(dm,NULL,"-dm_view");
368:   DMDestroy(&dm);
369:   PetscObjectSetName((PetscObject)field,"field");
370:   PetscObjectViewFromOptions((PetscObject)field,NULL,"-dmfield_view");
371:   if (numPoint) {TestEvaluate(field,numPoint,rand);}
372:   if (numFE) {TestEvaluateFE(field,numFE,cStart,cEnd,quad,rand);}
373:   if (numFV) {TestEvaluateFV(field,numFV,cStart,cEnd,rand);}
374:   PetscQuadratureDestroy(&quad);
375:   DMFieldDestroy(&field);
376:   PetscRandomDestroy(&rand);
377:   PetscFinalize();
378:   return ierr;
379: }

381: /*TEST

383:   test:
384:     suffix: da
385:     requires: !complex
386:     args: -dm_type da -dim 2 -num_components 2 -num_point_tests 2 -num_fe_tests 2 -num_fv_tests 2 -dmfield_view

388:   test:
389:     suffix: da_1
390:     requires: !complex
391:     args: -dm_type da -dim 1  -num_fe_tests 2

393:   test:
394:     suffix: da_2
395:     requires: !complex
396:     args: -dm_type da -dim 2  -num_fe_tests 2

398:   test:
399:     suffix: da_3
400:     requires: !complex
401:     args: -dm_type da -dim 3  -num_fe_tests 2

403:   test:
404:     suffix: ds
405:     requires: !complex triangle
406:     args: -dm_type plex -dim 2 -num_components 2 -num_point_tests 2 -num_fe_tests 2 -num_fv_tests 2 -dmfield_view -petscspace_degree 2 -num_quad_points 1

408:   test:
409:     suffix: ds_simplex_0
410:     requires: !complex triangle
411:     args: -dm_type plex -dim 2  -num_fe_tests 2  -petscspace_degree 0

413:   test:
414:     suffix: ds_simplex_1
415:     requires: !complex triangle
416:     args: -dm_type plex -dim 2  -num_fe_tests 2  -petscspace_degree 1

418:   test:
419:     suffix: ds_simplex_2
420:     requires: !complex triangle
421:     args: -dm_type plex -dim 2  -num_fe_tests 2  -petscspace_degree 2

423:   test:
424:     suffix: ds_tensor_2_0
425:     requires: !complex
426:     args: -dm_type plex -dim 2  -num_fe_tests 2  -petscspace_poly_tensor 1 -petscspace_degree 0 -simplex 0

428:   test:
429:     suffix: ds_tensor_2_1
430:     requires: !complex
431:     args: -dm_type plex -dim 2  -num_fe_tests 2  -petscspace_poly_tensor 1 -petscspace_degree 1 -simplex 0

433:   test:
434:     suffix: ds_tensor_2_2
435:     requires: !complex
436:     args: -dm_type plex -dim 2  -num_fe_tests 2  -petscspace_poly_tensor 1 -petscspace_degree 2 -simplex 0

438:   test:
439:     suffix: ds_tensor_3_0
440:     requires: !complex
441:     args: -dm_type plex -dim 3  -num_fe_tests 2  -petscspace_poly_tensor 1 -petscspace_degree 0 -simplex 0

443:   test:
444:     suffix: ds_tensor_3_1
445:     requires: !complex
446:     args: -dm_type plex -dim 3  -num_fe_tests 2  -petscspace_poly_tensor 1 -petscspace_degree 1 -simplex 0

448:   test:
449:     suffix: ds_tensor_3_2
450:     requires: !complex
451:     args: -dm_type plex -dim 3  -num_fe_tests 2  -petscspace_poly_tensor 1 -petscspace_degree 2 -simplex 0

453:   test:
454:     suffix: shell
455:     requires: !complex triangle
456:     args: -dm_type plex -dim 2 -num_components 2 -num_point_tests 2 -num_fe_tests 2 -num_fv_tests 2 -dmfield_view -num_quad_points 1 -test_shell

458: TEST*/