Actual source code: ex1.c
petsc-3.11.4 2019-09-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: 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*/