1 static char help[] = "Tests for mesh interpolation\n\n"; 2 3 #include <petscdmplex.h> 4 5 /* List of test meshes 6 7 Triangle 8 -------- 9 Test 0: 10 Two triangles sharing a face 11 12 4 13 / | \ 14 / | \ 15 / | \ 16 2 0 | 1 5 17 \ | / 18 \ | / 19 \ | / 20 3 21 22 should become 23 24 4 25 / | \ 26 8 | 9 27 / | \ 28 2 0 7 1 5 29 \ | / 30 6 | 10 31 \ | / 32 3 33 34 Tetrahedron 35 ----------- 36 Test 0: 37 Two tets sharing a face 38 39 cell 5 _______ cell 40 0 / | \ \ 1 41 / | \ \ 42 / | \ \ 43 2----|----4-----6 44 \ | / / 45 \ | / / 46 \ | / / 47 3------- 48 49 should become 50 51 cell 5 _______ cell 52 0 / | \ \ 1 53 / | \ \ 54 17 | 18 13 22 55 / 8 19 10 \ \ 56 / | \ \ 57 2---14-|------4--20--6 58 \ | / / 59 \ 9 | 7 / / 60 16 | 15 11 21 61 \ | / / 62 \ | / / 63 3------- 64 65 Quadrilateral 66 ------------- 67 Test 0: 68 Two quads sharing a face 69 70 5-------4-------7 71 | | | 72 | 0 | 1 | 73 | | | 74 2-------3-------6 75 76 should become 77 78 5--10---4--14---7 79 | | | 80 11 0 9 1 13 81 | | | 82 2---8---3--12---6 83 84 Test 1: 85 A quad and a triangle sharing a face 86 87 5-------4 88 | | \ 89 | 0 | \ 90 | | 1 \ 91 2-------3----6 92 93 should become 94 95 5---9---4 96 | | \ 97 10 0 8 12 98 | | 1 \ 99 2---7---3-11-6 100 101 Hexahedron 102 ---------- 103 Test 0: 104 Two hexes sharing a face 105 106 cell 9-------------8-------------13 cell 107 0 /| /| /| 1 108 / | 15 / | 21 / | 109 / | / | / | 110 6-------------7-------------12 | 111 | | 18 | | 24 | | 112 | | | | | | 113 |19 | |17 | |23 | 114 | | 16 | | 22 | | 115 | 5---------|---4---------|---11 116 | / | / | / 117 | / 14 | / 20 | / 118 |/ |/ |/ 119 2-------------3-------------10 120 121 should become 122 123 cell 9-----31------8-----42------13 cell 124 0 /| /| /| 1 125 32 | 15 30| 21 41| 126 / | / | / | 127 6-----29------7-----40------12 | 128 | | 17 | | 23 | | 129 | 35 | 36 | 44 130 |19 | |18 | |24 | 131 34 | 16 33 | 22 43 | 132 | 5-----26--|---4-----37--|---11 133 | / | / | / 134 | 25 14 | 27 20 | 38 135 |/ |/ |/ 136 2-----28------3-----39------10 137 138 */ 139 140 typedef struct { 141 PetscInt testNum; /* Indicates the mesh to create */ 142 PetscInt dim; /* The topological mesh dimension */ 143 PetscBool simplex; /* Use simplices or hexes */ 144 PetscBool useGenerator; /* Construct mesh with a mesh generator */ 145 PetscBool femGeometry; /* Flag to compute FEM geometry */ 146 } AppCtx; 147 148 PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 149 { 150 PetscErrorCode ierr; 151 152 PetscFunctionBegin; 153 options->testNum = 0; 154 options->dim = 2; 155 options->simplex = PETSC_TRUE; 156 options->useGenerator = PETSC_FALSE; 157 options->femGeometry = PETSC_TRUE; 158 159 ierr = PetscOptionsBegin(comm, "", "Meshing Interpolation Test Options", "DMPLEX");CHKERRQ(ierr); 160 ierr = PetscOptionsBoundedInt("-testnum", "The mesh to create", "ex7.c", options->testNum, &options->testNum, NULL,0);CHKERRQ(ierr); 161 ierr = PetscOptionsRangeInt("-dim", "The topological mesh dimension", "ex7.c", options->dim, &options->dim, NULL,1,3);CHKERRQ(ierr); 162 ierr = PetscOptionsBool("-simplex", "Use simplices if true, otherwise hexes", "ex7.c", options->simplex, &options->simplex, NULL);CHKERRQ(ierr); 163 ierr = PetscOptionsBool("-use_generator", "Use a mesh generator to build the mesh", "ex7.c", options->useGenerator, &options->useGenerator, NULL);CHKERRQ(ierr); 164 ierr = PetscOptionsBool("-fem_geometry", "Flag to compute FEM geometry", "ex7.c", options->femGeometry, &options->femGeometry, NULL);CHKERRQ(ierr); 165 ierr = PetscOptionsEnd(); 166 PetscFunctionReturn(0); 167 } 168 169 PetscErrorCode CreateSimplex_2D(MPI_Comm comm, DM dm) 170 { 171 PetscInt depth = 1, testNum = 0, p; 172 PetscMPIInt rank; 173 PetscErrorCode ierr; 174 175 PetscFunctionBegin; 176 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 177 if (!rank) { 178 switch (testNum) { 179 case 0: 180 { 181 PetscInt numPoints[2] = {4, 2}; 182 PetscInt coneSize[6] = {3, 3, 0, 0, 0, 0}; 183 PetscInt cones[6] = {2, 3, 4, 5, 4, 3}; 184 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 185 PetscScalar vertexCoords[8] = {-0.5, 0.5, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5}; 186 PetscInt markerPoints[8] = {2, 1, 3, 1, 4, 1, 5, 1}; 187 188 ierr = DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 189 for (p = 0; p < 4; ++p) { 190 ierr = DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr); 191 } 192 } 193 break; 194 default: 195 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %d", testNum); 196 } 197 } else { 198 PetscInt numPoints[2] = {0, 0}; 199 200 ierr = DMPlexCreateFromDAG(dm, depth, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 201 } 202 PetscFunctionReturn(0); 203 } 204 205 PetscErrorCode CreateSimplex_3D(MPI_Comm comm, DM dm) 206 { 207 PetscInt depth = 1, testNum = 0, p; 208 PetscMPIInt rank; 209 PetscErrorCode ierr; 210 211 PetscFunctionBegin; 212 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 213 if (!rank) { 214 switch (testNum) { 215 case 0: 216 { 217 PetscInt numPoints[2] = {5, 2}; 218 PetscInt coneSize[23] = {4, 4, 0, 0, 0, 0, 0}; 219 PetscInt cones[8] = {2, 4, 3, 5, 3, 4, 6, 5}; 220 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 221 PetscScalar vertexCoords[15] = {0.0, 0.0, -0.5, 0.0, -0.5, 0.0, 1.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5}; 222 PetscInt markerPoints[8] = {2, 1, 3, 1, 4, 1, 5, 1}; 223 224 ierr = DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 225 for (p = 0; p < 4; ++p) { 226 ierr = DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr); 227 } 228 } 229 break; 230 default: 231 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %d", testNum); 232 } 233 } else { 234 PetscInt numPoints[2] = {0, 0}; 235 236 ierr = DMPlexCreateFromDAG(dm, depth, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 237 } 238 PetscFunctionReturn(0); 239 } 240 241 PetscErrorCode CreateQuad_2D(MPI_Comm comm, PetscInt testNum, DM dm) 242 { 243 PetscInt depth = 1, p; 244 PetscMPIInt rank; 245 PetscErrorCode ierr; 246 247 PetscFunctionBegin; 248 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 249 if (!rank) { 250 switch (testNum) { 251 case 0: 252 { 253 PetscInt numPoints[2] = {6, 2}; 254 PetscInt coneSize[8] = {4, 4, 0, 0, 0, 0, 0, 0}; 255 PetscInt cones[8] = {2, 3, 4, 5, 3, 6, 7, 4}; 256 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 257 PetscScalar vertexCoords[12] = {-0.5, 0.0, 0.0, 0.0, 0.0, 1.0, -0.5, 1.0, 0.5, 0.0, 0.5, 1.0}; 258 PetscInt markerPoints[12] = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1, 7, 1}; 259 260 ierr = DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 261 for (p = 0; p < 6; ++p) { 262 ierr = DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr); 263 } 264 } 265 break; 266 case 1: 267 { 268 PetscInt numPoints[2] = {5, 2}; 269 PetscInt coneSize[7] = {4, 3, 0, 0, 0, 0, 0}; 270 PetscInt cones[7] = {2, 3, 4, 5, 3, 6, 4}; 271 PetscInt coneOrientations[7] = {0, 0, 0, 0, 0, 0, 0}; 272 PetscScalar vertexCoords[10] = {-0.5, 0.0, 0.0, 0.0, 0.0, 1.0, -0.5, 1.0, 0.5, 0.0}; 273 PetscInt markerPoints[10] = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1}; 274 275 ierr = DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 276 for (p = 0; p < 5; ++p) { 277 ierr = DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr); 278 } 279 } 280 break; 281 default: 282 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %d", testNum); 283 } 284 } else { 285 PetscInt numPoints[2] = {0, 0}; 286 287 ierr = DMPlexCreateFromDAG(dm, depth, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 288 } 289 PetscFunctionReturn(0); 290 } 291 292 PetscErrorCode CreateHex_3D(MPI_Comm comm, DM dm) 293 { 294 PetscInt depth = 1, testNum = 0, p; 295 PetscMPIInt rank; 296 PetscErrorCode ierr; 297 298 PetscFunctionBegin; 299 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 300 if (!rank) { 301 switch (testNum) { 302 case 0: 303 { 304 PetscInt numPoints[2] = {12, 2}; 305 PetscInt coneSize[14] = {8, 8, 0,0,0,0,0,0,0,0,0,0,0,0}; 306 PetscInt cones[16] = {2,5,4,3,6,7,8,9, 3,4,11,10,7,12,13,8}; 307 PetscInt coneOrientations[16] = {0,0,0,0,0,0,0,0, 0,0, 0, 0,0, 0, 0,0}; 308 PetscScalar vertexCoords[36] = {-0.5,0.0,0.0, 0.0,0.0,0.0, 0.0,1.0,0.0, -0.5,1.0,0.0, 309 -0.5,0.0,1.0, 0.0,0.0,1.0, 0.0,1.0,1.0, -0.5,1.0,1.0, 310 0.5,0.0,0.0, 0.5,1.0,0.0, 0.5,0.0,1.0, 0.5,1.0,1.0}; 311 PetscInt markerPoints[16] = {2,1,3,1,4,1,5,1,6,1,7,1,8,1,9,1}; 312 313 ierr = DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 314 for (p = 0; p < 8; ++p) { 315 ierr = DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr); 316 } 317 } 318 break; 319 default: 320 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %d", testNum); 321 } 322 } else { 323 PetscInt numPoints[2] = {0, 0}; 324 325 ierr = DMPlexCreateFromDAG(dm, depth, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 326 } 327 PetscFunctionReturn(0); 328 } 329 330 PetscErrorCode CheckMesh(DM dm, AppCtx *user) 331 { 332 PetscReal detJ, J[9], refVol = 1.0; 333 PetscReal vol; 334 PetscInt dim, depth, d, cStart, cEnd, c; 335 PetscErrorCode ierr; 336 337 PetscFunctionBegin; 338 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 339 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 340 for (d = 0; d < dim; ++d) { 341 refVol *= 2.0; 342 } 343 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 344 for (c = cStart; c < cEnd; ++c) { 345 if (user->femGeometry) { 346 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, NULL, J, NULL, &detJ);CHKERRQ(ierr); 347 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mesh cell %D is inverted, |J| = %g", c, (double) detJ); 348 ierr = PetscInfo1(dm, "FEM Volume: %g\n", (double) detJ*refVol);CHKERRQ(ierr); 349 } 350 if (depth > 1) { 351 ierr = DMPlexComputeCellGeometryFVM(dm, c, &vol, NULL, NULL);CHKERRQ(ierr); 352 if (vol <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mesh cell %D is inverted, vol = %g", c, (double) vol); 353 ierr = PetscInfo1(dm, "FVM Volume: %g\n", (double) vol);CHKERRQ(ierr); 354 } 355 } 356 PetscFunctionReturn(0); 357 } 358 359 PetscErrorCode CompareCones(DM dm, DM idm) 360 { 361 PetscInt cStart, cEnd, c, vStart, vEnd, v; 362 PetscErrorCode ierr; 363 364 PetscFunctionBegin; 365 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 366 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 367 for (c = cStart; c < cEnd; ++c) { 368 const PetscInt *cone; 369 PetscInt *points = NULL, numPoints, p, numVertices = 0, coneSize; 370 371 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 372 ierr = DMPlexGetConeSize(dm, c, &coneSize);CHKERRQ(ierr); 373 ierr = DMPlexGetTransitiveClosure(idm, c, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr); 374 for (p = 0; p < numPoints*2; p += 2) { 375 const PetscInt point = points[p]; 376 if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point; 377 } 378 if (numVertices != coneSize) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In cell %d, cone size %d != %d vertices in closure", c, coneSize, numVertices); 379 for (v = 0; v < numVertices; ++v) { 380 if (cone[v] != points[v]) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In cell %d, cone point %d is %d != %d vertex in closure", c, v, cone[v], points[v]); 381 } 382 ierr = DMPlexRestoreTransitiveClosure(idm, c, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr); 383 } 384 PetscFunctionReturn(0); 385 } 386 387 PetscErrorCode CreateMesh(MPI_Comm comm, PetscInt testNum, AppCtx *user, DM *dm) 388 { 389 PetscBool useGenerator = user->useGenerator; 390 PetscErrorCode ierr; 391 392 PetscFunctionBegin; 393 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 394 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 395 if (useGenerator) { 396 ierr = DMSetFromOptions(*dm);CHKERRQ(ierr); 397 } else { 398 PetscInt dim = user->dim; 399 PetscBool simplex = user->simplex; 400 401 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 402 switch (dim) { 403 case 2: 404 if (simplex) { 405 ierr = CreateSimplex_2D(comm, *dm);CHKERRQ(ierr); 406 } else { 407 ierr = CreateQuad_2D(comm, testNum, *dm);CHKERRQ(ierr); 408 } 409 break; 410 case 3: 411 if (simplex) { 412 ierr = CreateSimplex_3D(comm, *dm);CHKERRQ(ierr); 413 } else { 414 ierr = CreateHex_3D(comm, *dm);CHKERRQ(ierr); 415 } 416 break; 417 default: 418 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %d", dim); 419 } 420 } 421 { 422 DMPlexInterpolatedFlag interpolated; 423 424 ierr = CheckMesh(*dm, user);CHKERRQ(ierr); 425 ierr = DMPlexIsInterpolated(*dm, &interpolated);CHKERRQ(ierr); 426 if (interpolated == DMPLEX_INTERPOLATED_FULL) { 427 DM udm; 428 429 ierr = DMPlexUninterpolate(*dm, &udm);CHKERRQ(ierr); 430 ierr = CompareCones(udm, *dm);CHKERRQ(ierr); 431 ierr = DMDestroy(&udm);CHKERRQ(ierr); 432 } else { 433 DM idm; 434 435 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 436 ierr = CompareCones(*dm, idm);CHKERRQ(ierr); 437 ierr = DMDestroy(dm);CHKERRQ(ierr); 438 *dm = idm; 439 } 440 } 441 { 442 DM distributedMesh = NULL; 443 PetscPartitioner part; 444 445 ierr = DMPlexGetPartitioner(*dm, &part);CHKERRQ(ierr); 446 ierr = PetscPartitionerSetFromOptions(part);CHKERRQ(ierr); 447 /* Distribute mesh over processes */ 448 ierr = DMPlexDistribute(*dm, 0, NULL, &distributedMesh);CHKERRQ(ierr); 449 if (distributedMesh) { 450 ierr = DMViewFromOptions(distributedMesh, NULL, "-dm_view");CHKERRQ(ierr); 451 ierr = DMDestroy(dm);CHKERRQ(ierr); 452 *dm = distributedMesh; 453 } 454 } 455 ierr = PetscObjectSetName((PetscObject) *dm, "Interpolated Mesh");CHKERRQ(ierr); 456 ierr = DMSetFromOptions(*dm);CHKERRQ(ierr); 457 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 458 PetscFunctionReturn(0); 459 } 460 461 int main(int argc, char **argv) 462 { 463 DM dm; 464 AppCtx user; /* user-defined work context */ 465 PetscErrorCode ierr; 466 467 ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr; 468 ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr); 469 ierr = CreateMesh(PETSC_COMM_WORLD, user.testNum, &user, &dm);CHKERRQ(ierr); 470 ierr = CheckMesh(dm, &user);CHKERRQ(ierr); 471 ierr = DMDestroy(&dm);CHKERRQ(ierr); 472 ierr = PetscFinalize(); 473 return ierr; 474 } 475 476 /*TEST 477 478 # Two cell test meshes 0-7 479 test: 480 suffix: 0 481 args: -dim 2 -dm_view ascii::ascii_info_detail 482 test: 483 suffix: 1 484 nsize: 2 485 args: -petscpartitioner_type simple -dim 2 -dm_view ascii::ascii_info_detail 486 test: 487 suffix: 2 488 args: -dim 2 -simplex 0 -dm_view ascii::ascii_info_detail 489 test: 490 suffix: 3 491 nsize: 2 492 args: -petscpartitioner_type simple -dim 2 -simplex 0 -dm_view ascii::ascii_info_detail 493 test: 494 suffix: 4 495 args: -dim 3 -dm_view ascii::ascii_info_detail 496 test: 497 suffix: 5 498 nsize: 2 499 args: -petscpartitioner_type simple -dim 3 -dm_view ascii::ascii_info_detail 500 test: 501 suffix: 6 502 args: -dim 3 -simplex 0 -dm_view ascii::ascii_info_detail 503 test: 504 suffix: 7 505 nsize: 2 506 args: -petscpartitioner_type simple -dim 3 -simplex 0 -dm_view ascii::ascii_info_detail 507 # 2D Hybrid Mesh 8 508 test: 509 suffix: 8 510 args: -dim 2 -simplex 0 -testnum 1 -dm_view ascii::ascii_info_detail 511 # Reference cells 512 test: 513 suffix: 12 514 args: -use_generator -dm_plex_reference_cell_domain -dm_plex_cell pyramid -fem_geometry 0 -dm_plex_check_all 515 # TetGen meshes 9-10 516 test: 517 suffix: 9 518 requires: triangle 519 args: -use_generator -dm_view ascii::ascii_info_detail -dm_coord_space 0 520 test: 521 suffix: 10 522 requires: ctetgen 523 args: -use_generator -dm_plex_dim 3 -dm_view ascii::ascii_info_detail -dm_coord_space 0 524 # Cubit meshes 11 525 test: 526 suffix: 11 527 requires: exodusii 528 args: -use_generator -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/blockcylinder-50.exo -dm_view ascii::ascii_info_detail -dm_coord_space 0 529 530 TEST*/ 531