1 static char help[] = "Demonstration of creating and viewing DMFields objects.\n\n";
2
3 #include <petscdmfield.h>
4 #include <petscdmplex.h>
5 #include <petscdmda.h>
6
ViewResults(PetscViewer viewer,PetscInt N,PetscInt dim,PetscScalar * B,PetscScalar * D,PetscScalar * H,PetscReal * rB,PetscReal * rD,PetscReal * rH)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));
16
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 }
25
TestEvaluate(DMField field,PetscInt n,PetscRandom rand)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;
36
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);
51
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));
55
56 PetscCall(PetscFree6(B, rB, D, rD, H, rH));
57 PetscCall(VecDestroy(&points));
58 PetscFunctionReturn(PETSC_SUCCESS);
59 }
60
TestEvaluateFE(DMField field,PetscInt n,PetscInt cStart,PetscInt cEnd,PetscQuadrature quad,PetscRandom rand)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;
72
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;
82
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);
94
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));
99
100 PetscCall(PetscFree6(B, rB, D, rD, H, rH));
101 PetscCall(ISDestroy(&cellIS));
102 PetscFunctionReturn(PETSC_SUCCESS);
103 }
104
TestEvaluateFV(DMField field,PetscInt n,PetscInt cStart,PetscInt cEnd,PetscRandom rand)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;
116
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;
126
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);
137
138 PetscCall(ISView(cellIS, viewer));
139 PetscCall(ViewResults(viewer, N, dim, B, D, H, rB, rD, rH));
140
141 PetscCall(PetscFree6(B, rB, D, rD, H, rH));
142 PetscCall(ISDestroy(&cellIS));
143 PetscFunctionReturn(PETSC_SUCCESS);
144 }
145
radiusSquared(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nf,PetscScalar u[],PetscCtx ctx)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.;
150
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 }
156
TestShellEvaluate(DMField field,Vec points,PetscDataType type,void * B,void * D,void * H)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;
164
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.;
175
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 }
208
TestShellDestroy(DMField field)209 static PetscErrorCode TestShellDestroy(DMField field)
210 {
211 Vec ctxVec = NULL;
212
213 PetscFunctionBegin;
214 PetscCall(DMFieldShellGetContext(field, &ctxVec));
215 PetscCall(VecDestroy(&ctxVec));
216 PetscFunctionReturn(PETSC_SUCCESS);
217 }
218
main(int argc,char ** argv)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;
234
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();
248
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));
252
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;
261
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 PetscInt i;
281 PetscScalar *array;
282
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 PetscCtx ctxs[1];
301
302 func[0] = radiusSquared;
303 ctxs[0] = NULL;
304
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;
313
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;
330
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);
338
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 }
353
354 /*TEST
355
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
360
361 test:
362 suffix: da_1
363 requires: !complex
364 args: -dm_type da -dim 1 -num_fe_tests 2
365
366 test:
367 suffix: da_2
368 requires: !complex
369 args: -dm_type da -dim 2 -num_fe_tests 2
370
371 test:
372 suffix: da_3
373 requires: !complex
374 args: -dm_type da -dim 3 -num_fe_tests 2
375
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
380
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
385
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
390
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
395
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
400
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
405
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
410
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
415
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
420
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
425
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
430
431 TEST*/
432