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 3249566063dSJacob Faibussowitsch ierr = PetscOptionsBegin(comm, "", "Meshing Problem Options", "DMPLEX");PetscCall(ierr); 3259566063dSJacob Faibussowitsch PetscCall(PetscOptionsBoundedInt("-debug", "The debugging level", "ex5.c", options->debug, &options->debug, NULL,0)); 3269566063dSJacob Faibussowitsch PetscCall(PetscOptionsRangeInt("-dim", "The topological mesh dimension", "ex5.c", options->dim, &options->dim, NULL,1,3)); 3279566063dSJacob Faibussowitsch PetscCall(PetscOptionsBool("-cell_simplex", "Use simplices if true, otherwise hexes", "ex5.c", options->cellSimplex, &options->cellSimplex, NULL)); 3289566063dSJacob Faibussowitsch PetscCall(PetscOptionsBool("-test_partition", "Use a fixed partition for testing", "ex5.c", options->testPartition, &options->testPartition, NULL)); 3299566063dSJacob Faibussowitsch PetscCall(PetscOptionsBoundedInt("-test_num", "The particular mesh to test", "ex5.c", options->testNum, &options->testNum, NULL,0)); 3309566063dSJacob Faibussowitsch PetscCall(PetscOptionsBoundedInt("-cohesive_fields", "The number of cohesive fields", "ex5.c", options->cohesiveFields, &options->cohesiveFields, NULL, 0)); 3319566063dSJacob Faibussowitsch ierr = PetscOptionsEnd();PetscCall(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; 3429566063dSJacob Faibussowitsch PetscCallMPI(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 3559566063dSJacob Faibussowitsch PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3569566063dSJacob Faibussowitsch for (p = 0; p < 4; ++p) PetscCall(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1])); 3579566063dSJacob Faibussowitsch for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "fault", faultPoints[p], 1)); 3589566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 0, 1)); 3599566063dSJacob Faibussowitsch PetscCall(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 3729566063dSJacob Faibussowitsch PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 3739566063dSJacob Faibussowitsch for (p = 0; p < 3; ++p) PetscCall(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1])); 3749566063dSJacob Faibussowitsch for (p = 0; p < 3; ++p) PetscCall(DMSetLabelValue(*dm, "fault", faultPoints[p], 1)); 3759566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 0, 1));PetscCall(DMSetLabelValue(*dm, "material", 3, 1)); 3769566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 1, 2));PetscCall(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 3859566063dSJacob Faibussowitsch PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, NULL, NULL, NULL, NULL)); 3869566063dSJacob Faibussowitsch PetscCall(DMCreateLabel(*dm, "fault")); 387c4762a1bSJed Brown } 3889566063dSJacob Faibussowitsch PetscCall(DMPlexInterpolate(*dm, &idm)); 3899566063dSJacob Faibussowitsch PetscCall(DMViewFromOptions(idm, NULL, "-in_dm_view")); 3909566063dSJacob Faibussowitsch PetscCall(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; 4019566063dSJacob Faibussowitsch PetscCallMPI(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 4149566063dSJacob Faibussowitsch PetscCall(DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 415c4762a1bSJed Brown for (p = 0; p < 10; ++p) { 4169566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1])); 417c4762a1bSJed Brown } 418c4762a1bSJed Brown for (p = 0; p < 3; ++p) { 4199566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(dm, "fault", faultPoints[p], 1)); 420c4762a1bSJed Brown } 4219566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(dm, "material", 0, 1)); 4229566063dSJacob Faibussowitsch PetscCall(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 4359566063dSJacob Faibussowitsch PetscCall(DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 436c4762a1bSJed Brown for (p = 0; p < 7; ++p) { 4379566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1])); 438c4762a1bSJed Brown } 439c4762a1bSJed Brown for (p = 0; p < 4; ++p) { 4409566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(dm, "fault", faultPoints[p], 1)); 441c4762a1bSJed Brown } 4429566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(dm, "material", 0, 1)); 4439566063dSJacob Faibussowitsch PetscCall(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 4529566063dSJacob Faibussowitsch PetscCall(DMPlexCreateFromDAG(dm, depth, numPoints, NULL, NULL, NULL, NULL)); 4539566063dSJacob Faibussowitsch PetscCall(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; 4659566063dSJacob Faibussowitsch PetscCallMPI(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 4799566063dSJacob Faibussowitsch PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 4809566063dSJacob Faibussowitsch for (p = 0; p < 6; ++p) PetscCall(DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1])); 4819566063dSJacob Faibussowitsch if (testNum == 0) for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "fault", faultPoints[p], 1)); 4829566063dSJacob Faibussowitsch if (testNum == 2) for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "pfault", faultPoints[p], 1)); 4839566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 0, 1)); 4849566063dSJacob Faibussowitsch PetscCall(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 5069566063dSJacob Faibussowitsch PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 5079566063dSJacob Faibussowitsch for (p = 0; p < 3; ++p) PetscCall(DMSetLabelValue(*dm, "fault", faultPoints[p], 1)); 5089566063dSJacob Faibussowitsch for (p = 0; p < 2; ++p) PetscCall(DMSetLabelValue(*dm, "fault2", fault2Points[p], 1)); 5099566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 0, 1)); 5109566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 1, 1)); 5119566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 2, 1)); 5129566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 3, 1)); 5139566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 4, 1)); 5149566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 5, 2)); 5159566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 6, 2)); 5169566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 7, 2)); 5179566063dSJacob Faibussowitsch PetscCall(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 5269566063dSJacob Faibussowitsch PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, NULL, NULL, NULL, NULL)); 5279566063dSJacob Faibussowitsch if (testNum == 2) PetscCall(DMCreateLabel(*dm, "pfault")); 5289566063dSJacob Faibussowitsch else PetscCall(DMCreateLabel(*dm, "fault")); 529c4762a1bSJed Brown } 5309566063dSJacob Faibussowitsch PetscCall(DMPlexInterpolate(*dm, &idm)); 5319566063dSJacob Faibussowitsch PetscCall(DMViewFromOptions(idm, NULL, "-in_dm_view")); 5329566063dSJacob Faibussowitsch PetscCall(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; 5449566063dSJacob Faibussowitsch PetscCallMPI(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 5599566063dSJacob Faibussowitsch PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 5609566063dSJacob Faibussowitsch PetscCall(DMPlexInterpolate(*dm, &idm)); 5619566063dSJacob Faibussowitsch for (p = 0; p < 8; ++p) PetscCall(DMSetLabelValue(idm, "marker", markerPoints[p*2], markerPoints[p*2+1])); 5629566063dSJacob Faibussowitsch for (p = 0; p < 4; ++p) PetscCall(DMSetLabelValue(idm, "fault", faultPoints[p], 1)); 5639566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 0, 1)); 5649566063dSJacob Faibussowitsch PetscCall(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 5959566063dSJacob Faibussowitsch PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 5969566063dSJacob Faibussowitsch PetscCall(DMPlexInterpolate(*dm, &idm)); 5979566063dSJacob Faibussowitsch for (p = 0; p < 6; ++p) PetscCall(DMSetLabelValue(idm, "fault", faultPoints[p], 1)); 5989566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 0, 1)); 5999566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 1, 1)); 6009566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 2, 1)); 6019566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 3, 2)); 6029566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 4, 2)); 6039566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 5, 2)); 6049566063dSJacob Faibussowitsch PetscCall(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 6229566063dSJacob Faibussowitsch PetscCall(DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords)); 6239566063dSJacob Faibussowitsch PetscCall(DMPlexInterpolate(*dm, &idm)); 6249566063dSJacob Faibussowitsch for (p = 0; p < 4; ++p) PetscCall(DMSetLabelValue(idm, "fault", faultPoints[p], 1)); 6259566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 0, 1)); 6269566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 1, 1)); 6279566063dSJacob Faibussowitsch PetscCall(DMSetLabelValue(*dm, "material", 2, 2)); 6289566063dSJacob Faibussowitsch PetscCall(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 6369566063dSJacob Faibussowitsch PetscCall(DMPlexCreateFromDAG(*dm, depth, numPoints, NULL, NULL, NULL, NULL)); 6379566063dSJacob Faibussowitsch PetscCall(DMPlexInterpolate(*dm, &idm)); 6389566063dSJacob Faibussowitsch PetscCall(DMCreateLabel(idm, "fault")); 639c4762a1bSJed Brown } 6409566063dSJacob Faibussowitsch PetscCall(DMViewFromOptions(idm, NULL, "-in_dm_view")); 6419566063dSJacob Faibussowitsch PetscCall(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; 6529566063dSJacob Faibussowitsch PetscCall(DMGetDimension(dm, &dim)); 6539566063dSJacob Faibussowitsch PetscCall(DMCreateLabel(dm, "cohesive")); 6549566063dSJacob Faibussowitsch PetscCall(DMGetLabel(dm, "cohesive", &label)); 655b253942bSMatthew G. Knepley for (h = 0; h <= dim; ++h) { 6569566063dSJacob Faibussowitsch PetscCall(DMPlexGetSimplexOrBoxCells(dm, h, NULL, &pMax)); 6579566063dSJacob Faibussowitsch PetscCall(DMPlexGetHeightStratum(dm, h, &pStart, &pEnd)); 6589566063dSJacob Faibussowitsch for (p = pMax; p < pEnd; ++p) PetscCall(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; 6709566063dSJacob Faibussowitsch PetscCall(DMGetDimension(dm, &dim)); 6719566063dSJacob Faibussowitsch PetscCall(DMGetLabel(dm, "cohesive", &fault)); 6729566063dSJacob Faibussowitsch PetscCall(DMLabelView(fault, PETSC_VIEWER_STDOUT_WORLD)); 673b253942bSMatthew G. Knepley 6749566063dSJacob Faibussowitsch PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, dim, user->cellSimplex, "displacement_", PETSC_DETERMINE, &fe)); 6759566063dSJacob Faibussowitsch PetscCall(PetscFESetName(fe, "displacement")); 6769566063dSJacob Faibussowitsch PetscCall(DMAddField(dm, NULL, (PetscObject) fe)); 6779566063dSJacob Faibussowitsch PetscCall(PetscFEDestroy(&fe)); 678b253942bSMatthew G. Knepley 679b7519becSMatthew G. Knepley if (Ncf > 0) { 6809566063dSJacob Faibussowitsch PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim-1, dim, user->cellSimplex, "faulttraction_", PETSC_DETERMINE, &fe)); 6819566063dSJacob Faibussowitsch PetscCall(PetscFESetName(fe, "fault traction")); 6829566063dSJacob Faibussowitsch PetscCall(DMAddField(dm, fault, (PetscObject) fe)); 6839566063dSJacob Faibussowitsch PetscCall(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 6889566063dSJacob Faibussowitsch PetscCall(PetscSNPrintf(name, 256, "fault field %D", f)); 6899566063dSJacob Faibussowitsch PetscCall(PetscSNPrintf(opt, 256, "faultfield_%D_", f)); 6909566063dSJacob Faibussowitsch PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim-1, dim, user->cellSimplex, opt, PETSC_DETERMINE, &fe)); 6919566063dSJacob Faibussowitsch PetscCall(PetscFESetName(fe, name)); 6929566063dSJacob Faibussowitsch PetscCall(DMAddField(dm, fault, (PetscObject) fe)); 6939566063dSJacob Faibussowitsch PetscCall(PetscFEDestroy(&fe)); 694b7519becSMatthew G. Knepley } 695b253942bSMatthew G. Knepley 6969566063dSJacob Faibussowitsch PetscCall(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; 7089566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_rank(comm, &rank)); 7099566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_size(comm, &size)); 7109566063dSJacob Faibussowitsch PetscCall(DMCreate(comm, dm)); 7119566063dSJacob Faibussowitsch PetscCall(DMSetType(*dm, DMPLEX)); 7129566063dSJacob Faibussowitsch PetscCall(DMSetDimension(*dm, dim)); 713c4762a1bSJed Brown switch (dim) { 714c4762a1bSJed Brown case 2: 715c4762a1bSJed Brown if (cellSimplex) { 7169566063dSJacob Faibussowitsch PetscCall(CreateSimplex_2D(comm, user->testNum, dm)); 717c4762a1bSJed Brown } else { 7189566063dSJacob Faibussowitsch PetscCall(CreateQuad_2D(comm, user->testNum, dm)); 719c4762a1bSJed Brown } 720c4762a1bSJed Brown break; 721c4762a1bSJed Brown case 3: 722c4762a1bSJed Brown if (cellSimplex) { 7239566063dSJacob Faibussowitsch PetscCall(CreateSimplex_3D(comm, user, *dm)); 724c4762a1bSJed Brown } else { 7259566063dSJacob Faibussowitsch PetscCall(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 } 7319566063dSJacob Faibussowitsch PetscCall(PetscObjectSetOptionsPrefix((PetscObject) *dm, "orig_")); 7329566063dSJacob Faibussowitsch PetscCall(DMPlexDistributeSetDefault(*dm, PETSC_FALSE)); 7339566063dSJacob Faibussowitsch PetscCall(DMSetFromOptions(*dm)); 7349566063dSJacob Faibussowitsch PetscCall(DMGetLabel(*dm, "material", &matLabel)); 735dd96b1f9SToby Isaac if (matLabel) { 7369566063dSJacob Faibussowitsch PetscCall(DMPlexLabelComplete(*dm, matLabel)); 737dd96b1f9SToby Isaac } 7389566063dSJacob Faibussowitsch PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 7399566063dSJacob Faibussowitsch PetscCall(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 7449566063dSJacob Faibussowitsch PetscCall(DMGetLabel(*dm, "fault", &faultLabel)); 7459566063dSJacob Faibussowitsch PetscCall(DMGetLabel(*dm, "faultBd", &faultBdLabel)); 7469566063dSJacob Faibussowitsch PetscCall(DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, &hybridLabel, &splitLabel, &dmInterface, &dmHybrid)); 7479566063dSJacob Faibussowitsch PetscCall(DMLabelView(hybridLabel, PETSC_VIEWER_STDOUT_WORLD)); 7489566063dSJacob Faibussowitsch PetscCall(DMLabelDestroy(&hybridLabel)); 7499566063dSJacob Faibussowitsch PetscCall(DMLabelView(splitLabel, PETSC_VIEWER_STDOUT_WORLD)); 7509566063dSJacob Faibussowitsch PetscCall(DMLabelDestroy(&splitLabel)); 7519566063dSJacob Faibussowitsch PetscCall(DMViewFromOptions(dmInterface, NULL, "-dm_interface_view")); 7529566063dSJacob Faibussowitsch PetscCall(DMDestroy(&dmInterface)); 7539566063dSJacob Faibussowitsch PetscCall(DMDestroy(dm)); 754c4762a1bSJed Brown *dm = dmHybrid; 755c4762a1bSJed Brown } 7569566063dSJacob Faibussowitsch PetscCall(DMHasLabel(*dm, "fault2", &hasFault2)); 757c4762a1bSJed Brown if (hasFault2) { 758c4762a1bSJed Brown DM dmHybrid = NULL; 759c4762a1bSJed Brown DMLabel faultLabel, faultBdLabel, hybridLabel; 760c4762a1bSJed Brown 7619566063dSJacob Faibussowitsch PetscCall(PetscObjectSetOptionsPrefix((PetscObject) *dm, "faulted_")); 7629566063dSJacob Faibussowitsch PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view_pre")); 7639566063dSJacob Faibussowitsch PetscCall(DMPlexDistributeSetDefault(*dm, PETSC_FALSE)); 7649566063dSJacob Faibussowitsch PetscCall(DMSetFromOptions(*dm)); 7659566063dSJacob Faibussowitsch PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 7669566063dSJacob Faibussowitsch PetscCall(DMGetLabel(*dm, "fault2", &faultLabel)); 7679566063dSJacob Faibussowitsch PetscCall(DMGetLabel(*dm, "fault2Bd", &faultBdLabel)); 7689566063dSJacob Faibussowitsch PetscCall(DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, &hybridLabel, NULL, NULL, &dmHybrid)); 7699566063dSJacob Faibussowitsch PetscCall(DMLabelView(hybridLabel, PETSC_VIEWER_STDOUT_WORLD)); 7709566063dSJacob Faibussowitsch PetscCall(DMLabelDestroy(&hybridLabel)); 7719566063dSJacob Faibussowitsch PetscCall(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 7869566063dSJacob Faibussowitsch PetscCall(PetscMalloc2(2, &sizes, 3, &points)); 7879566063dSJacob Faibussowitsch PetscCall(PetscArraycpy(sizes, triSizes_p2, 2)); 7889566063dSJacob Faibussowitsch PetscCall(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 7989566063dSJacob Faibussowitsch PetscCall(PetscMalloc2(2, &sizes, 3, &points)); 7999566063dSJacob Faibussowitsch PetscCall(PetscArraycpy(sizes, quadSizes_p2, 2)); 8009566063dSJacob Faibussowitsch PetscCall(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 8059566063dSJacob Faibussowitsch PetscCall(PetscMalloc2(2, &sizes, 2, &points)); 8069566063dSJacob Faibussowitsch PetscCall(PetscArraycpy(sizes, quadSizes_p2, 2)); 8079566063dSJacob Faibussowitsch PetscCall(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 8179566063dSJacob Faibussowitsch PetscCall(PetscMalloc2(2, &sizes, 3, &points)); 8189566063dSJacob Faibussowitsch PetscCall(PetscArraycpy(sizes, tetSizes_p2, 2)); 8199566063dSJacob Faibussowitsch PetscCall(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 8299566063dSJacob Faibussowitsch PetscCall(PetscMalloc2(2, &sizes, 3, &points)); 8309566063dSJacob Faibussowitsch PetscCall(PetscArraycpy(sizes, hexSizes_p2, 2)); 8319566063dSJacob Faibussowitsch PetscCall(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 } 8379566063dSJacob Faibussowitsch PetscCall(DMPlexGetPartitioner(*dm, &part)); 8389566063dSJacob Faibussowitsch PetscCall(PetscPartitionerSetType(part, PETSCPARTITIONERSHELL)); 8399566063dSJacob Faibussowitsch PetscCall(PetscPartitionerShellSetPartition(part, size, sizes, points)); 8409566063dSJacob Faibussowitsch PetscCall(PetscFree2(sizes, points)); 841c4762a1bSJed Brown } 842c4762a1bSJed Brown { 843c4762a1bSJed Brown DM pdm = NULL; 844c4762a1bSJed Brown 845c4762a1bSJed Brown /* Distribute mesh over processes */ 8469566063dSJacob Faibussowitsch PetscCall(DMPlexDistribute(*dm, 0, NULL, &pdm)); 847c4762a1bSJed Brown if (pdm) { 8489566063dSJacob Faibussowitsch PetscCall(DMViewFromOptions(pdm, NULL, "-dm_view")); 8499566063dSJacob Faibussowitsch PetscCall(DMDestroy(dm)); 850c4762a1bSJed Brown *dm = pdm; 851c4762a1bSJed Brown } 852c4762a1bSJed Brown } 8539566063dSJacob Faibussowitsch PetscCall(DMHasLabel(*dm, "pfault", &hasParallelFault)); 854c4762a1bSJed Brown if (hasParallelFault) { 855c4762a1bSJed Brown DM dmHybrid = NULL; 856c4762a1bSJed Brown DMLabel faultLabel, faultBdLabel, hybridLabel; 857c4762a1bSJed Brown 8589566063dSJacob Faibussowitsch PetscCall(DMGetLabel(*dm, "pfault", &faultLabel)); 8599566063dSJacob Faibussowitsch PetscCall(DMGetLabel(*dm, "pfaultBd", &faultBdLabel)); 8609566063dSJacob Faibussowitsch PetscCall(DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, &hybridLabel, NULL, NULL, &dmHybrid)); 8619566063dSJacob Faibussowitsch PetscCall(DMLabelView(hybridLabel, PETSC_VIEWER_STDOUT_WORLD)); 8629566063dSJacob Faibussowitsch PetscCall(DMLabelDestroy(&hybridLabel)); 8639566063dSJacob Faibussowitsch PetscCall(DMDestroy(dm)); 864c4762a1bSJed Brown *dm = dmHybrid; 865c4762a1bSJed Brown } 8669566063dSJacob Faibussowitsch PetscCall(PetscObjectSetName((PetscObject) *dm, "Hybrid Mesh")); 8679566063dSJacob Faibussowitsch PetscCall(CreateFaultLabel(*dm)); 8689566063dSJacob Faibussowitsch PetscCall(CreateDiscretization(*dm, user)); 8699566063dSJacob Faibussowitsch PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view_pre")); 8709566063dSJacob Faibussowitsch PetscCall(DMPlexDistributeSetDefault(*dm, PETSC_FALSE)); 8719566063dSJacob Faibussowitsch PetscCall(DMSetFromOptions(*dm)); 8729566063dSJacob Faibussowitsch PetscCall(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; 8799566063dSJacob Faibussowitsch PetscCall(DMPlexCheckSymmetry(dm)); 8809566063dSJacob Faibussowitsch PetscCall(DMPlexCheckSkeleton(dm, 0)); 8819566063dSJacob Faibussowitsch PetscCall(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; 8909566063dSJacob Faibussowitsch PetscCall(DMGetSection(dm, &s)); 8919566063dSJacob Faibussowitsch PetscCall(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; 9889566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank)); 9899566063dSJacob Faibussowitsch PetscCall(DMGetDimension(dm, &dim)); 9909566063dSJacob Faibussowitsch PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, NULL, &cMax)); 9919566063dSJacob Faibussowitsch PetscCall(DMPlexGetHeightStratum(dm, 0, NULL, &cEnd)); 9929566063dSJacob Faibussowitsch PetscCall(ISCreateStride(PETSC_COMM_SELF, cEnd - cMax, cMax, 1, &cohesiveCells)); 9939566063dSJacob Faibussowitsch PetscCall(DMGetLabel(dm, "cohesive", &fault)); 9949566063dSJacob Faibussowitsch PetscCall(DMGetLocalVector(dm, &locX)); 9959566063dSJacob Faibussowitsch PetscCall(PetscObjectSetName((PetscObject) locX, "Local Solution")); 9969566063dSJacob Faibussowitsch PetscCall(DMGetLocalVector(dm, &locF)); 9979566063dSJacob Faibussowitsch PetscCall(PetscObjectSetName((PetscObject) locF, "Local Residual")); 9989566063dSJacob Faibussowitsch PetscCall(DMCreateMatrix(dm, &J)); 9999566063dSJacob Faibussowitsch PetscCall(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 */ 10029566063dSJacob Faibussowitsch PetscCall(DMGetLabel(dm, "material", &material)); 1003b253942bSMatthew G. Knepley id = 1; 1004b253942bSMatthew G. Knepley initialGuess[0] = r; 1005b253942bSMatthew G. Knepley initialGuess[1] = NULL; 10069566063dSJacob Faibussowitsch PetscCall(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; 10109566063dSJacob Faibussowitsch PetscCall(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; 10149566063dSJacob Faibussowitsch PetscCall(DMProjectFunctionLabelLocal(dm, 0.0, fault, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX)); 10159566063dSJacob Faibussowitsch PetscCall(VecViewFromOptions(locX, NULL, "-local_solution_view")); 1016b253942bSMatthew G. Knepley 10179566063dSJacob Faibussowitsch PetscCall(DMGetCellDS(dm, cMax, &probh)); 10189566063dSJacob Faibussowitsch PetscCall(PetscDSGetWeakForm(probh, &wf)); 10199566063dSJacob Faibussowitsch PetscCall(PetscDSGetNumFields(probh, &Nf)); 10209566063dSJacob Faibussowitsch PetscCall(PetscWeakFormSetIndexBdResidual(wf, material, 1, 0, 0, 0, f0_bd_u, 0, NULL)); 10219566063dSJacob Faibussowitsch PetscCall(PetscWeakFormSetIndexBdResidual(wf, material, 2, 0, 0, 0, f0_bd_u, 0, NULL)); 10229566063dSJacob Faibussowitsch PetscCall(PetscWeakFormSetIndexBdJacobian(wf, material, 1, 0, 1, 0, 0, g0_bd_ul, 0, NULL, 0, NULL, 0, NULL)); 10239566063dSJacob Faibussowitsch PetscCall(PetscWeakFormSetIndexBdJacobian(wf, material, 2, 0, 1, 0, 0, g0_bd_ul, 0, NULL, 0, NULL, 0, NULL)); 1024b7519becSMatthew G. Knepley if (Nf > 1) { 10259566063dSJacob Faibussowitsch PetscCall(PetscWeakFormSetIndexBdResidual(wf, fault, 1, 1, 0, 0, f0_bd_l, 0, NULL)); 10269566063dSJacob Faibussowitsch PetscCall(PetscWeakFormSetIndexBdJacobian(wf, fault, 1, 1, 0, 0, 0, g0_bd_lu, 0, NULL, 0, NULL, 0, NULL)); 1027b7519becSMatthew G. Knepley } 10289566063dSJacob Faibussowitsch if (!rank) PetscCall(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; 10429566063dSJacob Faibussowitsch PetscCall(VecSet(locF, 0.)); 10439566063dSJacob Faibussowitsch PetscCall(DMPlexComputeResidual_Hybrid_Internal(dm, keys, cohesiveCells, 0.0, locX, NULL, 0.0, locF, user)); 10449566063dSJacob Faibussowitsch PetscCall(VecViewFromOptions(locF, NULL, "-local_residual_view")); 10459566063dSJacob Faibussowitsch PetscCall(MatZeroEntries(J)); 10469566063dSJacob Faibussowitsch PetscCall(DMPlexComputeJacobian_Hybrid_Internal(dm, keys, cohesiveCells, 0.0, 0.0, locX, NULL, J, J, user)); 10479566063dSJacob Faibussowitsch PetscCall(MatViewFromOptions(J, NULL, "-local_jacobian_view")); 1048b253942bSMatthew G. Knepley 10499566063dSJacob Faibussowitsch PetscCall(DMRestoreLocalVector(dm, &locX)); 10509566063dSJacob Faibussowitsch PetscCall(DMRestoreLocalVector(dm, &locF)); 10519566063dSJacob Faibussowitsch PetscCall(MatDestroy(&J)); 10529566063dSJacob Faibussowitsch PetscCall(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 10619566063dSJacob Faibussowitsch PetscCall(PetscInitialize(&argc, &argv, NULL,help)); 10629566063dSJacob Faibussowitsch PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user)); 10639566063dSJacob Faibussowitsch PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm)); 10649566063dSJacob Faibussowitsch PetscCall(TestMesh(dm, &user)); 10659566063dSJacob Faibussowitsch PetscCall(TestDiscretization(dm, &user)); 10669566063dSJacob Faibussowitsch PetscCall(TestAssembly(dm, &user)); 10679566063dSJacob Faibussowitsch PetscCall(DMDestroy(&dm)); 10689566063dSJacob Faibussowitsch PetscCall(PetscFinalize()); 1069b122ec5aSJacob Faibussowitsch return 0; 1070c4762a1bSJed Brown } 1071c4762a1bSJed Brown 1072c4762a1bSJed Brown /*TEST 1073c4762a1bSJed Brown testset: 1074b253942bSMatthew G. Knepley args: -orig_dm_plex_check_all -dm_plex_check_all \ 1075b7519becSMatthew G. Knepley -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1 -local_section_view \ 1076b253942bSMatthew G. Knepley -local_solution_view -local_residual_view -local_jacobian_view 1077c4762a1bSJed Brown test: 1078c4762a1bSJed Brown suffix: tri_0 1079c4762a1bSJed Brown args: -dim 2 1080*1c6715b8SMatthew 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" 1081c4762a1bSJed Brown test: 1082c4762a1bSJed Brown suffix: tri_t1_0 1083c4762a1bSJed Brown args: -dim 2 -test_num 1 1084*1c6715b8SMatthew 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" 1085c4762a1bSJed Brown test: 1086c4762a1bSJed Brown suffix: tet_0 1087c4762a1bSJed Brown args: -dim 3 1088*1c6715b8SMatthew 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" 1089c4762a1bSJed Brown test: 1090b253942bSMatthew G. Knepley suffix: tet_t1_0 1091b253942bSMatthew G. Knepley args: -dim 3 -test_num 1 1092*1c6715b8SMatthew 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" 1093b253942bSMatthew G. Knepley 1094b253942bSMatthew G. Knepley testset: 1095b253942bSMatthew G. Knepley args: -orig_dm_plex_check_all -dm_plex_check_all \ 1096b7519becSMatthew G. Knepley -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1 1097b253942bSMatthew G. Knepley test: 1098c4762a1bSJed Brown suffix: tet_1 1099c4762a1bSJed Brown nsize: 2 1100c4762a1bSJed Brown args: -dim 3 1101*1c6715b8SMatthew 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" 1102c4762a1bSJed Brown test: 1103b253942bSMatthew G. Knepley suffix: tri_1 1104b253942bSMatthew G. Knepley nsize: 2 1105b253942bSMatthew G. Knepley args: -dim 2 1106*1c6715b8SMatthew 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" 1107c4762a1bSJed Brown 1108c4762a1bSJed Brown testset: 1109ecfb78b5SMatthew G. Knepley args: -orig_dm_plex_check_all -dm_plex_check_all \ 1110b7519becSMatthew G. Knepley -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1 1111c4762a1bSJed Brown # 2D Quads 1112c4762a1bSJed Brown test: 1113c4762a1bSJed Brown suffix: quad_0 1114c4762a1bSJed Brown args: -dim 2 -cell_simplex 0 1115*1c6715b8SMatthew 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" 1116c4762a1bSJed Brown test: 1117c4762a1bSJed Brown suffix: quad_1 1118c4762a1bSJed Brown nsize: 2 1119c4762a1bSJed Brown args: -dim 2 -cell_simplex 0 1120*1c6715b8SMatthew 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" 1121c4762a1bSJed Brown test: 1122c4762a1bSJed Brown suffix: quad_t1_0 112354fcfd0cSMatthew G. Knepley args: -dim 2 -cell_simplex 0 -test_num 1 -faulted_dm_plex_check_all 1124*1c6715b8SMatthew 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" 1125c4762a1bSJed Brown # 3D Hex 1126c4762a1bSJed Brown test: 1127c4762a1bSJed Brown suffix: hex_0 1128c4762a1bSJed Brown args: -dim 3 -cell_simplex 0 1129*1c6715b8SMatthew 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" 1130c4762a1bSJed Brown test: 1131c4762a1bSJed Brown suffix: hex_1 1132c4762a1bSJed Brown nsize: 2 1133c4762a1bSJed Brown args: -dim 3 -cell_simplex 0 1134*1c6715b8SMatthew 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" 1135c4762a1bSJed Brown test: 1136c4762a1bSJed Brown suffix: hex_t1_0 1137c4762a1bSJed Brown args: -dim 3 -cell_simplex 0 -test_num 1 1138*1c6715b8SMatthew 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" 1139c4762a1bSJed Brown test: 1140c4762a1bSJed Brown suffix: hex_t2_0 1141c4762a1bSJed Brown args: -dim 3 -cell_simplex 0 -test_num 2 1142*1c6715b8SMatthew 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" 1143c4762a1bSJed Brown 1144c4762a1bSJed Brown TEST*/ 1145