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:
437:     # Volume: 1.0
438:     suffix: gmsh_3d_q2
439:     requires: double
440:     args: -coord_space 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/cube_q2.msh -dm_plex_gmsh_project -volume 1.0 -tol 1e-6

442:   test:
443:     # Volume: 1.0
444:     suffix: gmsh_3d_q3
445:     requires: double
446:     nsize: {{1 2}}
447:     args: -coord_space 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/cube_q3.msh -dm_plex_gmsh_project -volume 1.0 -tol 1e-6

449: TEST*/