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 PetscFunctionBegin; 315c4762a1bSJed Brown options->debug = 0; 316c4762a1bSJed Brown options->dim = 2; 317c4762a1bSJed Brown options->cellSimplex = PETSC_TRUE; 318c4762a1bSJed Brown options->testPartition = PETSC_TRUE; 319c4762a1bSJed Brown options->testNum = 0; 320b7519becSMatthew G. Knepley options->cohesiveFields = 1; 321c4762a1bSJed Brown 322*d0609cedSBarry Smith PetscOptionsBegin(comm, "", "Meshing Problem Options", "DMPLEX"); 3239566063dSJacob Faibussowitsch PetscCall(PetscOptionsBoundedInt("-debug", "The debugging level", "ex5.c", options->debug, &options->debug, NULL,0)); 3249566063dSJacob Faibussowitsch PetscCall(PetscOptionsRangeInt("-dim", "The topological mesh dimension", "ex5.c", options->dim, &options->dim, NULL,1,3)); 3259566063dSJacob Faibussowitsch PetscCall(PetscOptionsBool("-cell_simplex", "Use simplices if true, otherwise hexes", "ex5.c", options->cellSimplex, &options->cellSimplex, NULL)); 3269566063dSJacob Faibussowitsch PetscCall(PetscOptionsBool("-test_partition", "Use a fixed partition for testing", "ex5.c", options->testPartition, &options->testPartition, NULL)); 3279566063dSJacob Faibussowitsch PetscCall(PetscOptionsBoundedInt("-test_num", "The particular mesh to test", "ex5.c", options->testNum, &options->testNum, NULL,0)); 3289566063dSJacob Faibussowitsch PetscCall(PetscOptionsBoundedInt("-cohesive_fields", "The number of cohesive fields", "ex5.c", options->cohesiveFields, &options->cohesiveFields, NULL, 0)); 329*d0609cedSBarry Smith PetscOptionsEnd(); 330c4762a1bSJed Brown PetscFunctionReturn(0); 331c4762a1bSJed Brown } 332c4762a1bSJed Brown 333b253942bSMatthew G. Knepley static PetscErrorCode CreateSimplex_2D(MPI_Comm comm, PetscInt testNum, DM *dm) 334c4762a1bSJed Brown { 335c4762a1bSJed Brown DM idm; 336c4762a1bSJed Brown PetscInt p; 337c4762a1bSJed Brown PetscMPIInt rank; 338c4762a1bSJed Brown 339c4762a1bSJed Brown PetscFunctionBegin; 3409566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_rank(comm, &rank)); 341dd400576SPatrick Sanan if (rank == 0) { 342c4762a1bSJed Brown switch (testNum) { 343c4762a1bSJed Brown case 0: 344c4762a1bSJed Brown { 345c4762a1bSJed Brown PetscInt numPoints[2] = {4, 2}; 346c4762a1bSJed Brown PetscInt coneSize[6] = {3, 3, 0, 0, 0, 0}; 347c4762a1bSJed Brown PetscInt cones[6] = {2, 3, 4, 5, 4, 3}; 348c4762a1bSJed Brown PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 349c4762a1bSJed Brown PetscScalar vertexCoords[8] = {-0.5, 0.5, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5}; 350c4762a1bSJed Brown PetscInt markerPoints[8] = {2, 1, 3, 1, 4, 1, 5, 1}; 351c4762a1bSJed Brown PetscInt faultPoints[2] = {3, 4}; 352c4762a1bSJed Brown 3539566063dSJacob Faibussowitsch PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3549566063dSJacob Faibussowitsch for (p = 0; p < 4; ++p) PetscCall(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1])); 3559566063dSJacob Faibussowitsch for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "fault", faultPoints[p], 1)); 3569566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 0, 1)); 3579566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 1, 2)); 358c4762a1bSJed Brown } 359c4762a1bSJed Brown break; 360c4762a1bSJed Brown case 1: 361c4762a1bSJed Brown { 362c4762a1bSJed Brown PetscInt numPoints[2] = {6, 4}; 363c4762a1bSJed Brown PetscInt coneSize[10] = {3, 3, 3, 3, 0, 0, 0, 0, 0, 0}; 364c4762a1bSJed Brown PetscInt cones[12] = {4, 6, 5, 5, 6, 7, 8, 5, 9, 8, 4, 5}; 365c4762a1bSJed Brown PetscInt coneOrientations[12] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 366c4762a1bSJed 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}; 367c4762a1bSJed Brown PetscInt markerPoints[6] = {4, 1, 6, 1, 8, 1}; 368c4762a1bSJed Brown PetscInt faultPoints[3] = {5, 6, 8}; 369c4762a1bSJed Brown 3709566063dSJacob Faibussowitsch PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3719566063dSJacob Faibussowitsch for (p = 0; p < 3; ++p) PetscCall(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1])); 3729566063dSJacob Faibussowitsch for (p = 0; p < 3; ++p) PetscCall(DMSetLabelValue(*dm, "fault", faultPoints[p], 1)); 3739566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 0, 1));PetscCall(DMSetLabelValue(*dm, "material", 3, 1)); 3749566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 1, 2));PetscCall(DMSetLabelValue(*dm, "material", 2, 2)); 375c4762a1bSJed Brown } 376c4762a1bSJed Brown break; 377c4762a1bSJed Brown default: 37898921bdaSJacob Faibussowitsch SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %d", testNum); 379c4762a1bSJed Brown } 380c4762a1bSJed Brown } else { 381c4762a1bSJed Brown PetscInt numPoints[3] = {0, 0, 0}; 382c4762a1bSJed Brown 3839566063dSJacob Faibussowitsch PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, NULL, NULL, NULL, NULL)); 3849566063dSJacob Faibussowitsch PetscCall(DMCreateLabel(*dm, "fault")); 385c4762a1bSJed Brown } 3869566063dSJacob Faibussowitsch PetscCall(DMPlexInterpolate(*dm, &idm)); 3879566063dSJacob Faibussowitsch PetscCall(DMViewFromOptions(idm, NULL, "-in_dm_view")); 3889566063dSJacob Faibussowitsch PetscCall(DMDestroy(dm)); 389c4762a1bSJed Brown *dm = idm; 390c4762a1bSJed Brown PetscFunctionReturn(0); 391c4762a1bSJed Brown } 392c4762a1bSJed Brown 393b253942bSMatthew G. Knepley static PetscErrorCode CreateSimplex_3D(MPI_Comm comm, AppCtx *user, DM dm) 394c4762a1bSJed Brown { 395c4762a1bSJed Brown PetscInt depth = 3, testNum = user->testNum, p; 396c4762a1bSJed Brown PetscMPIInt rank; 397c4762a1bSJed Brown 398c4762a1bSJed Brown PetscFunctionBegin; 3999566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_rank(comm, &rank)); 400dd400576SPatrick Sanan if (rank == 0) { 401c4762a1bSJed Brown switch (testNum) { 402c4762a1bSJed Brown case 0: 403c4762a1bSJed Brown { 404c4762a1bSJed Brown PetscInt numPoints[4] = {5, 7, 9, 2}; 405c4762a1bSJed 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}; 406c4762a1bSJed 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}; 407b5a892a1SMatthew 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}; 408c4762a1bSJed 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}; 409c4762a1bSJed Brown PetscInt markerPoints[20] = {2, 1, 3, 1, 4, 1, 5, 1, 14, 1, 15, 1, 16, 1, 17, 1, 18, 1, 19, 1}; 410c4762a1bSJed Brown PetscInt faultPoints[3] = {3, 4, 5}; 411c4762a1bSJed Brown 4129566063dSJacob Faibussowitsch PetscCall(DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 413c4762a1bSJed Brown for (p = 0; p < 10; ++p) { 4149566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1])); 415c4762a1bSJed Brown } 416c4762a1bSJed Brown for (p = 0; p < 3; ++p) { 4179566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(dm, "fault", faultPoints[p], 1)); 418c4762a1bSJed Brown } 4199566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(dm, "material", 0, 1)); 4209566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(dm, "material", 1, 2)); 421c4762a1bSJed Brown } 422c4762a1bSJed Brown break; 423c4762a1bSJed Brown case 1: 424c4762a1bSJed Brown { 425c4762a1bSJed Brown PetscInt numPoints[4] = {6, 13, 12, 4}; 426c4762a1bSJed 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}; 427c4762a1bSJed 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}; 428b5a892a1SMatthew 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}; 429c4762a1bSJed 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}; 430c4762a1bSJed Brown PetscInt markerPoints[14] = {5, 1, 6, 1, 7, 1, 10, 1, 22, 1, 23, 1, 24, 1}; 431c4762a1bSJed Brown PetscInt faultPoints[4] = {5, 6, 7, 8}; 432c4762a1bSJed Brown 4339566063dSJacob Faibussowitsch PetscCall(DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 434c4762a1bSJed Brown for (p = 0; p < 7; ++p) { 4359566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1])); 436c4762a1bSJed Brown } 437c4762a1bSJed Brown for (p = 0; p < 4; ++p) { 4389566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(dm, "fault", faultPoints[p], 1)); 439c4762a1bSJed Brown } 4409566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(dm, "material", 0, 1)); 4419566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(dm, "material", 1, 2)); 442c4762a1bSJed Brown } 443c4762a1bSJed Brown break; 444c4762a1bSJed Brown default: 44598921bdaSJacob Faibussowitsch SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %d", testNum); 446c4762a1bSJed Brown } 447c4762a1bSJed Brown } else { 448c4762a1bSJed Brown PetscInt numPoints[4] = {0, 0, 0, 0}; 449c4762a1bSJed Brown 4509566063dSJacob Faibussowitsch PetscCall(DMPlexCreateFromDAG(dm, depth, numPoints, NULL, NULL, NULL, NULL)); 4519566063dSJacob Faibussowitsch PetscCall(DMCreateLabel(dm, "fault")); 452c4762a1bSJed Brown } 453c4762a1bSJed Brown PetscFunctionReturn(0); 454c4762a1bSJed Brown } 455c4762a1bSJed Brown 456b253942bSMatthew G. Knepley static PetscErrorCode CreateQuad_2D(MPI_Comm comm, PetscInt testNum, DM *dm) 457c4762a1bSJed Brown { 458c4762a1bSJed Brown DM idm; 459c4762a1bSJed Brown PetscInt p; 460c4762a1bSJed Brown PetscMPIInt rank; 461c4762a1bSJed Brown 462c4762a1bSJed Brown PetscFunctionBegin; 4639566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_rank(comm, &rank)); 464dd400576SPatrick Sanan if (rank == 0) { 465c4762a1bSJed Brown switch (testNum) { 466c4762a1bSJed Brown case 0: 467c4762a1bSJed Brown case 2: 468c4762a1bSJed Brown { 469c4762a1bSJed Brown PetscInt numPoints[2] = {6, 2}; 470c4762a1bSJed Brown PetscInt coneSize[8] = {4, 4, 0, 0, 0, 0, 0, 0}; 471c4762a1bSJed Brown PetscInt cones[8] = {2, 3, 4, 5, 3, 6, 7, 4}; 472c4762a1bSJed Brown PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 473c4762a1bSJed 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}; 474c4762a1bSJed Brown PetscInt markerPoints[12] = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1, 7, 1}; 475c4762a1bSJed Brown PetscInt faultPoints[2] = {3, 4}; 476c4762a1bSJed Brown 4779566063dSJacob Faibussowitsch PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 4789566063dSJacob Faibussowitsch for (p = 0; p < 6; ++p) PetscCall(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1])); 4799566063dSJacob Faibussowitsch if (testNum == 0) for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "fault", faultPoints[p], 1)); 4809566063dSJacob Faibussowitsch if (testNum == 2) for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "pfault", faultPoints[p], 1)); 4819566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 0, 1)); 4829566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 1, 2)); 483c4762a1bSJed Brown } 484c4762a1bSJed Brown break; 485c4762a1bSJed Brown case 1: 486c4762a1bSJed Brown { 487c4762a1bSJed Brown PetscInt numPoints[2] = {16, 9}; 488c4762a1bSJed 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}; 489c4762a1bSJed Brown PetscInt cones[36] = {9, 13, 14, 10, 490c4762a1bSJed Brown 10, 14, 15, 11, 491c4762a1bSJed Brown 11, 15, 16, 12, 492c4762a1bSJed Brown 13, 17, 18, 14, 493c4762a1bSJed Brown 14, 18, 19, 15, 494c4762a1bSJed Brown 15, 19, 20, 16, 495c4762a1bSJed Brown 17, 21, 22, 18, 496c4762a1bSJed Brown 18, 22, 23, 19, 497c4762a1bSJed Brown 19, 23, 24, 20}; 498c4762a1bSJed 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}; 499c4762a1bSJed 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, 500c4762a1bSJed 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}; 501c4762a1bSJed Brown PetscInt faultPoints[3] = {11, 15, 19}; 502c4762a1bSJed Brown PetscInt fault2Points[2] = {17, 18}; 503c4762a1bSJed Brown 5049566063dSJacob Faibussowitsch PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 5059566063dSJacob Faibussowitsch for (p = 0; p < 3; ++p) PetscCall(DMSetLabelValue(*dm, "fault", faultPoints[p], 1)); 5069566063dSJacob Faibussowitsch for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "fault2", fault2Points[p], 1)); 5079566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 0, 1)); 5089566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 1, 1)); 5099566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 2, 1)); 5109566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 3, 1)); 5119566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 4, 1)); 5129566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 5, 2)); 5139566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 6, 2)); 5149566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 7, 2)); 5159566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 8, 2)); 516c4762a1bSJed Brown } 517c4762a1bSJed Brown break; 518c4762a1bSJed Brown default: 51998921bdaSJacob Faibussowitsch SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %d", testNum); 520c4762a1bSJed Brown } 521c4762a1bSJed Brown } else { 522c4762a1bSJed Brown PetscInt numPoints[3] = {0, 0, 0}; 523c4762a1bSJed Brown 5249566063dSJacob Faibussowitsch PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, NULL, NULL, NULL, NULL)); 5259566063dSJacob Faibussowitsch if (testNum == 2) PetscCall(DMCreateLabel(*dm, "pfault")); 5269566063dSJacob Faibussowitsch else PetscCall(DMCreateLabel(*dm, "fault")); 527c4762a1bSJed Brown } 5289566063dSJacob Faibussowitsch PetscCall(DMPlexInterpolate(*dm, &idm)); 5299566063dSJacob Faibussowitsch PetscCall(DMViewFromOptions(idm, NULL, "-in_dm_view")); 5309566063dSJacob Faibussowitsch PetscCall(DMDestroy(dm)); 531c4762a1bSJed Brown *dm = idm; 532c4762a1bSJed Brown PetscFunctionReturn(0); 533c4762a1bSJed Brown } 534c4762a1bSJed Brown 535b253942bSMatthew G. Knepley static PetscErrorCode CreateHex_3D(MPI_Comm comm, PetscInt testNum, DM *dm) 536c4762a1bSJed Brown { 537c4762a1bSJed Brown DM idm; 538c4762a1bSJed Brown PetscInt depth = 3, p; 539c4762a1bSJed Brown PetscMPIInt rank; 540c4762a1bSJed Brown 541c4762a1bSJed Brown PetscFunctionBegin; 5429566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_rank(comm, &rank)); 543dd400576SPatrick Sanan if (rank == 0) { 544c4762a1bSJed Brown switch (testNum) { 545c4762a1bSJed Brown case 0: 546c4762a1bSJed Brown { 547c4762a1bSJed Brown PetscInt numPoints[2] = {12, 2}; 548c4762a1bSJed Brown PetscInt coneSize[14] = {8,8, 0,0,0,0,0,0,0,0,0,0,0,0}; 549c4762a1bSJed Brown PetscInt cones[16] = {2,5,4,3,6,7,8,9, 3,4,11,10,7,12,13,8}; 550c4762a1bSJed Brown PetscInt coneOrientations[16] = {0,0,0,0,0,0,0,0, 0,0, 0,0 ,0, 0, 0,0}; 551c4762a1bSJed 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, 552c4762a1bSJed 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, 553c4762a1bSJed 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}; 554c4762a1bSJed Brown PetscInt markerPoints[52] = {2,1,3,1,4,1,5,1,6,1,7,1,8,1,9,1}; 555c4762a1bSJed Brown PetscInt faultPoints[4] = {3, 4, 7, 8}; 556c4762a1bSJed Brown 5579566063dSJacob Faibussowitsch PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 5589566063dSJacob Faibussowitsch PetscCall(DMPlexInterpolate(*dm, &idm)); 5599566063dSJacob Faibussowitsch for (p = 0; p < 8; ++p) PetscCall(DMSetLabelValue(idm, "marker", markerPoints[p*2], markerPoints[p*2+1])); 5609566063dSJacob Faibussowitsch for (p = 0; p < 4; ++p) PetscCall(DMSetLabelValue(idm, "fault", faultPoints[p], 1)); 5619566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 0, 1)); 5629566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 1, 2)); 563c4762a1bSJed Brown } 564c4762a1bSJed Brown break; 565c4762a1bSJed Brown case 1: 566c4762a1bSJed Brown { 567c4762a1bSJed Brown /* Cell Adjacency Graph: 568c4762a1bSJed Brown 0 -- { 8, 13, 21, 24} --> 1 569c4762a1bSJed Brown 0 -- {20, 21, 23, 24} --> 5 F 570c4762a1bSJed Brown 1 -- {10, 15, 21, 24} --> 2 571c4762a1bSJed Brown 1 -- {13, 14, 15, 24} --> 6 572c4762a1bSJed Brown 2 -- {21, 22, 24, 25} --> 4 F 573c4762a1bSJed Brown 3 -- {21, 24, 30, 35} --> 4 574c4762a1bSJed Brown 3 -- {21, 24, 28, 33} --> 5 575c4762a1bSJed Brown */ 576c4762a1bSJed Brown PetscInt numPoints[2] = {30, 7}; 577c4762a1bSJed 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}; 578c4762a1bSJed Brown PetscInt cones[56] = { 8, 21, 20, 7, 13, 12, 23, 24, 579c4762a1bSJed Brown 14, 15, 10, 9, 13, 8, 21, 24, 580c4762a1bSJed Brown 15, 16, 11, 10, 24, 21, 22, 25, 581c4762a1bSJed Brown 30, 29, 28, 21, 35, 24, 33, 34, 582c4762a1bSJed Brown 24, 21, 30, 35, 25, 36, 31, 22, 583c4762a1bSJed Brown 27, 20, 21, 28, 32, 33, 24, 23, 584c4762a1bSJed Brown 15, 24, 13, 14, 19, 18, 17, 26}; 585c4762a1bSJed 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}; 586c4762a1bSJed 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, 587c4762a1bSJed 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, 588c4762a1bSJed 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, 589c4762a1bSJed 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, 590c4762a1bSJed 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}; 591c4762a1bSJed Brown PetscInt faultPoints[6] = {20, 21, 22, 23, 24, 25}; 592c4762a1bSJed Brown 5939566063dSJacob Faibussowitsch PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 5949566063dSJacob Faibussowitsch PetscCall(DMPlexInterpolate(*dm, &idm)); 5959566063dSJacob Faibussowitsch for (p = 0; p < 6; ++p) PetscCall(DMSetLabelValue(idm, "fault", faultPoints[p], 1)); 5969566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 0, 1)); 5979566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 1, 1)); 5989566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 2, 1)); 5999566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 3, 2)); 6009566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 4, 2)); 6019566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 5, 2)); 6029566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 6, 2)); 603c4762a1bSJed Brown } 604c4762a1bSJed Brown break; 605c4762a1bSJed Brown case 2: 606c4762a1bSJed Brown { 607c4762a1bSJed Brown /* Buried fault edge */ 608c4762a1bSJed Brown PetscInt numPoints[2] = {18, 4}; 609c4762a1bSJed 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}; 610c4762a1bSJed Brown PetscInt cones[32] = { 4, 5, 8, 7, 13, 16, 17, 14, 611c4762a1bSJed Brown 5, 6, 9, 8, 14, 17, 18, 15, 612c4762a1bSJed Brown 7, 8, 11, 10, 16, 19, 20, 17, 613c4762a1bSJed Brown 8, 9, 12, 11, 17, 20, 21, 18}; 614c4762a1bSJed 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}; 615c4762a1bSJed 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, 616c4762a1bSJed 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, 617c4762a1bSJed 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}; 618c4762a1bSJed Brown PetscInt faultPoints[4] = {7, 8, 16, 17}; 619c4762a1bSJed Brown 6209566063dSJacob Faibussowitsch PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 6219566063dSJacob Faibussowitsch PetscCall(DMPlexInterpolate(*dm, &idm)); 6229566063dSJacob Faibussowitsch for (p = 0; p < 4; ++p) PetscCall(DMSetLabelValue(idm, "fault", faultPoints[p], 1)); 6239566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 0, 1)); 6249566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 1, 1)); 6259566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 2, 2)); 6269566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 3, 2)); 627c4762a1bSJed Brown } 628c4762a1bSJed Brown break; 62998921bdaSJacob Faibussowitsch default: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %d", testNum); 630c4762a1bSJed Brown } 631c4762a1bSJed Brown } else { 632c4762a1bSJed Brown PetscInt numPoints[4] = {0, 0, 0, 0}; 633c4762a1bSJed Brown 6349566063dSJacob Faibussowitsch PetscCall(DMPlexCreateFromDAG(*dm, depth, numPoints, NULL, NULL, NULL, NULL)); 6359566063dSJacob Faibussowitsch PetscCall(DMPlexInterpolate(*dm, &idm)); 6369566063dSJacob Faibussowitsch PetscCall(DMCreateLabel(idm, "fault")); 637c4762a1bSJed Brown } 6389566063dSJacob Faibussowitsch PetscCall(DMViewFromOptions(idm, NULL, "-in_dm_view")); 6399566063dSJacob Faibussowitsch PetscCall(DMDestroy(dm)); 640c4762a1bSJed Brown *dm = idm; 641c4762a1bSJed Brown PetscFunctionReturn(0); 642c4762a1bSJed Brown } 643c4762a1bSJed Brown 644b253942bSMatthew G. Knepley static PetscErrorCode CreateFaultLabel(DM dm) 645b253942bSMatthew G. Knepley { 646b253942bSMatthew G. Knepley DMLabel label; 647b253942bSMatthew G. Knepley PetscInt dim, h, pStart, pEnd, pMax, p; 648b253942bSMatthew G. Knepley 649b253942bSMatthew G. Knepley PetscFunctionBegin; 6509566063dSJacob Faibussowitsch PetscCall(DMGetDimension(dm, &dim)); 6519566063dSJacob Faibussowitsch PetscCall(DMCreateLabel(dm, "cohesive")); 6529566063dSJacob Faibussowitsch PetscCall(DMGetLabel(dm, "cohesive", &label)); 653b253942bSMatthew G. Knepley for (h = 0; h <= dim; ++h) { 6549566063dSJacob Faibussowitsch PetscCall(DMPlexGetSimplexOrBoxCells(dm, h, NULL, &pMax)); 6559566063dSJacob Faibussowitsch PetscCall(DMPlexGetHeightStratum(dm, h, &pStart, &pEnd)); 6569566063dSJacob Faibussowitsch for (p = pMax; p < pEnd; ++p) PetscCall(DMLabelSetValue(label, p, 1)); 657b253942bSMatthew G. Knepley } 658b253942bSMatthew G. Knepley PetscFunctionReturn(0); 659b253942bSMatthew G. Knepley } 660b253942bSMatthew G. Knepley 661b253942bSMatthew G. Knepley static PetscErrorCode CreateDiscretization(DM dm, AppCtx *user) 662b253942bSMatthew G. Knepley { 663b253942bSMatthew G. Knepley PetscFE fe; 664b253942bSMatthew G. Knepley DMLabel fault; 665b7519becSMatthew G. Knepley PetscInt dim, Ncf = user->cohesiveFields, f; 666b253942bSMatthew G. Knepley 667b253942bSMatthew G. Knepley PetscFunctionBegin; 6689566063dSJacob Faibussowitsch PetscCall(DMGetDimension(dm, &dim)); 6699566063dSJacob Faibussowitsch PetscCall(DMGetLabel(dm, "cohesive", &fault)); 6709566063dSJacob Faibussowitsch PetscCall(DMLabelView(fault, PETSC_VIEWER_STDOUT_WORLD)); 671b253942bSMatthew G. Knepley 6729566063dSJacob Faibussowitsch PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, dim, user->cellSimplex, "displacement_", PETSC_DETERMINE, &fe)); 6739566063dSJacob Faibussowitsch PetscCall(PetscFESetName(fe, "displacement")); 6749566063dSJacob Faibussowitsch PetscCall(DMAddField(dm, NULL, (PetscObject) fe)); 6759566063dSJacob Faibussowitsch PetscCall(PetscFEDestroy(&fe)); 676b253942bSMatthew G. Knepley 677b7519becSMatthew G. Knepley if (Ncf > 0) { 6789566063dSJacob Faibussowitsch PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim-1, dim, user->cellSimplex, "faulttraction_", PETSC_DETERMINE, &fe)); 6799566063dSJacob Faibussowitsch PetscCall(PetscFESetName(fe, "fault traction")); 6809566063dSJacob Faibussowitsch PetscCall(DMAddField(dm, fault, (PetscObject) fe)); 6819566063dSJacob Faibussowitsch PetscCall(PetscFEDestroy(&fe)); 682b7519becSMatthew G. Knepley } 683b7519becSMatthew G. Knepley for (f = 1; f < Ncf; ++f) { 684b7519becSMatthew G. Knepley char name[256], opt[256]; 685b7519becSMatthew G. Knepley 6869566063dSJacob Faibussowitsch PetscCall(PetscSNPrintf(name, 256, "fault field %D", f)); 6879566063dSJacob Faibussowitsch PetscCall(PetscSNPrintf(opt, 256, "faultfield_%D_", f)); 6889566063dSJacob Faibussowitsch PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim-1, dim, user->cellSimplex, opt, PETSC_DETERMINE, &fe)); 6899566063dSJacob Faibussowitsch PetscCall(PetscFESetName(fe, name)); 6909566063dSJacob Faibussowitsch PetscCall(DMAddField(dm, fault, (PetscObject) fe)); 6919566063dSJacob Faibussowitsch PetscCall(PetscFEDestroy(&fe)); 692b7519becSMatthew G. Knepley } 693b253942bSMatthew G. Knepley 6949566063dSJacob Faibussowitsch PetscCall(DMCreateDS(dm)); 695b253942bSMatthew G. Knepley PetscFunctionReturn(0); 696b253942bSMatthew G. Knepley } 697b253942bSMatthew G. Knepley 698b253942bSMatthew G. Knepley static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 699c4762a1bSJed Brown { 700c4762a1bSJed Brown PetscInt dim = user->dim; 701c4762a1bSJed Brown PetscBool cellSimplex = user->cellSimplex, hasFault, hasFault2, hasParallelFault; 702c4762a1bSJed Brown PetscMPIInt rank, size; 703dd96b1f9SToby Isaac DMLabel matLabel; 704c4762a1bSJed Brown 705c4762a1bSJed Brown PetscFunctionBegin; 7069566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_rank(comm, &rank)); 7079566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_size(comm, &size)); 7089566063dSJacob Faibussowitsch PetscCall(DMCreate(comm, dm)); 7099566063dSJacob Faibussowitsch PetscCall(DMSetType(*dm, DMPLEX)); 7109566063dSJacob Faibussowitsch PetscCall(DMSetDimension(*dm, dim)); 711c4762a1bSJed Brown switch (dim) { 712c4762a1bSJed Brown case 2: 713c4762a1bSJed Brown if (cellSimplex) { 7149566063dSJacob Faibussowitsch PetscCall(CreateSimplex_2D(comm, user->testNum, dm)); 715c4762a1bSJed Brown } else { 7169566063dSJacob Faibussowitsch PetscCall(CreateQuad_2D(comm, user->testNum, dm)); 717c4762a1bSJed Brown } 718c4762a1bSJed Brown break; 719c4762a1bSJed Brown case 3: 720c4762a1bSJed Brown if (cellSimplex) { 7219566063dSJacob Faibussowitsch PetscCall(CreateSimplex_3D(comm, user, *dm)); 722c4762a1bSJed Brown } else { 7239566063dSJacob Faibussowitsch PetscCall(CreateHex_3D(comm, user->testNum, dm)); 724c4762a1bSJed Brown } 725c4762a1bSJed Brown break; 726c4762a1bSJed Brown default: 72798921bdaSJacob Faibussowitsch SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make hybrid meshes for dimension %d", dim); 728c4762a1bSJed Brown } 7299566063dSJacob Faibussowitsch PetscCall(PetscObjectSetOptionsPrefix((PetscObject) *dm, "orig_")); 7309566063dSJacob Faibussowitsch PetscCall(DMPlexDistributeSetDefault(*dm, PETSC_FALSE)); 7319566063dSJacob Faibussowitsch PetscCall(DMSetFromOptions(*dm)); 7329566063dSJacob Faibussowitsch PetscCall(DMGetLabel(*dm, "material", &matLabel)); 733dd96b1f9SToby Isaac if (matLabel) { 7349566063dSJacob Faibussowitsch PetscCall(DMPlexLabelComplete(*dm, matLabel)); 735dd96b1f9SToby Isaac } 7369566063dSJacob Faibussowitsch PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 7379566063dSJacob Faibussowitsch PetscCall(DMHasLabel(*dm, "fault", &hasFault)); 738c4762a1bSJed Brown if (hasFault) { 739b253942bSMatthew G. Knepley DM dmHybrid = NULL, dmInterface = NULL; 740b253942bSMatthew G. Knepley DMLabel faultLabel, faultBdLabel, hybridLabel, splitLabel; 741c4762a1bSJed Brown 7429566063dSJacob Faibussowitsch PetscCall(DMGetLabel(*dm, "fault", &faultLabel)); 7439566063dSJacob Faibussowitsch PetscCall(DMGetLabel(*dm, "faultBd", &faultBdLabel)); 7449566063dSJacob Faibussowitsch PetscCall(DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, &hybridLabel, &splitLabel, &dmInterface, &dmHybrid)); 7459566063dSJacob Faibussowitsch PetscCall(DMLabelView(hybridLabel, PETSC_VIEWER_STDOUT_WORLD)); 7469566063dSJacob Faibussowitsch PetscCall(DMLabelDestroy(&hybridLabel)); 7479566063dSJacob Faibussowitsch PetscCall(DMLabelView(splitLabel, PETSC_VIEWER_STDOUT_WORLD)); 7489566063dSJacob Faibussowitsch PetscCall(DMLabelDestroy(&splitLabel)); 7499566063dSJacob Faibussowitsch PetscCall(DMViewFromOptions(dmInterface, NULL, "-dm_interface_view")); 7509566063dSJacob Faibussowitsch PetscCall(DMDestroy(&dmInterface)); 7519566063dSJacob Faibussowitsch PetscCall(DMDestroy(dm)); 752c4762a1bSJed Brown *dm = dmHybrid; 753c4762a1bSJed Brown } 7549566063dSJacob Faibussowitsch PetscCall(DMHasLabel(*dm, "fault2", &hasFault2)); 755c4762a1bSJed Brown if (hasFault2) { 756c4762a1bSJed Brown DM dmHybrid = NULL; 757c4762a1bSJed Brown DMLabel faultLabel, faultBdLabel, hybridLabel; 758c4762a1bSJed Brown 7599566063dSJacob Faibussowitsch PetscCall(PetscObjectSetOptionsPrefix((PetscObject) *dm, "faulted_")); 7609566063dSJacob Faibussowitsch PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view_pre")); 7619566063dSJacob Faibussowitsch PetscCall(DMPlexDistributeSetDefault(*dm, PETSC_FALSE)); 7629566063dSJacob Faibussowitsch PetscCall(DMSetFromOptions(*dm)); 7639566063dSJacob Faibussowitsch PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 7649566063dSJacob Faibussowitsch PetscCall(DMGetLabel(*dm, "fault2", &faultLabel)); 7659566063dSJacob Faibussowitsch PetscCall(DMGetLabel(*dm, "fault2Bd", &faultBdLabel)); 7669566063dSJacob Faibussowitsch PetscCall(DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, &hybridLabel, NULL, NULL, &dmHybrid)); 7679566063dSJacob Faibussowitsch PetscCall(DMLabelView(hybridLabel, PETSC_VIEWER_STDOUT_WORLD)); 7689566063dSJacob Faibussowitsch PetscCall(DMLabelDestroy(&hybridLabel)); 7699566063dSJacob Faibussowitsch PetscCall(DMDestroy(dm)); 770c4762a1bSJed Brown *dm = dmHybrid; 771c4762a1bSJed Brown } 772c4762a1bSJed Brown if (user->testPartition && size > 1) { 773c4762a1bSJed Brown PetscPartitioner part; 774c4762a1bSJed Brown PetscInt *sizes = NULL; 775c4762a1bSJed Brown PetscInt *points = NULL; 776c4762a1bSJed Brown 777dd400576SPatrick Sanan if (rank == 0) { 778c4762a1bSJed Brown if (dim == 2 && cellSimplex && size == 2) { 779c4762a1bSJed Brown switch (user->testNum) { 780c4762a1bSJed Brown case 0: { 781c4762a1bSJed Brown PetscInt triSizes_p2[2] = {1, 2}; 782c4762a1bSJed Brown PetscInt triPoints_p2[3] = {0, 1, 2}; 783c4762a1bSJed Brown 7849566063dSJacob Faibussowitsch PetscCall(PetscMalloc2(2, &sizes, 3, &points)); 7859566063dSJacob Faibussowitsch PetscCall(PetscArraycpy(sizes, triSizes_p2, 2)); 7869566063dSJacob Faibussowitsch PetscCall(PetscArraycpy(points, triPoints_p2, 3));break;} 787c4762a1bSJed Brown default: 78898921bdaSJacob Faibussowitsch SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %d for triangular mesh on 2 procs", user->testNum); 789c4762a1bSJed Brown } 790c4762a1bSJed Brown } else if (dim == 2 && !cellSimplex && size == 2) { 791c4762a1bSJed Brown switch (user->testNum) { 792c4762a1bSJed Brown case 0: { 793c4762a1bSJed Brown PetscInt quadSizes_p2[2] = {1, 2}; 794c4762a1bSJed Brown PetscInt quadPoints_p2[3] = {0, 1, 2}; 795c4762a1bSJed Brown 7969566063dSJacob Faibussowitsch PetscCall(PetscMalloc2(2, &sizes, 3, &points)); 7979566063dSJacob Faibussowitsch PetscCall(PetscArraycpy(sizes, quadSizes_p2, 2)); 7989566063dSJacob Faibussowitsch PetscCall(PetscArraycpy(points, quadPoints_p2, 3));break;} 799c4762a1bSJed Brown case 2: { 800c4762a1bSJed Brown PetscInt quadSizes_p2[2] = {1, 1}; 801c4762a1bSJed Brown PetscInt quadPoints_p2[2] = {0, 1}; 802c4762a1bSJed Brown 8039566063dSJacob Faibussowitsch PetscCall(PetscMalloc2(2, &sizes, 2, &points)); 8049566063dSJacob Faibussowitsch PetscCall(PetscArraycpy(sizes, quadSizes_p2, 2)); 8059566063dSJacob Faibussowitsch PetscCall(PetscArraycpy(points, quadPoints_p2, 2));break;} 806c4762a1bSJed Brown default: 80798921bdaSJacob Faibussowitsch SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %d for quadrilateral mesh on 2 procs", user->testNum); 808c4762a1bSJed Brown } 809c4762a1bSJed Brown } else if (dim == 3 && cellSimplex && size == 2) { 810c4762a1bSJed Brown switch (user->testNum) { 811c4762a1bSJed Brown case 0: { 812c4762a1bSJed Brown PetscInt tetSizes_p2[2] = {1, 2}; 813c4762a1bSJed Brown PetscInt tetPoints_p2[3] = {0, 1, 2}; 814c4762a1bSJed Brown 8159566063dSJacob Faibussowitsch PetscCall(PetscMalloc2(2, &sizes, 3, &points)); 8169566063dSJacob Faibussowitsch PetscCall(PetscArraycpy(sizes, tetSizes_p2, 2)); 8179566063dSJacob Faibussowitsch PetscCall(PetscArraycpy(points, tetPoints_p2, 3));break;} 818c4762a1bSJed Brown default: 81998921bdaSJacob Faibussowitsch SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %d for triangular mesh on 2 procs", user->testNum); 820c4762a1bSJed Brown } 821c4762a1bSJed Brown } else if (dim == 3 && !cellSimplex && size == 2) { 822c4762a1bSJed Brown switch (user->testNum) { 823c4762a1bSJed Brown case 0: { 824c4762a1bSJed Brown PetscInt hexSizes_p2[2] = {1, 2}; 825c4762a1bSJed Brown PetscInt hexPoints_p2[3] = {0, 1, 2}; 826c4762a1bSJed Brown 8279566063dSJacob Faibussowitsch PetscCall(PetscMalloc2(2, &sizes, 3, &points)); 8289566063dSJacob Faibussowitsch PetscCall(PetscArraycpy(sizes, hexSizes_p2, 2)); 8299566063dSJacob Faibussowitsch PetscCall(PetscArraycpy(points, hexPoints_p2, 3));break;} 830c4762a1bSJed Brown default: 83198921bdaSJacob Faibussowitsch SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %d for triangular mesh on 2 procs", user->testNum); 832c4762a1bSJed Brown } 833c4762a1bSJed Brown } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test partition"); 834c4762a1bSJed Brown } 8359566063dSJacob Faibussowitsch PetscCall(DMPlexGetPartitioner(*dm, &part)); 8369566063dSJacob Faibussowitsch PetscCall(PetscPartitionerSetType(part, PETSCPARTITIONERSHELL)); 8379566063dSJacob Faibussowitsch PetscCall(PetscPartitionerShellSetPartition(part, size, sizes, points)); 8389566063dSJacob Faibussowitsch PetscCall(PetscFree2(sizes, points)); 839c4762a1bSJed Brown } 840c4762a1bSJed Brown { 841c4762a1bSJed Brown DM pdm = NULL; 842c4762a1bSJed Brown 843c4762a1bSJed Brown /* Distribute mesh over processes */ 8449566063dSJacob Faibussowitsch PetscCall(DMPlexDistribute(*dm, 0, NULL, &pdm)); 845c4762a1bSJed Brown if (pdm) { 8469566063dSJacob Faibussowitsch PetscCall(DMViewFromOptions(pdm, NULL, "-dm_view")); 8479566063dSJacob Faibussowitsch PetscCall(DMDestroy(dm)); 848c4762a1bSJed Brown *dm = pdm; 849c4762a1bSJed Brown } 850c4762a1bSJed Brown } 8519566063dSJacob Faibussowitsch PetscCall(DMHasLabel(*dm, "pfault", &hasParallelFault)); 852c4762a1bSJed Brown if (hasParallelFault) { 853c4762a1bSJed Brown DM dmHybrid = NULL; 854c4762a1bSJed Brown DMLabel faultLabel, faultBdLabel, hybridLabel; 855c4762a1bSJed Brown 8569566063dSJacob Faibussowitsch PetscCall(DMGetLabel(*dm, "pfault", &faultLabel)); 8579566063dSJacob Faibussowitsch PetscCall(DMGetLabel(*dm, "pfaultBd", &faultBdLabel)); 8589566063dSJacob Faibussowitsch PetscCall(DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, &hybridLabel, NULL, NULL, &dmHybrid)); 8599566063dSJacob Faibussowitsch PetscCall(DMLabelView(hybridLabel, PETSC_VIEWER_STDOUT_WORLD)); 8609566063dSJacob Faibussowitsch PetscCall(DMLabelDestroy(&hybridLabel)); 8619566063dSJacob Faibussowitsch PetscCall(DMDestroy(dm)); 862c4762a1bSJed Brown *dm = dmHybrid; 863c4762a1bSJed Brown } 8649566063dSJacob Faibussowitsch PetscCall(PetscObjectSetName((PetscObject) *dm, "Hybrid Mesh")); 8659566063dSJacob Faibussowitsch PetscCall(CreateFaultLabel(*dm)); 8669566063dSJacob Faibussowitsch PetscCall(CreateDiscretization(*dm, user)); 8679566063dSJacob Faibussowitsch PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view_pre")); 8689566063dSJacob Faibussowitsch PetscCall(DMPlexDistributeSetDefault(*dm, PETSC_FALSE)); 8699566063dSJacob Faibussowitsch PetscCall(DMSetFromOptions(*dm)); 8709566063dSJacob Faibussowitsch PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 871c4762a1bSJed Brown PetscFunctionReturn(0); 872c4762a1bSJed Brown } 873c4762a1bSJed Brown 874b253942bSMatthew G. Knepley static PetscErrorCode TestMesh(DM dm, AppCtx *user) 875b253942bSMatthew G. Knepley { 876b253942bSMatthew G. Knepley PetscFunctionBegin; 8779566063dSJacob Faibussowitsch PetscCall(DMPlexCheckSymmetry(dm)); 8789566063dSJacob Faibussowitsch PetscCall(DMPlexCheckSkeleton(dm, 0)); 8799566063dSJacob Faibussowitsch PetscCall(DMPlexCheckFaces(dm, 0)); 880b253942bSMatthew G. Knepley PetscFunctionReturn(0); 881b253942bSMatthew G. Knepley } 882b253942bSMatthew G. Knepley 883b253942bSMatthew G. Knepley static PetscErrorCode TestDiscretization(DM dm, AppCtx *user) 884b253942bSMatthew G. Knepley { 885b253942bSMatthew G. Knepley PetscSection s; 886b253942bSMatthew G. Knepley 887b253942bSMatthew G. Knepley PetscFunctionBegin; 8889566063dSJacob Faibussowitsch PetscCall(DMGetSection(dm, &s)); 8899566063dSJacob Faibussowitsch PetscCall(PetscObjectViewFromOptions((PetscObject) s, NULL, "-local_section_view")); 890b253942bSMatthew G. Knepley PetscFunctionReturn(0); 891b253942bSMatthew G. Knepley } 892b253942bSMatthew G. Knepley 893b253942bSMatthew G. Knepley static PetscErrorCode r(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 894b253942bSMatthew G. Knepley { 895b253942bSMatthew G. Knepley PetscInt d; 896b253942bSMatthew G. Knepley for (d = 0; d < dim; ++d) u[d] = x[d]; 897b253942bSMatthew G. Knepley return 0; 898b253942bSMatthew G. Knepley } 899b253942bSMatthew G. Knepley 900b253942bSMatthew G. Knepley static PetscErrorCode rp1(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 901b253942bSMatthew G. Knepley { 902b253942bSMatthew G. Knepley PetscInt d; 903b253942bSMatthew G. Knepley for (d = 0; d < dim; ++d) u[d] = x[d] + (d > 0 ? 1.0 : 0.0); 904b253942bSMatthew G. Knepley return 0; 905b253942bSMatthew G. Knepley } 906b253942bSMatthew G. Knepley 907b253942bSMatthew G. Knepley static PetscErrorCode phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 908b253942bSMatthew G. Knepley { 909b253942bSMatthew G. Knepley PetscInt d; 910b253942bSMatthew G. Knepley u[0] = -x[1]; 911b253942bSMatthew G. Knepley u[1] = x[0]; 912b253942bSMatthew G. Knepley for (d = 2; d < dim; ++d) u[d] = x[d]; 913b253942bSMatthew G. Knepley return 0; 914b253942bSMatthew G. Knepley } 915b253942bSMatthew G. Knepley 916b253942bSMatthew G. Knepley /* \lambda \cdot (\psi_u^- - \psi_u^+) */ 917b253942bSMatthew G. Knepley static void f0_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 918b253942bSMatthew G. Knepley const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 919b253942bSMatthew G. Knepley const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 920b253942bSMatthew G. Knepley PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 921b253942bSMatthew G. Knepley { 922b253942bSMatthew G. Knepley const PetscInt Nc = uOff[1]-uOff[0]; 923b253942bSMatthew G. Knepley PetscInt c; 924b253942bSMatthew G. Knepley for (c = 0; c < Nc; ++c) { 925b253942bSMatthew G. Knepley f0[c] = u[Nc*2+c] + x[Nc-c-1]; 926b253942bSMatthew G. Knepley f0[Nc+c] = -(u[Nc*2+c] + x[Nc-c-1]); 927b253942bSMatthew G. Knepley } 928b253942bSMatthew G. Knepley } 929b253942bSMatthew G. Knepley 930b253942bSMatthew G. Knepley /* (d - u^+ + u^-) \cdot \psi_\lambda */ 931b253942bSMatthew G. Knepley static void f0_bd_l(PetscInt dim, PetscInt Nf, PetscInt NfAux, 932b253942bSMatthew G. Knepley const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 933b253942bSMatthew G. Knepley const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 934b253942bSMatthew G. Knepley PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 935b253942bSMatthew G. Knepley { 936b253942bSMatthew G. Knepley const PetscInt Nc = uOff[2]-uOff[1]; 937b253942bSMatthew G. Knepley PetscInt c; 938b253942bSMatthew G. Knepley 939b253942bSMatthew G. Knepley for (c = 0; c < Nc; ++c) f0[c] = (c > 0 ? 1.0 : 0.0) + u[c] - u[Nc+c]; 940b253942bSMatthew G. Knepley } 941b253942bSMatthew G. Knepley 942b253942bSMatthew G. Knepley /* \psi_lambda \cdot (\psi_u^- - \psi_u^+) */ 943b253942bSMatthew G. Knepley static void g0_bd_ul(PetscInt dim, PetscInt Nf, PetscInt NfAux, 944b253942bSMatthew G. Knepley const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 945b253942bSMatthew G. Knepley const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 946b253942bSMatthew G. Knepley PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 947b253942bSMatthew G. Knepley { 948b253942bSMatthew G. Knepley const PetscInt Nc = uOff[1]-uOff[0]; 949b253942bSMatthew G. Knepley PetscInt c; 950b253942bSMatthew G. Knepley 951b253942bSMatthew G. Knepley for (c = 0; c < Nc; ++c) { 952b253942bSMatthew G. Knepley g0[(0 +c)*Nc+c] = 1.0; 953b253942bSMatthew G. Knepley g0[(Nc+c)*Nc+c] = -1.0; 954b253942bSMatthew G. Knepley } 955b253942bSMatthew G. Knepley } 956b253942bSMatthew G. Knepley 957b253942bSMatthew G. Knepley /* (-\psi_u^+ + \psi_u^-) \cdot \psi_\lambda */ 958b253942bSMatthew G. Knepley static void g0_bd_lu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 959b253942bSMatthew G. Knepley const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 960b253942bSMatthew G. Knepley const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 961b253942bSMatthew G. Knepley PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 962b253942bSMatthew G. Knepley { 963b253942bSMatthew G. Knepley const PetscInt Nc = uOff[2]-uOff[1]; 964b253942bSMatthew G. Knepley PetscInt c; 965b253942bSMatthew G. Knepley 966b253942bSMatthew G. Knepley for (c = 0; c < Nc; ++c) { 967b253942bSMatthew G. Knepley g0[c*Nc*2+c] = 1.0; 968b253942bSMatthew G. Knepley g0[c*Nc*2+Nc+c] = -1.0; 969b253942bSMatthew G. Knepley } 970b253942bSMatthew G. Knepley } 971b253942bSMatthew G. Knepley 972b253942bSMatthew G. Knepley static PetscErrorCode TestAssembly(DM dm, AppCtx *user) 973b253942bSMatthew G. Knepley { 974b253942bSMatthew G. Knepley Mat J; 975b253942bSMatthew G. Knepley Vec locX, locF; 976b253942bSMatthew G. Knepley PetscDS probh; 977b253942bSMatthew G. Knepley DMLabel fault, material; 978b253942bSMatthew G. Knepley IS cohesiveCells; 979b7519becSMatthew G. Knepley PetscWeakForm wf; 98006ad1575SMatthew G. Knepley PetscFormKey keys[3]; 981b253942bSMatthew G. Knepley PetscErrorCode (*initialGuess[2])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx); 982b253942bSMatthew G. Knepley PetscInt dim, Nf, cMax, cEnd, id; 983b7519becSMatthew G. Knepley PetscMPIInt rank; 984b253942bSMatthew G. Knepley 985b253942bSMatthew G. Knepley PetscFunctionBegin; 9869566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank)); 9879566063dSJacob Faibussowitsch PetscCall(DMGetDimension(dm, &dim)); 9889566063dSJacob Faibussowitsch PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, NULL, &cMax)); 9899566063dSJacob Faibussowitsch PetscCall(DMPlexGetHeightStratum(dm, 0, NULL, &cEnd)); 9909566063dSJacob Faibussowitsch PetscCall(ISCreateStride(PETSC_COMM_SELF, cEnd - cMax, cMax, 1, &cohesiveCells)); 9919566063dSJacob Faibussowitsch PetscCall(DMGetLabel(dm, "cohesive", &fault)); 9929566063dSJacob Faibussowitsch PetscCall(DMGetLocalVector(dm, &locX)); 9939566063dSJacob Faibussowitsch PetscCall(PetscObjectSetName((PetscObject) locX, "Local Solution")); 9949566063dSJacob Faibussowitsch PetscCall(DMGetLocalVector(dm, &locF)); 9959566063dSJacob Faibussowitsch PetscCall(PetscObjectSetName((PetscObject) locF, "Local Residual")); 9969566063dSJacob Faibussowitsch PetscCall(DMCreateMatrix(dm, &J)); 9979566063dSJacob Faibussowitsch PetscCall(PetscObjectSetName((PetscObject) J, "Jacobian")); 998b253942bSMatthew G. Knepley 999b253942bSMatthew G. Knepley /* The initial guess has displacement shifted by one unit in each fault parallel direction across the fault */ 10009566063dSJacob Faibussowitsch PetscCall(DMGetLabel(dm, "material", &material)); 1001b253942bSMatthew G. Knepley id = 1; 1002b253942bSMatthew G. Knepley initialGuess[0] = r; 1003b253942bSMatthew G. Knepley initialGuess[1] = NULL; 10049566063dSJacob Faibussowitsch PetscCall(DMProjectFunctionLabelLocal(dm, 0.0, material, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX)); 1005b253942bSMatthew G. Knepley id = 2; 1006b253942bSMatthew G. Knepley initialGuess[0] = rp1; 1007b253942bSMatthew G. Knepley initialGuess[1] = NULL; 10089566063dSJacob Faibussowitsch PetscCall(DMProjectFunctionLabelLocal(dm, 0.0, material, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX)); 1009b253942bSMatthew G. Knepley id = 1; 1010b253942bSMatthew G. Knepley initialGuess[0] = NULL; 1011b253942bSMatthew G. Knepley initialGuess[1] = phi; 10129566063dSJacob Faibussowitsch PetscCall(DMProjectFunctionLabelLocal(dm, 0.0, fault, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX)); 10139566063dSJacob Faibussowitsch PetscCall(VecViewFromOptions(locX, NULL, "-local_solution_view")); 1014b253942bSMatthew G. Knepley 10159566063dSJacob Faibussowitsch PetscCall(DMGetCellDS(dm, cMax, &probh)); 10169566063dSJacob Faibussowitsch PetscCall(PetscDSGetWeakForm(probh, &wf)); 10179566063dSJacob Faibussowitsch PetscCall(PetscDSGetNumFields(probh, &Nf)); 10189566063dSJacob Faibussowitsch PetscCall(PetscWeakFormSetIndexBdResidual(wf, material, 1, 0, 0, 0, f0_bd_u, 0, NULL)); 10199566063dSJacob Faibussowitsch PetscCall(PetscWeakFormSetIndexBdResidual(wf, material, 2, 0, 0, 0, f0_bd_u, 0, NULL)); 10209566063dSJacob Faibussowitsch PetscCall(PetscWeakFormSetIndexBdJacobian(wf, material, 1, 0, 1, 0, 0, g0_bd_ul, 0, NULL, 0, NULL, 0, NULL)); 10219566063dSJacob Faibussowitsch PetscCall(PetscWeakFormSetIndexBdJacobian(wf, material, 2, 0, 1, 0, 0, g0_bd_ul, 0, NULL, 0, NULL, 0, NULL)); 1022b7519becSMatthew G. Knepley if (Nf > 1) { 10239566063dSJacob Faibussowitsch PetscCall(PetscWeakFormSetIndexBdResidual(wf, fault, 1, 1, 0, 0, f0_bd_l, 0, NULL)); 10249566063dSJacob Faibussowitsch PetscCall(PetscWeakFormSetIndexBdJacobian(wf, fault, 1, 1, 0, 0, 0, g0_bd_lu, 0, NULL, 0, NULL, 0, NULL)); 1025b7519becSMatthew G. Knepley } 10269566063dSJacob Faibussowitsch if (!rank) PetscCall(PetscDSView(probh, NULL)); 1027b253942bSMatthew G. Knepley 1028dd96b1f9SToby Isaac keys[0].label = NULL; 1029dd96b1f9SToby Isaac keys[0].value = 0; 10306528b96dSMatthew G. Knepley keys[0].field = 0; 1031dd96b1f9SToby Isaac keys[0].part = 0; 1032b7519becSMatthew G. Knepley keys[1].label = material; 1033b7519becSMatthew G. Knepley keys[1].value = 2; 10346528b96dSMatthew G. Knepley keys[1].field = 0; 1035dd96b1f9SToby Isaac keys[1].part = 0; 1036b7519becSMatthew G. Knepley keys[2].label = fault; 1037b7519becSMatthew G. Knepley keys[2].value = 1; 1038b7519becSMatthew G. Knepley keys[2].field = 1; 1039dd96b1f9SToby Isaac keys[2].part = 0; 10409566063dSJacob Faibussowitsch PetscCall(VecSet(locF, 0.)); 10419566063dSJacob Faibussowitsch PetscCall(DMPlexComputeResidual_Hybrid_Internal(dm, keys, cohesiveCells, 0.0, locX, NULL, 0.0, locF, user)); 10429566063dSJacob Faibussowitsch PetscCall(VecViewFromOptions(locF, NULL, "-local_residual_view")); 10439566063dSJacob Faibussowitsch PetscCall(MatZeroEntries(J)); 10449566063dSJacob Faibussowitsch PetscCall(DMPlexComputeJacobian_Hybrid_Internal(dm, keys, cohesiveCells, 0.0, 0.0, locX, NULL, J, J, user)); 10459566063dSJacob Faibussowitsch PetscCall(MatViewFromOptions(J, NULL, "-local_jacobian_view")); 1046b253942bSMatthew G. Knepley 10479566063dSJacob Faibussowitsch PetscCall(DMRestoreLocalVector(dm, &locX)); 10489566063dSJacob Faibussowitsch PetscCall(DMRestoreLocalVector(dm, &locF)); 10499566063dSJacob Faibussowitsch PetscCall(MatDestroy(&J)); 10509566063dSJacob Faibussowitsch PetscCall(ISDestroy(&cohesiveCells)); 1051b253942bSMatthew G. Knepley PetscFunctionReturn(0); 1052b253942bSMatthew G. Knepley } 1053b253942bSMatthew G. Knepley 1054c4762a1bSJed Brown int main(int argc, char **argv) 1055c4762a1bSJed Brown { 1056c4762a1bSJed Brown DM dm; 1057c4762a1bSJed Brown AppCtx user; /* user-defined work context */ 1058c4762a1bSJed Brown 10599566063dSJacob Faibussowitsch PetscCall(PetscInitialize(&argc, &argv, NULL,help)); 10609566063dSJacob Faibussowitsch PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user)); 10619566063dSJacob Faibussowitsch PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm)); 10629566063dSJacob Faibussowitsch PetscCall(TestMesh(dm, &user)); 10639566063dSJacob Faibussowitsch PetscCall(TestDiscretization(dm, &user)); 10649566063dSJacob Faibussowitsch PetscCall(TestAssembly(dm, &user)); 10659566063dSJacob Faibussowitsch PetscCall(DMDestroy(&dm)); 10669566063dSJacob Faibussowitsch PetscCall(PetscFinalize()); 1067b122ec5aSJacob Faibussowitsch return 0; 1068c4762a1bSJed Brown } 1069c4762a1bSJed Brown 1070c4762a1bSJed Brown /*TEST 1071c4762a1bSJed Brown testset: 1072b253942bSMatthew G. Knepley args: -orig_dm_plex_check_all -dm_plex_check_all \ 1073b7519becSMatthew G. Knepley -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1 -local_section_view \ 1074b253942bSMatthew G. Knepley -local_solution_view -local_residual_view -local_jacobian_view 1075c4762a1bSJed Brown test: 1076c4762a1bSJed Brown suffix: tri_0 1077c4762a1bSJed Brown args: -dim 2 10781c6715b8SMatthew G. Knepley filter: sed -e "s/_start//g" -e "s/f0_bd_u//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul//g" -e "s/g0_bd_lu//g" 1079c4762a1bSJed Brown test: 1080c4762a1bSJed Brown suffix: tri_t1_0 1081c4762a1bSJed Brown args: -dim 2 -test_num 1 10821c6715b8SMatthew G. Knepley filter: sed -e "s/_start//g" -e "s/f0_bd_u//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul//g" -e "s/g0_bd_lu//g" 1083c4762a1bSJed Brown test: 1084c4762a1bSJed Brown suffix: tet_0 1085c4762a1bSJed Brown args: -dim 3 10861c6715b8SMatthew G. Knepley filter: sed -e "s/_start//g" -e "s/f0_bd_u//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul//g" -e "s/g0_bd_lu//g" 1087c4762a1bSJed Brown test: 1088b253942bSMatthew G. Knepley suffix: tet_t1_0 1089b253942bSMatthew G. Knepley args: -dim 3 -test_num 1 10901c6715b8SMatthew G. Knepley filter: sed -e "s/_start//g" -e "s/f0_bd_u//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul//g" -e "s/g0_bd_lu//g" 1091b253942bSMatthew G. Knepley 1092b253942bSMatthew G. Knepley testset: 1093b253942bSMatthew G. Knepley args: -orig_dm_plex_check_all -dm_plex_check_all \ 1094b7519becSMatthew G. Knepley -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1 1095b253942bSMatthew G. Knepley test: 1096c4762a1bSJed Brown suffix: tet_1 1097c4762a1bSJed Brown nsize: 2 1098c4762a1bSJed Brown args: -dim 3 10991c6715b8SMatthew G. Knepley filter: sed -e "s/_start//g" -e "s/f0_bd_u//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul//g" -e "s/g0_bd_lu//g" 1100c4762a1bSJed Brown test: 1101b253942bSMatthew G. Knepley suffix: tri_1 1102b253942bSMatthew G. Knepley nsize: 2 1103b253942bSMatthew G. Knepley args: -dim 2 11041c6715b8SMatthew G. Knepley filter: sed -e "s/_start//g" -e "s/f0_bd_u//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul//g" -e "s/g0_bd_lu//g" 1105c4762a1bSJed Brown 1106c4762a1bSJed Brown testset: 1107ecfb78b5SMatthew G. Knepley args: -orig_dm_plex_check_all -dm_plex_check_all \ 1108b7519becSMatthew G. Knepley -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1 1109c4762a1bSJed Brown # 2D Quads 1110c4762a1bSJed Brown test: 1111c4762a1bSJed Brown suffix: quad_0 1112c4762a1bSJed Brown args: -dim 2 -cell_simplex 0 11131c6715b8SMatthew G. Knepley filter: sed -e "s/_start//g" -e "s/f0_bd_u//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul//g" -e "s/g0_bd_lu//g" 1114c4762a1bSJed Brown test: 1115c4762a1bSJed Brown suffix: quad_1 1116c4762a1bSJed Brown nsize: 2 1117c4762a1bSJed Brown args: -dim 2 -cell_simplex 0 11181c6715b8SMatthew G. Knepley filter: sed -e "s/_start//g" -e "s/f0_bd_u//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul//g" -e "s/g0_bd_lu//g" 1119c4762a1bSJed Brown test: 1120c4762a1bSJed Brown suffix: quad_t1_0 112154fcfd0cSMatthew G. Knepley args: -dim 2 -cell_simplex 0 -test_num 1 -faulted_dm_plex_check_all 11221c6715b8SMatthew G. Knepley filter: sed -e "s/_start//g" -e "s/f0_bd_u//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul//g" -e "s/g0_bd_lu//g" 1123c4762a1bSJed Brown # 3D Hex 1124c4762a1bSJed Brown test: 1125c4762a1bSJed Brown suffix: hex_0 1126c4762a1bSJed Brown args: -dim 3 -cell_simplex 0 11271c6715b8SMatthew G. Knepley filter: sed -e "s/_start//g" -e "s/f0_bd_u//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul//g" -e "s/g0_bd_lu//g" 1128c4762a1bSJed Brown test: 1129c4762a1bSJed Brown suffix: hex_1 1130c4762a1bSJed Brown nsize: 2 1131c4762a1bSJed Brown args: -dim 3 -cell_simplex 0 11321c6715b8SMatthew G. Knepley filter: sed -e "s/_start//g" -e "s/f0_bd_u//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul//g" -e "s/g0_bd_lu//g" 1133c4762a1bSJed Brown test: 1134c4762a1bSJed Brown suffix: hex_t1_0 1135c4762a1bSJed Brown args: -dim 3 -cell_simplex 0 -test_num 1 11361c6715b8SMatthew G. Knepley filter: sed -e "s/_start//g" -e "s/f0_bd_u//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul//g" -e "s/g0_bd_lu//g" 1137c4762a1bSJed Brown test: 1138c4762a1bSJed Brown suffix: hex_t2_0 1139c4762a1bSJed Brown args: -dim 3 -cell_simplex 0 -test_num 2 11401c6715b8SMatthew G. Knepley filter: sed -e "s/_start//g" -e "s/f0_bd_u//g" -e "s/f0_bd_l//g" -e "s/g0_bd_ul//g" -e "s/g0_bd_lu//g" 1141c4762a1bSJed Brown 1142c4762a1bSJed Brown TEST*/ 1143