Actual source code: ex33.c
1: static char help[] = "Tests for high order geometry\n\n";
3: #include <petscdmplex.h>
4: #include <petscds.h>
6: typedef enum {
7: TRANSFORM_NONE,
8: TRANSFORM_SHEAR,
9: TRANSFORM_FLARE,
10: TRANSFORM_ANNULUS,
11: TRANSFORM_SHELL
12: } Transform;
13: const char *const TransformTypes[] = {"none", "shear", "flare", "annulus", "shell", "Mesh Transform", "TRANSFORM_", NULL};
15: typedef struct {
16: PetscBool coordSpace; /* Flag to create coordinate space */
17: Transform meshTransform; /* Transform for initial box mesh */
18: PetscReal *transformDataReal; /* Parameters for mesh transform */
19: PetscScalar *transformData; /* Parameters for mesh transform */
20: PetscReal volume; /* Analytical volume of the mesh */
21: PetscReal tol; /* Tolerance for volume check */
22: } AppCtx;
24: PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
25: {
26: PetscInt n = 0, i;
28: options->coordSpace = PETSC_TRUE;
29: options->meshTransform = TRANSFORM_NONE;
30: options->transformDataReal = NULL;
31: options->transformData = NULL;
32: options->volume = -1.0;
33: options->tol = PETSC_SMALL;
35: PetscOptionsBegin(comm, "", "Meshing Interpolation Test Options", "DMPLEX");
36: PetscOptionsBool("-coord_space", "Flag to create a coordinate space", "ex33.c", options->coordSpace, &options->coordSpace, NULL);
37: PetscOptionsEnum("-mesh_transform", "Method to transform initial box mesh <none, shear, annulus, shell>", "ex33.c", TransformTypes, (PetscEnum)options->meshTransform, (PetscEnum *)&options->meshTransform, NULL);
38: switch (options->meshTransform) {
39: case TRANSFORM_NONE:
40: break;
41: case TRANSFORM_SHEAR:
42: n = 2;
43: PetscMalloc1(n, &options->transformDataReal);
44: for (i = 0; i < n; ++i) options->transformDataReal[i] = 1.0;
45: PetscOptionsRealArray("-transform_data", "Parameters for mesh transforms", "ex33.c", options->transformDataReal, &n, NULL);
46: break;
47: case TRANSFORM_FLARE:
48: n = 4;
49: PetscMalloc1(n, &options->transformData);
50: for (i = 0; i < n; ++i) options->transformData[i] = 1.0;
51: options->transformData[0] = (PetscScalar)0;
52: PetscOptionsScalarArray("-transform_data", "Parameters for mesh transforms", "ex33.c", options->transformData, &n, NULL);
53: break;
54: case TRANSFORM_ANNULUS:
55: n = 2;
56: PetscMalloc1(n, &options->transformData);
57: options->transformData[0] = 1.0;
58: options->transformData[1] = 2.0;
59: PetscOptionsScalarArray("-transform_data", "Parameters for mesh transforms", "ex33.c", options->transformData, &n, NULL);
60: break;
61: case TRANSFORM_SHELL:
62: n = 2;
63: PetscMalloc1(n, &options->transformData);
64: options->transformData[0] = 1.0;
65: options->transformData[1] = 2.0;
66: PetscOptionsScalarArray("-transform_data", "Parameters for mesh transforms", "ex33.c", options->transformData, &n, NULL);
67: break;
68: default:
69: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Unknown mesh transform %d", options->meshTransform);
70: }
71: PetscOptionsReal("-volume", "The analytical volume of the mesh", "ex33.c", options->volume, &options->volume, NULL);
72: PetscOptionsReal("-tol", "The tolerance for the volume check", "ex33.c", options->tol, &options->tol, NULL);
73: PetscOptionsEnd();
74: return 0;
75: }
77: static void identity(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
78: {
79: const PetscInt Nc = uOff[1] - uOff[0];
80: PetscInt c;
82: for (c = 0; c < Nc; ++c) f0[c] = u[c];
83: }
85: /* Flare applies the transformation, assuming we fix x_f,
87: x_i = x_i * alpha_i x_f
88: */
89: static void f0_flare(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar coords[])
90: {
91: const PetscInt Nc = uOff[1] - uOff[0];
92: const PetscInt cf = (PetscInt)PetscRealPart(constants[0]);
93: PetscInt c;
95: for (c = 0; c < Nc; ++c) coords[c] = u[c] * (c == cf ? 1.0 : constants[c + 1] * u[cf]);
96: }
98: /*
99: We would like to map the unit square to a quarter of the annulus between circles of radius 1 and 2. We start by mapping the straight sections, which
100: will correspond to the top and bottom of our square. So
102: (0,0)--(1,0) ==> (1,0)--(2,0) Just a shift of (1,0)
103: (0,1)--(1,1) ==> (0,1)--(0,2) Switch x and y
105: So it looks like we want to map each layer in y to a ray, so x is the radius and y is the angle:
107: (x, y) ==> (x+1, \pi/2 y) in (r', \theta') space
108: ==> ((x+1) cos(\pi/2 y), (x+1) sin(\pi/2 y)) in (x', y') space
109: */
110: static void f0_annulus(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar xp[])
111: {
112: const PetscReal ri = PetscRealPart(constants[0]);
113: const PetscReal ro = PetscRealPart(constants[1]);
115: xp[0] = (x[0] * (ro - ri) + ri) * PetscCosReal(0.5 * PETSC_PI * x[1]);
116: xp[1] = (x[0] * (ro - ri) + ri) * PetscSinReal(0.5 * PETSC_PI * x[1]);
117: }
119: /*
120: We would like to map the unit cube to a hemisphere of the spherical shell between balls of radius 1 and 2. We want to map the bottom surface onto the
121: lower hemisphere and the upper surface onto the top, letting z be the radius.
123: (x, y) ==> ((z+3)/2, \pi/2 (|x| or |y|), arctan y/x) in (r', \theta', \phi') space
124: ==> ((z+3)/2 \cos(\theta') cos(\phi'), (z+3)/2 \cos(\theta') sin(\phi'), (z+3)/2 sin(\theta')) in (x', y', z') space
125: */
126: static void f0_shell(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar xp[])
127: {
128: const PetscReal pi4 = PETSC_PI / 4.0;
129: const PetscReal ri = PetscRealPart(constants[0]);
130: const PetscReal ro = PetscRealPart(constants[1]);
131: const PetscReal rp = (x[2] + 1) * 0.5 * (ro - ri) + ri;
132: const PetscReal phip = PetscAtan2Real(x[1], x[0]);
133: const PetscReal thetap = 0.5 * PETSC_PI * (1.0 - ((((phip <= pi4) && (phip >= -pi4)) || ((phip >= 3.0 * pi4) || (phip <= -3.0 * pi4))) ? PetscAbsReal(x[0]) : PetscAbsReal(x[1])));
135: xp[0] = rp * PetscCosReal(thetap) * PetscCosReal(phip);
136: xp[1] = rp * PetscCosReal(thetap) * PetscSinReal(phip);
137: xp[2] = rp * PetscSinReal(thetap);
138: }
140: static PetscErrorCode DMCreateCoordinateDisc(DM dm)
141: {
142: DM cdm;
143: PetscFE fe;
144: DMPolytopeType ct;
145: PetscInt dim, dE, cStart;
146: PetscBool simplex;
148: DMGetCoordinateDM(dm, &cdm);
149: DMGetDimension(dm, &dim);
150: DMGetCoordinateDim(dm, &dE);
151: DMPlexGetHeightStratum(cdm, 0, &cStart, NULL);
152: DMPlexGetCellType(dm, cStart, &ct);
153: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
154: PetscFECreateDefault(PETSC_COMM_SELF, dim, dE, simplex, "dm_coord_", -1, &fe);
155: DMProjectCoordinates(dm, fe);
156: PetscFEDestroy(&fe);
157: return 0;
158: }
160: PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *ctx, DM *dm)
161: {
162: DM cdm;
163: PetscDS cds;
165: DMCreate(comm, dm);
166: DMSetType(*dm, DMPLEX);
167: DMSetFromOptions(*dm);
169: if (ctx->coordSpace) DMCreateCoordinateDisc(*dm);
170: switch (ctx->meshTransform) {
171: case TRANSFORM_NONE:
172: DMPlexRemapGeometry(*dm, 0.0, identity);
173: break;
174: case TRANSFORM_SHEAR:
175: DMPlexShearGeometry(*dm, DM_X, ctx->transformDataReal);
176: break;
177: case TRANSFORM_FLARE:
178: DMGetCoordinateDM(*dm, &cdm);
179: DMGetDS(cdm, &cds);
180: PetscDSSetConstants(cds, 4, ctx->transformData);
181: DMPlexRemapGeometry(*dm, 0.0, f0_flare);
182: break;
183: case TRANSFORM_ANNULUS:
184: DMGetCoordinateDM(*dm, &cdm);
185: DMGetDS(cdm, &cds);
186: PetscDSSetConstants(cds, 2, ctx->transformData);
187: DMPlexRemapGeometry(*dm, 0.0, f0_annulus);
188: break;
189: case TRANSFORM_SHELL:
190: DMGetCoordinateDM(*dm, &cdm);
191: DMGetDS(cdm, &cds);
192: PetscDSSetConstants(cds, 2, ctx->transformData);
193: DMPlexRemapGeometry(*dm, 0.0, f0_shell);
194: break;
195: default:
196: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Unknown mesh transform %d", ctx->meshTransform);
197: }
198: DMViewFromOptions(*dm, NULL, "-dm_view");
199: return 0;
200: }
202: static void volume(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar vol[])
203: {
204: vol[0] = 1.;
205: }
207: static PetscErrorCode CreateDiscretization(DM dm, AppCtx *ctx)
208: {
209: PetscDS ds;
210: PetscFE fe;
211: DMPolytopeType ct;
212: PetscInt dim, cStart;
213: PetscBool simplex;
216: DMGetDimension(dm, &dim);
217: DMPlexGetHeightStratum(dm, 0, &cStart, NULL);
218: DMPlexGetCellType(dm, cStart, &ct);
219: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
220: PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, NULL, PETSC_DETERMINE, &fe);
221: PetscFESetName(fe, "scalar");
222: DMAddField(dm, NULL, (PetscObject)fe);
223: PetscFEDestroy(&fe);
224: DMCreateDS(dm);
225: DMGetDS(dm, &ds);
226: PetscDSSetObjective(ds, 0, volume);
227: return 0;
228: }
230: static PetscErrorCode CheckVolume(DM dm, AppCtx *ctx)
231: {
232: Vec u;
233: PetscScalar result;
234: PetscReal vol, tol = ctx->tol;
237: DMGetGlobalVector(dm, &u);
238: DMPlexComputeIntegralFEM(dm, u, &result, ctx);
239: vol = PetscRealPart(result);
240: DMRestoreGlobalVector(dm, &u);
241: PetscPrintf(PetscObjectComm((PetscObject)dm), "Volume: %g\n", (double)vol);
242: if (ctx->volume > 0.0 && PetscAbsReal(ctx->volume - vol) > tol) {
243: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Calculated volume %g != %g actual volume (error %g > %g tol)", (double)vol, (double)ctx->volume, (double)PetscAbsReal(ctx->volume - vol), (double)tol);
244: }
245: return 0;
246: }
248: int main(int argc, char **argv)
249: {
250: DM dm;
251: AppCtx user;
254: PetscInitialize(&argc, &argv, NULL, help);
255: ProcessOptions(PETSC_COMM_WORLD, &user);
256: CreateMesh(PETSC_COMM_WORLD, &user, &dm);
257: CreateDiscretization(dm, &user);
258: CheckVolume(dm, &user);
259: DMDestroy(&dm);
260: PetscFree(user.transformDataReal);
261: PetscFree(user.transformData);
262: PetscFinalize();
263: return 0;
264: }
266: /*TEST
268: testset:
269: args: -dm_plex_simplex 0 -dm_plex_box_faces 1,1 -dm_plex_box_lower -1.,-1. -dm_plex_box_upper 1.,1. -volume 4. -dm_coord_space 0
271: test:
272: suffix: square_0
273: args: -dm_coord_petscspace_degree 1
275: test:
276: suffix: square_1
277: args: -dm_coord_petscspace_degree 2
279: test:
280: suffix: square_2
281: args: -dm_refine 1 -dm_coord_petscspace_degree 1
283: test:
284: suffix: square_3
285: args: -dm_refine 1 -dm_coord_petscspace_degree 2
287: testset:
288: args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1 -dm_plex_box_lower -1.,-1.,-1. -dm_plex_box_upper 1.,1.,1. -volume 8. -dm_coord_space 0
290: test:
291: suffix: cube_0
292: args: -dm_coord_petscspace_degree 1
294: test:
295: suffix: cube_1
296: args: -dm_coord_petscspace_degree 2
298: test:
299: suffix: cube_2
300: args: -dm_refine 1 -dm_coord_petscspace_degree 1
302: test:
303: suffix: cube_3
304: args: -dm_refine 1 -dm_coord_petscspace_degree 2
306: testset:
307: args: -dm_plex_simplex 0 -dm_plex_box_faces 1,1 -dm_plex_box_lower -1.,-1. -dm_plex_box_upper 1.,1. -volume 4. -dm_coord_space 0
309: test:
310: suffix: shear_0
311: args: -dm_coord_petscspace_degree 1 -mesh_transform shear -transform_data 3.0
313: test:
314: suffix: shear_1
315: args: -dm_coord_petscspace_degree 2 -mesh_transform shear -transform_data 3.0
317: test:
318: suffix: shear_2
319: args: -dm_refine 1 -dm_coord_petscspace_degree 1 -mesh_transform shear -transform_data 3.0
321: test:
322: suffix: shear_3
323: args: -dm_refine 1 -dm_coord_petscspace_degree 2 -mesh_transform shear -transform_data 3.0
325: testset:
326: args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1 -dm_plex_box_lower -1.,-1.,-1. -dm_plex_box_upper 1.,1.,1. -volume 8. -dm_coord_space 0
328: test:
329: suffix: shear_4
330: args: -dm_coord_petscspace_degree 1 -mesh_transform shear -transform_data 3.0
332: test:
333: suffix: shear_5
334: args: -dm_coord_petscspace_degree 2 -mesh_transform shear -transform_data 3.0
336: test:
337: suffix: shear_6
338: args: -dm_refine 1 -dm_coord_petscspace_degree 1 -mesh_transform shear -transform_data 3.0,4.0
340: test:
341: suffix: shear_7
342: args: -dm_refine 1 -dm_coord_petscspace_degree 2 -mesh_transform shear -transform_data 3.0,4.0
344: testset:
345: args: -dm_coord_space 0 -dm_plex_simplex 0 -dm_plex_box_faces 1,1 -dm_plex_box_lower 1.,-1. -dm_plex_box_upper 3.,1. \
346: -dm_coord_petscspace_degree 1 -mesh_transform flare -volume 8.
348: test:
349: suffix: flare_0
350: args:
352: test:
353: suffix: flare_1
354: args: -dm_refine 2
356: testset:
357: # Area: 3/4 \pi = 2.3562
358: args: -dm_plex_simplex 0 -dm_plex_box_faces 1,1 -mesh_transform annulus -volume 2.35619449019235 -dm_coord_space 0
360: test:
361: # Area: (a+b)/2 h = 3/\sqrt{2} (sqrt{2} - 1/\sqrt{2}) = 3/2
362: suffix: annulus_0
363: requires: double
364: args: -dm_coord_petscspace_degree 1 -volume 1.5
366: test:
367: suffix: annulus_1
368: requires: double
369: args: -dm_refine 3 -dm_coord_petscspace_degree 1 -tol .016
371: test:
372: suffix: annulus_2
373: requires: double
374: args: -dm_refine 3 -dm_coord_petscspace_degree 2 -tol .0038
376: test:
377: suffix: annulus_3
378: requires: double
379: args: -dm_refine 3 -dm_coord_petscspace_degree 3 -tol 2.2e-6
381: test:
382: suffix: annulus_4
383: requires: double
384: args: -dm_refine 2 -dm_coord_petscspace_degree 2 -petscfe_default_quadrature_order 2 -tol .00012
386: test:
387: suffix: annulus_5
388: requires: double
389: args: -dm_refine 2 -dm_coord_petscspace_degree 3 -petscfe_default_quadrature_order 3 -tol 1.2e-7
391: testset:
392: # Volume: 4/3 \pi (8 - 1)/2 = 14/3 \pi = 14.66076571675238
393: args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1 -dm_plex_box_lower -1.,-1.,-1. -dm_plex_box_upper 1.,1.,1. -mesh_transform shell -volume 14.66076571675238 -dm_coord_space 0
395: test:
396: suffix: shell_0
397: requires: double
398: args: -dm_refine 1 -dm_coord_petscspace_degree 1 -petscfe_default_quadrature_order 1 -volume 5.633164922 -tol 1.0e-7
400: test:
401: suffix: shell_1
402: requires: double
403: args: -dm_refine 2 -dm_coord_petscspace_degree 1 -petscfe_default_quadrature_order 1 -tol 3.1
405: test:
406: suffix: shell_2
407: requires: double
408: args: -dm_refine 2 -dm_coord_petscspace_degree 2 -petscfe_default_quadrature_order 2 -tol .1
410: test:
411: suffix: shell_3
412: requires: double
413: args: -dm_refine 2 -dm_coord_petscspace_degree 3 -petscfe_default_quadrature_order 3 -tol .02
415: test:
416: suffix: shell_4
417: requires: double
418: args: -dm_refine 2 -dm_coord_petscspace_degree 4 -petscfe_default_quadrature_order 4 -tol .006
420: test:
421: # Volume: 1.0
422: suffix: gmsh_q2
423: requires: double
424: args: -coord_space 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/quads-q2.msh -dm_plex_gmsh_project -volume 1.0 -tol 1e-6
426: test:
427: # Volume: 1.0
428: suffix: gmsh_q3
429: requires: double
430: nsize: {{1 2}}
431: args: -coord_space 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/quads-q3.msh -dm_plex_gmsh_project -volume 1.0 -tol 1e-6
433: TEST*/