1c4762a1bSJed Brown static char help[] = "Tests for creation of hybrid meshes\n\n"; 2c4762a1bSJed Brown 3c4762a1bSJed Brown /* TODO 4c4762a1bSJed Brown - Propagate hybridSize with distribution 5c4762a1bSJed Brown - Test with multiple fault segments 6c4762a1bSJed Brown - Test with embedded fault 7c4762a1bSJed Brown - Test with multiple faults 8c4762a1bSJed Brown - Move over all PyLith tests? 9c4762a1bSJed Brown */ 10c4762a1bSJed Brown 11c4762a1bSJed Brown #include <petscdmplex.h> 12b253942bSMatthew G. Knepley #include <petscds.h> 13b253942bSMatthew G. Knepley #include <petsc/private/dmpleximpl.h> 14c4762a1bSJed Brown /* List of test meshes 15c4762a1bSJed Brown 16c4762a1bSJed Brown Triangle 17c4762a1bSJed Brown -------- 18c4762a1bSJed Brown Test 0: 19c4762a1bSJed Brown Two triangles sharing a face 20c4762a1bSJed Brown 21c4762a1bSJed Brown 4 22c4762a1bSJed Brown / | \ 23c4762a1bSJed Brown 8 | 9 24c4762a1bSJed Brown / | \ 25c4762a1bSJed Brown 2 0 7 1 5 26c4762a1bSJed Brown \ | / 27c4762a1bSJed Brown 6 | 10 28c4762a1bSJed Brown \ | / 29c4762a1bSJed Brown 3 30c4762a1bSJed Brown 31c4762a1bSJed Brown should become two triangles separated by a zero-volume cell with 4 vertices 32c4762a1bSJed Brown 33c4762a1bSJed Brown 5--16--8 4--12--6 3 34c4762a1bSJed Brown / | | \ / | | | \ 35c4762a1bSJed Brown 11 | | 12 9 | | | 4 36c4762a1bSJed Brown / | | \ / | | | \ 37c4762a1bSJed Brown 3 0 10 2 14 1 6 2 0 8 1 10 6 0 1 38c4762a1bSJed Brown \ | | / \ | | | / 39c4762a1bSJed Brown 9 | | 13 7 | | | 5 40c4762a1bSJed Brown \ | | / \ | | | / 41c4762a1bSJed Brown 4--15--7 3--11--5 2 42c4762a1bSJed Brown 43c4762a1bSJed Brown Test 1: 44c4762a1bSJed Brown Four triangles sharing two faces which are oriented against each other 45c4762a1bSJed Brown 46c4762a1bSJed Brown 9 47c4762a1bSJed Brown / \ 48c4762a1bSJed Brown / \ 49c4762a1bSJed Brown 17 2 16 50c4762a1bSJed Brown / \ 51c4762a1bSJed Brown / \ 52c4762a1bSJed Brown 8-----15----5 53c4762a1bSJed Brown \ /|\ 54c4762a1bSJed Brown \ / | \ 55c4762a1bSJed Brown 18 3 12 | 14 56c4762a1bSJed Brown \ / | \ 57c4762a1bSJed Brown \ / | \ 58c4762a1bSJed Brown 4 0 11 1 7 59c4762a1bSJed Brown \ | / 60c4762a1bSJed Brown \ | / 61c4762a1bSJed Brown 10 | 13 62c4762a1bSJed Brown \ | / 63c4762a1bSJed Brown \|/ 64c4762a1bSJed Brown 6 65c4762a1bSJed Brown 66c4762a1bSJed Brown Fault mesh 67c4762a1bSJed Brown 68c4762a1bSJed Brown 0 --> 0 69c4762a1bSJed Brown 1 --> 1 70c4762a1bSJed Brown 2 --> 2 71c4762a1bSJed Brown 3 --> 3 72c4762a1bSJed Brown 4 --> 5 73c4762a1bSJed Brown 5 --> 6 74c4762a1bSJed Brown 6 --> 8 75c4762a1bSJed Brown 7 --> 11 76c4762a1bSJed Brown 8 --> 15 77c4762a1bSJed Brown 78c4762a1bSJed Brown 2 79c4762a1bSJed Brown | 80c4762a1bSJed Brown 6----8----4 81c4762a1bSJed Brown | | 82c4762a1bSJed Brown 3 | 83c4762a1bSJed Brown 0-7-1 84c4762a1bSJed Brown | 85c4762a1bSJed Brown | 86c4762a1bSJed Brown 5 87c4762a1bSJed Brown 88c4762a1bSJed Brown should become four triangles separated by two zero-volume cells with 4 vertices 89c4762a1bSJed Brown 90c4762a1bSJed Brown 11 91c4762a1bSJed Brown / \ 92c4762a1bSJed Brown / \ 93c4762a1bSJed Brown / \ 94c4762a1bSJed Brown 22 2 21 95c4762a1bSJed Brown / \ 96c4762a1bSJed Brown / \ 97c4762a1bSJed Brown 10-----20------7 98c4762a1bSJed Brown 28 | 5 26/ \ 99c4762a1bSJed Brown 14----25----12 \ 100c4762a1bSJed Brown \ /| |\ 101c4762a1bSJed Brown \ / | | \ 102c4762a1bSJed Brown 23 3 17 | | 19 103c4762a1bSJed Brown \ / | | \ 104c4762a1bSJed Brown \ / | | \ 105c4762a1bSJed Brown 6 0 24 4 16 1 9 106c4762a1bSJed Brown \ | | / 107c4762a1bSJed Brown \ | | / 108c4762a1bSJed Brown 15 | | 18 109c4762a1bSJed Brown \ | | / 110c4762a1bSJed Brown \| |/ 111c4762a1bSJed Brown 13---8 112c4762a1bSJed Brown 27 113c4762a1bSJed Brown 114c4762a1bSJed Brown Tetrahedron 115c4762a1bSJed Brown ----------- 116c4762a1bSJed Brown Test 0: 117c4762a1bSJed Brown Two tets sharing a face 118c4762a1bSJed Brown 119c4762a1bSJed Brown cell 5 _______ cell 120c4762a1bSJed Brown 0 / | \ \ 1 121c4762a1bSJed Brown 16 | 18 22 122c4762a1bSJed Brown /8 19 10\ \ 123c4762a1bSJed Brown 2-15-|----4--21--6 124c4762a1bSJed Brown \ 9| 7 / / 125c4762a1bSJed Brown 14 | 17 20 126c4762a1bSJed Brown \ | / / 127c4762a1bSJed Brown 3------- 128c4762a1bSJed Brown 129c4762a1bSJed Brown should become two tetrahedrons separated by a zero-volume cell with 3 faces/3 edges/6 vertices 130c4762a1bSJed Brown 131c4762a1bSJed Brown cell 6 ___36___10______ cell 132c4762a1bSJed Brown 0 / | \ |\ \ 1 133c4762a1bSJed Brown 24 | 26 | 32 30 134c4762a1bSJed Brown /12 27 14\ 33 \ \ 135c4762a1bSJed Brown 3-23-|----5--35-|---9--29--7 136c4762a1bSJed Brown \ 13| 11/ |18 / / 137c4762a1bSJed Brown 22 | 25 | 31 28 138c4762a1bSJed Brown \ | / |/ / 139c4762a1bSJed Brown 4----34----8------ 140c4762a1bSJed Brown cell 2 141c4762a1bSJed Brown 142c4762a1bSJed Brown In parallel, 143c4762a1bSJed Brown 144c4762a1bSJed Brown cell 5 ___28____8 4______ cell 145c4762a1bSJed Brown 0 / | \ |\ |\ \ 0 146c4762a1bSJed Brown 19 | 21 | 24 | 13 6 11 147c4762a1bSJed Brown /10 22 12\ 25 \ |8 \ \ 148c4762a1bSJed Brown 2-18-|----4--27-|---7 14 3--10--1 149c4762a1bSJed Brown \ 11| 9 / |13 / | / / 150c4762a1bSJed Brown 17 | 20 | 23 | 12 5 9 151c4762a1bSJed Brown \ | / |/ |/ / 152c4762a1bSJed Brown 3----26----6 2------ 153c4762a1bSJed Brown cell 1 154c4762a1bSJed Brown 155c4762a1bSJed Brown Test 1: 156c4762a1bSJed Brown Four tets sharing two faces 157c4762a1bSJed Brown 158c4762a1bSJed Brown Cells: 0-3,4-5 159c4762a1bSJed Brown Vertices: 6-15 160c4762a1bSJed Brown Faces: 16-29,30-34 161c4762a1bSJed Brown Edges: 35-52,53-56 162c4762a1bSJed Brown 163c4762a1bSJed Brown Quadrilateral 164c4762a1bSJed Brown ------------- 165c4762a1bSJed Brown Test 0: 166c4762a1bSJed Brown Two quads sharing a face 167c4762a1bSJed Brown 168c4762a1bSJed Brown 5--10---4--14---7 169c4762a1bSJed Brown | | | 170c4762a1bSJed Brown 11 0 9 1 13 171c4762a1bSJed Brown | | | 172c4762a1bSJed Brown 2---8---3--12---6 173c4762a1bSJed Brown 174c4762a1bSJed Brown should become two quads separated by a zero-volume cell with 4 vertices 175c4762a1bSJed Brown 176c4762a1bSJed Brown 6--13---5-20-10--17---8 5--10---4-14--7 4---7---2 177c4762a1bSJed Brown | | | | | | | | | 178c4762a1bSJed Brown 14 0 12 2 18 1 16 11 0 9 1 12 8 0 6 179c4762a1bSJed Brown | | | | | | | | | 180c4762a1bSJed Brown 3--11---4-19--9--15---7 2---8---3-13--6 3---5---1 181c4762a1bSJed Brown 182c4762a1bSJed Brown Test 1: 183c4762a1bSJed Brown 184c4762a1bSJed Brown Original mesh with 9 cells, 185c4762a1bSJed Brown 186c4762a1bSJed Brown 9 ----10 ----11 ----12 187c4762a1bSJed Brown | | | | 188c4762a1bSJed Brown | | | | 189c4762a1bSJed Brown | | | | 190c4762a1bSJed Brown | | | | 191c4762a1bSJed Brown 13 ----14 ----15 ----16 192c4762a1bSJed Brown | | | | 193c4762a1bSJed Brown | | | | 194c4762a1bSJed Brown | | | | 195c4762a1bSJed Brown | | | | 196c4762a1bSJed Brown 17 ----18 ----19 ----20 197c4762a1bSJed Brown | | | | 198c4762a1bSJed Brown | | | | 199c4762a1bSJed Brown | | | | 200c4762a1bSJed Brown | | | | 201c4762a1bSJed Brown 21 ----22 ----23 ----24 202c4762a1bSJed Brown 203c4762a1bSJed Brown After first fault, 204c4762a1bSJed Brown 205c4762a1bSJed Brown 12 ----13 ----14-28 ----15 206c4762a1bSJed Brown | | | | | 207c4762a1bSJed Brown | 0 | 1 | 9| 2 | 208c4762a1bSJed Brown | | | | | 209c4762a1bSJed Brown | | | | | 210c4762a1bSJed Brown 16 ----17 ----18-29 ----19 211c4762a1bSJed Brown | | | | | 212c4762a1bSJed Brown | 3 | 4 |10| 5 | 213c4762a1bSJed Brown | | | | | 214c4762a1bSJed Brown | | | | | 215c4762a1bSJed Brown 20 ----21-----22-30 ----23 216c4762a1bSJed Brown | | | \--11- | 217c4762a1bSJed Brown | 6 | 7 | 8 | 218c4762a1bSJed Brown | | | | 219c4762a1bSJed Brown | | | | 220c4762a1bSJed Brown 24 ----25 ----26--------27 221c4762a1bSJed Brown 222c4762a1bSJed Brown After second fault, 223c4762a1bSJed Brown 224c4762a1bSJed Brown 14 ----15 ----16-30 ----17 225c4762a1bSJed Brown | | | | | 226c4762a1bSJed Brown | 0 | 1 | 9| 2 | 227c4762a1bSJed Brown | | | | | 228c4762a1bSJed Brown | | | | | 229c4762a1bSJed Brown 18 ----19 ----20-31 ----21 230c4762a1bSJed Brown | | | | | 231c4762a1bSJed Brown | 3 | 4 |10| 5 | 232c4762a1bSJed Brown | | | | | 233c4762a1bSJed Brown | | | | | 234c4762a1bSJed Brown 33 ----34-----24-32 ----25 235c4762a1bSJed Brown | 12 | 13 / | \-11-- | 236c4762a1bSJed Brown 22 ----23---/ | | 237c4762a1bSJed Brown | | 7 | 8 | 238c4762a1bSJed Brown | 6 | | | 239c4762a1bSJed Brown | | | | 240c4762a1bSJed Brown | | | | 241c4762a1bSJed Brown 26 ----27 ----28--------29 242c4762a1bSJed Brown 243c4762a1bSJed Brown Hexahedron 244c4762a1bSJed Brown ---------- 245c4762a1bSJed Brown Test 0: 246c4762a1bSJed Brown Two hexes sharing a face 247c4762a1bSJed Brown 248c4762a1bSJed Brown cell 9-----31------8-----42------13 cell 249c4762a1bSJed Brown 0 /| /| /| 1 250c4762a1bSJed Brown 32 | 15 30| 21 41| 251c4762a1bSJed Brown / | / | / | 252c4762a1bSJed Brown 6-----29------7-----40------12 | 253c4762a1bSJed Brown | | 18 | | 24 | | 254c4762a1bSJed Brown | 36 | 35 | 44 255c4762a1bSJed Brown |19 | |17 | |23 | 256c4762a1bSJed Brown 33 | 16 34 | 22 43 | 257c4762a1bSJed Brown | 5-----27--|---4-----39--|---11 258c4762a1bSJed Brown | / | / | / 259c4762a1bSJed Brown | 28 14 | 26 20 | 38 260c4762a1bSJed Brown |/ |/ |/ 261c4762a1bSJed Brown 2-----25------3-----37------10 262c4762a1bSJed Brown 263c4762a1bSJed Brown should become two hexes separated by a zero-volume cell with 8 vertices 264c4762a1bSJed Brown 265c4762a1bSJed Brown cell 2 266c4762a1bSJed Brown cell 10-----41------9-----62------18----52------14 cell 267c4762a1bSJed Brown 0 /| /| /| /| 1 268c4762a1bSJed Brown 42 | 20 40| 32 56| 26 51| 269c4762a1bSJed Brown / | / | / | / | 270c4762a1bSJed Brown 7-----39------8-----61------17--|-50------13 | 271c4762a1bSJed Brown | | 23 | | | | 29 | | 272c4762a1bSJed Brown | 46 | 45 | 58 | 54 273c4762a1bSJed Brown |24 | |22 | |30 | |28 | 274c4762a1bSJed Brown 43 | 21 44 | 57 | 27 53 | 275c4762a1bSJed Brown | 6-----37--|---5-----60--|---16----49--|---12 276c4762a1bSJed Brown | / | / | / | / 277c4762a1bSJed Brown | 38 19 | 36 31 | 55 25 | 48 278c4762a1bSJed Brown |/ |/ |/ |/ 279c4762a1bSJed Brown 3-----35------4-----59------15----47------11 280c4762a1bSJed Brown 281c4762a1bSJed Brown In parallel, 282c4762a1bSJed Brown 283c4762a1bSJed Brown cell 2 284c4762a1bSJed Brown cell 9-----31------8-----44------13 8----20------4 cell 285c4762a1bSJed Brown 0 /| /| /| /| /| 1 286c4762a1bSJed Brown 32 | 15 30| 22 38| 24 | 10 19| 287c4762a1bSJed Brown / | / | / | / | / | 288c4762a1bSJed Brown 6-----29------7-----43------12 | 7----18------3 | 289c4762a1bSJed Brown | | 18 | | | | | | 13 | | 290c4762a1bSJed Brown | 36 | 35 | 40 | 26 | 22 291c4762a1bSJed Brown |19 | |17 | |20 | |14 | |12 | 292c4762a1bSJed Brown 33 | 16 34 | 39 | 25 | 11 21 | 293c4762a1bSJed Brown | 5-----27--|---4-----42--|---11 | 6----17--|---2 294c4762a1bSJed Brown | / | / | / | / | / 295c4762a1bSJed Brown | 28 14 | 26 21 | 37 |23 9 | 16 296c4762a1bSJed Brown |/ |/ |/ |/ |/ 297c4762a1bSJed Brown 2-----25------3-----41------10 5----15------1 298c4762a1bSJed Brown 299c4762a1bSJed Brown Test 1: 300c4762a1bSJed Brown 301c4762a1bSJed Brown */ 302c4762a1bSJed Brown 303c4762a1bSJed Brown typedef struct { 304c4762a1bSJed Brown PetscInt debug; /* The debugging level */ 305c4762a1bSJed Brown PetscInt dim; /* The topological mesh dimension */ 306c4762a1bSJed Brown PetscBool cellSimplex; /* Use simplices or hexes */ 307c4762a1bSJed Brown PetscBool testPartition; /* Use a fixed partitioning for testing */ 308c4762a1bSJed Brown PetscInt testNum; /* The particular mesh to test */ 309b7519becSMatthew G. Knepley PetscInt cohesiveFields; /* The number of cohesive fields */ 310c4762a1bSJed Brown } AppCtx; 311c4762a1bSJed Brown 312b253942bSMatthew G. Knepley static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 313c4762a1bSJed Brown { 314c4762a1bSJed Brown PetscErrorCode ierr; 315c4762a1bSJed Brown 316c4762a1bSJed Brown PetscFunctionBegin; 317c4762a1bSJed Brown options->debug = 0; 318c4762a1bSJed Brown options->dim = 2; 319c4762a1bSJed Brown options->cellSimplex = PETSC_TRUE; 320c4762a1bSJed Brown options->testPartition = PETSC_TRUE; 321c4762a1bSJed Brown options->testNum = 0; 322b7519becSMatthew G. Knepley options->cohesiveFields = 1; 323c4762a1bSJed Brown 324c4762a1bSJed Brown ierr = PetscOptionsBegin(comm, "", "Meshing Problem Options", "DMPLEX");CHKERRQ(ierr); 325*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscOptionsBoundedInt("-debug", "The debugging level", "ex5.c", options->debug, &options->debug, NULL,0)); 326*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscOptionsRangeInt("-dim", "The topological mesh dimension", "ex5.c", options->dim, &options->dim, NULL,1,3)); 327*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscOptionsBool("-cell_simplex", "Use simplices if true, otherwise hexes", "ex5.c", options->cellSimplex, &options->cellSimplex, NULL)); 328*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscOptionsBool("-test_partition", "Use a fixed partition for testing", "ex5.c", options->testPartition, &options->testPartition, NULL)); 329*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscOptionsBoundedInt("-test_num", "The particular mesh to test", "ex5.c", options->testNum, &options->testNum, NULL,0)); 330*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscOptionsBoundedInt("-cohesive_fields", "The number of cohesive fields", "ex5.c", options->cohesiveFields, &options->cohesiveFields, NULL, 0)); 3311e1ea65dSPierre Jolivet ierr = PetscOptionsEnd();CHKERRQ(ierr); 332c4762a1bSJed Brown PetscFunctionReturn(0); 333c4762a1bSJed Brown } 334c4762a1bSJed Brown 335b253942bSMatthew G. Knepley static PetscErrorCode CreateSimplex_2D(MPI_Comm comm, PetscInt testNum, DM *dm) 336c4762a1bSJed Brown { 337c4762a1bSJed Brown DM idm; 338c4762a1bSJed Brown PetscInt p; 339c4762a1bSJed Brown PetscMPIInt rank; 340c4762a1bSJed Brown 341c4762a1bSJed Brown PetscFunctionBegin; 342*5f80ce2aSJacob Faibussowitsch CHKERRMPI(MPI_Comm_rank(comm, &rank)); 343dd400576SPatrick Sanan if (rank == 0) { 344c4762a1bSJed Brown switch (testNum) { 345c4762a1bSJed Brown case 0: 346c4762a1bSJed Brown { 347c4762a1bSJed Brown PetscInt numPoints[2] = {4, 2}; 348c4762a1bSJed Brown PetscInt coneSize[6] = {3, 3, 0, 0, 0, 0}; 349c4762a1bSJed Brown PetscInt cones[6] = {2, 3, 4, 5, 4, 3}; 350c4762a1bSJed Brown PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 351c4762a1bSJed Brown PetscScalar vertexCoords[8] = {-0.5, 0.5, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5}; 352c4762a1bSJed Brown PetscInt markerPoints[8] = {2, 1, 3, 1, 4, 1, 5, 1}; 353c4762a1bSJed Brown PetscInt faultPoints[2] = {3, 4}; 354c4762a1bSJed Brown 355*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 356*5f80ce2aSJacob Faibussowitsch for (p = 0; p < 4; ++p) CHKERRQ(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1])); 357*5f80ce2aSJacob Faibussowitsch for (p = 0; p < 2; ++p) CHKERRQ(DMSetLabelValue(*dm, "fault", faultPoints[p], 1)); 358*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetLabelValue(*dm, "material", 0, 1)); 359*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetLabelValue(*dm, "material", 1, 2)); 360c4762a1bSJed Brown } 361c4762a1bSJed Brown break; 362c4762a1bSJed Brown case 1: 363c4762a1bSJed Brown { 364c4762a1bSJed Brown PetscInt numPoints[2] = {6, 4}; 365c4762a1bSJed Brown PetscInt coneSize[10] = {3, 3, 3, 3, 0, 0, 0, 0, 0, 0}; 366c4762a1bSJed Brown PetscInt cones[12] = {4, 6, 5, 5, 6, 7, 8, 5, 9, 8, 4, 5}; 367c4762a1bSJed Brown PetscInt coneOrientations[12] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 368c4762a1bSJed Brown PetscScalar vertexCoords[12] = {-1.0, 0.0, 0.0, 1.0, 0.0, -1.0, 1.0, 0.0, -2.0, 1.0, -1.0, 2.0}; 369c4762a1bSJed Brown PetscInt markerPoints[6] = {4, 1, 6, 1, 8, 1}; 370c4762a1bSJed Brown PetscInt faultPoints[3] = {5, 6, 8}; 371c4762a1bSJed Brown 372*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 373*5f80ce2aSJacob Faibussowitsch for (p = 0; p < 3; ++p) CHKERRQ(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1])); 374*5f80ce2aSJacob Faibussowitsch for (p = 0; p < 3; ++p) CHKERRQ(DMSetLabelValue(*dm, "fault", faultPoints[p], 1)); 375*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetLabelValue(*dm, "material", 0, 1));CHKERRQ(DMSetLabelValue(*dm, "material", 3, 1)); 376*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetLabelValue(*dm, "material", 1, 2));CHKERRQ(DMSetLabelValue(*dm, "material", 2, 2)); 377c4762a1bSJed Brown } 378c4762a1bSJed Brown break; 379c4762a1bSJed Brown default: 38098921bdaSJacob Faibussowitsch SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %d", testNum); 381c4762a1bSJed Brown } 382c4762a1bSJed Brown } else { 383c4762a1bSJed Brown PetscInt numPoints[3] = {0, 0, 0}; 384c4762a1bSJed Brown 385*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexCreateFromDAG(*dm, 1, numPoints, NULL, NULL, NULL, NULL)); 386*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMCreateLabel(*dm, "fault")); 387c4762a1bSJed Brown } 388*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexInterpolate(*dm, &idm)); 389*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMViewFromOptions(idm, NULL, "-in_dm_view")); 390*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMDestroy(dm)); 391c4762a1bSJed Brown *dm = idm; 392c4762a1bSJed Brown PetscFunctionReturn(0); 393c4762a1bSJed Brown } 394c4762a1bSJed Brown 395b253942bSMatthew G. Knepley static PetscErrorCode CreateSimplex_3D(MPI_Comm comm, AppCtx *user, DM dm) 396c4762a1bSJed Brown { 397c4762a1bSJed Brown PetscInt depth = 3, testNum = user->testNum, p; 398c4762a1bSJed Brown PetscMPIInt rank; 399c4762a1bSJed Brown 400c4762a1bSJed Brown PetscFunctionBegin; 401*5f80ce2aSJacob Faibussowitsch CHKERRMPI(MPI_Comm_rank(comm, &rank)); 402dd400576SPatrick Sanan if (rank == 0) { 403c4762a1bSJed Brown switch (testNum) { 404c4762a1bSJed Brown case 0: 405c4762a1bSJed Brown { 406c4762a1bSJed Brown PetscInt numPoints[4] = {5, 7, 9, 2}; 407c4762a1bSJed Brown PetscInt coneSize[23] = {4, 4, 0, 0, 0, 0, 0, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2}; 408c4762a1bSJed Brown PetscInt cones[47] = {7, 8, 9, 10, 11, 10, 13, 12, 15, 17, 14, 16, 18, 15, 14, 19, 16, 17, 18, 19, 17, 21, 20, 18, 22, 21, 22, 19, 20, 2, 3, 2, 4, 2, 5, 3, 4, 4, 5, 5, 3, 3, 6, 4, 6, 5, 6}; 409b5a892a1SMatthew G. Knepley PetscInt coneOrientations[47] = {0, 0, 0, 0, 0, -2, 2, 2, 0, -1, -1, 0, -1, -1, 0, -1, -1, 0, 0, 0, 0, 0, -1, 0, 0, -1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 410c4762a1bSJed 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}; 411c4762a1bSJed Brown PetscInt markerPoints[20] = {2, 1, 3, 1, 4, 1, 5, 1, 14, 1, 15, 1, 16, 1, 17, 1, 18, 1, 19, 1}; 412c4762a1bSJed Brown PetscInt faultPoints[3] = {3, 4, 5}; 413c4762a1bSJed Brown 414*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 415c4762a1bSJed Brown for (p = 0; p < 10; ++p) { 416*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1])); 417c4762a1bSJed Brown } 418c4762a1bSJed Brown for (p = 0; p < 3; ++p) { 419*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetLabelValue(dm, "fault", faultPoints[p], 1)); 420c4762a1bSJed Brown } 421*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetLabelValue(dm, "material", 0, 1)); 422*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetLabelValue(dm, "material", 1, 2)); 423c4762a1bSJed Brown } 424c4762a1bSJed Brown break; 425c4762a1bSJed Brown case 1: 426c4762a1bSJed Brown { 427c4762a1bSJed Brown PetscInt numPoints[4] = {6, 13, 12, 4}; 428c4762a1bSJed Brown PetscInt coneSize[35] = {4, 4, 4, 4, 0, 0, 0, 0, 0, 0, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2}; 429c4762a1bSJed Brown PetscInt cones[78] = {10, 11, 12, 13, 10, 15, 16, 14, 17, 18, 14, 19, 20, 13, 19, 21, 22, 23, 24, 25, 26, 22, 24, 27, 25, 23, 26, 27, 28, 29, 23, 24, 30, 28, 22, 29, 30, 31, 32, 28, 29, 33, 31, 32, 33, 23, 26, 34, 33, 34, 27, 32, 6, 5, 5, 7, 7, 6, 6, 4, 4, 5, 7, 4, 7, 9, 9, 5, 6, 9, 9, 8, 8, 7, 5, 8, 4, 8}; 430b5a892a1SMatthew G. Knepley PetscInt coneOrientations[78] = { 0, 0, 0, 0, -2, 1, 0, 2, 0, 0, -3, 0, 0, -3, -1, 0, 0, 0, 0, 0, 0, -1, -1, 0, -1, -1, -1, -1, 0, 0, 0, 0, 0, -1, 0, -1, -1, 0, 0, 0, 0, 0, -1, -1, -1, 0, -1, 0, -1, -1, -1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 431c4762a1bSJed Brown PetscScalar vertexCoords[18] = {-1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, -1.0, 1.0, 0.0, 0.0}; 432c4762a1bSJed Brown PetscInt markerPoints[14] = {5, 1, 6, 1, 7, 1, 10, 1, 22, 1, 23, 1, 24, 1}; 433c4762a1bSJed Brown PetscInt faultPoints[4] = {5, 6, 7, 8}; 434c4762a1bSJed Brown 435*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 436c4762a1bSJed Brown for (p = 0; p < 7; ++p) { 437*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1])); 438c4762a1bSJed Brown } 439c4762a1bSJed Brown for (p = 0; p < 4; ++p) { 440*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetLabelValue(dm, "fault", faultPoints[p], 1)); 441c4762a1bSJed Brown } 442*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetLabelValue(dm, "material", 0, 1)); 443*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetLabelValue(dm, "material", 1, 2)); 444c4762a1bSJed Brown } 445c4762a1bSJed Brown break; 446c4762a1bSJed Brown default: 44798921bdaSJacob Faibussowitsch SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %d", testNum); 448c4762a1bSJed Brown } 449c4762a1bSJed Brown } else { 450c4762a1bSJed Brown PetscInt numPoints[4] = {0, 0, 0, 0}; 451c4762a1bSJed Brown 452*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexCreateFromDAG(dm, depth, numPoints, NULL, NULL, NULL, NULL)); 453*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMCreateLabel(dm, "fault")); 454c4762a1bSJed Brown } 455c4762a1bSJed Brown PetscFunctionReturn(0); 456c4762a1bSJed Brown } 457c4762a1bSJed Brown 458b253942bSMatthew G. Knepley static PetscErrorCode CreateQuad_2D(MPI_Comm comm, PetscInt testNum, DM *dm) 459c4762a1bSJed Brown { 460c4762a1bSJed Brown DM idm; 461c4762a1bSJed Brown PetscInt p; 462c4762a1bSJed Brown PetscMPIInt rank; 463c4762a1bSJed Brown 464c4762a1bSJed Brown PetscFunctionBegin; 465*5f80ce2aSJacob Faibussowitsch CHKERRMPI(MPI_Comm_rank(comm, &rank)); 466dd400576SPatrick Sanan if (rank == 0) { 467c4762a1bSJed Brown switch (testNum) { 468c4762a1bSJed Brown case 0: 469c4762a1bSJed Brown case 2: 470c4762a1bSJed Brown { 471c4762a1bSJed Brown PetscInt numPoints[2] = {6, 2}; 472c4762a1bSJed Brown PetscInt coneSize[8] = {4, 4, 0, 0, 0, 0, 0, 0}; 473c4762a1bSJed Brown PetscInt cones[8] = {2, 3, 4, 5, 3, 6, 7, 4}; 474c4762a1bSJed Brown PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 475c4762a1bSJed 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}; 476c4762a1bSJed Brown PetscInt markerPoints[12] = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1, 7, 1}; 477c4762a1bSJed Brown PetscInt faultPoints[2] = {3, 4}; 478c4762a1bSJed Brown 479*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 480*5f80ce2aSJacob Faibussowitsch for (p = 0; p < 6; ++p) CHKERRQ(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1])); 481*5f80ce2aSJacob Faibussowitsch if (testNum == 0) for (p = 0; p < 2; ++p) CHKERRQ(DMSetLabelValue(*dm, "fault", faultPoints[p], 1)); 482*5f80ce2aSJacob Faibussowitsch if (testNum == 2) for (p = 0; p < 2; ++p) CHKERRQ(DMSetLabelValue(*dm, "pfault", faultPoints[p], 1)); 483*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetLabelValue(*dm, "material", 0, 1)); 484*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetLabelValue(*dm, "material", 1, 2)); 485c4762a1bSJed Brown } 486c4762a1bSJed Brown break; 487c4762a1bSJed Brown case 1: 488c4762a1bSJed Brown { 489c4762a1bSJed Brown PetscInt numPoints[2] = {16, 9}; 490c4762a1bSJed Brown PetscInt coneSize[25] = {4, 4, 4, 4, 4, 4, 4, 4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 491c4762a1bSJed Brown PetscInt cones[36] = {9, 13, 14, 10, 492c4762a1bSJed Brown 10, 14, 15, 11, 493c4762a1bSJed Brown 11, 15, 16, 12, 494c4762a1bSJed Brown 13, 17, 18, 14, 495c4762a1bSJed Brown 14, 18, 19, 15, 496c4762a1bSJed Brown 15, 19, 20, 16, 497c4762a1bSJed Brown 17, 21, 22, 18, 498c4762a1bSJed Brown 18, 22, 23, 19, 499c4762a1bSJed Brown 19, 23, 24, 20}; 500c4762a1bSJed Brown PetscInt coneOrientations[36] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; 501c4762a1bSJed Brown PetscScalar vertexCoords[32] = {-3.0, 3.0, -1.0, 3.0, 1.0, 3.0, 3.0, 3.0, -3.0, 1.0, -1.0, 1.0, 1.0, 1.0, 3.0, 1.0, 502c4762a1bSJed Brown -3.0, -1.0, -1.0, -1.0, 1.0, -1.0, 3.0, -1.0, -3.0, -3.0, -1.0, -3.0, 1.0, -3.0, 3.0, -3.0}; 503c4762a1bSJed Brown PetscInt faultPoints[3] = {11, 15, 19}; 504c4762a1bSJed Brown PetscInt fault2Points[2] = {17, 18}; 505c4762a1bSJed Brown 506*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 507*5f80ce2aSJacob Faibussowitsch for (p = 0; p < 3; ++p) CHKERRQ(DMSetLabelValue(*dm, "fault", faultPoints[p], 1)); 508*5f80ce2aSJacob Faibussowitsch for (p = 0; p < 2; ++p) CHKERRQ(DMSetLabelValue(*dm, "fault2", fault2Points[p], 1)); 509*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetLabelValue(*dm, "material", 0, 1)); 510*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetLabelValue(*dm, "material", 1, 1)); 511*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetLabelValue(*dm, "material", 2, 1)); 512*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetLabelValue(*dm, "material", 3, 1)); 513*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetLabelValue(*dm, "material", 4, 1)); 514*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetLabelValue(*dm, "material", 5, 2)); 515*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetLabelValue(*dm, "material", 6, 2)); 516*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetLabelValue(*dm, "material", 7, 2)); 517*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetLabelValue(*dm, "material", 8, 2)); 518c4762a1bSJed Brown } 519c4762a1bSJed Brown break; 520c4762a1bSJed Brown default: 52198921bdaSJacob Faibussowitsch SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %d", testNum); 522c4762a1bSJed Brown } 523c4762a1bSJed Brown } else { 524c4762a1bSJed Brown PetscInt numPoints[3] = {0, 0, 0}; 525c4762a1bSJed Brown 526*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexCreateFromDAG(*dm, 1, numPoints, NULL, NULL, NULL, NULL)); 527*5f80ce2aSJacob Faibussowitsch if (testNum == 2) CHKERRQ(DMCreateLabel(*dm, "pfault")); 528*5f80ce2aSJacob Faibussowitsch else CHKERRQ(DMCreateLabel(*dm, "fault")); 529c4762a1bSJed Brown } 530*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexInterpolate(*dm, &idm)); 531*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMViewFromOptions(idm, NULL, "-in_dm_view")); 532*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMDestroy(dm)); 533c4762a1bSJed Brown *dm = idm; 534c4762a1bSJed Brown PetscFunctionReturn(0); 535c4762a1bSJed Brown } 536c4762a1bSJed Brown 537b253942bSMatthew G. Knepley static PetscErrorCode CreateHex_3D(MPI_Comm comm, PetscInt testNum, DM *dm) 538c4762a1bSJed Brown { 539c4762a1bSJed Brown DM idm; 540c4762a1bSJed Brown PetscInt depth = 3, p; 541c4762a1bSJed Brown PetscMPIInt rank; 542c4762a1bSJed Brown 543c4762a1bSJed Brown PetscFunctionBegin; 544*5f80ce2aSJacob Faibussowitsch CHKERRMPI(MPI_Comm_rank(comm, &rank)); 545dd400576SPatrick Sanan if (rank == 0) { 546c4762a1bSJed Brown switch (testNum) { 547c4762a1bSJed Brown case 0: 548c4762a1bSJed Brown { 549c4762a1bSJed Brown PetscInt numPoints[2] = {12, 2}; 550c4762a1bSJed Brown PetscInt coneSize[14] = {8,8, 0,0,0,0,0,0,0,0,0,0,0,0}; 551c4762a1bSJed Brown PetscInt cones[16] = {2,5,4,3,6,7,8,9, 3,4,11,10,7,12,13,8}; 552c4762a1bSJed Brown PetscInt coneOrientations[16] = {0,0,0,0,0,0,0,0, 0,0, 0,0 ,0, 0, 0,0}; 553c4762a1bSJed 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, 554c4762a1bSJed 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, 555c4762a1bSJed 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}; 556c4762a1bSJed Brown PetscInt markerPoints[52] = {2,1,3,1,4,1,5,1,6,1,7,1,8,1,9,1}; 557c4762a1bSJed Brown PetscInt faultPoints[4] = {3, 4, 7, 8}; 558c4762a1bSJed Brown 559*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 560*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexInterpolate(*dm, &idm)); 561*5f80ce2aSJacob Faibussowitsch for (p = 0; p < 8; ++p) CHKERRQ(DMSetLabelValue(idm, "marker", markerPoints[p*2], markerPoints[p*2+1])); 562*5f80ce2aSJacob Faibussowitsch for (p = 0; p < 4; ++p) CHKERRQ(DMSetLabelValue(idm, "fault", faultPoints[p], 1)); 563*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetLabelValue(*dm, "material", 0, 1)); 564*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetLabelValue(*dm, "material", 1, 2)); 565c4762a1bSJed Brown } 566c4762a1bSJed Brown break; 567c4762a1bSJed Brown case 1: 568c4762a1bSJed Brown { 569c4762a1bSJed Brown /* Cell Adjacency Graph: 570c4762a1bSJed Brown 0 -- { 8, 13, 21, 24} --> 1 571c4762a1bSJed Brown 0 -- {20, 21, 23, 24} --> 5 F 572c4762a1bSJed Brown 1 -- {10, 15, 21, 24} --> 2 573c4762a1bSJed Brown 1 -- {13, 14, 15, 24} --> 6 574c4762a1bSJed Brown 2 -- {21, 22, 24, 25} --> 4 F 575c4762a1bSJed Brown 3 -- {21, 24, 30, 35} --> 4 576c4762a1bSJed Brown 3 -- {21, 24, 28, 33} --> 5 577c4762a1bSJed Brown */ 578c4762a1bSJed Brown PetscInt numPoints[2] = {30, 7}; 579c4762a1bSJed Brown PetscInt coneSize[37] = {8,8,8,8,8,8,8, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; 580c4762a1bSJed Brown PetscInt cones[56] = { 8, 21, 20, 7, 13, 12, 23, 24, 581c4762a1bSJed Brown 14, 15, 10, 9, 13, 8, 21, 24, 582c4762a1bSJed Brown 15, 16, 11, 10, 24, 21, 22, 25, 583c4762a1bSJed Brown 30, 29, 28, 21, 35, 24, 33, 34, 584c4762a1bSJed Brown 24, 21, 30, 35, 25, 36, 31, 22, 585c4762a1bSJed Brown 27, 20, 21, 28, 32, 33, 24, 23, 586c4762a1bSJed Brown 15, 24, 13, 14, 19, 18, 17, 26}; 587c4762a1bSJed Brown PetscInt coneOrientations[56] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; 588c4762a1bSJed Brown PetscScalar vertexCoords[90] = {-2.0, -2.0, -2.0, -2.0, -1.0, -2.0, -3.0, 0.0, -2.0, -2.0, 1.0, -2.0, -2.0, 2.0, -2.0, -2.0, -2.0, 0.0, 589c4762a1bSJed Brown -2.0, -1.0, 0.0, -3.0, 0.0, 0.0, -2.0, 1.0, 0.0, -2.0, 2.0, 0.0, -2.0, -1.0, 2.0, -3.0, 0.0, 2.0, 590c4762a1bSJed Brown -2.0, 1.0, 2.0, 0.0, -2.0, -2.0, 0.0, 0.0, -2.0, 0.0, 2.0, -2.0, 0.0, -2.0, 0.0, 0.0, 0.0, 0.0, 591c4762a1bSJed Brown 0.0, 2.0, 0.0, 0.0, 0.0, 2.0, 2.0, -2.0, -2.0, 2.0, -1.0, -2.0, 3.0, 0.0, -2.0, 2.0, 1.0, -2.0, 592c4762a1bSJed Brown 2.0, 2.0, -2.0, 2.0, -2.0, 0.0, 2.0, -1.0, 0.0, 3.0, 0.0, 0.0, 2.0, 1.0, 0.0, 2.0, 2.0, 0.0}; 593c4762a1bSJed Brown PetscInt faultPoints[6] = {20, 21, 22, 23, 24, 25}; 594c4762a1bSJed Brown 595*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 596*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexInterpolate(*dm, &idm)); 597*5f80ce2aSJacob Faibussowitsch for (p = 0; p < 6; ++p) CHKERRQ(DMSetLabelValue(idm, "fault", faultPoints[p], 1)); 598*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetLabelValue(*dm, "material", 0, 1)); 599*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetLabelValue(*dm, "material", 1, 1)); 600*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetLabelValue(*dm, "material", 2, 1)); 601*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetLabelValue(*dm, "material", 3, 2)); 602*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetLabelValue(*dm, "material", 4, 2)); 603*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetLabelValue(*dm, "material", 5, 2)); 604*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetLabelValue(*dm, "material", 6, 2)); 605c4762a1bSJed Brown } 606c4762a1bSJed Brown break; 607c4762a1bSJed Brown case 2: 608c4762a1bSJed Brown { 609c4762a1bSJed Brown /* Buried fault edge */ 610c4762a1bSJed Brown PetscInt numPoints[2] = {18, 4}; 611c4762a1bSJed Brown PetscInt coneSize[22] = {8,8,8,8, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; 612c4762a1bSJed Brown PetscInt cones[32] = { 4, 5, 8, 7, 13, 16, 17, 14, 613c4762a1bSJed Brown 5, 6, 9, 8, 14, 17, 18, 15, 614c4762a1bSJed Brown 7, 8, 11, 10, 16, 19, 20, 17, 615c4762a1bSJed Brown 8, 9, 12, 11, 17, 20, 21, 18}; 616c4762a1bSJed Brown PetscInt coneOrientations[32] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; 617c4762a1bSJed Brown PetscScalar vertexCoords[54] = {-2.0, -2.0, 0.0, -2.0, 0.0, 0.0, -2.0, 2.0, 0.0, 0.0, -2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0, 0.0, 618c4762a1bSJed Brown 2.0, -2.0, 0.0, 2.0, 0.0, 0.0, 2.0, 2.0, 0.0, -2.0, -2.0, 2.0, -2.0, 0.0, 2.0, -2.0, 2.0, 2.0, 619c4762a1bSJed Brown 0.0, -2.0, 2.0, 0.0, 0.0, 2.0, 0.0, 2.0, 2.0, 2.0, -2.0, 2.0, 2.0, 0.0, 2.0, 2.0, 2.0, 2.0}; 620c4762a1bSJed Brown PetscInt faultPoints[4] = {7, 8, 16, 17}; 621c4762a1bSJed Brown 622*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 623*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexInterpolate(*dm, &idm)); 624*5f80ce2aSJacob Faibussowitsch for (p = 0; p < 4; ++p) CHKERRQ(DMSetLabelValue(idm, "fault", faultPoints[p], 1)); 625*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetLabelValue(*dm, "material", 0, 1)); 626*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetLabelValue(*dm, "material", 1, 1)); 627*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetLabelValue(*dm, "material", 2, 2)); 628*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetLabelValue(*dm, "material", 3, 2)); 629c4762a1bSJed Brown } 630c4762a1bSJed Brown break; 63198921bdaSJacob Faibussowitsch default: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %d", testNum); 632c4762a1bSJed Brown } 633c4762a1bSJed Brown } else { 634c4762a1bSJed Brown PetscInt numPoints[4] = {0, 0, 0, 0}; 635c4762a1bSJed Brown 636*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexCreateFromDAG(*dm, depth, numPoints, NULL, NULL, NULL, NULL)); 637*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexInterpolate(*dm, &idm)); 638*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMCreateLabel(idm, "fault")); 639c4762a1bSJed Brown } 640*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMViewFromOptions(idm, NULL, "-in_dm_view")); 641*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMDestroy(dm)); 642c4762a1bSJed Brown *dm = idm; 643c4762a1bSJed Brown PetscFunctionReturn(0); 644c4762a1bSJed Brown } 645c4762a1bSJed Brown 646b253942bSMatthew G. Knepley static PetscErrorCode CreateFaultLabel(DM dm) 647b253942bSMatthew G. Knepley { 648b253942bSMatthew G. Knepley DMLabel label; 649b253942bSMatthew G. Knepley PetscInt dim, h, pStart, pEnd, pMax, p; 650b253942bSMatthew G. Knepley 651b253942bSMatthew G. Knepley PetscFunctionBegin; 652*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMGetDimension(dm, &dim)); 653*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMCreateLabel(dm, "cohesive")); 654*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMGetLabel(dm, "cohesive", &label)); 655b253942bSMatthew G. Knepley for (h = 0; h <= dim; ++h) { 656*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexGetSimplexOrBoxCells(dm, h, NULL, &pMax)); 657*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexGetHeightStratum(dm, h, &pStart, &pEnd)); 658*5f80ce2aSJacob Faibussowitsch for (p = pMax; p < pEnd; ++p) CHKERRQ(DMLabelSetValue(label, p, 1)); 659b253942bSMatthew G. Knepley } 660b253942bSMatthew G. Knepley PetscFunctionReturn(0); 661b253942bSMatthew G. Knepley } 662b253942bSMatthew G. Knepley 663b253942bSMatthew G. Knepley static PetscErrorCode CreateDiscretization(DM dm, AppCtx *user) 664b253942bSMatthew G. Knepley { 665b253942bSMatthew G. Knepley PetscFE fe; 666b253942bSMatthew G. Knepley DMLabel fault; 667b7519becSMatthew G. Knepley PetscInt dim, Ncf = user->cohesiveFields, f; 668b253942bSMatthew G. Knepley 669b253942bSMatthew G. Knepley PetscFunctionBegin; 670*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMGetDimension(dm, &dim)); 671*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMGetLabel(dm, "cohesive", &fault)); 672*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMLabelView(fault, PETSC_VIEWER_STDOUT_WORLD)); 673b253942bSMatthew G. Knepley 674*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscFECreateDefault(PETSC_COMM_SELF, dim, dim, user->cellSimplex, "displacement_", PETSC_DETERMINE, &fe)); 675*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscFESetName(fe, "displacement")); 676*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMAddField(dm, NULL, (PetscObject) fe)); 677*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscFEDestroy(&fe)); 678b253942bSMatthew G. Knepley 679b7519becSMatthew G. Knepley if (Ncf > 0) { 680*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscFECreateDefault(PETSC_COMM_SELF, dim-1, dim, user->cellSimplex, "faulttraction_", PETSC_DETERMINE, &fe)); 681*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscFESetName(fe, "fault traction")); 682*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMAddField(dm, fault, (PetscObject) fe)); 683*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscFEDestroy(&fe)); 684b7519becSMatthew G. Knepley } 685b7519becSMatthew G. Knepley for (f = 1; f < Ncf; ++f) { 686b7519becSMatthew G. Knepley char name[256], opt[256]; 687b7519becSMatthew G. Knepley 688*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscSNPrintf(name, 256, "fault field %D", f)); 689*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscSNPrintf(opt, 256, "faultfield_%D_", f)); 690*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscFECreateDefault(PETSC_COMM_SELF, dim-1, dim, user->cellSimplex, opt, PETSC_DETERMINE, &fe)); 691*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscFESetName(fe, name)); 692*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMAddField(dm, fault, (PetscObject) fe)); 693*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscFEDestroy(&fe)); 694b7519becSMatthew G. Knepley } 695b253942bSMatthew G. Knepley 696*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMCreateDS(dm)); 697b253942bSMatthew G. Knepley PetscFunctionReturn(0); 698b253942bSMatthew G. Knepley } 699b253942bSMatthew G. Knepley 700b253942bSMatthew G. Knepley static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 701c4762a1bSJed Brown { 702c4762a1bSJed Brown PetscInt dim = user->dim; 703c4762a1bSJed Brown PetscBool cellSimplex = user->cellSimplex, hasFault, hasFault2, hasParallelFault; 704c4762a1bSJed Brown PetscMPIInt rank, size; 705dd96b1f9SToby Isaac DMLabel matLabel; 706c4762a1bSJed Brown 707c4762a1bSJed Brown PetscFunctionBegin; 708*5f80ce2aSJacob Faibussowitsch CHKERRMPI(MPI_Comm_rank(comm, &rank)); 709*5f80ce2aSJacob Faibussowitsch CHKERRMPI(MPI_Comm_size(comm, &size)); 710*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMCreate(comm, dm)); 711*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetType(*dm, DMPLEX)); 712*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetDimension(*dm, dim)); 713c4762a1bSJed Brown switch (dim) { 714c4762a1bSJed Brown case 2: 715c4762a1bSJed Brown if (cellSimplex) { 716*5f80ce2aSJacob Faibussowitsch CHKERRQ(CreateSimplex_2D(comm, user->testNum, dm)); 717c4762a1bSJed Brown } else { 718*5f80ce2aSJacob Faibussowitsch CHKERRQ(CreateQuad_2D(comm, user->testNum, dm)); 719c4762a1bSJed Brown } 720c4762a1bSJed Brown break; 721c4762a1bSJed Brown case 3: 722c4762a1bSJed Brown if (cellSimplex) { 723*5f80ce2aSJacob Faibussowitsch CHKERRQ(CreateSimplex_3D(comm, user, *dm)); 724c4762a1bSJed Brown } else { 725*5f80ce2aSJacob Faibussowitsch CHKERRQ(CreateHex_3D(comm, user->testNum, dm)); 726c4762a1bSJed Brown } 727c4762a1bSJed Brown break; 728c4762a1bSJed Brown default: 72998921bdaSJacob Faibussowitsch SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make hybrid meshes for dimension %d", dim); 730c4762a1bSJed Brown } 731*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscObjectSetOptionsPrefix((PetscObject) *dm, "orig_")); 732*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexDistributeSetDefault(*dm, PETSC_FALSE)); 733*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetFromOptions(*dm)); 734*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMGetLabel(*dm, "material", &matLabel)); 735dd96b1f9SToby Isaac if (matLabel) { 736*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexLabelComplete(*dm, matLabel)); 737dd96b1f9SToby Isaac } 738*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMViewFromOptions(*dm, NULL, "-dm_view")); 739*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMHasLabel(*dm, "fault", &hasFault)); 740c4762a1bSJed Brown if (hasFault) { 741b253942bSMatthew G. Knepley DM dmHybrid = NULL, dmInterface = NULL; 742b253942bSMatthew G. Knepley DMLabel faultLabel, faultBdLabel, hybridLabel, splitLabel; 743c4762a1bSJed Brown 744*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMGetLabel(*dm, "fault", &faultLabel)); 745*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMGetLabel(*dm, "faultBd", &faultBdLabel)); 746*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, &hybridLabel, &splitLabel, &dmInterface, &dmHybrid)); 747*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMLabelView(hybridLabel, PETSC_VIEWER_STDOUT_WORLD)); 748*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMLabelDestroy(&hybridLabel)); 749*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMLabelView(splitLabel, PETSC_VIEWER_STDOUT_WORLD)); 750*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMLabelDestroy(&splitLabel)); 751*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMViewFromOptions(dmInterface, NULL, "-dm_interface_view")); 752*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMDestroy(&dmInterface)); 753*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMDestroy(dm)); 754c4762a1bSJed Brown *dm = dmHybrid; 755c4762a1bSJed Brown } 756*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMHasLabel(*dm, "fault2", &hasFault2)); 757c4762a1bSJed Brown if (hasFault2) { 758c4762a1bSJed Brown DM dmHybrid = NULL; 759c4762a1bSJed Brown DMLabel faultLabel, faultBdLabel, hybridLabel; 760c4762a1bSJed Brown 761*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscObjectSetOptionsPrefix((PetscObject) *dm, "faulted_")); 762*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMViewFromOptions(*dm, NULL, "-dm_view_pre")); 763*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexDistributeSetDefault(*dm, PETSC_FALSE)); 764*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetFromOptions(*dm)); 765*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMViewFromOptions(*dm, NULL, "-dm_view")); 766*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMGetLabel(*dm, "fault2", &faultLabel)); 767*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMGetLabel(*dm, "fault2Bd", &faultBdLabel)); 768*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, &hybridLabel, NULL, NULL, &dmHybrid)); 769*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMLabelView(hybridLabel, PETSC_VIEWER_STDOUT_WORLD)); 770*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMLabelDestroy(&hybridLabel)); 771*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMDestroy(dm)); 772c4762a1bSJed Brown *dm = dmHybrid; 773c4762a1bSJed Brown } 774c4762a1bSJed Brown if (user->testPartition && size > 1) { 775c4762a1bSJed Brown PetscPartitioner part; 776c4762a1bSJed Brown PetscInt *sizes = NULL; 777c4762a1bSJed Brown PetscInt *points = NULL; 778c4762a1bSJed Brown 779dd400576SPatrick Sanan if (rank == 0) { 780c4762a1bSJed Brown if (dim == 2 && cellSimplex && size == 2) { 781c4762a1bSJed Brown switch (user->testNum) { 782c4762a1bSJed Brown case 0: { 783c4762a1bSJed Brown PetscInt triSizes_p2[2] = {1, 2}; 784c4762a1bSJed Brown PetscInt triPoints_p2[3] = {0, 1, 2}; 785c4762a1bSJed Brown 786*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscMalloc2(2, &sizes, 3, &points)); 787*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscArraycpy(sizes, triSizes_p2, 2)); 788*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscArraycpy(points, triPoints_p2, 3));break;} 789c4762a1bSJed Brown default: 79098921bdaSJacob Faibussowitsch SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %d for triangular mesh on 2 procs", user->testNum); 791c4762a1bSJed Brown } 792c4762a1bSJed Brown } else if (dim == 2 && !cellSimplex && size == 2) { 793c4762a1bSJed Brown switch (user->testNum) { 794c4762a1bSJed Brown case 0: { 795c4762a1bSJed Brown PetscInt quadSizes_p2[2] = {1, 2}; 796c4762a1bSJed Brown PetscInt quadPoints_p2[3] = {0, 1, 2}; 797c4762a1bSJed Brown 798*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscMalloc2(2, &sizes, 3, &points)); 799*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscArraycpy(sizes, quadSizes_p2, 2)); 800*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscArraycpy(points, quadPoints_p2, 3));break;} 801c4762a1bSJed Brown case 2: { 802c4762a1bSJed Brown PetscInt quadSizes_p2[2] = {1, 1}; 803c4762a1bSJed Brown PetscInt quadPoints_p2[2] = {0, 1}; 804c4762a1bSJed Brown 805*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscMalloc2(2, &sizes, 2, &points)); 806*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscArraycpy(sizes, quadSizes_p2, 2)); 807*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscArraycpy(points, quadPoints_p2, 2));break;} 808c4762a1bSJed Brown default: 80998921bdaSJacob Faibussowitsch SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %d for quadrilateral mesh on 2 procs", user->testNum); 810c4762a1bSJed Brown } 811c4762a1bSJed Brown } else if (dim == 3 && cellSimplex && size == 2) { 812c4762a1bSJed Brown switch (user->testNum) { 813c4762a1bSJed Brown case 0: { 814c4762a1bSJed Brown PetscInt tetSizes_p2[2] = {1, 2}; 815c4762a1bSJed Brown PetscInt tetPoints_p2[3] = {0, 1, 2}; 816c4762a1bSJed Brown 817*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscMalloc2(2, &sizes, 3, &points)); 818*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscArraycpy(sizes, tetSizes_p2, 2)); 819*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscArraycpy(points, tetPoints_p2, 3));break;} 820c4762a1bSJed Brown default: 82198921bdaSJacob Faibussowitsch SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %d for triangular mesh on 2 procs", user->testNum); 822c4762a1bSJed Brown } 823c4762a1bSJed Brown } else if (dim == 3 && !cellSimplex && size == 2) { 824c4762a1bSJed Brown switch (user->testNum) { 825c4762a1bSJed Brown case 0: { 826c4762a1bSJed Brown PetscInt hexSizes_p2[2] = {1, 2}; 827c4762a1bSJed Brown PetscInt hexPoints_p2[3] = {0, 1, 2}; 828c4762a1bSJed Brown 829*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscMalloc2(2, &sizes, 3, &points)); 830*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscArraycpy(sizes, hexSizes_p2, 2)); 831*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscArraycpy(points, hexPoints_p2, 3));break;} 832c4762a1bSJed Brown default: 83398921bdaSJacob Faibussowitsch SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %d for triangular mesh on 2 procs", user->testNum); 834c4762a1bSJed Brown } 835c4762a1bSJed Brown } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test partition"); 836c4762a1bSJed Brown } 837*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexGetPartitioner(*dm, &part)); 838*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscPartitionerSetType(part, PETSCPARTITIONERSHELL)); 839*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscPartitionerShellSetPartition(part, size, sizes, points)); 840*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscFree2(sizes, points)); 841c4762a1bSJed Brown } 842c4762a1bSJed Brown { 843c4762a1bSJed Brown DM pdm = NULL; 844c4762a1bSJed Brown 845c4762a1bSJed Brown /* Distribute mesh over processes */ 846*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexDistribute(*dm, 0, NULL, &pdm)); 847c4762a1bSJed Brown if (pdm) { 848*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMViewFromOptions(pdm, NULL, "-dm_view")); 849*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMDestroy(dm)); 850c4762a1bSJed Brown *dm = pdm; 851c4762a1bSJed Brown } 852c4762a1bSJed Brown } 853*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMHasLabel(*dm, "pfault", &hasParallelFault)); 854c4762a1bSJed Brown if (hasParallelFault) { 855c4762a1bSJed Brown DM dmHybrid = NULL; 856c4762a1bSJed Brown DMLabel faultLabel, faultBdLabel, hybridLabel; 857c4762a1bSJed Brown 858*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMGetLabel(*dm, "pfault", &faultLabel)); 859*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMGetLabel(*dm, "pfaultBd", &faultBdLabel)); 860*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, &hybridLabel, NULL, NULL, &dmHybrid)); 861*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMLabelView(hybridLabel, PETSC_VIEWER_STDOUT_WORLD)); 862*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMLabelDestroy(&hybridLabel)); 863*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMDestroy(dm)); 864c4762a1bSJed Brown *dm = dmHybrid; 865c4762a1bSJed Brown } 866*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscObjectSetName((PetscObject) *dm, "Hybrid Mesh")); 867*5f80ce2aSJacob Faibussowitsch CHKERRQ(CreateFaultLabel(*dm)); 868*5f80ce2aSJacob Faibussowitsch CHKERRQ(CreateDiscretization(*dm, user)); 869*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMViewFromOptions(*dm, NULL, "-dm_view_pre")); 870*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexDistributeSetDefault(*dm, PETSC_FALSE)); 871*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetFromOptions(*dm)); 872*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMViewFromOptions(*dm, NULL, "-dm_view")); 873c4762a1bSJed Brown PetscFunctionReturn(0); 874c4762a1bSJed Brown } 875c4762a1bSJed Brown 876b253942bSMatthew G. Knepley static PetscErrorCode TestMesh(DM dm, AppCtx *user) 877b253942bSMatthew G. Knepley { 878b253942bSMatthew G. Knepley PetscFunctionBegin; 879*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexCheckSymmetry(dm)); 880*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexCheckSkeleton(dm, 0)); 881*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexCheckFaces(dm, 0)); 882b253942bSMatthew G. Knepley PetscFunctionReturn(0); 883b253942bSMatthew G. Knepley } 884b253942bSMatthew G. Knepley 885b253942bSMatthew G. Knepley static PetscErrorCode TestDiscretization(DM dm, AppCtx *user) 886b253942bSMatthew G. Knepley { 887b253942bSMatthew G. Knepley PetscSection s; 888b253942bSMatthew G. Knepley 889b253942bSMatthew G. Knepley PetscFunctionBegin; 890*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMGetSection(dm, &s)); 891*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscObjectViewFromOptions((PetscObject) s, NULL, "-local_section_view")); 892b253942bSMatthew G. Knepley PetscFunctionReturn(0); 893b253942bSMatthew G. Knepley } 894b253942bSMatthew G. Knepley 895b253942bSMatthew G. Knepley static PetscErrorCode r(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 896b253942bSMatthew G. Knepley { 897b253942bSMatthew G. Knepley PetscInt d; 898b253942bSMatthew G. Knepley for (d = 0; d < dim; ++d) u[d] = x[d]; 899b253942bSMatthew G. Knepley return 0; 900b253942bSMatthew G. Knepley } 901b253942bSMatthew G. Knepley 902b253942bSMatthew G. Knepley static PetscErrorCode rp1(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 903b253942bSMatthew G. Knepley { 904b253942bSMatthew G. Knepley PetscInt d; 905b253942bSMatthew G. Knepley for (d = 0; d < dim; ++d) u[d] = x[d] + (d > 0 ? 1.0 : 0.0); 906b253942bSMatthew G. Knepley return 0; 907b253942bSMatthew G. Knepley } 908b253942bSMatthew G. Knepley 909b253942bSMatthew G. Knepley static PetscErrorCode phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 910b253942bSMatthew G. Knepley { 911b253942bSMatthew G. Knepley PetscInt d; 912b253942bSMatthew G. Knepley u[0] = -x[1]; 913b253942bSMatthew G. Knepley u[1] = x[0]; 914b253942bSMatthew G. Knepley for (d = 2; d < dim; ++d) u[d] = x[d]; 915b253942bSMatthew G. Knepley return 0; 916b253942bSMatthew G. Knepley } 917b253942bSMatthew G. Knepley 918b253942bSMatthew G. Knepley /* \lambda \cdot (\psi_u^- - \psi_u^+) */ 919b253942bSMatthew G. Knepley static void f0_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 920b253942bSMatthew G. Knepley const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 921b253942bSMatthew G. Knepley const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 922b253942bSMatthew G. Knepley PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 923b253942bSMatthew G. Knepley { 924b253942bSMatthew G. Knepley const PetscInt Nc = uOff[1]-uOff[0]; 925b253942bSMatthew G. Knepley PetscInt c; 926b253942bSMatthew G. Knepley for (c = 0; c < Nc; ++c) { 927b253942bSMatthew G. Knepley f0[c] = u[Nc*2+c] + x[Nc-c-1]; 928b253942bSMatthew G. Knepley f0[Nc+c] = -(u[Nc*2+c] + x[Nc-c-1]); 929b253942bSMatthew G. Knepley } 930b253942bSMatthew G. Knepley } 931b253942bSMatthew G. Knepley 932b253942bSMatthew G. Knepley /* (d - u^+ + u^-) \cdot \psi_\lambda */ 933b253942bSMatthew G. Knepley static void f0_bd_l(PetscInt dim, PetscInt Nf, PetscInt NfAux, 934b253942bSMatthew G. Knepley const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 935b253942bSMatthew G. Knepley const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 936b253942bSMatthew G. Knepley PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 937b253942bSMatthew G. Knepley { 938b253942bSMatthew G. Knepley const PetscInt Nc = uOff[2]-uOff[1]; 939b253942bSMatthew G. Knepley PetscInt c; 940b253942bSMatthew G. Knepley 941b253942bSMatthew G. Knepley for (c = 0; c < Nc; ++c) f0[c] = (c > 0 ? 1.0 : 0.0) + u[c] - u[Nc+c]; 942b253942bSMatthew G. Knepley } 943b253942bSMatthew G. Knepley 944b253942bSMatthew G. Knepley /* \psi_lambda \cdot (\psi_u^- - \psi_u^+) */ 945b253942bSMatthew G. Knepley static void g0_bd_ul(PetscInt dim, PetscInt Nf, PetscInt NfAux, 946b253942bSMatthew G. Knepley const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 947b253942bSMatthew G. Knepley const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 948b253942bSMatthew G. Knepley PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 949b253942bSMatthew G. Knepley { 950b253942bSMatthew G. Knepley const PetscInt Nc = uOff[1]-uOff[0]; 951b253942bSMatthew G. Knepley PetscInt c; 952b253942bSMatthew G. Knepley 953b253942bSMatthew G. Knepley for (c = 0; c < Nc; ++c) { 954b253942bSMatthew G. Knepley g0[(0 +c)*Nc+c] = 1.0; 955b253942bSMatthew G. Knepley g0[(Nc+c)*Nc+c] = -1.0; 956b253942bSMatthew G. Knepley } 957b253942bSMatthew G. Knepley } 958b253942bSMatthew G. Knepley 959b253942bSMatthew G. Knepley /* (-\psi_u^+ + \psi_u^-) \cdot \psi_\lambda */ 960b253942bSMatthew G. Knepley static void g0_bd_lu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 961b253942bSMatthew G. Knepley const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 962b253942bSMatthew G. Knepley const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 963b253942bSMatthew G. Knepley PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 964b253942bSMatthew G. Knepley { 965b253942bSMatthew G. Knepley const PetscInt Nc = uOff[2]-uOff[1]; 966b253942bSMatthew G. Knepley PetscInt c; 967b253942bSMatthew G. Knepley 968b253942bSMatthew G. Knepley for (c = 0; c < Nc; ++c) { 969b253942bSMatthew G. Knepley g0[c*Nc*2+c] = 1.0; 970b253942bSMatthew G. Knepley g0[c*Nc*2+Nc+c] = -1.0; 971b253942bSMatthew G. Knepley } 972b253942bSMatthew G. Knepley } 973b253942bSMatthew G. Knepley 974b253942bSMatthew G. Knepley static PetscErrorCode TestAssembly(DM dm, AppCtx *user) 975b253942bSMatthew G. Knepley { 976b253942bSMatthew G. Knepley Mat J; 977b253942bSMatthew G. Knepley Vec locX, locF; 978b253942bSMatthew G. Knepley PetscDS probh; 979b253942bSMatthew G. Knepley DMLabel fault, material; 980b253942bSMatthew G. Knepley IS cohesiveCells; 981b7519becSMatthew G. Knepley PetscWeakForm wf; 98206ad1575SMatthew G. Knepley PetscFormKey keys[3]; 983b253942bSMatthew G. Knepley PetscErrorCode (*initialGuess[2])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx); 984b253942bSMatthew G. Knepley PetscInt dim, Nf, cMax, cEnd, id; 985b7519becSMatthew G. Knepley PetscMPIInt rank; 986b253942bSMatthew G. Knepley 987b253942bSMatthew G. Knepley PetscFunctionBegin; 988*5f80ce2aSJacob Faibussowitsch CHKERRMPI(MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank)); 989*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMGetDimension(dm, &dim)); 990*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexGetSimplexOrBoxCells(dm, 0, NULL, &cMax)); 991*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexGetHeightStratum(dm, 0, NULL, &cEnd)); 992*5f80ce2aSJacob Faibussowitsch CHKERRQ(ISCreateStride(PETSC_COMM_SELF, cEnd - cMax, cMax, 1, &cohesiveCells)); 993*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMGetLabel(dm, "cohesive", &fault)); 994*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMGetLocalVector(dm, &locX)); 995*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscObjectSetName((PetscObject) locX, "Local Solution")); 996*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMGetLocalVector(dm, &locF)); 997*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscObjectSetName((PetscObject) locF, "Local Residual")); 998*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMCreateMatrix(dm, &J)); 999*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscObjectSetName((PetscObject) J, "Jacobian")); 1000b253942bSMatthew G. Knepley 1001b253942bSMatthew G. Knepley /* The initial guess has displacement shifted by one unit in each fault parallel direction across the fault */ 1002*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMGetLabel(dm, "material", &material)); 1003b253942bSMatthew G. Knepley id = 1; 1004b253942bSMatthew G. Knepley initialGuess[0] = r; 1005b253942bSMatthew G. Knepley initialGuess[1] = NULL; 1006*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMProjectFunctionLabelLocal(dm, 0.0, material, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX)); 1007b253942bSMatthew G. Knepley id = 2; 1008b253942bSMatthew G. Knepley initialGuess[0] = rp1; 1009b253942bSMatthew G. Knepley initialGuess[1] = NULL; 1010*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMProjectFunctionLabelLocal(dm, 0.0, material, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX)); 1011b253942bSMatthew G. Knepley id = 1; 1012b253942bSMatthew G. Knepley initialGuess[0] = NULL; 1013b253942bSMatthew G. Knepley initialGuess[1] = phi; 1014*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMProjectFunctionLabelLocal(dm, 0.0, fault, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX)); 1015*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecViewFromOptions(locX, NULL, "-local_solution_view")); 1016b253942bSMatthew G. Knepley 1017*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMGetCellDS(dm, cMax, &probh)); 1018*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscDSGetWeakForm(probh, &wf)); 1019*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscDSGetNumFields(probh, &Nf)); 1020*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscWeakFormSetIndexBdResidual(wf, material, 1, 0, 0, 0, f0_bd_u, 0, NULL)); 1021*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscWeakFormSetIndexBdResidual(wf, material, 2, 0, 0, 0, f0_bd_u, 0, NULL)); 1022*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscWeakFormSetIndexBdJacobian(wf, material, 1, 0, 1, 0, 0, g0_bd_ul, 0, NULL, 0, NULL, 0, NULL)); 1023*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscWeakFormSetIndexBdJacobian(wf, material, 2, 0, 1, 0, 0, g0_bd_ul, 0, NULL, 0, NULL, 0, NULL)); 1024b7519becSMatthew G. Knepley if (Nf > 1) { 1025*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscWeakFormSetIndexBdResidual(wf, fault, 1, 1, 0, 0, f0_bd_l, 0, NULL)); 1026*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscWeakFormSetIndexBdJacobian(wf, fault, 1, 1, 0, 0, 0, g0_bd_lu, 0, NULL, 0, NULL, 0, NULL)); 1027b7519becSMatthew G. Knepley } 1028*5f80ce2aSJacob Faibussowitsch if (!rank) CHKERRQ(PetscDSView(probh, NULL)); 1029b253942bSMatthew G. Knepley 1030dd96b1f9SToby Isaac keys[0].label = NULL; 1031dd96b1f9SToby Isaac keys[0].value = 0; 10326528b96dSMatthew G. Knepley keys[0].field = 0; 1033dd96b1f9SToby Isaac keys[0].part = 0; 1034b7519becSMatthew G. Knepley keys[1].label = material; 1035b7519becSMatthew G. Knepley keys[1].value = 2; 10366528b96dSMatthew G. Knepley keys[1].field = 0; 1037dd96b1f9SToby Isaac keys[1].part = 0; 1038b7519becSMatthew G. Knepley keys[2].label = fault; 1039b7519becSMatthew G. Knepley keys[2].value = 1; 1040b7519becSMatthew G. Knepley keys[2].field = 1; 1041dd96b1f9SToby Isaac keys[2].part = 0; 1042*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecSet(locF, 0.)); 1043*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexComputeResidual_Hybrid_Internal(dm, keys, cohesiveCells, 0.0, locX, NULL, 0.0, locF, user)); 1044*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecViewFromOptions(locF, NULL, "-local_residual_view")); 1045*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatZeroEntries(J)); 1046*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexComputeJacobian_Hybrid_Internal(dm, keys, cohesiveCells, 0.0, 0.0, locX, NULL, J, J, user)); 1047*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatViewFromOptions(J, NULL, "-local_jacobian_view")); 1048b253942bSMatthew G. Knepley 1049*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMRestoreLocalVector(dm, &locX)); 1050*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMRestoreLocalVector(dm, &locF)); 1051*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatDestroy(&J)); 1052*5f80ce2aSJacob Faibussowitsch CHKERRQ(ISDestroy(&cohesiveCells)); 1053b253942bSMatthew G. Knepley PetscFunctionReturn(0); 1054b253942bSMatthew G. Knepley } 1055b253942bSMatthew G. Knepley 1056c4762a1bSJed Brown int main(int argc, char **argv) 1057c4762a1bSJed Brown { 1058c4762a1bSJed Brown DM dm; 1059c4762a1bSJed Brown AppCtx user; /* user-defined work context */ 1060c4762a1bSJed Brown PetscErrorCode ierr; 1061c4762a1bSJed Brown 1062c4762a1bSJed Brown ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr; 1063*5f80ce2aSJacob Faibussowitsch CHKERRQ(ProcessOptions(PETSC_COMM_WORLD, &user)); 1064*5f80ce2aSJacob Faibussowitsch CHKERRQ(CreateMesh(PETSC_COMM_WORLD, &user, &dm)); 1065*5f80ce2aSJacob Faibussowitsch CHKERRQ(TestMesh(dm, &user)); 1066*5f80ce2aSJacob Faibussowitsch CHKERRQ(TestDiscretization(dm, &user)); 1067*5f80ce2aSJacob Faibussowitsch CHKERRQ(TestAssembly(dm, &user)); 1068*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMDestroy(&dm)); 1069c4762a1bSJed Brown ierr = PetscFinalize(); 1070c4762a1bSJed Brown return ierr; 1071c4762a1bSJed Brown } 1072c4762a1bSJed Brown 1073c4762a1bSJed Brown /*TEST 1074c4762a1bSJed Brown testset: 1075b253942bSMatthew G. Knepley args: -orig_dm_plex_check_all -dm_plex_check_all \ 1076b7519becSMatthew G. Knepley -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1 -local_section_view \ 1077b253942bSMatthew G. Knepley -local_solution_view -local_residual_view -local_jacobian_view 1078c4762a1bSJed Brown test: 1079c4762a1bSJed Brown suffix: tri_0 1080c4762a1bSJed Brown args: -dim 2 1081c4762a1bSJed Brown test: 1082c4762a1bSJed Brown suffix: tri_t1_0 1083c4762a1bSJed Brown args: -dim 2 -test_num 1 1084c4762a1bSJed Brown test: 1085c4762a1bSJed Brown suffix: tet_0 1086c4762a1bSJed Brown args: -dim 3 1087c4762a1bSJed Brown test: 1088b253942bSMatthew G. Knepley suffix: tet_t1_0 1089b253942bSMatthew G. Knepley args: -dim 3 -test_num 1 1090b253942bSMatthew G. Knepley 1091b253942bSMatthew G. Knepley testset: 1092b253942bSMatthew G. Knepley args: -orig_dm_plex_check_all -dm_plex_check_all \ 1093b7519becSMatthew G. Knepley -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1 1094b253942bSMatthew G. Knepley test: 1095c4762a1bSJed Brown suffix: tet_1 1096c4762a1bSJed Brown nsize: 2 1097c4762a1bSJed Brown args: -dim 3 1098c4762a1bSJed Brown test: 1099b253942bSMatthew G. Knepley suffix: tri_1 1100b253942bSMatthew G. Knepley nsize: 2 1101b253942bSMatthew G. Knepley args: -dim 2 1102c4762a1bSJed Brown 1103c4762a1bSJed Brown testset: 1104ecfb78b5SMatthew G. Knepley args: -orig_dm_plex_check_all -dm_plex_check_all \ 1105b7519becSMatthew G. Knepley -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1 1106c4762a1bSJed Brown # 2D Quads 1107c4762a1bSJed Brown test: 1108c4762a1bSJed Brown suffix: quad_0 1109c4762a1bSJed Brown args: -dim 2 -cell_simplex 0 1110c4762a1bSJed Brown test: 1111c4762a1bSJed Brown suffix: quad_1 1112c4762a1bSJed Brown nsize: 2 1113c4762a1bSJed Brown args: -dim 2 -cell_simplex 0 1114c4762a1bSJed Brown test: 1115c4762a1bSJed Brown suffix: quad_t1_0 111654fcfd0cSMatthew G. Knepley args: -dim 2 -cell_simplex 0 -test_num 1 -faulted_dm_plex_check_all 1117c4762a1bSJed Brown # 3D Hex 1118c4762a1bSJed Brown test: 1119c4762a1bSJed Brown suffix: hex_0 1120c4762a1bSJed Brown args: -dim 3 -cell_simplex 0 1121c4762a1bSJed Brown test: 1122c4762a1bSJed Brown suffix: hex_1 1123c4762a1bSJed Brown nsize: 2 1124c4762a1bSJed Brown args: -dim 3 -cell_simplex 0 1125c4762a1bSJed Brown test: 1126c4762a1bSJed Brown suffix: hex_t1_0 1127c4762a1bSJed Brown args: -dim 3 -cell_simplex 0 -test_num 1 1128c4762a1bSJed Brown test: 1129c4762a1bSJed Brown suffix: hex_t2_0 1130c4762a1bSJed Brown args: -dim 3 -cell_simplex 0 -test_num 2 1131c4762a1bSJed Brown 1132c4762a1bSJed Brown TEST*/ 1133