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*/