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