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 { 26 PetscInt n = 0, i; 27 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; 35 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(0); 76 } 77 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; 82 83 for (c = 0; c < Nc; ++c) f0[c] = u[c]; 84 } 85 86 /* Flare applies the transformation, assuming we fix x_f, 87 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; 95 96 for (c = 0; c < Nc; ++c) coords[c] = u[c] * (c == cf ? 1.0 : constants[c + 1] * u[cf]); 97 } 98 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 102 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 105 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: 107 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]); 115 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 } 119 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. 123 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]))); 135 136 xp[0] = rp * PetscCosReal(thetap) * PetscCosReal(phip); 137 xp[1] = rp * PetscCosReal(thetap) * PetscSinReal(phip); 138 xp[2] = rp * PetscSinReal(thetap); 139 } 140 141 static PetscErrorCode DMCreateCoordinateDisc(DM dm) 142 { 143 DM cdm; 144 PetscFE fe; 145 DMPolytopeType ct; 146 PetscInt dim, dE, cStart; 147 PetscBool simplex; 148 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(0); 160 } 161 162 PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *ctx, DM *dm) 163 { 164 DM cdm; 165 PetscDS cds; 166 167 PetscFunctionBegin; 168 PetscCall(DMCreate(comm, dm)); 169 PetscCall(DMSetType(*dm, DMPLEX)); 170 PetscCall(DMSetFromOptions(*dm)); 171 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(0); 203 } 204 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 } 209 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; 217 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(0); 231 } 232 233 static PetscErrorCode CheckVolume(DM dm, AppCtx *ctx) 234 { 235 Vec u; 236 PetscScalar result; 237 PetscReal vol, tol = ctx->tol; 238 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(0); 249 } 250 251 int main(int argc, char **argv) 252 { 253 DM dm; 254 AppCtx user; 255 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 } 268 269 /*TEST 270 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 273 274 test: 275 suffix: square_0 276 args: -dm_coord_petscspace_degree 1 277 278 test: 279 suffix: square_1 280 args: -dm_coord_petscspace_degree 2 281 282 test: 283 suffix: square_2 284 args: -dm_refine 1 -dm_coord_petscspace_degree 1 285 286 test: 287 suffix: square_3 288 args: -dm_refine 1 -dm_coord_petscspace_degree 2 289 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 292 293 test: 294 suffix: cube_0 295 args: -dm_coord_petscspace_degree 1 296 297 test: 298 suffix: cube_1 299 args: -dm_coord_petscspace_degree 2 300 301 test: 302 suffix: cube_2 303 args: -dm_refine 1 -dm_coord_petscspace_degree 1 304 305 test: 306 suffix: cube_3 307 args: -dm_refine 1 -dm_coord_petscspace_degree 2 308 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 311 312 test: 313 suffix: shear_0 314 args: -dm_coord_petscspace_degree 1 -mesh_transform shear -transform_data 3.0 315 316 test: 317 suffix: shear_1 318 args: -dm_coord_petscspace_degree 2 -mesh_transform shear -transform_data 3.0 319 320 test: 321 suffix: shear_2 322 args: -dm_refine 1 -dm_coord_petscspace_degree 1 -mesh_transform shear -transform_data 3.0 323 324 test: 325 suffix: shear_3 326 args: -dm_refine 1 -dm_coord_petscspace_degree 2 -mesh_transform shear -transform_data 3.0 327 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 330 331 test: 332 suffix: shear_4 333 args: -dm_coord_petscspace_degree 1 -mesh_transform shear -transform_data 3.0 334 335 test: 336 suffix: shear_5 337 args: -dm_coord_petscspace_degree 2 -mesh_transform shear -transform_data 3.0 338 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 342 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 346 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. 350 351 test: 352 suffix: flare_0 353 args: 354 355 test: 356 suffix: flare_1 357 args: -dm_refine 2 358 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 362 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 368 369 test: 370 suffix: annulus_1 371 requires: double 372 args: -dm_refine 3 -dm_coord_petscspace_degree 1 -tol .016 373 374 test: 375 suffix: annulus_2 376 requires: double 377 args: -dm_refine 3 -dm_coord_petscspace_degree 2 -tol .0038 378 379 test: 380 suffix: annulus_3 381 requires: double 382 args: -dm_refine 3 -dm_coord_petscspace_degree 3 -tol 2.2e-6 383 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 388 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 393 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 397 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 402 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 407 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 412 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 417 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 422 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 428 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 435 436 TEST*/ 437