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