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