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;
227: PetscErrorCode ierr;
229: PetscInitialize(&argc, &argv, NULL, help);
230: comm = PETSC_COMM_WORLD;
231: PetscOptionsBegin(comm, "", "DMField Tutorial Options", "DM");
232: PetscOptionsFList("-dm_type","DM implementation on which to define field","ex1.c",DMList,type,type,256,NULL);
233: PetscOptionsRangeInt("-dim","DM intrinsic dimension", "ex1.c", dim, &dim, NULL,1,3);
234: PetscOptionsBoundedInt("-num_components","Number of components in field", "ex1.c", nc, &nc, NULL,1);
235: PetscOptionsBoundedInt("-num_quad_points","Number of quadrature points per dimension", "ex1.c", pointsPerEdge, &pointsPerEdge, NULL,1);
236: PetscOptionsBoundedInt("-num_point_tests", "Number of test points for DMFieldEvaluate()", "ex1.c", numPoint, &numPoint, NULL,0);
237: PetscOptionsBoundedInt("-num_fe_tests", "Number of test cells for DMFieldEvaluateFE()", "ex1.c", numFE, &numFE, NULL,0);
238: PetscOptionsBoundedInt("-num_fv_tests", "Number of test cells for DMFieldEvaluateFV()", "ex1.c", numFV, &numFV, NULL,0);
239: PetscOptionsBool("-test_shell", "Test the DMFIELDSHELL implementation of DMField", "ex1.c", testShell, &testShell, NULL);
240: PetscOptionsEnd();
243: PetscStrncmp(type,DMPLEX,256,&isplex);
244: PetscStrncmp(type,DMDA,256,&isda);
246: PetscRandomCreate(PETSC_COMM_SELF,&rand);
247: PetscRandomSetFromOptions(rand);
248: if (isplex) {
249: PetscInt overlap = 0;
250: Vec fieldvec;
251: PetscInt cells[3] = {3,3,3};
252: PetscBool simplex;
253: PetscFE fe;
255: PetscOptionsBegin(comm, "", "DMField DMPlex Options", "DM");
256: PetscOptionsBoundedInt("-overlap","DMPlex parallel overlap","ex1.c",overlap,&overlap,NULL,0);
257: PetscOptionsEnd();
258: if (0) {
259: DMPlexCreateBoxMesh(comm,2,PETSC_TRUE,cells,NULL,NULL,NULL,PETSC_TRUE,&dm);
260: } else {
261: DMCreate(comm, &dm);
262: DMSetType(dm, DMPLEX);
263: DMSetFromOptions(dm);
264: CHKMEMQ;
265: }
266: DMGetDimension(dm, &dim);
267: DMPlexIsSimplex(dm, &simplex);
268: if (simplex) {
269: PetscDTStroudConicalQuadrature(dim, 1, pointsPerEdge, -1.0, 1.0, &quad);
270: } else {
271: PetscDTGaussTensorQuadrature(dim, 1, pointsPerEdge, -1.0, 1.0, &quad);
272: }
273: DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
274: if (testShell) {
275: Vec ctxVec;
276: PetscInt i;
277: PetscScalar *array;
279: VecCreateSeq(PETSC_COMM_SELF, nc, &ctxVec);
280: VecSetUp(ctxVec);
281: VecGetArray(ctxVec,&array);
282: for (i = 0; i < nc; i++) array[i] = i + 1.;
283: VecRestoreArray(ctxVec,&array);
284: DMFieldCreateShell(dm, nc, DMFIELD_VERTEX, (void *) ctxVec, &field);
285: DMFieldShellSetEvaluate(field, TestShellEvaluate);
286: DMFieldShellSetDestroy(field, TestShellDestroy);
287: } else {
288: PetscFECreateDefault(PETSC_COMM_SELF,dim,nc,simplex,NULL,PETSC_DEFAULT,&fe);
289: PetscFESetName(fe,"MyPetscFE");
290: DMSetField(dm,0,NULL,(PetscObject)fe);
291: PetscFEDestroy(&fe);
292: DMCreateDS(dm);
293: DMCreateLocalVector(dm,&fieldvec);
294: {
295: PetscErrorCode (*func[1]) (PetscInt,PetscReal,const PetscReal [],PetscInt, PetscScalar *,void *);
296: void *ctxs[1];
298: func[0] = radiusSquared;
299: ctxs[0] = NULL;
301: DMProjectFunctionLocal(dm,0.0,func,ctxs,INSERT_ALL_VALUES,fieldvec);
302: }
303: DMFieldCreateDS(dm,0,fieldvec,&field);
304: VecDestroy(&fieldvec);
305: }
306: } else if (isda) {
307: PetscInt i;
308: PetscScalar *cv;
310: switch (dim) {
311: case 1:
312: DMDACreate1d(comm, DM_BOUNDARY_NONE, 3, 1, 1, NULL, &dm);
313: break;
314: case 2:
315: DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, 3, 3, PETSC_DETERMINE, PETSC_DETERMINE, 1, 1, NULL, NULL, &dm);
316: break;
317: default:
318: 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);
319: break;
320: }
321: DMSetUp(dm);
322: DMDAGetHeightStratum(dm,0,&cStart,&cEnd);
323: PetscMalloc1(nc * (1 << dim),&cv);
324: for (i = 0; i < nc * (1 << dim); i++) {
325: PetscReal rv;
327: PetscRandomGetValueReal(rand,&rv);
328: cv[i] = rv;
329: }
330: DMFieldCreateDA(dm,nc,cv,&field);
331: PetscFree(cv);
332: PetscDTGaussTensorQuadrature(dim, 1, pointsPerEdge, -1.0, 1.0, &quad);
333: } else SETERRQ(comm,PETSC_ERR_SUP,"This test does not run for DM type %s",type);
335: PetscObjectSetName((PetscObject)dm,"mesh");
336: DMViewFromOptions(dm,NULL,"-dm_view");
337: DMDestroy(&dm);
338: PetscObjectSetName((PetscObject)field,"field");
339: PetscObjectViewFromOptions((PetscObject)field,NULL,"-dmfield_view");
340: if (numPoint) TestEvaluate(field,numPoint,rand);
341: if (numFE) TestEvaluateFE(field,numFE,cStart,cEnd,quad,rand);
342: if (numFV) TestEvaluateFV(field,numFV,cStart,cEnd,rand);
343: DMFieldDestroy(&field);
344: PetscQuadratureDestroy(&quad);
345: PetscRandomDestroy(&rand);
346: PetscFinalize();
347: return 0;
348: }
350: /*TEST
352: test:
353: suffix: da
354: requires: !complex
355: args: -dm_type da -dim 2 -num_components 2 -num_point_tests 2 -num_fe_tests 2 -num_fv_tests 2 -dmfield_view
357: test:
358: suffix: da_1
359: requires: !complex
360: args: -dm_type da -dim 1 -num_fe_tests 2
362: test:
363: suffix: da_2
364: requires: !complex
365: args: -dm_type da -dim 2 -num_fe_tests 2
367: test:
368: suffix: da_3
369: requires: !complex
370: args: -dm_type da -dim 3 -num_fe_tests 2
372: test:
373: suffix: ds
374: requires: !complex triangle
375: 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
377: test:
378: suffix: ds_simplex_0
379: requires: !complex triangle
380: args: -dm_type plex -dm_plex_dim 2 -dm_plex_box_faces 3,3 -num_fe_tests 2 -petscspace_degree 0
382: test:
383: suffix: ds_simplex_1
384: requires: !complex triangle
385: args: -dm_type plex -dm_plex_dim 2 -dm_plex_box_faces 3,3 -num_fe_tests 2 -petscspace_degree 1
387: test:
388: suffix: ds_simplex_2
389: requires: !complex triangle
390: args: -dm_type plex -dm_plex_dim 2 -dm_plex_box_faces 3,3 -num_fe_tests 2 -petscspace_degree 2
392: test:
393: suffix: ds_tensor_2_0
394: requires: !complex
395: 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
397: test:
398: suffix: ds_tensor_2_1
399: requires: !complex
400: 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
402: test:
403: suffix: ds_tensor_2_2
404: requires: !complex
405: 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
407: test:
408: suffix: ds_tensor_3_0
409: requires: !complex
410: 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
412: test:
413: suffix: ds_tensor_3_1
414: requires: !complex
415: 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
417: test:
418: suffix: ds_tensor_3_2
419: requires: !complex
420: 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
422: test:
423: suffix: shell
424: requires: !complex triangle
425: 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
427: TEST*/