Actual source code: ex1.c
petsc-3.10.5 2019-03-28
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: PetscOptionsInt("-dim","DM intrinsic dimension", "ex1.c", dim, &dim, NULL);
248: PetscOptionsInt("-num_components","Number of components in field", "ex1.c", nc, &nc, NULL);
249: PetscOptionsInt("-num_quad_points","Number of quadrature points per dimension", "ex1.c", pointsPerEdge, &pointsPerEdge, NULL);
250: PetscOptionsInt("-num_point_tests", "Number of test points for DMFieldEvaluate()", "ex1.c", numPoint, &numPoint, NULL);
251: PetscOptionsInt("-num_fe_tests", "Number of test cells for DMFieldEvaluateFE()", "ex1.c", numFE, &numFE, NULL);
252: PetscOptionsInt("-num_fv_tests", "Number of test cells for DMFieldEvaluateFV()", "ex1.c", numFV, &numFV, NULL);
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: PetscOptionsInt("-overlap","DMPlex parallel overlap","ex1.c",overlap,&overlap,NULL);
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: PetscDTGaussJacobiQuadrature(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: DMSetField(dm,0,(PetscObject)fe);
321: PetscFEDestroy(&fe);
322: DMCreateLocalVector(dm,&fieldvec);
323: {
324: PetscErrorCode (*func[1]) (PetscInt,PetscReal,const PetscReal [],PetscInt, PetscScalar *,void *);
325: void *ctxs[1];
327: func[0] = radiusSquared;
328: ctxs[0] = NULL;
330: DMProjectFunctionLocal(dm,0.0,func,ctxs,INSERT_ALL_VALUES,fieldvec);
331: }
332: DMFieldCreateDS(dm,0,fieldvec,&field);
333: VecDestroy(&fieldvec);
334: }
335: } else if (isda) {
336: PetscInt i;
337: PetscScalar *cv;
339: switch (dim) {
340: case 1:
341: DMDACreate1d(comm, DM_BOUNDARY_NONE, 3, 1, 1, NULL, &dm);
342: break;
343: case 2:
344: DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, 3, 3, PETSC_DETERMINE, PETSC_DETERMINE, 1, 1, NULL, NULL, &dm);
345: break;
346: default:
347: 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);
348: break;
349: }
350: DMSetUp(dm);
351: DMDAGetHeightStratum(dm,0,&cStart,&cEnd);
352: PetscMalloc1(nc * (1 << dim),&cv);
353: for (i = 0; i < nc * (1 << dim); i++) {
354: PetscReal rv;
356: PetscRandomGetValueReal(rand,&rv);
357: cv[i] = rv;
358: }
359: DMFieldCreateDA(dm,nc,cv,&field);
360: PetscFree(cv);
361: PetscDTGaussTensorQuadrature(dim, 1, pointsPerEdge, -1.0, 1.0, &quad);
362: } else SETERRQ1(comm,PETSC_ERR_SUP,"This test does not run for DM type %s",type);
364: PetscObjectSetName((PetscObject)dm,"mesh");
365: DMViewFromOptions(dm,NULL,"-dm_view");
366: DMDestroy(&dm);
367: PetscObjectSetName((PetscObject)field,"field");
368: PetscObjectViewFromOptions((PetscObject)field,NULL,"-dmfield_view");
369: if (numPoint) {TestEvaluate(field,numPoint,rand);}
370: if (numFE) {TestEvaluateFE(field,numFE,cStart,cEnd,quad,rand);}
371: if (numFV) {TestEvaluateFV(field,numFV,cStart,cEnd,rand);}
372: PetscQuadratureDestroy(&quad);
373: DMFieldDestroy(&field);
374: PetscRandomDestroy(&rand);
375: PetscFinalize();
376: return ierr;
377: }
379: /*TEST
381: test:
382: suffix: da
383: requires: !complex
384: args: -dm_type da -dim 2 -num_components 2 -num_point_tests 2 -num_fe_tests 2 -num_fv_tests 2 -dmfield_view
386: test:
387: suffix: da_1
388: requires: !complex
389: args: -dm_type da -dim 1 -num_fe_tests 2
391: test:
392: suffix: da_2
393: requires: !complex
394: args: -dm_type da -dim 2 -num_fe_tests 2
396: test:
397: suffix: da_3
398: requires: !complex
399: args: -dm_type da -dim 3 -num_fe_tests 2
401: test:
402: suffix: ds
403: requires: !complex triangle
404: 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
406: test:
407: suffix: ds_simplex_0
408: requires: !complex triangle
409: args: -dm_type plex -dim 2 -num_fe_tests 2 -petscspace_degree 0
411: test:
412: suffix: ds_simplex_1
413: requires: !complex triangle
414: args: -dm_type plex -dim 2 -num_fe_tests 2 -petscspace_degree 1
416: test:
417: suffix: ds_simplex_2
418: requires: !complex triangle
419: args: -dm_type plex -dim 2 -num_fe_tests 2 -petscspace_degree 2
421: test:
422: suffix: ds_tensor_2_0
423: requires: !complex
424: args: -dm_type plex -dim 2 -num_fe_tests 2 -petscspace_poly_tensor 1 -petscspace_degree 0 -simplex 0
426: test:
427: suffix: ds_tensor_2_1
428: requires: !complex
429: args: -dm_type plex -dim 2 -num_fe_tests 2 -petscspace_poly_tensor 1 -petscspace_degree 1 -simplex 0
431: test:
432: suffix: ds_tensor_2_2
433: requires: !complex
434: args: -dm_type plex -dim 2 -num_fe_tests 2 -petscspace_poly_tensor 1 -petscspace_degree 2 -simplex 0
436: test:
437: suffix: ds_tensor_3_0
438: requires: !complex
439: args: -dm_type plex -dim 3 -num_fe_tests 2 -petscspace_poly_tensor 1 -petscspace_degree 0 -simplex 0
441: test:
442: suffix: ds_tensor_3_1
443: requires: !complex
444: args: -dm_type plex -dim 3 -num_fe_tests 2 -petscspace_poly_tensor 1 -petscspace_degree 1 -simplex 0
446: test:
447: suffix: ds_tensor_3_2
448: requires: !complex
449: args: -dm_type plex -dim 3 -num_fe_tests 2 -petscspace_poly_tensor 1 -petscspace_degree 2 -simplex 0
451: test:
452: suffix: shell
453: requires: !complex triangle
454: 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
456: TEST*/