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: {
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, &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, &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: PetscInt overlap = 0;
264: Vec fieldvec;
265: PetscInt cells[3] = {3,3,3};
266: PetscBool simplex;
267: PetscFE fe;
269: PetscOptionsBegin(comm, "", "DMField DMPlex Options", "DM");
270: PetscOptionsBoundedInt("-overlap","DMPlex parallel overlap","ex1.c",overlap,&overlap,NULL,0);
271: PetscOptionsEnd();
272: if (0) {
273: DMPlexCreateBoxMesh(comm,2,PETSC_TRUE,cells,NULL,NULL,NULL,PETSC_TRUE,&dm);
274: } else {
275: DMCreate(comm, &dm);
276: DMSetType(dm, DMPLEX);
277: DMSetFromOptions(dm);
278: CHKMEMQ;
279: }
280: DMGetDimension(dm, &dim);
281: DMPlexIsSimplex(dm, &simplex);
282: if (simplex) {
283: PetscDTStroudConicalQuadrature(dim, 1, pointsPerEdge, -1.0, 1.0, &quad);
284: } else {
285: PetscDTGaussTensorQuadrature(dim, 1, pointsPerEdge, -1.0, 1.0, &quad);
286: }
287: DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
288: if (testShell) {
289: Vec ctxVec;
290: PetscInt i;
291: PetscScalar *array;
293: VecCreateSeq(PETSC_COMM_SELF, nc, &ctxVec);
294: VecSetUp(ctxVec);
295: VecGetArray(ctxVec,&array);
296: for (i = 0; i < nc; i++) array[i] = i + 1.;
297: VecRestoreArray(ctxVec,&array);
298: DMFieldCreateShell(dm, nc, DMFIELD_VERTEX, (void *) ctxVec, &field);
299: DMFieldShellSetEvaluate(field, TestShellEvaluate);
300: DMFieldShellSetDestroy(field, TestShellDestroy);
301: } else {
302: PetscFECreateDefault(PETSC_COMM_SELF,dim,nc,simplex,NULL,PETSC_DEFAULT,&fe);
303: PetscFESetName(fe,"MyPetscFE");
304: DMSetField(dm,0,NULL,(PetscObject)fe);
305: PetscFEDestroy(&fe);
306: DMCreateDS(dm);
307: DMCreateLocalVector(dm,&fieldvec);
308: {
309: PetscErrorCode (*func[1]) (PetscInt,PetscReal,const PetscReal [],PetscInt, PetscScalar *,void *);
310: void *ctxs[1];
312: func[0] = radiusSquared;
313: ctxs[0] = NULL;
315: DMProjectFunctionLocal(dm,0.0,func,ctxs,INSERT_ALL_VALUES,fieldvec);
316: }
317: DMFieldCreateDS(dm,0,fieldvec,&field);
318: VecDestroy(&fieldvec);
319: }
320: } else if (isda) {
321: PetscInt i;
322: PetscScalar *cv;
324: switch (dim) {
325: case 1:
326: DMDACreate1d(comm, DM_BOUNDARY_NONE, 3, 1, 1, NULL, &dm);
327: break;
328: case 2:
329: DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, 3, 3, PETSC_DETERMINE, PETSC_DETERMINE, 1, 1, NULL, NULL, &dm);
330: break;
331: default:
332: 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);
333: break;
334: }
335: DMSetUp(dm);
336: DMDAGetHeightStratum(dm,0,&cStart,&cEnd);
337: PetscMalloc1(nc * (1 << dim),&cv);
338: for (i = 0; i < nc * (1 << dim); i++) {
339: PetscReal rv;
341: PetscRandomGetValueReal(rand,&rv);
342: cv[i] = rv;
343: }
344: DMFieldCreateDA(dm,nc,cv,&field);
345: PetscFree(cv);
346: PetscDTGaussTensorQuadrature(dim, 1, pointsPerEdge, -1.0, 1.0, &quad);
347: } else SETERRQ1(comm,PETSC_ERR_SUP,"This test does not run for DM type %s",type);
349: PetscObjectSetName((PetscObject)dm,"mesh");
350: DMViewFromOptions(dm,NULL,"-dm_view");
351: DMDestroy(&dm);
352: PetscObjectSetName((PetscObject)field,"field");
353: PetscObjectViewFromOptions((PetscObject)field,NULL,"-dmfield_view");
354: if (numPoint) {TestEvaluate(field,numPoint,rand);}
355: if (numFE) {TestEvaluateFE(field,numFE,cStart,cEnd,quad,rand);}
356: if (numFV) {TestEvaluateFV(field,numFV,cStart,cEnd,rand);}
357: DMFieldDestroy(&field);
358: PetscQuadratureDestroy(&quad);
359: PetscRandomDestroy(&rand);
360: PetscFinalize();
361: return ierr;
362: }
364: /*TEST
366: test:
367: suffix: da
368: requires: !complex
369: args: -dm_type da -dim 2 -num_components 2 -num_point_tests 2 -num_fe_tests 2 -num_fv_tests 2 -dmfield_view
371: test:
372: suffix: da_1
373: requires: !complex
374: args: -dm_type da -dim 1 -num_fe_tests 2
376: test:
377: suffix: da_2
378: requires: !complex
379: args: -dm_type da -dim 2 -num_fe_tests 2
381: test:
382: suffix: da_3
383: requires: !complex
384: args: -dm_type da -dim 3 -num_fe_tests 2
386: test:
387: suffix: ds
388: requires: !complex triangle
389: 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
391: test:
392: suffix: ds_simplex_0
393: requires: !complex triangle
394: args: -dm_type plex -dm_plex_dim 2 -dm_plex_box_faces 3,3 -num_fe_tests 2 -petscspace_degree 0
396: test:
397: suffix: ds_simplex_1
398: requires: !complex triangle
399: args: -dm_type plex -dm_plex_dim 2 -dm_plex_box_faces 3,3 -num_fe_tests 2 -petscspace_degree 1
401: test:
402: suffix: ds_simplex_2
403: requires: !complex triangle
404: args: -dm_type plex -dm_plex_dim 2 -dm_plex_box_faces 3,3 -num_fe_tests 2 -petscspace_degree 2
406: test:
407: suffix: ds_tensor_2_0
408: requires: !complex
409: 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
411: test:
412: suffix: ds_tensor_2_1
413: requires: !complex
414: 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
416: test:
417: suffix: ds_tensor_2_2
418: requires: !complex
419: 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
421: test:
422: suffix: ds_tensor_3_0
423: requires: !complex
424: 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
426: test:
427: suffix: ds_tensor_3_1
428: requires: !complex
429: 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
431: test:
432: suffix: ds_tensor_3_2
433: requires: !complex
434: 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
436: test:
437: suffix: shell
438: requires: !complex triangle
439: 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
441: TEST*/