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: PetscFunctionBegin;
10: PetscCall(PetscViewerASCIIPrintf(viewer, "B:\n"));
11: PetscCall(PetscScalarView(N, B, viewer));
12: PetscCall(PetscViewerASCIIPrintf(viewer, "D:\n"));
13: PetscCall(PetscScalarView(N * dim, D, viewer));
14: PetscCall(PetscViewerASCIIPrintf(viewer, "H:\n"));
15: PetscCall(PetscScalarView(N * dim * dim, H, viewer));
17: PetscCall(PetscViewerASCIIPrintf(viewer, "rB:\n"));
18: PetscCall(PetscRealView(N, rB, viewer));
19: PetscCall(PetscViewerASCIIPrintf(viewer, "rD:\n"));
20: PetscCall(PetscRealView(N * dim, rD, viewer));
21: PetscCall(PetscViewerASCIIPrintf(viewer, "rH:\n"));
22: PetscCall(PetscRealView(N * dim * dim, rH, viewer));
23: PetscFunctionReturn(PETSC_SUCCESS);
24: }
26: static PetscErrorCode TestEvaluate(DMField field, PetscInt n, PetscRandom rand)
27: {
28: DM dm;
29: PetscInt dim, i, nc;
30: PetscScalar *B, *D, *H;
31: PetscReal *rB, *rD, *rH;
32: Vec points;
33: PetscScalar *array;
34: PetscViewer viewer;
35: MPI_Comm comm;
37: PetscFunctionBegin;
38: comm = PetscObjectComm((PetscObject)field);
39: PetscCall(DMFieldGetNumComponents(field, &nc));
40: PetscCall(DMFieldGetDM(field, &dm));
41: PetscCall(DMGetDimension(dm, &dim));
42: PetscCall(VecCreateFromOptions(PetscObjectComm((PetscObject)field), NULL, 1, n * dim, PETSC_DETERMINE, &points));
43: PetscCall(VecSetBlockSize(points, dim));
44: PetscCall(VecGetArray(points, &array));
45: for (i = 0; i < n * dim; i++) PetscCall(PetscRandomGetValue(rand, &array[i]));
46: PetscCall(VecRestoreArray(points, &array));
47: PetscCall(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));
48: PetscCall(DMFieldEvaluate(field, points, PETSC_SCALAR, B, D, H));
49: PetscCall(DMFieldEvaluate(field, points, PETSC_REAL, rB, rD, rH));
50: viewer = PETSC_VIEWER_STDOUT_(comm);
52: PetscCall(PetscObjectSetName((PetscObject)points, "Test Points"));
53: PetscCall(VecView(points, viewer));
54: PetscCall(ViewResults(viewer, n * nc, dim, B, D, H, rB, rD, rH));
56: PetscCall(PetscFree6(B, rB, D, rD, H, rH));
57: PetscCall(VecDestroy(&points));
58: PetscFunctionReturn(PETSC_SUCCESS);
59: }
61: static PetscErrorCode TestEvaluateFE(DMField field, PetscInt n, PetscInt cStart, PetscInt cEnd, PetscQuadrature quad, PetscRandom rand)
62: {
63: DM dm;
64: PetscInt dim, i, nc, nq;
65: PetscInt N;
66: PetscScalar *B, *D, *H;
67: PetscReal *rB, *rD, *rH;
68: PetscInt *cells;
69: IS cellIS;
70: PetscViewer viewer;
71: MPI_Comm comm;
73: PetscFunctionBegin;
74: comm = PetscObjectComm((PetscObject)field);
75: PetscCall(DMFieldGetNumComponents(field, &nc));
76: PetscCall(DMFieldGetDM(field, &dm));
77: PetscCall(DMGetDimension(dm, &dim));
78: PetscCall(PetscRandomSetInterval(rand, (PetscScalar)cStart, (PetscScalar)cEnd));
79: PetscCall(PetscMalloc1(n, &cells));
80: for (i = 0; i < n; i++) {
81: PetscReal rc;
83: PetscCall(PetscRandomGetValueReal(rand, &rc));
84: cells[i] = PetscFloorReal(rc);
85: }
86: PetscCall(ISCreateGeneral(comm, n, cells, PETSC_OWN_POINTER, &cellIS));
87: PetscCall(PetscObjectSetName((PetscObject)cellIS, "FE Test Cells"));
88: PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &nq, NULL, NULL));
89: N = n * nq * nc;
90: PetscCall(PetscMalloc6(N, &B, N, &rB, N * dim, &D, N * dim, &rD, N * dim * dim, &H, N * dim * dim, &rH));
91: PetscCall(DMFieldEvaluateFE(field, cellIS, quad, PETSC_SCALAR, B, D, H));
92: PetscCall(DMFieldEvaluateFE(field, cellIS, quad, PETSC_REAL, rB, rD, rH));
93: viewer = PETSC_VIEWER_STDOUT_(comm);
95: PetscCall(PetscObjectSetName((PetscObject)quad, "Test quadrature"));
96: PetscCall(PetscQuadratureView(quad, viewer));
97: PetscCall(ISView(cellIS, viewer));
98: PetscCall(ViewResults(viewer, N, dim, B, D, H, rB, rD, rH));
100: PetscCall(PetscFree6(B, rB, D, rD, H, rH));
101: PetscCall(ISDestroy(&cellIS));
102: PetscFunctionReturn(PETSC_SUCCESS);
103: }
105: static PetscErrorCode TestEvaluateFV(DMField field, PetscInt n, PetscInt cStart, PetscInt cEnd, PetscRandom rand)
106: {
107: DM dm;
108: PetscInt dim, i, nc;
109: PetscInt N;
110: PetscScalar *B, *D, *H;
111: PetscReal *rB, *rD, *rH;
112: PetscInt *cells;
113: IS cellIS;
114: PetscViewer viewer;
115: MPI_Comm comm;
117: PetscFunctionBegin;
118: comm = PetscObjectComm((PetscObject)field);
119: PetscCall(DMFieldGetNumComponents(field, &nc));
120: PetscCall(DMFieldGetDM(field, &dm));
121: PetscCall(DMGetDimension(dm, &dim));
122: PetscCall(PetscRandomSetInterval(rand, (PetscScalar)cStart, (PetscScalar)cEnd));
123: PetscCall(PetscMalloc1(n, &cells));
124: for (i = 0; i < n; i++) {
125: PetscReal rc;
127: PetscCall(PetscRandomGetValueReal(rand, &rc));
128: cells[i] = PetscFloorReal(rc);
129: }
130: PetscCall(ISCreateGeneral(comm, n, cells, PETSC_OWN_POINTER, &cellIS));
131: PetscCall(PetscObjectSetName((PetscObject)cellIS, "FV Test Cells"));
132: N = n * nc;
133: PetscCall(PetscMalloc6(N, &B, N, &rB, N * dim, &D, N * dim, &rD, N * dim * dim, &H, N * dim * dim, &rH));
134: PetscCall(DMFieldEvaluateFV(field, cellIS, PETSC_SCALAR, B, D, H));
135: PetscCall(DMFieldEvaluateFV(field, cellIS, PETSC_REAL, rB, rD, rH));
136: viewer = PETSC_VIEWER_STDOUT_(comm);
138: PetscCall(ISView(cellIS, viewer));
139: PetscCall(ViewResults(viewer, N, dim, B, D, H, rB, rD, rH));
141: PetscCall(PetscFree6(B, rB, D, rD, H, rH));
142: PetscCall(ISDestroy(&cellIS));
143: PetscFunctionReturn(PETSC_SUCCESS);
144: }
146: static PetscErrorCode radiusSquared(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx)
147: {
148: PetscInt i;
149: PetscReal r2 = 0.;
151: PetscFunctionBegin;
152: for (i = 0; i < dim; i++) r2 += PetscSqr(x[i]);
153: for (i = 0; i < Nf; i++) u[i] = (i + 1) * r2;
154: PetscFunctionReturn(PETSC_SUCCESS);
155: }
157: static PetscErrorCode TestShellEvaluate(DMField field, Vec points, PetscDataType type, void *B, void *D, void *H)
158: {
159: Vec ctxVec = NULL;
160: const PetscScalar *mult;
161: PetscInt dim;
162: const PetscScalar *x;
163: PetscInt Nc, n, i, j, k, l;
165: PetscFunctionBegin;
166: PetscCall(DMFieldGetNumComponents(field, &Nc));
167: PetscCall(DMFieldShellGetContext(field, &ctxVec));
168: PetscCall(VecGetBlockSize(points, &dim));
169: PetscCall(VecGetLocalSize(points, &n));
170: n /= Nc;
171: PetscCall(VecGetArrayRead(ctxVec, &mult));
172: PetscCall(VecGetArrayRead(points, &x));
173: for (i = 0; i < n; i++) {
174: PetscReal r2 = 0.;
176: for (j = 0; j < dim; j++) r2 += PetscSqr(PetscRealPart(x[i * dim + j]));
177: for (j = 0; j < Nc; j++) {
178: PetscReal m = PetscRealPart(mult[j]);
179: if (B) {
180: if (type == PETSC_SCALAR) {
181: ((PetscScalar *)B)[i * Nc + j] = m * r2;
182: } else {
183: ((PetscReal *)B)[i * Nc + j] = m * r2;
184: }
185: }
186: if (D) {
187: if (type == PETSC_SCALAR) {
188: for (k = 0; k < dim; k++) ((PetscScalar *)D)[(i * Nc + j) * dim + k] = 2. * m * x[i * dim + k];
189: } else {
190: for (k = 0; k < dim; k++) ((PetscReal *)D)[(i * Nc + j) * dim + k] = 2. * m * PetscRealPart(x[i * dim + k]);
191: }
192: }
193: if (H) {
194: if (type == PETSC_SCALAR) {
195: for (k = 0; k < dim; k++)
196: for (l = 0; l < dim; l++) ((PetscScalar *)H)[((i * Nc + j) * dim + k) * dim + l] = (k == l) ? 2. * m : 0.;
197: } else {
198: for (k = 0; k < dim; k++)
199: for (l = 0; l < dim; l++) ((PetscReal *)H)[((i * Nc + j) * dim + k) * dim + l] = (k == l) ? 2. * m : 0.;
200: }
201: }
202: }
203: }
204: PetscCall(VecRestoreArrayRead(points, &x));
205: PetscCall(VecRestoreArrayRead(ctxVec, &mult));
206: PetscFunctionReturn(PETSC_SUCCESS);
207: }
209: static PetscErrorCode TestShellDestroy(DMField field)
210: {
211: Vec ctxVec = NULL;
213: PetscFunctionBegin;
214: PetscCall(DMFieldShellGetContext(field, &ctxVec));
215: PetscCall(VecDestroy(&ctxVec));
216: PetscFunctionReturn(PETSC_SUCCESS);
217: }
219: int main(int argc, char **argv)
220: {
221: DM dm = NULL;
222: MPI_Comm comm;
223: char type[256] = DMPLEX;
224: PetscBool isda, isplex;
225: PetscInt dim = 2;
226: DMField field = NULL;
227: PetscInt nc = 1;
228: PetscInt cStart = -1, cEnd = -1;
229: PetscRandom rand;
230: PetscQuadrature quad = NULL;
231: PetscInt pointsPerEdge = 2;
232: PetscInt numPoint = 0, numFE = 0, numFV = 0;
233: PetscBool testShell = PETSC_FALSE;
235: PetscFunctionBeginUser;
236: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
237: comm = PETSC_COMM_WORLD;
238: PetscOptionsBegin(comm, "", "DMField Tutorial Options", "DM");
239: PetscCall(PetscOptionsFList("-dm_type", "DM implementation on which to define field", "ex1.c", DMList, type, type, 256, NULL));
240: PetscCall(PetscOptionsRangeInt("-dim", "DM intrinsic dimension", "ex1.c", dim, &dim, NULL, 1, 3));
241: PetscCall(PetscOptionsBoundedInt("-num_components", "Number of components in field", "ex1.c", nc, &nc, NULL, 1));
242: PetscCall(PetscOptionsBoundedInt("-num_quad_points", "Number of quadrature points per dimension", "ex1.c", pointsPerEdge, &pointsPerEdge, NULL, 1));
243: PetscCall(PetscOptionsBoundedInt("-num_point_tests", "Number of test points for DMFieldEvaluate()", "ex1.c", numPoint, &numPoint, NULL, 0));
244: PetscCall(PetscOptionsBoundedInt("-num_fe_tests", "Number of test cells for DMFieldEvaluateFE()", "ex1.c", numFE, &numFE, NULL, 0));
245: PetscCall(PetscOptionsBoundedInt("-num_fv_tests", "Number of test cells for DMFieldEvaluateFV()", "ex1.c", numFV, &numFV, NULL, 0));
246: PetscCall(PetscOptionsBool("-test_shell", "Test the DMFIELDSHELL implementation of DMField", "ex1.c", testShell, &testShell, NULL));
247: PetscOptionsEnd();
249: PetscCheck(dim <= 3, comm, PETSC_ERR_ARG_OUTOFRANGE, "This examples works for dim <= 3, not %" PetscInt_FMT, dim);
250: PetscCall(PetscStrncmp(type, DMPLEX, 256, &isplex));
251: PetscCall(PetscStrncmp(type, DMDA, 256, &isda));
253: PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rand));
254: PetscCall(PetscRandomSetFromOptions(rand));
255: if (isplex) {
256: PetscInt overlap = 0;
257: Vec fieldvec;
258: PetscInt cells[3] = {3, 3, 3};
259: PetscBool simplex;
260: PetscFE fe;
262: PetscOptionsBegin(comm, "", "DMField DMPlex Options", "DM");
263: PetscCall(PetscOptionsBoundedInt("-overlap", "DMPlex parallel overlap", "ex1.c", overlap, &overlap, NULL, 0));
264: PetscOptionsEnd();
265: if (0) {
266: PetscCall(DMPlexCreateBoxMesh(comm, 2, PETSC_TRUE, cells, NULL, NULL, NULL, PETSC_TRUE, &dm));
267: } else {
268: PetscCall(DMCreate(comm, &dm));
269: PetscCall(DMSetType(dm, DMPLEX));
270: PetscCall(DMSetFromOptions(dm));
271: CHKMEMQ;
272: }
273: PetscCall(DMGetDimension(dm, &dim));
274: PetscCall(DMPlexIsSimplex(dm, &simplex));
275: if (simplex) PetscCall(PetscDTStroudConicalQuadrature(dim, 1, pointsPerEdge, -1.0, 1.0, &quad));
276: else PetscCall(PetscDTGaussTensorQuadrature(dim, 1, pointsPerEdge, -1.0, 1.0, &quad));
277: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
278: if (testShell) {
279: Vec ctxVec;
280: PetscInt i;
281: PetscScalar *array;
283: PetscCall(VecCreateSeq(PETSC_COMM_SELF, nc, &ctxVec));
284: PetscCall(VecSetUp(ctxVec));
285: PetscCall(VecGetArray(ctxVec, &array));
286: for (i = 0; i < nc; i++) array[i] = i + 1.;
287: PetscCall(VecRestoreArray(ctxVec, &array));
288: PetscCall(DMFieldCreateShell(dm, nc, DMFIELD_VERTEX, (void *)ctxVec, &field));
289: PetscCall(DMFieldShellSetEvaluate(field, TestShellEvaluate));
290: PetscCall(DMFieldShellSetDestroy(field, TestShellDestroy));
291: } else {
292: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, nc, simplex, NULL, PETSC_DEFAULT, &fe));
293: PetscCall(PetscFESetName(fe, "MyPetscFE"));
294: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
295: PetscCall(PetscFEDestroy(&fe));
296: PetscCall(DMCreateDS(dm));
297: PetscCall(DMCreateLocalVector(dm, &fieldvec));
298: {
299: PetscErrorCode (*func[1])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
300: void *ctxs[1];
302: func[0] = radiusSquared;
303: ctxs[0] = NULL;
305: PetscCall(DMProjectFunctionLocal(dm, 0.0, func, ctxs, INSERT_ALL_VALUES, fieldvec));
306: }
307: PetscCall(DMFieldCreateDS(dm, 0, fieldvec, &field));
308: PetscCall(VecDestroy(&fieldvec));
309: }
310: } else if (isda) {
311: PetscInt i;
312: PetscScalar *cv;
314: switch (dim) {
315: case 1:
316: PetscCall(DMDACreate1d(comm, DM_BOUNDARY_NONE, 3, 1, 1, NULL, &dm));
317: break;
318: case 2:
319: PetscCall(DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, 3, 3, PETSC_DETERMINE, PETSC_DETERMINE, 1, 1, NULL, NULL, &dm));
320: break;
321: default:
322: PetscCall(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));
323: break;
324: }
325: PetscCall(DMSetUp(dm));
326: PetscCall(DMDAGetHeightStratum(dm, 0, &cStart, &cEnd));
327: PetscCall(PetscMalloc1(nc * (1 << dim), &cv));
328: for (i = 0; i < nc * (1 << dim); i++) {
329: PetscReal rv;
331: PetscCall(PetscRandomGetValueReal(rand, &rv));
332: cv[i] = rv;
333: }
334: PetscCall(DMFieldCreateDA(dm, nc, cv, &field));
335: PetscCall(PetscFree(cv));
336: PetscCall(PetscDTGaussTensorQuadrature(dim, 1, pointsPerEdge, -1.0, 1.0, &quad));
337: } else SETERRQ(comm, PETSC_ERR_SUP, "This test does not run for DM type %s", type);
339: PetscCall(PetscObjectSetName((PetscObject)dm, "mesh"));
340: PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
341: PetscCall(DMDestroy(&dm));
342: PetscCall(PetscObjectSetName((PetscObject)field, "field"));
343: PetscCall(PetscObjectViewFromOptions((PetscObject)field, NULL, "-dmfield_view"));
344: if (numPoint) PetscCall(TestEvaluate(field, numPoint, rand));
345: if (numFE) PetscCall(TestEvaluateFE(field, numFE, cStart, cEnd, quad, rand));
346: if (numFV) PetscCall(TestEvaluateFV(field, numFV, cStart, cEnd, rand));
347: PetscCall(DMFieldDestroy(&field));
348: PetscCall(PetscQuadratureDestroy(&quad));
349: PetscCall(PetscRandomDestroy(&rand));
350: PetscCall(PetscFinalize());
351: return 0;
352: }
354: /*TEST
356: test:
357: suffix: da
358: requires: !complex
359: args: -dm_type da -dim 2 -num_components 2 -num_point_tests 2 -num_fe_tests 2 -num_fv_tests 2 -dmfield_view
361: test:
362: suffix: da_1
363: requires: !complex
364: args: -dm_type da -dim 1 -num_fe_tests 2
366: test:
367: suffix: da_2
368: requires: !complex
369: args: -dm_type da -dim 2 -num_fe_tests 2
371: test:
372: suffix: da_3
373: requires: !complex
374: args: -dm_type da -dim 3 -num_fe_tests 2
376: test:
377: suffix: ds
378: requires: !complex triangle
379: 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
381: test:
382: suffix: ds_simplex_0
383: requires: !complex triangle
384: args: -dm_type plex -dm_plex_dim 2 -dm_plex_box_faces 3,3 -num_fe_tests 2 -petscspace_degree 0
386: test:
387: suffix: ds_simplex_1
388: requires: !complex triangle
389: args: -dm_type plex -dm_plex_dim 2 -dm_plex_box_faces 3,3 -num_fe_tests 2 -petscspace_degree 1
391: test:
392: suffix: ds_simplex_2
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 2
396: test:
397: suffix: ds_tensor_2_0
398: requires: !complex
399: 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
401: test:
402: suffix: ds_tensor_2_1
403: requires: !complex
404: 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
406: test:
407: suffix: ds_tensor_2_2
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 2 -dm_plex_simplex 0
411: test:
412: suffix: ds_tensor_3_0
413: requires: !complex
414: 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
416: test:
417: suffix: ds_tensor_3_1
418: requires: !complex
419: 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
421: test:
422: suffix: ds_tensor_3_2
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 2 -dm_plex_simplex 0
426: test:
427: suffix: shell
428: requires: !complex triangle
429: 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
431: TEST*/