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 struct {
  7:   PetscReal volume; /* Analytical volume of the mesh */
  8:   PetscReal tol;    /* Tolerance for volume check */
  9: } AppCtx;

 11: PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
 12: {
 13:   PetscFunctionBegin;
 14:   options->volume = -1.0;
 15:   options->tol    = PETSC_SMALL;

 17:   PetscOptionsBegin(comm, "", "Meshing Interpolation Test Options", "DMPLEX");
 18:   PetscCall(PetscOptionsReal("-volume", "The analytical volume of the mesh", "ex33.c", options->volume, &options->volume, NULL));
 19:   PetscCall(PetscOptionsReal("-tol", "The tolerance for the volume check", "ex33.c", options->tol, &options->tol, NULL));
 20:   PetscOptionsEnd();
 21:   PetscFunctionReturn(PETSC_SUCCESS);
 22: }

 24: PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *ctx, DM *dm)
 25: {
 26:   PetscFunctionBegin;
 27:   PetscCall(DMCreate(comm, dm));
 28:   PetscCall(DMSetType(*dm, DMPLEX));
 29:   PetscCall(DMSetFromOptions(*dm));
 30:   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
 31:   PetscFunctionReturn(PETSC_SUCCESS);
 32: }

 34: 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[])
 35: {
 36:   vol[0] = 1.;
 37: }

 39: static PetscErrorCode CreateDiscretization(DM dm, AppCtx *ctx)
 40: {
 41:   PetscDS        ds;
 42:   PetscFE        fe;
 43:   DMPolytopeType ct;
 44:   PetscInt       dim, cStart;

 46:   PetscFunctionBeginUser;
 47:   PetscCall(DMGetDimension(dm, &dim));
 48:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
 49:   PetscCall(DMPlexGetCellType(dm, cStart, &ct));
 50:   PetscCall(PetscFECreateByCell(PETSC_COMM_SELF, dim, 1, ct, NULL, PETSC_DETERMINE, &fe));
 51:   PetscCall(PetscFESetName(fe, "scalar"));
 52:   PetscCall(DMAddField(dm, NULL, (PetscObject)fe));
 53:   PetscCall(PetscFEDestroy(&fe));
 54:   PetscCall(DMCreateDS(dm));
 55:   PetscCall(DMGetDS(dm, &ds));
 56:   PetscCall(PetscDSSetObjective(ds, 0, volume));
 57:   PetscFunctionReturn(PETSC_SUCCESS);
 58: }

 60: static PetscErrorCode CheckVolume(DM dm, AppCtx *ctx)
 61: {
 62:   Vec         u;
 63:   PetscScalar result;
 64:   PetscReal   vol, tol = ctx->tol;

 66:   PetscFunctionBeginUser;
 67:   PetscCall(DMGetGlobalVector(dm, &u));
 68:   PetscCall(DMPlexComputeIntegralFEM(dm, u, &result, ctx));
 69:   vol = PetscRealPart(result);
 70:   PetscCall(DMRestoreGlobalVector(dm, &u));
 71:   PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Volume: %g\n", (double)vol));
 72:   if (ctx->volume > 0.0 && PetscAbsReal(ctx->volume - vol) > tol) {
 73:     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);
 74:   }
 75:   PetscFunctionReturn(PETSC_SUCCESS);
 76: }

 78: int main(int argc, char **argv)
 79: {
 80:   DM     dm;
 81:   AppCtx user;

 83:   PetscFunctionBeginUser;
 84:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 85:   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
 86:   PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
 87:   PetscCall(CreateDiscretization(dm, &user));
 88:   PetscCall(CheckVolume(dm, &user));
 89:   PetscCall(DMDestroy(&dm));
 90:   PetscCall(PetscFinalize());
 91:   return 0;
 92: }

 94: /*TEST

 96:   testset:
 97:     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.

 99:     test:
100:       suffix: square_0
101:       args: -dm_coord_petscspace_degree 1

103:     test:
104:       suffix: square_1
105:       args: -dm_coord_petscspace_degree 2

107:     test:
108:       suffix: square_2
109:       args: -dm_refine 1 -dm_coord_petscspace_degree 1

111:     test:
112:       suffix: square_3
113:       args: -dm_refine 1 -dm_coord_petscspace_degree 2

115:   testset:
116:     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.

118:     test:
119:       suffix: cube_0
120:       args: -dm_coord_petscspace_degree 1

122:     test:
123:       suffix: cube_1
124:       args: -dm_coord_petscspace_degree 2

126:     test:
127:       suffix: cube_2
128:       args: -dm_refine 1 -dm_coord_petscspace_degree 1

130:     test:
131:       suffix: cube_3
132:       args: -dm_refine 1 -dm_coord_petscspace_degree 2

134:   testset:
135:     args: -dm_plex_simplex 0 -dm_plex_box_faces 1,1 -dm_plex_box_lower -1.,-1. -dm_plex_box_upper 1.,1. \
136:           -volume 4. -dm_coord_remap -dm_coord_map shear -dm_coord_map_params 0.0,0.0,3.0

138:     test:
139:       suffix: shear_0
140:       args: -dm_coord_petscspace_degree 1

142:     test:
143:       suffix: shear_1
144:       args: -dm_coord_petscspace_degree 2

146:     test:
147:       suffix: shear_2
148:       args: -dm_refine 1 -dm_coord_petscspace_degree 1

150:     test:
151:       suffix: shear_3
152:       args: -dm_refine 1 -dm_coord_petscspace_degree 2

154:   testset:
155:     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. \
156:           -volume 8. -dm_coord_remap -dm_coord_map shear -dm_coord_map_params 0.0,0.0,3.0,4.0

158:     test:
159:       suffix: shear_4
160:       args: -dm_coord_petscspace_degree 1

162:     test:
163:       suffix: shear_5
164:       args: -dm_coord_petscspace_degree 2

166:     test:
167:       suffix: shear_6
168:       args: -dm_refine 1 -dm_coord_petscspace_degree 1

170:     test:
171:       suffix: shear_7
172:       args: -dm_refine 1 -dm_coord_petscspace_degree 2

174:   testset:
175:     args: -dm_plex_simplex 0 -dm_plex_box_faces 1,1 -dm_plex_box_lower 1.,-1. -dm_plex_box_upper 3.,1. \
176:           -dm_coord_petscspace_degree 1 -volume 8. -dm_coord_remap -dm_coord_map flare

178:     test:
179:       suffix: flare_0
180:       args:

182:     test:
183:       suffix: flare_1
184:       args: -dm_refine 2

186:   testset:
187:     # Area: 3/4 \pi = 2.3562
188:     args: -dm_plex_simplex 0 -dm_plex_box_faces 1,1 -petscfe_default_quadrature_order 2 \
189:           -volume 2.35619449019235 -dm_coord_remap -dm_coord_map annulus

191:     test:
192:       # Area: (a+b)/2 h = 3/\sqrt{2} (sqrt{2} - 1/\sqrt{2}) = 3/2
193:       suffix: annulus_0
194:       requires: double
195:       args: -dm_coord_petscspace_degree 1 -volume 1.5

197:     test:
198:       suffix: annulus_1
199:       requires: double
200:       args: -dm_refine 3 -dm_coord_petscspace_degree 1 -tol .016

202:     test:
203:       suffix: annulus_2
204:       requires: double
205:       args: -dm_refine 3 -dm_coord_petscspace_degree 2 -tol .0038

207:     test:
208:       suffix: annulus_3
209:       requires: double
210:       args: -dm_refine 3 -dm_coord_petscspace_degree 3 -tol 2.2e-6

212:     test:
213:       suffix: annulus_4
214:       requires: double
215:       args: -dm_refine 2 -dm_coord_petscspace_degree 2 -tol .00012

217:     test:
218:       suffix: annulus_5
219:       requires: double
220:       args: -dm_refine 2 -dm_coord_petscspace_degree 3 -petscfe_default_quadrature_order 3 -tol 1.2e-7

222:   testset:
223:     # Volume: 4/3 \pi (8 - 1)/2 = 14/3 \pi = 14.66076571675238
224:     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. \
225:           -volume 14.66076571675238 -dm_coord_remap -dm_coord_map shell

227:     test:
228:       suffix: shell_0
229:       requires: double
230:       args: -dm_refine 1 -dm_coord_petscspace_degree 1 -petscfe_default_quadrature_order 1 -volume 5.633164922 -tol 1.0e-7

232:     test:
233:       suffix: shell_1
234:       requires: double
235:       args: -dm_refine 2 -dm_coord_petscspace_degree 1 -petscfe_default_quadrature_order 1 -tol 3.1

237:     test:
238:       suffix: shell_2
239:       requires: double
240:       args: -dm_refine 2 -dm_coord_petscspace_degree 2 -petscfe_default_quadrature_order 2 -tol .1

242:     test:
243:       suffix: shell_3
244:       requires: double
245:       args: -dm_refine 2 -dm_coord_petscspace_degree 3 -petscfe_default_quadrature_order 3 -tol .02

247:     test:
248:       suffix: shell_4
249:       requires: double
250:       args: -dm_refine 2 -dm_coord_petscspace_degree 4 -petscfe_default_quadrature_order 4 -tol .006

252:   test:
253:     # Volume: 1.0
254:     suffix: gmsh_q2
255:     requires: double
256:     args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/quads-q2.msh -dm_plex_gmsh_project -volume 1.0 -tol 1e-6

258:   test:
259:     # Volume: 1.0
260:     suffix: gmsh_q3
261:     requires: double
262:     nsize: {{1 2}}
263:     args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/quads-q3.msh -dm_plex_gmsh_project -volume 1.0 -tol 1e-6

265:   test:
266:     # Volume: 1.0
267:     suffix: gmsh_3d_q2
268:     requires: double
269:     args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/cube_q2.msh -dm_plex_gmsh_project -volume 1.0 -tol 1e-6

271:   test:
272:     # Volume: 1.0
273:     suffix: gmsh_3d_q3
274:     requires: double
275:     nsize: {{1 2}}
276:     args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/cube_q3.msh -dm_plex_gmsh_project -volume 1.0 -tol 1e-6

278: TEST*/