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