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: PetscFunctionBegin;
29: options->coordSpace = PETSC_TRUE;
30: options->meshTransform = TRANSFORM_NONE;
31: options->transformDataReal = NULL;
32: options->transformData = NULL;
33: options->volume = -1.0;
34: options->tol = PETSC_SMALL;
36: PetscOptionsBegin(comm, "", "Meshing Interpolation Test Options", "DMPLEX");
37: PetscCall(PetscOptionsBool("-coord_space", "Flag to create a coordinate space", "ex33.c", options->coordSpace, &options->coordSpace, NULL));
38: PetscCall(PetscOptionsEnum("-mesh_transform", "Method to transform initial box mesh <none, shear, annulus, shell>", "ex33.c", TransformTypes, (PetscEnum)options->meshTransform, (PetscEnum *)&options->meshTransform, NULL));
39: switch (options->meshTransform) {
40: case TRANSFORM_NONE:
41: break;
42: case TRANSFORM_SHEAR:
43: n = 2;
44: PetscCall(PetscMalloc1(n, &options->transformDataReal));
45: for (i = 0; i < n; ++i) options->transformDataReal[i] = 1.0;
46: PetscCall(PetscOptionsRealArray("-transform_data", "Parameters for mesh transforms", "ex33.c", options->transformDataReal, &n, NULL));
47: break;
48: case TRANSFORM_FLARE:
49: n = 4;
50: PetscCall(PetscMalloc1(n, &options->transformData));
51: for (i = 0; i < n; ++i) options->transformData[i] = 1.0;
52: options->transformData[0] = (PetscScalar)0;
53: PetscCall(PetscOptionsScalarArray("-transform_data", "Parameters for mesh transforms", "ex33.c", options->transformData, &n, NULL));
54: break;
55: case TRANSFORM_ANNULUS:
56: n = 2;
57: PetscCall(PetscMalloc1(n, &options->transformData));
58: options->transformData[0] = 1.0;
59: options->transformData[1] = 2.0;
60: PetscCall(PetscOptionsScalarArray("-transform_data", "Parameters for mesh transforms", "ex33.c", options->transformData, &n, NULL));
61: break;
62: case TRANSFORM_SHELL:
63: n = 2;
64: PetscCall(PetscMalloc1(n, &options->transformData));
65: options->transformData[0] = 1.0;
66: options->transformData[1] = 2.0;
67: PetscCall(PetscOptionsScalarArray("-transform_data", "Parameters for mesh transforms", "ex33.c", options->transformData, &n, NULL));
68: break;
69: default:
70: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Unknown mesh transform %d", options->meshTransform);
71: }
72: PetscCall(PetscOptionsReal("-volume", "The analytical volume of the mesh", "ex33.c", options->volume, &options->volume, NULL));
73: PetscCall(PetscOptionsReal("-tol", "The tolerance for the volume check", "ex33.c", options->tol, &options->tol, NULL));
74: PetscOptionsEnd();
75: PetscFunctionReturn(PETSC_SUCCESS);
76: }
78: 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[])
79: {
80: const PetscInt Nc = uOff[1] - uOff[0];
81: PetscInt c;
83: for (c = 0; c < Nc; ++c) f0[c] = u[c];
84: }
86: /* Flare applies the transformation, assuming we fix x_f,
88: x_i = x_i * alpha_i x_f
89: */
90: 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[])
91: {
92: const PetscInt Nc = uOff[1] - uOff[0];
93: const PetscInt cf = (PetscInt)PetscRealPart(constants[0]);
94: PetscInt c;
96: for (c = 0; c < Nc; ++c) coords[c] = u[c] * (c == cf ? 1.0 : constants[c + 1] * u[cf]);
97: }
99: /*
100: 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
101: will correspond to the top and bottom of our square. So
103: (0,0)--(1,0) ==> (1,0)--(2,0) Just a shift of (1,0)
104: (0,1)--(1,1) ==> (0,1)--(0,2) Switch x and y
106: 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:
108: (x, y) ==> (x+1, \pi/2 y) in (r', \theta') space
109: ==> ((x+1) cos(\pi/2 y), (x+1) sin(\pi/2 y)) in (x', y') space
110: */
111: 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[])
112: {
113: const PetscReal ri = PetscRealPart(constants[0]);
114: const PetscReal ro = PetscRealPart(constants[1]);
116: xp[0] = (x[0] * (ro - ri) + ri) * PetscCosReal(0.5 * PETSC_PI * x[1]);
117: xp[1] = (x[0] * (ro - ri) + ri) * PetscSinReal(0.5 * PETSC_PI * x[1]);
118: }
120: /*
121: 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
122: lower hemisphere and the upper surface onto the top, letting z be the radius.
124: (x, y) ==> ((z+3)/2, \pi/2 (|x| or |y|), arctan y/x) in (r', \theta', \phi') space
125: ==> ((z+3)/2 \cos(\theta') cos(\phi'), (z+3)/2 \cos(\theta') sin(\phi'), (z+3)/2 sin(\theta')) in (x', y', z') space
126: */
127: 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[])
128: {
129: const PetscReal pi4 = PETSC_PI / 4.0;
130: const PetscReal ri = PetscRealPart(constants[0]);
131: const PetscReal ro = PetscRealPart(constants[1]);
132: const PetscReal rp = (x[2] + 1) * 0.5 * (ro - ri) + ri;
133: const PetscReal phip = PetscAtan2Real(x[1], x[0]);
134: 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])));
136: xp[0] = rp * PetscCosReal(thetap) * PetscCosReal(phip);
137: xp[1] = rp * PetscCosReal(thetap) * PetscSinReal(phip);
138: xp[2] = rp * PetscSinReal(thetap);
139: }
141: static PetscErrorCode DMCreateCoordinateDisc(DM dm)
142: {
143: DM cdm;
144: PetscFE fe;
145: DMPolytopeType ct;
146: PetscInt dim, dE, cStart;
147: PetscBool simplex;
149: PetscFunctionBegin;
150: PetscCall(DMGetCoordinateDM(dm, &cdm));
151: PetscCall(DMGetDimension(dm, &dim));
152: PetscCall(DMGetCoordinateDim(dm, &dE));
153: PetscCall(DMPlexGetHeightStratum(cdm, 0, &cStart, NULL));
154: PetscCall(DMPlexGetCellType(dm, cStart, &ct));
155: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
156: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, dE, simplex, "dm_coord_", -1, &fe));
157: PetscCall(DMProjectCoordinates(dm, fe));
158: PetscCall(PetscFEDestroy(&fe));
159: PetscFunctionReturn(PETSC_SUCCESS);
160: }
162: PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *ctx, DM *dm)
163: {
164: DM cdm;
165: PetscDS cds;
167: PetscFunctionBegin;
168: PetscCall(DMCreate(comm, dm));
169: PetscCall(DMSetType(*dm, DMPLEX));
170: PetscCall(DMSetFromOptions(*dm));
172: if (ctx->coordSpace) PetscCall(DMCreateCoordinateDisc(*dm));
173: switch (ctx->meshTransform) {
174: case TRANSFORM_NONE:
175: PetscCall(DMPlexRemapGeometry(*dm, 0.0, identity));
176: break;
177: case TRANSFORM_SHEAR:
178: PetscCall(DMPlexShearGeometry(*dm, DM_X, ctx->transformDataReal));
179: break;
180: case TRANSFORM_FLARE:
181: PetscCall(DMGetCoordinateDM(*dm, &cdm));
182: PetscCall(DMGetDS(cdm, &cds));
183: PetscCall(PetscDSSetConstants(cds, 4, ctx->transformData));
184: PetscCall(DMPlexRemapGeometry(*dm, 0.0, f0_flare));
185: break;
186: case TRANSFORM_ANNULUS:
187: PetscCall(DMGetCoordinateDM(*dm, &cdm));
188: PetscCall(DMGetDS(cdm, &cds));
189: PetscCall(PetscDSSetConstants(cds, 2, ctx->transformData));
190: PetscCall(DMPlexRemapGeometry(*dm, 0.0, f0_annulus));
191: break;
192: case TRANSFORM_SHELL:
193: PetscCall(DMGetCoordinateDM(*dm, &cdm));
194: PetscCall(DMGetDS(cdm, &cds));
195: PetscCall(PetscDSSetConstants(cds, 2, ctx->transformData));
196: PetscCall(DMPlexRemapGeometry(*dm, 0.0, f0_shell));
197: break;
198: default:
199: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Unknown mesh transform %d", ctx->meshTransform);
200: }
201: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
202: PetscFunctionReturn(PETSC_SUCCESS);
203: }
205: 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[])
206: {
207: vol[0] = 1.;
208: }
210: static PetscErrorCode CreateDiscretization(DM dm, AppCtx *ctx)
211: {
212: PetscDS ds;
213: PetscFE fe;
214: DMPolytopeType ct;
215: PetscInt dim, cStart;
216: PetscBool simplex;
218: PetscFunctionBeginUser;
219: PetscCall(DMGetDimension(dm, &dim));
220: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
221: PetscCall(DMPlexGetCellType(dm, cStart, &ct));
222: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
223: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, NULL, PETSC_DETERMINE, &fe));
224: PetscCall(PetscFESetName(fe, "scalar"));
225: PetscCall(DMAddField(dm, NULL, (PetscObject)fe));
226: PetscCall(PetscFEDestroy(&fe));
227: PetscCall(DMCreateDS(dm));
228: PetscCall(DMGetDS(dm, &ds));
229: PetscCall(PetscDSSetObjective(ds, 0, volume));
230: PetscFunctionReturn(PETSC_SUCCESS);
231: }
233: static PetscErrorCode CheckVolume(DM dm, AppCtx *ctx)
234: {
235: Vec u;
236: PetscScalar result;
237: PetscReal vol, tol = ctx->tol;
239: PetscFunctionBeginUser;
240: PetscCall(DMGetGlobalVector(dm, &u));
241: PetscCall(DMPlexComputeIntegralFEM(dm, u, &result, ctx));
242: vol = PetscRealPart(result);
243: PetscCall(DMRestoreGlobalVector(dm, &u));
244: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Volume: %g\n", (double)vol));
245: if (ctx->volume > 0.0 && PetscAbsReal(ctx->volume - vol) > tol) {
246: 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);
247: }
248: PetscFunctionReturn(PETSC_SUCCESS);
249: }
251: int main(int argc, char **argv)
252: {
253: DM dm;
254: AppCtx user;
256: PetscFunctionBeginUser;
257: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
258: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
259: PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
260: PetscCall(CreateDiscretization(dm, &user));
261: PetscCall(CheckVolume(dm, &user));
262: PetscCall(DMDestroy(&dm));
263: PetscCall(PetscFree(user.transformDataReal));
264: PetscCall(PetscFree(user.transformData));
265: PetscCall(PetscFinalize());
266: return 0;
267: }
269: /*TEST
271: testset:
272: 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
274: test:
275: suffix: square_0
276: args: -dm_coord_petscspace_degree 1
278: test:
279: suffix: square_1
280: args: -dm_coord_petscspace_degree 2
282: test:
283: suffix: square_2
284: args: -dm_refine 1 -dm_coord_petscspace_degree 1
286: test:
287: suffix: square_3
288: args: -dm_refine 1 -dm_coord_petscspace_degree 2
290: testset:
291: 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
293: test:
294: suffix: cube_0
295: args: -dm_coord_petscspace_degree 1
297: test:
298: suffix: cube_1
299: args: -dm_coord_petscspace_degree 2
301: test:
302: suffix: cube_2
303: args: -dm_refine 1 -dm_coord_petscspace_degree 1
305: test:
306: suffix: cube_3
307: args: -dm_refine 1 -dm_coord_petscspace_degree 2
309: testset:
310: 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
312: test:
313: suffix: shear_0
314: args: -dm_coord_petscspace_degree 1 -mesh_transform shear -transform_data 3.0
316: test:
317: suffix: shear_1
318: args: -dm_coord_petscspace_degree 2 -mesh_transform shear -transform_data 3.0
320: test:
321: suffix: shear_2
322: args: -dm_refine 1 -dm_coord_petscspace_degree 1 -mesh_transform shear -transform_data 3.0
324: test:
325: suffix: shear_3
326: args: -dm_refine 1 -dm_coord_petscspace_degree 2 -mesh_transform shear -transform_data 3.0
328: testset:
329: 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
331: test:
332: suffix: shear_4
333: args: -dm_coord_petscspace_degree 1 -mesh_transform shear -transform_data 3.0
335: test:
336: suffix: shear_5
337: args: -dm_coord_petscspace_degree 2 -mesh_transform shear -transform_data 3.0
339: test:
340: suffix: shear_6
341: args: -dm_refine 1 -dm_coord_petscspace_degree 1 -mesh_transform shear -transform_data 3.0,4.0
343: test:
344: suffix: shear_7
345: args: -dm_refine 1 -dm_coord_petscspace_degree 2 -mesh_transform shear -transform_data 3.0,4.0
347: testset:
348: 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. \
349: -dm_coord_petscspace_degree 1 -mesh_transform flare -volume 8.
351: test:
352: suffix: flare_0
353: args:
355: test:
356: suffix: flare_1
357: args: -dm_refine 2
359: testset:
360: # Area: 3/4 \pi = 2.3562
361: args: -dm_plex_simplex 0 -dm_plex_box_faces 1,1 -mesh_transform annulus -volume 2.35619449019235 -dm_coord_space 0
363: test:
364: # Area: (a+b)/2 h = 3/\sqrt{2} (sqrt{2} - 1/\sqrt{2}) = 3/2
365: suffix: annulus_0
366: requires: double
367: args: -dm_coord_petscspace_degree 1 -volume 1.5
369: test:
370: suffix: annulus_1
371: requires: double
372: args: -dm_refine 3 -dm_coord_petscspace_degree 1 -tol .016
374: test:
375: suffix: annulus_2
376: requires: double
377: args: -dm_refine 3 -dm_coord_petscspace_degree 2 -tol .0038
379: test:
380: suffix: annulus_3
381: requires: double
382: args: -dm_refine 3 -dm_coord_petscspace_degree 3 -tol 2.2e-6
384: test:
385: suffix: annulus_4
386: requires: double
387: args: -dm_refine 2 -dm_coord_petscspace_degree 2 -petscfe_default_quadrature_order 2 -tol .00012
389: test:
390: suffix: annulus_5
391: requires: double
392: args: -dm_refine 2 -dm_coord_petscspace_degree 3 -petscfe_default_quadrature_order 3 -tol 1.2e-7
394: testset:
395: # Volume: 4/3 \pi (8 - 1)/2 = 14/3 \pi = 14.66076571675238
396: 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
398: test:
399: suffix: shell_0
400: requires: double
401: args: -dm_refine 1 -dm_coord_petscspace_degree 1 -petscfe_default_quadrature_order 1 -volume 5.633164922 -tol 1.0e-7
403: test:
404: suffix: shell_1
405: requires: double
406: args: -dm_refine 2 -dm_coord_petscspace_degree 1 -petscfe_default_quadrature_order 1 -tol 3.1
408: test:
409: suffix: shell_2
410: requires: double
411: args: -dm_refine 2 -dm_coord_petscspace_degree 2 -petscfe_default_quadrature_order 2 -tol .1
413: test:
414: suffix: shell_3
415: requires: double
416: args: -dm_refine 2 -dm_coord_petscspace_degree 3 -petscfe_default_quadrature_order 3 -tol .02
418: test:
419: suffix: shell_4
420: requires: double
421: args: -dm_refine 2 -dm_coord_petscspace_degree 4 -petscfe_default_quadrature_order 4 -tol .006
423: test:
424: # Volume: 1.0
425: suffix: gmsh_q2
426: requires: double
427: 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
429: test:
430: # Volume: 1.0
431: suffix: gmsh_q3
432: requires: double
433: nsize: {{1 2}}
434: 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
436: TEST*/