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[], PetscCtx 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, 0, 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: PetscScalar *array;
282: PetscCall(VecCreateSeq(PETSC_COMM_SELF, nc, &ctxVec));
283: PetscCall(VecSetUp(ctxVec));
284: PetscCall(VecGetArray(ctxVec, &array));
285: for (PetscInt i = 0; i < nc; i++) array[i] = i + 1.;
286: PetscCall(VecRestoreArray(ctxVec, &array));
287: PetscCall(DMFieldCreateShell(dm, nc, DMFIELD_VERTEX, (void *)ctxVec, &field));
288: PetscCall(DMFieldShellSetEvaluate(field, TestShellEvaluate));
289: PetscCall(DMFieldShellSetDestroy(field, TestShellDestroy));
290: } else {
291: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, nc, simplex, NULL, PETSC_DEFAULT, &fe));
292: PetscCall(PetscFESetName(fe, "MyPetscFE"));
293: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
294: PetscCall(PetscFEDestroy(&fe));
295: PetscCall(DMCreateDS(dm));
296: PetscCall(DMCreateLocalVector(dm, &fieldvec));
297: {
298: PetscErrorCode (*func[1])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
299: PetscCtx ctxs[1];
301: func[0] = radiusSquared;
302: ctxs[0] = NULL;
304: PetscCall(DMProjectFunctionLocal(dm, 0.0, func, ctxs, INSERT_ALL_VALUES, fieldvec));
305: }
306: PetscCall(DMFieldCreateDS(dm, 0, fieldvec, &field));
307: PetscCall(VecDestroy(&fieldvec));
308: }
309: } else if (isda) {
310: PetscScalar *cv;
312: switch (dim) {
313: case 1:
314: PetscCall(DMDACreate1d(comm, DM_BOUNDARY_NONE, 3, 1, 1, NULL, &dm));
315: break;
316: case 2:
317: PetscCall(DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, 3, 3, PETSC_DETERMINE, PETSC_DETERMINE, 1, 1, NULL, NULL, &dm));
318: break;
319: default:
320: 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));
321: break;
322: }
323: PetscCall(DMSetUp(dm));
324: PetscCall(DMDAGetHeightStratum(dm, 0, &cStart, &cEnd));
325: PetscCall(PetscMalloc1(nc * (1 << dim), &cv));
326: for (PetscInt i = 0; i < nc * (1 << dim); i++) {
327: PetscReal rv;
329: PetscCall(PetscRandomGetValueReal(rand, &rv));
330: cv[i] = rv;
331: }
332: PetscCall(DMFieldCreateDA(dm, nc, cv, &field));
333: PetscCall(PetscFree(cv));
334: PetscCall(PetscDTGaussTensorQuadrature(dim, 1, pointsPerEdge, -1.0, 1.0, &quad));
335: } else SETERRQ(comm, PETSC_ERR_SUP, "This test does not run for DM type %s", type);
337: PetscCall(PetscObjectSetName((PetscObject)dm, "mesh"));
338: PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
339: PetscCall(DMDestroy(&dm));
340: PetscCall(PetscObjectSetName((PetscObject)field, "field"));
341: PetscCall(PetscObjectViewFromOptions((PetscObject)field, NULL, "-dmfield_view"));
342: if (numPoint) PetscCall(TestEvaluate(field, numPoint, rand));
343: if (numFE) PetscCall(TestEvaluateFE(field, numFE, cStart, cEnd, quad, rand));
344: if (numFV) PetscCall(TestEvaluateFV(field, numFV, cStart, cEnd, rand));
345: PetscCall(DMFieldDestroy(&field));
346: PetscCall(PetscQuadratureDestroy(&quad));
347: PetscCall(PetscRandomDestroy(&rand));
348: PetscCall(PetscFinalize());
349: return 0;
350: }
352: /*TEST
354: test:
355: suffix: da
356: requires: !complex
357: args: -dm_type da -dim 2 -num_components 2 -num_point_tests 2 -num_fe_tests 2 -num_fv_tests 2 -dmfield_view
359: test:
360: suffix: da_1
361: requires: !complex
362: args: -dm_type da -dim 1 -num_fe_tests 2
364: test:
365: suffix: da_2
366: requires: !complex
367: args: -dm_type da -dim 2 -num_fe_tests 2
369: test:
370: suffix: da_3
371: requires: !complex
372: args: -dm_type da -dim 3 -num_fe_tests 2
374: test:
375: suffix: ds
376: requires: !complex triangle
377: 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
379: test:
380: suffix: ds_simplex_0
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 0
384: test:
385: suffix: ds_simplex_1
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 1
389: test:
390: suffix: ds_simplex_2
391: requires: !complex triangle
392: args: -dm_type plex -dm_plex_dim 2 -dm_plex_box_faces 3,3 -num_fe_tests 2 -petscspace_degree 2
394: test:
395: suffix: ds_tensor_2_0
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 0 -dm_plex_simplex 0
399: test:
400: suffix: ds_tensor_2_1
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 1 -dm_plex_simplex 0
404: test:
405: suffix: ds_tensor_2_2
406: requires: !complex
407: 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
409: test:
410: suffix: ds_tensor_3_0
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 0 -dm_plex_simplex 0
414: test:
415: suffix: ds_tensor_3_1
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 1 -dm_plex_simplex 0
419: test:
420: suffix: ds_tensor_3_2
421: requires: !complex
422: 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
424: test:
425: suffix: shell
426: requires: !complex triangle
427: 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
429: TEST*/