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