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