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