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++) u[i] = (i + 1) * r2;
149: return 0;
150: }
152: static PetscErrorCode TestShellEvaluate(DMField field, Vec points, PetscDataType type, void *B, void *D, void *H)
153: {
154: Vec ctxVec = NULL;
155: const PetscScalar *mult;
156: PetscInt dim;
157: const PetscScalar *x;
158: PetscInt Nc, n, i, j, k, l;
160: DMFieldGetNumComponents(field, &Nc);
161: DMFieldShellGetContext(field, &ctxVec);
162: VecGetBlockSize(points, &dim);
163: VecGetLocalSize(points, &n);
164: n /= Nc;
165: VecGetArrayRead(ctxVec, &mult);
166: VecGetArrayRead(points, &x);
167: for (i = 0; i < n; i++) {
168: PetscReal r2 = 0.;
170: for (j = 0; j < dim; j++) r2 += PetscSqr(PetscRealPart(x[i * dim + j]));
171: for (j = 0; j < Nc; j++) {
172: PetscReal m = PetscRealPart(mult[j]);
173: if (B) {
174: if (type == PETSC_SCALAR) {
175: ((PetscScalar *)B)[i * Nc + j] = m * r2;
176: } else {
177: ((PetscReal *)B)[i * Nc + j] = m * r2;
178: }
179: }
180: if (D) {
181: if (type == PETSC_SCALAR) {
182: for (k = 0; k < dim; k++) ((PetscScalar *)D)[(i * Nc + j) * dim + k] = 2. * m * x[i * dim + k];
183: } else {
184: for (k = 0; k < dim; k++) ((PetscReal *)D)[(i * Nc + j) * dim + k] = 2. * m * PetscRealPart(x[i * dim + k]);
185: }
186: }
187: if (H) {
188: if (type == PETSC_SCALAR) {
189: for (k = 0; k < dim; k++)
190: for (l = 0; l < dim; l++) ((PetscScalar *)H)[((i * Nc + j) * dim + k) * dim + l] = (k == l) ? 2. * m : 0.;
191: } else {
192: for (k = 0; k < dim; k++)
193: 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;
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) PetscDTStroudConicalQuadrature(dim, 1, pointsPerEdge, -1.0, 1.0, &quad);
269: else PetscDTGaussTensorQuadrature(dim, 1, pointsPerEdge, -1.0, 1.0, &quad);
270: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
271: if (testShell) {
272: Vec ctxVec;
273: PetscInt i;
274: PetscScalar *array;
276: VecCreateSeq(PETSC_COMM_SELF, nc, &ctxVec);
277: VecSetUp(ctxVec);
278: VecGetArray(ctxVec, &array);
279: for (i = 0; i < nc; i++) array[i] = i + 1.;
280: VecRestoreArray(ctxVec, &array);
281: DMFieldCreateShell(dm, nc, DMFIELD_VERTEX, (void *)ctxVec, &field);
282: DMFieldShellSetEvaluate(field, TestShellEvaluate);
283: DMFieldShellSetDestroy(field, TestShellDestroy);
284: } else {
285: PetscFECreateDefault(PETSC_COMM_SELF, dim, nc, simplex, NULL, PETSC_DEFAULT, &fe);
286: PetscFESetName(fe, "MyPetscFE");
287: DMSetField(dm, 0, NULL, (PetscObject)fe);
288: PetscFEDestroy(&fe);
289: DMCreateDS(dm);
290: DMCreateLocalVector(dm, &fieldvec);
291: {
292: PetscErrorCode (*func[1])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
293: void *ctxs[1];
295: func[0] = radiusSquared;
296: ctxs[0] = NULL;
298: DMProjectFunctionLocal(dm, 0.0, func, ctxs, INSERT_ALL_VALUES, fieldvec);
299: }
300: DMFieldCreateDS(dm, 0, fieldvec, &field);
301: VecDestroy(&fieldvec);
302: }
303: } else if (isda) {
304: PetscInt i;
305: PetscScalar *cv;
307: switch (dim) {
308: case 1:
309: DMDACreate1d(comm, DM_BOUNDARY_NONE, 3, 1, 1, NULL, &dm);
310: break;
311: case 2:
312: DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, 3, 3, PETSC_DETERMINE, PETSC_DETERMINE, 1, 1, NULL, NULL, &dm);
313: break;
314: default:
315: 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);
316: break;
317: }
318: DMSetUp(dm);
319: DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);
320: PetscMalloc1(nc * (1 << dim), &cv);
321: for (i = 0; i < nc * (1 << dim); i++) {
322: PetscReal rv;
324: PetscRandomGetValueReal(rand, &rv);
325: cv[i] = rv;
326: }
327: DMFieldCreateDA(dm, nc, cv, &field);
328: PetscFree(cv);
329: PetscDTGaussTensorQuadrature(dim, 1, pointsPerEdge, -1.0, 1.0, &quad);
330: } else SETERRQ(comm, PETSC_ERR_SUP, "This test does not run for DM type %s", type);
332: PetscObjectSetName((PetscObject)dm, "mesh");
333: DMViewFromOptions(dm, NULL, "-dm_view");
334: DMDestroy(&dm);
335: PetscObjectSetName((PetscObject)field, "field");
336: PetscObjectViewFromOptions((PetscObject)field, NULL, "-dmfield_view");
337: if (numPoint) TestEvaluate(field, numPoint, rand);
338: if (numFE) TestEvaluateFE(field, numFE, cStart, cEnd, quad, rand);
339: if (numFV) TestEvaluateFV(field, numFV, cStart, cEnd, rand);
340: DMFieldDestroy(&field);
341: PetscQuadratureDestroy(&quad);
342: PetscRandomDestroy(&rand);
343: PetscFinalize();
344: return 0;
345: }
347: /*TEST
349: test:
350: suffix: da
351: requires: !complex
352: args: -dm_type da -dim 2 -num_components 2 -num_point_tests 2 -num_fe_tests 2 -num_fv_tests 2 -dmfield_view
354: test:
355: suffix: da_1
356: requires: !complex
357: args: -dm_type da -dim 1 -num_fe_tests 2
359: test:
360: suffix: da_2
361: requires: !complex
362: args: -dm_type da -dim 2 -num_fe_tests 2
364: test:
365: suffix: da_3
366: requires: !complex
367: args: -dm_type da -dim 3 -num_fe_tests 2
369: test:
370: suffix: ds
371: requires: !complex triangle
372: 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
374: test:
375: suffix: ds_simplex_0
376: requires: !complex triangle
377: args: -dm_type plex -dm_plex_dim 2 -dm_plex_box_faces 3,3 -num_fe_tests 2 -petscspace_degree 0
379: test:
380: suffix: ds_simplex_1
381: requires: !complex triangle
382: args: -dm_type plex -dm_plex_dim 2 -dm_plex_box_faces 3,3 -num_fe_tests 2 -petscspace_degree 1
384: test:
385: suffix: ds_simplex_2
386: requires: !complex triangle
387: args: -dm_type plex -dm_plex_dim 2 -dm_plex_box_faces 3,3 -num_fe_tests 2 -petscspace_degree 2
389: test:
390: suffix: ds_tensor_2_0
391: requires: !complex
392: 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
394: test:
395: suffix: ds_tensor_2_1
396: requires: !complex
397: 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
399: test:
400: suffix: ds_tensor_2_2
401: requires: !complex
402: 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
404: test:
405: suffix: ds_tensor_3_0
406: requires: !complex
407: 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
409: test:
410: suffix: ds_tensor_3_1
411: requires: !complex
412: 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
414: test:
415: suffix: ds_tensor_3_2
416: requires: !complex
417: 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
419: test:
420: suffix: shell
421: requires: !complex triangle
422: 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
424: TEST*/