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 */ 309c4762a1bSJed Brown } AppCtx; 310c4762a1bSJed Brown 311b253942bSMatthew G. Knepley static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 312c4762a1bSJed Brown { 313c4762a1bSJed Brown PetscErrorCode ierr; 314c4762a1bSJed Brown 315c4762a1bSJed Brown PetscFunctionBegin; 316c4762a1bSJed Brown options->debug = 0; 317c4762a1bSJed Brown options->dim = 2; 318c4762a1bSJed Brown options->cellSimplex = PETSC_TRUE; 319c4762a1bSJed Brown options->testPartition = PETSC_TRUE; 320c4762a1bSJed Brown options->testNum = 0; 321c4762a1bSJed Brown 322c4762a1bSJed Brown ierr = PetscOptionsBegin(comm, "", "Meshing Problem Options", "DMPLEX");CHKERRQ(ierr); 323c4762a1bSJed Brown ierr = PetscOptionsBoundedInt("-debug", "The debugging level", "ex5.c", options->debug, &options->debug, NULL,0);CHKERRQ(ierr); 324c4762a1bSJed Brown ierr = PetscOptionsRangeInt("-dim", "The topological mesh dimension", "ex5.c", options->dim, &options->dim, NULL,1,3);CHKERRQ(ierr); 325c4762a1bSJed Brown ierr = PetscOptionsBool("-cell_simplex", "Use simplices if true, otherwise hexes", "ex5.c", options->cellSimplex, &options->cellSimplex, NULL);CHKERRQ(ierr); 326c4762a1bSJed Brown ierr = PetscOptionsBool("-test_partition", "Use a fixed partition for testing", "ex5.c", options->testPartition, &options->testPartition, NULL);CHKERRQ(ierr); 327c4762a1bSJed Brown ierr = PetscOptionsBoundedInt("-test_num", "The particular mesh to test", "ex5.c", options->testNum, &options->testNum, NULL,0);CHKERRQ(ierr); 3281e1ea65dSPierre Jolivet ierr = PetscOptionsEnd();CHKERRQ(ierr); 329c4762a1bSJed Brown PetscFunctionReturn(0); 330c4762a1bSJed Brown } 331c4762a1bSJed Brown 332b253942bSMatthew G. Knepley static PetscErrorCode CreateSimplex_2D(MPI_Comm comm, PetscInt testNum, DM *dm) 333c4762a1bSJed Brown { 334c4762a1bSJed Brown DM idm; 335c4762a1bSJed Brown PetscInt p; 336c4762a1bSJed Brown PetscMPIInt rank; 337c4762a1bSJed Brown PetscErrorCode ierr; 338c4762a1bSJed Brown 339c4762a1bSJed Brown PetscFunctionBegin; 340ffc4695bSBarry Smith ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 341c4762a1bSJed Brown if (!rank) { 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 353c4762a1bSJed Brown ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 354c4762a1bSJed Brown for (p = 0; p < 4; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 355c4762a1bSJed Brown for (p = 0; p < 2; ++p) {ierr = DMSetLabelValue(*dm, "fault", faultPoints[p], 1);CHKERRQ(ierr);} 356b253942bSMatthew G. Knepley ierr = DMSetLabelValue(*dm, "material", 0, 1);CHKERRQ(ierr); 357b253942bSMatthew G. Knepley ierr = DMSetLabelValue(*dm, "material", 1, 2);CHKERRQ(ierr); 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 370c4762a1bSJed Brown ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 371c4762a1bSJed Brown for (p = 0; p < 3; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 372c4762a1bSJed Brown for (p = 0; p < 3; ++p) {ierr = DMSetLabelValue(*dm, "fault", faultPoints[p], 1);CHKERRQ(ierr);} 373b253942bSMatthew G. Knepley ierr = DMSetLabelValue(*dm, "material", 0, 1);CHKERRQ(ierr);ierr = DMSetLabelValue(*dm, "material", 3, 1);CHKERRQ(ierr); 374b253942bSMatthew G. Knepley ierr = DMSetLabelValue(*dm, "material", 1, 2);CHKERRQ(ierr);ierr = DMSetLabelValue(*dm, "material", 2, 2);CHKERRQ(ierr); 375c4762a1bSJed Brown } 376c4762a1bSJed Brown break; 377c4762a1bSJed Brown default: 378c4762a1bSJed Brown SETERRQ1(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 383c4762a1bSJed Brown ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 384c4762a1bSJed Brown ierr = DMCreateLabel(*dm, "fault");CHKERRQ(ierr); 385c4762a1bSJed Brown } 386c4762a1bSJed Brown ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 387c4762a1bSJed Brown ierr = DMViewFromOptions(idm, NULL, "-in_dm_view");CHKERRQ(ierr); 388c4762a1bSJed Brown ierr = DMDestroy(dm);CHKERRQ(ierr); 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 PetscErrorCode ierr; 398c4762a1bSJed Brown 399c4762a1bSJed Brown PetscFunctionBegin; 400ffc4695bSBarry Smith ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 401c4762a1bSJed Brown if (!rank) { 402c4762a1bSJed Brown switch (testNum) { 403c4762a1bSJed Brown case 0: 404c4762a1bSJed Brown { 405c4762a1bSJed Brown PetscInt numPoints[4] = {5, 7, 9, 2}; 406c4762a1bSJed 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}; 407c4762a1bSJed 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}; 408b5a892a1SMatthew 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}; 409c4762a1bSJed 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}; 410c4762a1bSJed Brown PetscInt markerPoints[20] = {2, 1, 3, 1, 4, 1, 5, 1, 14, 1, 15, 1, 16, 1, 17, 1, 18, 1, 19, 1}; 411c4762a1bSJed Brown PetscInt faultPoints[3] = {3, 4, 5}; 412c4762a1bSJed Brown 413c4762a1bSJed Brown ierr = DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 414c4762a1bSJed Brown for (p = 0; p < 10; ++p) { 415c4762a1bSJed Brown ierr = DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr); 416c4762a1bSJed Brown } 417c4762a1bSJed Brown for (p = 0; p < 3; ++p) { 418c4762a1bSJed Brown ierr = DMSetLabelValue(dm, "fault", faultPoints[p], 1);CHKERRQ(ierr); 419c4762a1bSJed Brown } 420c4762a1bSJed Brown } 421c4762a1bSJed Brown break; 422c4762a1bSJed Brown case 1: 423c4762a1bSJed Brown { 424c4762a1bSJed Brown PetscInt numPoints[4] = {6, 13, 12, 4}; 425c4762a1bSJed 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}; 426c4762a1bSJed 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}; 427b5a892a1SMatthew 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}; 428c4762a1bSJed 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}; 429c4762a1bSJed Brown PetscInt markerPoints[14] = {5, 1, 6, 1, 7, 1, 10, 1, 22, 1, 23, 1, 24, 1}; 430c4762a1bSJed Brown PetscInt faultPoints[4] = {5, 6, 7, 8}; 431c4762a1bSJed Brown 432c4762a1bSJed Brown ierr = DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 433c4762a1bSJed Brown for (p = 0; p < 7; ++p) { 434c4762a1bSJed Brown ierr = DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr); 435c4762a1bSJed Brown } 436c4762a1bSJed Brown for (p = 0; p < 4; ++p) { 437c4762a1bSJed Brown ierr = DMSetLabelValue(dm, "fault", faultPoints[p], 1);CHKERRQ(ierr); 438c4762a1bSJed Brown } 439c4762a1bSJed Brown } 440c4762a1bSJed Brown break; 441c4762a1bSJed Brown default: 442c4762a1bSJed Brown SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %d", testNum); 443c4762a1bSJed Brown } 444c4762a1bSJed Brown } else { 445c4762a1bSJed Brown PetscInt numPoints[4] = {0, 0, 0, 0}; 446c4762a1bSJed Brown 447c4762a1bSJed Brown ierr = DMPlexCreateFromDAG(dm, depth, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 448c4762a1bSJed Brown ierr = DMCreateLabel(dm, "fault");CHKERRQ(ierr); 449c4762a1bSJed Brown } 450c4762a1bSJed Brown PetscFunctionReturn(0); 451c4762a1bSJed Brown } 452c4762a1bSJed Brown 453b253942bSMatthew G. Knepley static PetscErrorCode CreateQuad_2D(MPI_Comm comm, PetscInt testNum, DM *dm) 454c4762a1bSJed Brown { 455c4762a1bSJed Brown DM idm; 456c4762a1bSJed Brown PetscInt p; 457c4762a1bSJed Brown PetscMPIInt rank; 458c4762a1bSJed Brown PetscErrorCode ierr; 459c4762a1bSJed Brown 460c4762a1bSJed Brown PetscFunctionBegin; 461ffc4695bSBarry Smith ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 462c4762a1bSJed Brown if (!rank) { 463c4762a1bSJed Brown switch (testNum) { 464c4762a1bSJed Brown case 0: 465c4762a1bSJed Brown case 2: 466c4762a1bSJed Brown { 467c4762a1bSJed Brown PetscInt numPoints[2] = {6, 2}; 468c4762a1bSJed Brown PetscInt coneSize[8] = {4, 4, 0, 0, 0, 0, 0, 0}; 469c4762a1bSJed Brown PetscInt cones[8] = {2, 3, 4, 5, 3, 6, 7, 4}; 470c4762a1bSJed Brown PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 471c4762a1bSJed 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}; 472c4762a1bSJed Brown PetscInt markerPoints[12] = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1, 7, 1}; 473c4762a1bSJed Brown PetscInt faultPoints[2] = {3, 4}; 474c4762a1bSJed Brown 475c4762a1bSJed Brown ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 476c4762a1bSJed Brown for (p = 0; p < 6; ++p) {ierr = DMSetLabelValue(*dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 477c4762a1bSJed Brown if (testNum == 0) for (p = 0; p < 2; ++p) {ierr = DMSetLabelValue(*dm, "fault", faultPoints[p], 1);CHKERRQ(ierr);} 478c4762a1bSJed Brown if (testNum == 2) for (p = 0; p < 2; ++p) {ierr = DMSetLabelValue(*dm, "pfault", faultPoints[p], 1);CHKERRQ(ierr);} 479c4762a1bSJed Brown } 480c4762a1bSJed Brown break; 481c4762a1bSJed Brown case 1: 482c4762a1bSJed Brown { 483c4762a1bSJed Brown PetscInt numPoints[2] = {16, 9}; 484c4762a1bSJed 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}; 485c4762a1bSJed Brown PetscInt cones[36] = {9, 13, 14, 10, 486c4762a1bSJed Brown 10, 14, 15, 11, 487c4762a1bSJed Brown 11, 15, 16, 12, 488c4762a1bSJed Brown 13, 17, 18, 14, 489c4762a1bSJed Brown 14, 18, 19, 15, 490c4762a1bSJed Brown 15, 19, 20, 16, 491c4762a1bSJed Brown 17, 21, 22, 18, 492c4762a1bSJed Brown 18, 22, 23, 19, 493c4762a1bSJed Brown 19, 23, 24, 20}; 494c4762a1bSJed 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}; 495c4762a1bSJed 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, 496c4762a1bSJed 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}; 497c4762a1bSJed Brown PetscInt faultPoints[3] = {11, 15, 19}; 498c4762a1bSJed Brown PetscInt fault2Points[2] = {17, 18}; 499c4762a1bSJed Brown 500c4762a1bSJed Brown ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 501c4762a1bSJed Brown for (p = 0; p < 3; ++p) {ierr = DMSetLabelValue(*dm, "fault", faultPoints[p], 1);CHKERRQ(ierr);} 502c4762a1bSJed Brown for (p = 0; p < 2; ++p) {ierr = DMSetLabelValue(*dm, "fault2", fault2Points[p], 1);CHKERRQ(ierr);} 503c4762a1bSJed Brown } 504c4762a1bSJed Brown break; 505c4762a1bSJed Brown default: 506c4762a1bSJed Brown SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %d", testNum); 507c4762a1bSJed Brown } 508c4762a1bSJed Brown } else { 509c4762a1bSJed Brown PetscInt numPoints[3] = {0, 0, 0}; 510c4762a1bSJed Brown 511c4762a1bSJed Brown ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 512c4762a1bSJed Brown if (testNum == 2) {ierr = DMCreateLabel(*dm, "pfault");CHKERRQ(ierr);} 513c4762a1bSJed Brown else {ierr = DMCreateLabel(*dm, "fault");CHKERRQ(ierr);} 514c4762a1bSJed Brown } 515c4762a1bSJed Brown ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 516c4762a1bSJed Brown ierr = DMViewFromOptions(idm, NULL, "-in_dm_view");CHKERRQ(ierr); 517c4762a1bSJed Brown ierr = DMDestroy(dm);CHKERRQ(ierr); 518c4762a1bSJed Brown *dm = idm; 519c4762a1bSJed Brown PetscFunctionReturn(0); 520c4762a1bSJed Brown } 521c4762a1bSJed Brown 522b253942bSMatthew G. Knepley static PetscErrorCode CreateHex_3D(MPI_Comm comm, PetscInt testNum, DM *dm) 523c4762a1bSJed Brown { 524c4762a1bSJed Brown DM idm; 525c4762a1bSJed Brown PetscInt depth = 3, p; 526c4762a1bSJed Brown PetscMPIInt rank; 527c4762a1bSJed Brown PetscErrorCode ierr; 528c4762a1bSJed Brown 529c4762a1bSJed Brown PetscFunctionBegin; 530ffc4695bSBarry Smith ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 531c4762a1bSJed Brown if (!rank) { 532c4762a1bSJed Brown switch (testNum) { 533c4762a1bSJed Brown case 0: 534c4762a1bSJed Brown { 535c4762a1bSJed Brown PetscInt numPoints[2] = {12, 2}; 536c4762a1bSJed Brown PetscInt coneSize[14] = {8,8, 0,0,0,0,0,0,0,0,0,0,0,0}; 537c4762a1bSJed Brown PetscInt cones[16] = {2,5,4,3,6,7,8,9, 3,4,11,10,7,12,13,8}; 538c4762a1bSJed Brown PetscInt coneOrientations[16] = {0,0,0,0,0,0,0,0, 0,0, 0,0 ,0, 0, 0,0}; 539c4762a1bSJed 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, 540c4762a1bSJed 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, 541c4762a1bSJed 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}; 542c4762a1bSJed Brown PetscInt markerPoints[52] = {2,1,3,1,4,1,5,1,6,1,7,1,8,1,9,1}; 543c4762a1bSJed Brown PetscInt faultPoints[4] = {3, 4, 7, 8}; 544c4762a1bSJed Brown 545c4762a1bSJed Brown ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 546c4762a1bSJed Brown ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 547c4762a1bSJed Brown for (p = 0; p < 8; ++p) {ierr = DMSetLabelValue(idm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 548c4762a1bSJed Brown for (p = 0; p < 4; ++p) {ierr = DMSetLabelValue(idm, "fault", faultPoints[p], 1);CHKERRQ(ierr);} 549c4762a1bSJed Brown } 550c4762a1bSJed Brown break; 551c4762a1bSJed Brown case 1: 552c4762a1bSJed Brown { 553c4762a1bSJed Brown /* Cell Adjacency Graph: 554c4762a1bSJed Brown 0 -- { 8, 13, 21, 24} --> 1 555c4762a1bSJed Brown 0 -- {20, 21, 23, 24} --> 5 F 556c4762a1bSJed Brown 1 -- {10, 15, 21, 24} --> 2 557c4762a1bSJed Brown 1 -- {13, 14, 15, 24} --> 6 558c4762a1bSJed Brown 2 -- {21, 22, 24, 25} --> 4 F 559c4762a1bSJed Brown 3 -- {21, 24, 30, 35} --> 4 560c4762a1bSJed Brown 3 -- {21, 24, 28, 33} --> 5 561c4762a1bSJed Brown */ 562c4762a1bSJed Brown PetscInt numPoints[2] = {30, 7}; 563c4762a1bSJed 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}; 564c4762a1bSJed Brown PetscInt cones[56] = { 8, 21, 20, 7, 13, 12, 23, 24, 565c4762a1bSJed Brown 14, 15, 10, 9, 13, 8, 21, 24, 566c4762a1bSJed Brown 15, 16, 11, 10, 24, 21, 22, 25, 567c4762a1bSJed Brown 30, 29, 28, 21, 35, 24, 33, 34, 568c4762a1bSJed Brown 24, 21, 30, 35, 25, 36, 31, 22, 569c4762a1bSJed Brown 27, 20, 21, 28, 32, 33, 24, 23, 570c4762a1bSJed Brown 15, 24, 13, 14, 19, 18, 17, 26}; 571c4762a1bSJed 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}; 572c4762a1bSJed 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, 573c4762a1bSJed 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, 574c4762a1bSJed 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, 575c4762a1bSJed 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, 576c4762a1bSJed 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}; 577c4762a1bSJed Brown PetscInt faultPoints[6] = {20, 21, 22, 23, 24, 25}; 578c4762a1bSJed Brown 579c4762a1bSJed Brown ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 580c4762a1bSJed Brown ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 581c4762a1bSJed Brown for (p = 0; p < 6; ++p) {ierr = DMSetLabelValue(idm, "fault", faultPoints[p], 1);CHKERRQ(ierr);} 582c4762a1bSJed Brown } 583c4762a1bSJed Brown break; 584c4762a1bSJed Brown case 2: 585c4762a1bSJed Brown { 586c4762a1bSJed Brown /* Buried fault edge */ 587c4762a1bSJed Brown PetscInt numPoints[2] = {18, 4}; 588c4762a1bSJed 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}; 589c4762a1bSJed Brown PetscInt cones[32] = { 4, 5, 8, 7, 13, 16, 17, 14, 590c4762a1bSJed Brown 5, 6, 9, 8, 14, 17, 18, 15, 591c4762a1bSJed Brown 7, 8, 11, 10, 16, 19, 20, 17, 592c4762a1bSJed Brown 8, 9, 12, 11, 17, 20, 21, 18}; 593c4762a1bSJed 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}; 594c4762a1bSJed 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, 595c4762a1bSJed 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, 596c4762a1bSJed 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}; 597c4762a1bSJed Brown PetscInt faultPoints[4] = {7, 8, 16, 17}; 598c4762a1bSJed Brown 599c4762a1bSJed Brown ierr = DMPlexCreateFromDAG(*dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 600c4762a1bSJed Brown ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 601c4762a1bSJed Brown for (p = 0; p < 4; ++p) {ierr = DMSetLabelValue(idm, "fault", faultPoints[p], 1);CHKERRQ(ierr);} 602c4762a1bSJed Brown } 603c4762a1bSJed Brown break; 604c4762a1bSJed Brown default: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No test mesh %d", testNum); 605c4762a1bSJed Brown } 606c4762a1bSJed Brown } else { 607c4762a1bSJed Brown PetscInt numPoints[4] = {0, 0, 0, 0}; 608c4762a1bSJed Brown 609c4762a1bSJed Brown ierr = DMPlexCreateFromDAG(*dm, depth, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 610c4762a1bSJed Brown ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 611c4762a1bSJed Brown ierr = DMCreateLabel(idm, "fault");CHKERRQ(ierr); 612c4762a1bSJed Brown } 613c4762a1bSJed Brown ierr = DMViewFromOptions(idm, NULL, "-in_dm_view");CHKERRQ(ierr); 614c4762a1bSJed Brown ierr = DMDestroy(dm);CHKERRQ(ierr); 615c4762a1bSJed Brown *dm = idm; 616c4762a1bSJed Brown PetscFunctionReturn(0); 617c4762a1bSJed Brown } 618c4762a1bSJed Brown 619b253942bSMatthew G. Knepley static PetscErrorCode CreateFaultLabel(DM dm) 620b253942bSMatthew G. Knepley { 621b253942bSMatthew G. Knepley DMLabel label; 622b253942bSMatthew G. Knepley PetscInt dim, h, pStart, pEnd, pMax, p; 623b253942bSMatthew G. Knepley PetscErrorCode ierr; 624b253942bSMatthew G. Knepley 625b253942bSMatthew G. Knepley PetscFunctionBegin; 626b253942bSMatthew G. Knepley ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 627b253942bSMatthew G. Knepley ierr = DMCreateLabel(dm, "cohesive");CHKERRQ(ierr); 628b253942bSMatthew G. Knepley ierr = DMGetLabel(dm, "cohesive", &label);CHKERRQ(ierr); 629b253942bSMatthew G. Knepley for (h = 0; h <= dim; ++h) { 630b253942bSMatthew G. Knepley ierr = DMPlexGetSimplexOrBoxCells(dm, h, NULL, &pMax);CHKERRQ(ierr); 631b253942bSMatthew G. Knepley ierr = DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);CHKERRQ(ierr); 632b253942bSMatthew G. Knepley for (p = pMax; p < pEnd; ++p) {ierr = DMLabelSetValue(label, p, 1);CHKERRQ(ierr);} 633b253942bSMatthew G. Knepley } 634b253942bSMatthew G. Knepley PetscFunctionReturn(0); 635b253942bSMatthew G. Knepley } 636b253942bSMatthew G. Knepley 637b253942bSMatthew G. Knepley static PetscErrorCode CreateDiscretization(DM dm, AppCtx *user) 638b253942bSMatthew G. Knepley { 639b253942bSMatthew G. Knepley PetscFE fe; 640b253942bSMatthew G. Knepley DMLabel fault; 641b253942bSMatthew G. Knepley PetscInt dim; 642b253942bSMatthew G. Knepley PetscErrorCode ierr; 643b253942bSMatthew G. Knepley 644b253942bSMatthew G. Knepley PetscFunctionBegin; 645b253942bSMatthew G. Knepley ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 646b253942bSMatthew G. Knepley ierr = DMGetLabel(dm, "cohesive", &fault);CHKERRQ(ierr); 647b253942bSMatthew G. Knepley ierr = DMLabelView(fault, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 648b253942bSMatthew G. Knepley 649b253942bSMatthew G. Knepley ierr = PetscFECreateDefault(PETSC_COMM_SELF, dim, dim, user->cellSimplex, "displacement_", PETSC_DETERMINE, &fe);CHKERRQ(ierr); 650b253942bSMatthew G. Knepley ierr = PetscFESetName(fe, "displacement");CHKERRQ(ierr); 6511e1ea65dSPierre Jolivet ierr = DMAddField(dm, NULL, (PetscObject) fe);CHKERRQ(ierr); 652b253942bSMatthew G. Knepley ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 653b253942bSMatthew G. Knepley 654b253942bSMatthew G. Knepley ierr = PetscFECreateDefault(PETSC_COMM_SELF, dim-1, dim, user->cellSimplex, "faulttraction_", PETSC_DETERMINE, &fe);CHKERRQ(ierr); 655b253942bSMatthew G. Knepley ierr = PetscFESetName(fe, "fault traction");CHKERRQ(ierr); 6561e1ea65dSPierre Jolivet ierr = DMAddField(dm, fault, (PetscObject) fe);CHKERRQ(ierr); 657b253942bSMatthew G. Knepley ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 658b253942bSMatthew G. Knepley 659b253942bSMatthew G. Knepley ierr = DMCreateDS(dm);CHKERRQ(ierr); 660b253942bSMatthew G. Knepley PetscFunctionReturn(0); 661b253942bSMatthew G. Knepley } 662b253942bSMatthew G. Knepley 663b253942bSMatthew G. Knepley static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 664c4762a1bSJed Brown { 665c4762a1bSJed Brown PetscInt dim = user->dim; 666c4762a1bSJed Brown PetscBool cellSimplex = user->cellSimplex, hasFault, hasFault2, hasParallelFault; 667c4762a1bSJed Brown PetscMPIInt rank, size; 668c4762a1bSJed Brown PetscErrorCode ierr; 669c4762a1bSJed Brown 670c4762a1bSJed Brown PetscFunctionBegin; 671ffc4695bSBarry Smith ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 672ffc4695bSBarry Smith ierr = MPI_Comm_size(comm, &size);CHKERRMPI(ierr); 673c4762a1bSJed Brown ierr = DMCreate(comm, dm);CHKERRQ(ierr); 674c4762a1bSJed Brown ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 675c4762a1bSJed Brown ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 676c4762a1bSJed Brown switch (dim) { 677c4762a1bSJed Brown case 2: 678c4762a1bSJed Brown if (cellSimplex) { 679c4762a1bSJed Brown ierr = CreateSimplex_2D(comm, user->testNum, dm);CHKERRQ(ierr); 680c4762a1bSJed Brown } else { 681c4762a1bSJed Brown ierr = CreateQuad_2D(comm, user->testNum, dm);CHKERRQ(ierr); 682c4762a1bSJed Brown } 683c4762a1bSJed Brown break; 684c4762a1bSJed Brown case 3: 685c4762a1bSJed Brown if (cellSimplex) { 686c4762a1bSJed Brown ierr = CreateSimplex_3D(comm, user, *dm);CHKERRQ(ierr); 687c4762a1bSJed Brown } else { 688c4762a1bSJed Brown ierr = CreateHex_3D(comm, user->testNum, dm);CHKERRQ(ierr); 689c4762a1bSJed Brown } 690c4762a1bSJed Brown break; 691c4762a1bSJed Brown default: 692c4762a1bSJed Brown SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make hybrid meshes for dimension %d", dim); 693c4762a1bSJed Brown } 694c4762a1bSJed Brown ierr = PetscObjectSetOptionsPrefix((PetscObject) *dm, "orig_");CHKERRQ(ierr); 695c4762a1bSJed Brown ierr = DMSetFromOptions(*dm);CHKERRQ(ierr); 696*012bc364SMatthew G. Knepley ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 697c4762a1bSJed Brown ierr = DMHasLabel(*dm, "fault", &hasFault);CHKERRQ(ierr); 698c4762a1bSJed Brown if (hasFault) { 699b253942bSMatthew G. Knepley DM dmHybrid = NULL, dmInterface = NULL; 700b253942bSMatthew G. Knepley DMLabel faultLabel, faultBdLabel, hybridLabel, splitLabel; 701c4762a1bSJed Brown 702c4762a1bSJed Brown ierr = DMGetLabel(*dm, "fault", &faultLabel);CHKERRQ(ierr); 703c4762a1bSJed Brown ierr = DMGetLabel(*dm, "faultBd", &faultBdLabel);CHKERRQ(ierr); 704b253942bSMatthew G. Knepley ierr = DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, &hybridLabel, &splitLabel, &dmInterface, &dmHybrid);CHKERRQ(ierr); 705c4762a1bSJed Brown ierr = DMLabelView(hybridLabel, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 706c4762a1bSJed Brown ierr = DMLabelDestroy(&hybridLabel);CHKERRQ(ierr); 707b253942bSMatthew G. Knepley ierr = DMLabelView(splitLabel, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 708b253942bSMatthew G. Knepley ierr = DMLabelDestroy(&splitLabel);CHKERRQ(ierr); 709b253942bSMatthew G. Knepley ierr = DMViewFromOptions(dmInterface, NULL, "-dm_interface_view");CHKERRQ(ierr); 710b253942bSMatthew G. Knepley ierr = DMDestroy(&dmInterface);CHKERRQ(ierr); 711c4762a1bSJed Brown ierr = DMDestroy(dm);CHKERRQ(ierr); 712c4762a1bSJed Brown *dm = dmHybrid; 713c4762a1bSJed Brown } 714c4762a1bSJed Brown ierr = DMHasLabel(*dm, "fault2", &hasFault2);CHKERRQ(ierr); 715c4762a1bSJed Brown if (hasFault2) { 716c4762a1bSJed Brown DM dmHybrid = NULL; 717c4762a1bSJed Brown DMLabel faultLabel, faultBdLabel, hybridLabel; 718c4762a1bSJed Brown 719c4762a1bSJed Brown ierr = PetscObjectSetOptionsPrefix((PetscObject) *dm, "faulted_");CHKERRQ(ierr); 720*012bc364SMatthew G. Knepley ierr = DMViewFromOptions(*dm, NULL, "-dm_view_pre");CHKERRQ(ierr); 721c4762a1bSJed Brown ierr = DMSetFromOptions(*dm);CHKERRQ(ierr); 722*012bc364SMatthew G. Knepley ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 723c4762a1bSJed Brown ierr = DMGetLabel(*dm, "fault2", &faultLabel);CHKERRQ(ierr); 724c4762a1bSJed Brown ierr = DMGetLabel(*dm, "fault2Bd", &faultBdLabel);CHKERRQ(ierr); 725c4762a1bSJed Brown ierr = DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, &hybridLabel, NULL, NULL, &dmHybrid);CHKERRQ(ierr); 726c4762a1bSJed Brown ierr = DMLabelView(hybridLabel, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 727c4762a1bSJed Brown ierr = DMLabelDestroy(&hybridLabel);CHKERRQ(ierr); 728c4762a1bSJed Brown ierr = DMDestroy(dm);CHKERRQ(ierr); 729c4762a1bSJed Brown *dm = dmHybrid; 730c4762a1bSJed Brown } 731c4762a1bSJed Brown if (user->testPartition && size > 1) { 732c4762a1bSJed Brown PetscPartitioner part; 733c4762a1bSJed Brown PetscInt *sizes = NULL; 734c4762a1bSJed Brown PetscInt *points = NULL; 735c4762a1bSJed Brown 736c4762a1bSJed Brown if (!rank) { 737c4762a1bSJed Brown if (dim == 2 && cellSimplex && size == 2) { 738c4762a1bSJed Brown switch (user->testNum) { 739c4762a1bSJed Brown case 0: { 740c4762a1bSJed Brown PetscInt triSizes_p2[2] = {1, 2}; 741c4762a1bSJed Brown PetscInt triPoints_p2[3] = {0, 1, 2}; 742c4762a1bSJed Brown 743c4762a1bSJed Brown ierr = PetscMalloc2(2, &sizes, 3, &points);CHKERRQ(ierr); 744c4762a1bSJed Brown ierr = PetscArraycpy(sizes, triSizes_p2, 2);CHKERRQ(ierr); 745c4762a1bSJed Brown ierr = PetscArraycpy(points, triPoints_p2, 3);CHKERRQ(ierr);break;} 746c4762a1bSJed Brown default: 747c4762a1bSJed Brown SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %d for triangular mesh on 2 procs", user->testNum); 748c4762a1bSJed Brown } 749c4762a1bSJed Brown } else if (dim == 2 && !cellSimplex && size == 2) { 750c4762a1bSJed Brown switch (user->testNum) { 751c4762a1bSJed Brown case 0: { 752c4762a1bSJed Brown PetscInt quadSizes_p2[2] = {1, 2}; 753c4762a1bSJed Brown PetscInt quadPoints_p2[3] = {0, 1, 2}; 754c4762a1bSJed Brown 755c4762a1bSJed Brown ierr = PetscMalloc2(2, &sizes, 3, &points);CHKERRQ(ierr); 756c4762a1bSJed Brown ierr = PetscArraycpy(sizes, quadSizes_p2, 2);CHKERRQ(ierr); 757c4762a1bSJed Brown ierr = PetscArraycpy(points, quadPoints_p2, 3);CHKERRQ(ierr);break;} 758c4762a1bSJed Brown case 2: { 759c4762a1bSJed Brown PetscInt quadSizes_p2[2] = {1, 1}; 760c4762a1bSJed Brown PetscInt quadPoints_p2[2] = {0, 1}; 761c4762a1bSJed Brown 762c4762a1bSJed Brown ierr = PetscMalloc2(2, &sizes, 2, &points);CHKERRQ(ierr); 763c4762a1bSJed Brown ierr = PetscArraycpy(sizes, quadSizes_p2, 2);CHKERRQ(ierr); 764c4762a1bSJed Brown ierr = PetscArraycpy(points, quadPoints_p2, 2);CHKERRQ(ierr);break;} 765c4762a1bSJed Brown default: 766c4762a1bSJed Brown SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %d for quadrilateral mesh on 2 procs", user->testNum); 767c4762a1bSJed Brown } 768c4762a1bSJed Brown } else if (dim == 3 && cellSimplex && size == 2) { 769c4762a1bSJed Brown switch (user->testNum) { 770c4762a1bSJed Brown case 0: { 771c4762a1bSJed Brown PetscInt tetSizes_p2[2] = {1, 2}; 772c4762a1bSJed Brown PetscInt tetPoints_p2[3] = {0, 1, 2}; 773c4762a1bSJed Brown 774c4762a1bSJed Brown ierr = PetscMalloc2(2, &sizes, 3, &points);CHKERRQ(ierr); 775c4762a1bSJed Brown ierr = PetscArraycpy(sizes, tetSizes_p2, 2);CHKERRQ(ierr); 776c4762a1bSJed Brown ierr = PetscArraycpy(points, tetPoints_p2, 3);CHKERRQ(ierr);break;} 777c4762a1bSJed Brown default: 778c4762a1bSJed Brown SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test number %d for triangular mesh on 2 procs", user->testNum); 779c4762a1bSJed Brown } 780c4762a1bSJed Brown } else if (dim == 3 && !cellSimplex && size == 2) { 781c4762a1bSJed Brown switch (user->testNum) { 782c4762a1bSJed Brown case 0: { 783c4762a1bSJed Brown PetscInt hexSizes_p2[2] = {1, 2}; 784c4762a1bSJed Brown PetscInt hexPoints_p2[3] = {0, 1, 2}; 785c4762a1bSJed Brown 786c4762a1bSJed Brown ierr = PetscMalloc2(2, &sizes, 3, &points);CHKERRQ(ierr); 787c4762a1bSJed Brown ierr = PetscArraycpy(sizes, hexSizes_p2, 2);CHKERRQ(ierr); 788c4762a1bSJed Brown ierr = PetscArraycpy(points, hexPoints_p2, 3);CHKERRQ(ierr);break;} 789c4762a1bSJed Brown default: 790c4762a1bSJed Brown SETERRQ1(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 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Could not find matching test partition"); 793c4762a1bSJed Brown } 794c4762a1bSJed Brown ierr = DMPlexGetPartitioner(*dm, &part);CHKERRQ(ierr); 795c4762a1bSJed Brown ierr = PetscPartitionerSetType(part, PETSCPARTITIONERSHELL);CHKERRQ(ierr); 796c4762a1bSJed Brown ierr = PetscPartitionerShellSetPartition(part, size, sizes, points);CHKERRQ(ierr); 797c4762a1bSJed Brown ierr = PetscFree2(sizes, points);CHKERRQ(ierr); 798c4762a1bSJed Brown } 799c4762a1bSJed Brown { 800c4762a1bSJed Brown DM pdm = NULL; 801c4762a1bSJed Brown 802c4762a1bSJed Brown /* Distribute mesh over processes */ 803c4762a1bSJed Brown ierr = DMPlexDistribute(*dm, 0, NULL, &pdm);CHKERRQ(ierr); 804c4762a1bSJed Brown if (pdm) { 805c4762a1bSJed Brown ierr = DMViewFromOptions(pdm, NULL, "-dm_view");CHKERRQ(ierr); 806c4762a1bSJed Brown ierr = DMDestroy(dm);CHKERRQ(ierr); 807c4762a1bSJed Brown *dm = pdm; 808c4762a1bSJed Brown } 809c4762a1bSJed Brown } 810c4762a1bSJed Brown ierr = DMHasLabel(*dm, "pfault", &hasParallelFault);CHKERRQ(ierr); 811c4762a1bSJed Brown if (hasParallelFault) { 812c4762a1bSJed Brown DM dmHybrid = NULL; 813c4762a1bSJed Brown DMLabel faultLabel, faultBdLabel, hybridLabel; 814c4762a1bSJed Brown 815c4762a1bSJed Brown ierr = DMGetLabel(*dm, "pfault", &faultLabel);CHKERRQ(ierr); 816c4762a1bSJed Brown ierr = DMGetLabel(*dm, "pfaultBd", &faultBdLabel);CHKERRQ(ierr); 817c4762a1bSJed Brown ierr = DMPlexCreateHybridMesh(*dm, faultLabel, faultBdLabel, &hybridLabel, NULL, NULL, &dmHybrid);CHKERRQ(ierr); 818c4762a1bSJed Brown ierr = DMLabelView(hybridLabel, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 819c4762a1bSJed Brown ierr = DMLabelDestroy(&hybridLabel);CHKERRQ(ierr); 820c4762a1bSJed Brown ierr = DMDestroy(dm);CHKERRQ(ierr); 821c4762a1bSJed Brown *dm = dmHybrid; 822c4762a1bSJed Brown } 823c4762a1bSJed Brown ierr = PetscObjectSetName((PetscObject) *dm, "Hybrid Mesh");CHKERRQ(ierr); 824b253942bSMatthew G. Knepley ierr = CreateFaultLabel(*dm);CHKERRQ(ierr); 825b253942bSMatthew G. Knepley ierr = CreateDiscretization(*dm, user);CHKERRQ(ierr); 826*012bc364SMatthew G. Knepley ierr = DMViewFromOptions(*dm, NULL, "-dm_view_pre");CHKERRQ(ierr); 827c4762a1bSJed Brown ierr = DMSetFromOptions(*dm);CHKERRQ(ierr); 828c4762a1bSJed Brown ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 829c4762a1bSJed Brown PetscFunctionReturn(0); 830c4762a1bSJed Brown } 831c4762a1bSJed Brown 832b253942bSMatthew G. Knepley static PetscErrorCode TestMesh(DM dm, AppCtx *user) 833b253942bSMatthew G. Knepley { 834b253942bSMatthew G. Knepley PetscErrorCode ierr; 835b253942bSMatthew G. Knepley 836b253942bSMatthew G. Knepley PetscFunctionBegin; 837b253942bSMatthew G. Knepley ierr = DMPlexCheckSymmetry(dm);CHKERRQ(ierr); 838b253942bSMatthew G. Knepley ierr = DMPlexCheckSkeleton(dm, 0);CHKERRQ(ierr); 839b253942bSMatthew G. Knepley ierr = DMPlexCheckFaces(dm, 0);CHKERRQ(ierr); 840b253942bSMatthew G. Knepley PetscFunctionReturn(0); 841b253942bSMatthew G. Knepley } 842b253942bSMatthew G. Knepley 843b253942bSMatthew G. Knepley static PetscErrorCode TestDiscretization(DM dm, AppCtx *user) 844b253942bSMatthew G. Knepley { 845b253942bSMatthew G. Knepley PetscSection s; 846b253942bSMatthew G. Knepley PetscErrorCode ierr; 847b253942bSMatthew G. Knepley 848b253942bSMatthew G. Knepley PetscFunctionBegin; 849b253942bSMatthew G. Knepley ierr = DMGetSection(dm, &s);CHKERRQ(ierr); 850b253942bSMatthew G. Knepley ierr = PetscObjectViewFromOptions((PetscObject) s, NULL, "-local_section_view");CHKERRQ(ierr); 851b253942bSMatthew G. Knepley PetscFunctionReturn(0); 852b253942bSMatthew G. Knepley } 853b253942bSMatthew G. Knepley 854b253942bSMatthew G. Knepley static PetscErrorCode r(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 855b253942bSMatthew G. Knepley { 856b253942bSMatthew G. Knepley PetscInt d; 857b253942bSMatthew G. Knepley for (d = 0; d < dim; ++d) u[d] = x[d]; 858b253942bSMatthew G. Knepley return 0; 859b253942bSMatthew G. Knepley } 860b253942bSMatthew G. Knepley 861b253942bSMatthew G. Knepley static PetscErrorCode rp1(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 862b253942bSMatthew G. Knepley { 863b253942bSMatthew G. Knepley PetscInt d; 864b253942bSMatthew G. Knepley for (d = 0; d < dim; ++d) u[d] = x[d] + (d > 0 ? 1.0 : 0.0); 865b253942bSMatthew G. Knepley return 0; 866b253942bSMatthew G. Knepley } 867b253942bSMatthew G. Knepley 868b253942bSMatthew G. Knepley static PetscErrorCode phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 869b253942bSMatthew G. Knepley { 870b253942bSMatthew G. Knepley PetscInt d; 871b253942bSMatthew G. Knepley u[0] = -x[1]; 872b253942bSMatthew G. Knepley u[1] = x[0]; 873b253942bSMatthew G. Knepley for (d = 2; d < dim; ++d) u[d] = x[d]; 874b253942bSMatthew G. Knepley return 0; 875b253942bSMatthew G. Knepley } 876b253942bSMatthew G. Knepley 877b253942bSMatthew G. Knepley /* \lambda \cdot (\psi_u^- - \psi_u^+) */ 878b253942bSMatthew G. Knepley static void f0_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 879b253942bSMatthew G. Knepley const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 880b253942bSMatthew G. Knepley const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 881b253942bSMatthew G. Knepley PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 882b253942bSMatthew G. Knepley { 883b253942bSMatthew G. Knepley const PetscInt Nc = uOff[1]-uOff[0]; 884b253942bSMatthew G. Knepley PetscInt c; 885b253942bSMatthew G. Knepley for (c = 0; c < Nc; ++c) { 886b253942bSMatthew G. Knepley f0[c] = u[Nc*2+c] + x[Nc-c-1]; 887b253942bSMatthew G. Knepley f0[Nc+c] = -(u[Nc*2+c] + x[Nc-c-1]); 888b253942bSMatthew G. Knepley } 889b253942bSMatthew G. Knepley } 890b253942bSMatthew G. Knepley 891b253942bSMatthew G. Knepley /* (d - u^+ + u^-) \cdot \psi_\lambda */ 892b253942bSMatthew G. Knepley static void f0_bd_l(PetscInt dim, PetscInt Nf, PetscInt NfAux, 893b253942bSMatthew G. Knepley const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 894b253942bSMatthew G. Knepley const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 895b253942bSMatthew G. Knepley PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 896b253942bSMatthew G. Knepley { 897b253942bSMatthew G. Knepley const PetscInt Nc = uOff[2]-uOff[1]; 898b253942bSMatthew G. Knepley PetscInt c; 899b253942bSMatthew G. Knepley 900b253942bSMatthew G. Knepley for (c = 0; c < Nc; ++c) f0[c] = (c > 0 ? 1.0 : 0.0) + u[c] - u[Nc+c]; 901b253942bSMatthew G. Knepley } 902b253942bSMatthew G. Knepley 903b253942bSMatthew G. Knepley /* \psi_lambda \cdot (\psi_u^- - \psi_u^+) */ 904b253942bSMatthew G. Knepley static void g0_bd_ul(PetscInt dim, PetscInt Nf, PetscInt NfAux, 905b253942bSMatthew G. Knepley const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 906b253942bSMatthew G. Knepley const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 907b253942bSMatthew G. Knepley PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 908b253942bSMatthew G. Knepley { 909b253942bSMatthew G. Knepley const PetscInt Nc = uOff[1]-uOff[0]; 910b253942bSMatthew G. Knepley PetscInt c; 911b253942bSMatthew G. Knepley 912b253942bSMatthew G. Knepley for (c = 0; c < Nc; ++c) { 913b253942bSMatthew G. Knepley g0[(0 +c)*Nc+c] = 1.0; 914b253942bSMatthew G. Knepley g0[(Nc+c)*Nc+c] = -1.0; 915b253942bSMatthew G. Knepley } 916b253942bSMatthew G. Knepley } 917b253942bSMatthew G. Knepley 918b253942bSMatthew G. Knepley /* (-\psi_u^+ + \psi_u^-) \cdot \psi_\lambda */ 919b253942bSMatthew G. Knepley static void g0_bd_lu(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, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 923b253942bSMatthew G. Knepley { 924b253942bSMatthew G. Knepley const PetscInt Nc = uOff[2]-uOff[1]; 925b253942bSMatthew G. Knepley PetscInt c; 926b253942bSMatthew G. Knepley 927b253942bSMatthew G. Knepley for (c = 0; c < Nc; ++c) { 928b253942bSMatthew G. Knepley g0[c*Nc*2+c] = 1.0; 929b253942bSMatthew G. Knepley g0[c*Nc*2+Nc+c] = -1.0; 930b253942bSMatthew G. Knepley } 931b253942bSMatthew G. Knepley } 932b253942bSMatthew G. Knepley 933b253942bSMatthew G. Knepley static PetscErrorCode TestAssembly(DM dm, AppCtx *user) 934b253942bSMatthew G. Knepley { 935b253942bSMatthew G. Knepley Mat J; 936b253942bSMatthew G. Knepley Vec locX, locF; 937b253942bSMatthew G. Knepley PetscDS probh; 938b253942bSMatthew G. Knepley DMLabel fault, material; 939b253942bSMatthew G. Knepley IS cohesiveCells; 94006ad1575SMatthew G. Knepley PetscFormKey keys[3]; 941b253942bSMatthew G. Knepley PetscErrorCode (*initialGuess[2])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx); 942b253942bSMatthew G. Knepley PetscInt dim, Nf, cMax, cEnd, id; 943b253942bSMatthew G. Knepley PetscErrorCode ierr; 944b253942bSMatthew G. Knepley 945b253942bSMatthew G. Knepley PetscFunctionBegin; 946b253942bSMatthew G. Knepley ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 947b253942bSMatthew G. Knepley ierr = DMPlexGetSimplexOrBoxCells(dm, 0, NULL, &cMax);CHKERRQ(ierr); 948b253942bSMatthew G. Knepley ierr = DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);CHKERRQ(ierr); 949b253942bSMatthew G. Knepley ierr = ISCreateStride(PETSC_COMM_SELF, cEnd - cMax, cMax, 1, &cohesiveCells);CHKERRQ(ierr); 950b253942bSMatthew G. Knepley ierr = DMGetLabel(dm, "cohesive", &fault);CHKERRQ(ierr); 951b253942bSMatthew G. Knepley ierr = DMGetLocalVector(dm, &locX);CHKERRQ(ierr); 952b253942bSMatthew G. Knepley ierr = PetscObjectSetName((PetscObject) locX, "Local Solution");CHKERRQ(ierr); 953b253942bSMatthew G. Knepley ierr = DMGetLocalVector(dm, &locF);CHKERRQ(ierr); 954b253942bSMatthew G. Knepley ierr = PetscObjectSetName((PetscObject) locF, "Local Residual");CHKERRQ(ierr); 955b253942bSMatthew G. Knepley ierr = DMCreateMatrix(dm, &J);CHKERRQ(ierr); 956b253942bSMatthew G. Knepley ierr = PetscObjectSetName((PetscObject) J, "Jacobian");CHKERRQ(ierr); 957b253942bSMatthew G. Knepley 958b253942bSMatthew G. Knepley /* The initial guess has displacement shifted by one unit in each fault parallel direction across the fault */ 959b253942bSMatthew G. Knepley ierr = DMGetLabel(dm, "material", &material);CHKERRQ(ierr); 960b253942bSMatthew G. Knepley id = 1; 961b253942bSMatthew G. Knepley initialGuess[0] = r; 962b253942bSMatthew G. Knepley initialGuess[1] = NULL; 963b253942bSMatthew G. Knepley ierr = DMProjectFunctionLabelLocal(dm, 0.0, material, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX);CHKERRQ(ierr); 964b253942bSMatthew G. Knepley id = 2; 965b253942bSMatthew G. Knepley initialGuess[0] = rp1; 966b253942bSMatthew G. Knepley initialGuess[1] = NULL; 967b253942bSMatthew G. Knepley ierr = DMProjectFunctionLabelLocal(dm, 0.0, material, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX);CHKERRQ(ierr); 968b253942bSMatthew G. Knepley id = 1; 969b253942bSMatthew G. Knepley initialGuess[0] = NULL; 970b253942bSMatthew G. Knepley initialGuess[1] = phi; 971b253942bSMatthew G. Knepley ierr = DMProjectFunctionLabelLocal(dm, 0.0, fault, 1, &id, PETSC_DETERMINE, NULL, initialGuess, NULL, INSERT_VALUES, locX);CHKERRQ(ierr); 972b253942bSMatthew G. Knepley ierr = VecViewFromOptions(locX, NULL, "-local_solution_view");CHKERRQ(ierr); 973b253942bSMatthew G. Knepley 974b253942bSMatthew G. Knepley ierr = DMGetCellDS(dm, cMax, &probh);CHKERRQ(ierr); 975b253942bSMatthew G. Knepley ierr = PetscDSGetNumFields(probh, &Nf);CHKERRQ(ierr); 976b253942bSMatthew G. Knepley ierr = PetscDSSetBdResidual(probh, 0, f0_bd_u, NULL);CHKERRQ(ierr); 977b253942bSMatthew G. Knepley if (Nf > 1) {ierr = PetscDSSetBdResidual(probh, 1, f0_bd_l, NULL);CHKERRQ(ierr);} 978b253942bSMatthew G. Knepley ierr = PetscDSSetBdJacobian(probh, 0, 1, g0_bd_ul, NULL, NULL, NULL);CHKERRQ(ierr); 979b253942bSMatthew G. Knepley if (Nf > 1) {ierr = PetscDSSetBdJacobian(probh, 1, 0, g0_bd_lu, NULL, NULL, NULL);CHKERRQ(ierr);} 980b253942bSMatthew G. Knepley 9816528b96dSMatthew G. Knepley keys[0].label = material; 9826528b96dSMatthew G. Knepley keys[0].value = 1; 9836528b96dSMatthew G. Knepley keys[0].field = 0; 9846528b96dSMatthew G. Knepley keys[1].label = material; 9856528b96dSMatthew G. Knepley keys[1].value = 2; 9866528b96dSMatthew G. Knepley keys[1].field = 0; 9876528b96dSMatthew G. Knepley keys[2].label = fault; 9886528b96dSMatthew G. Knepley keys[2].value = 1; 9896528b96dSMatthew G. Knepley keys[2].field = 0; 990b253942bSMatthew G. Knepley ierr = VecSet(locF, 0.);CHKERRQ(ierr); 9916528b96dSMatthew G. Knepley ierr = DMPlexComputeResidual_Hybrid_Internal(dm, keys, cohesiveCells, 0.0, locX, NULL, 0.0, locF, user);CHKERRQ(ierr); 992b253942bSMatthew G. Knepley ierr = VecViewFromOptions(locF, NULL, "-local_residual_view");CHKERRQ(ierr); 993b253942bSMatthew G. Knepley ierr = MatZeroEntries(J);CHKERRQ(ierr); 994148442b3SMatthew G. Knepley ierr = DMPlexComputeJacobian_Hybrid_Internal(dm, keys, cohesiveCells, 0.0, 0.0, locX, NULL, J, J, user);CHKERRQ(ierr); 995b253942bSMatthew G. Knepley ierr = MatViewFromOptions(J, NULL, "-local_jacobian_view");CHKERRQ(ierr); 996b253942bSMatthew G. Knepley 997b253942bSMatthew G. Knepley ierr = DMRestoreLocalVector(dm, &locX);CHKERRQ(ierr); 998b253942bSMatthew G. Knepley ierr = DMRestoreLocalVector(dm, &locF);CHKERRQ(ierr); 999b253942bSMatthew G. Knepley ierr = MatDestroy(&J);CHKERRQ(ierr); 1000b253942bSMatthew G. Knepley ierr = ISDestroy(&cohesiveCells);CHKERRQ(ierr); 1001b253942bSMatthew G. Knepley PetscFunctionReturn(0); 1002b253942bSMatthew G. Knepley } 1003b253942bSMatthew G. Knepley 1004c4762a1bSJed Brown int main(int argc, char **argv) 1005c4762a1bSJed Brown { 1006c4762a1bSJed Brown DM dm; 1007c4762a1bSJed Brown AppCtx user; /* user-defined work context */ 1008c4762a1bSJed Brown PetscErrorCode ierr; 1009c4762a1bSJed Brown 1010c4762a1bSJed Brown ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr; 1011c4762a1bSJed Brown ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr); 1012c4762a1bSJed Brown ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr); 1013b253942bSMatthew G. Knepley ierr = TestMesh(dm, &user);CHKERRQ(ierr); 1014b253942bSMatthew G. Knepley ierr = TestDiscretization(dm, &user);CHKERRQ(ierr); 1015b253942bSMatthew G. Knepley ierr = TestAssembly(dm, &user);CHKERRQ(ierr); 1016c4762a1bSJed Brown ierr = DMDestroy(&dm);CHKERRQ(ierr); 1017c4762a1bSJed Brown ierr = PetscFinalize(); 1018c4762a1bSJed Brown return ierr; 1019c4762a1bSJed Brown } 1020c4762a1bSJed Brown 1021c4762a1bSJed Brown /*TEST 1022c4762a1bSJed Brown testset: 1023b253942bSMatthew G. Knepley args: -orig_dm_plex_check_all -dm_plex_check_all \ 102445480ffeSMatthew G. Knepley -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1 -dm_petscds_view -local_section_view \ 1025b253942bSMatthew G. Knepley -local_solution_view -local_residual_view -local_jacobian_view 1026c4762a1bSJed Brown test: 1027c4762a1bSJed Brown suffix: tri_0 1028c4762a1bSJed Brown args: -dim 2 1029c4762a1bSJed Brown test: 1030c4762a1bSJed Brown suffix: tri_t1_0 1031c4762a1bSJed Brown args: -dim 2 -test_num 1 1032c4762a1bSJed Brown test: 1033c4762a1bSJed Brown suffix: tet_0 1034c4762a1bSJed Brown args: -dim 3 1035c4762a1bSJed Brown test: 1036b253942bSMatthew G. Knepley suffix: tet_t1_0 1037b253942bSMatthew G. Knepley args: -dim 3 -test_num 1 1038b253942bSMatthew G. Knepley 1039b253942bSMatthew G. Knepley testset: 1040b253942bSMatthew G. Knepley args: -orig_dm_plex_check_all -dm_plex_check_all \ 104145480ffeSMatthew G. Knepley -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1 -dm_petscds_view 1042b253942bSMatthew G. Knepley test: 1043c4762a1bSJed Brown suffix: tet_1 1044c4762a1bSJed Brown nsize: 2 1045c4762a1bSJed Brown args: -dim 3 1046c4762a1bSJed Brown test: 1047b253942bSMatthew G. Knepley suffix: tri_1 1048b253942bSMatthew G. Knepley nsize: 2 1049b253942bSMatthew G. Knepley args: -dim 2 1050c4762a1bSJed Brown 1051c4762a1bSJed Brown testset: 1052ecfb78b5SMatthew G. Knepley args: -orig_dm_plex_check_all -dm_plex_check_all \ 105345480ffeSMatthew G. Knepley -displacement_petscspace_degree 1 -faulttraction_petscspace_degree 1 -dm_petscds_view 1054c4762a1bSJed Brown # 2D Quads 1055c4762a1bSJed Brown test: 1056c4762a1bSJed Brown suffix: quad_0 1057c4762a1bSJed Brown args: -dim 2 -cell_simplex 0 1058c4762a1bSJed Brown test: 1059c4762a1bSJed Brown suffix: quad_1 1060c4762a1bSJed Brown nsize: 2 1061c4762a1bSJed Brown args: -dim 2 -cell_simplex 0 1062c4762a1bSJed Brown test: 1063c4762a1bSJed Brown suffix: quad_t1_0 106454fcfd0cSMatthew G. Knepley args: -dim 2 -cell_simplex 0 -test_num 1 -faulted_dm_plex_check_all 1065c4762a1bSJed Brown # 3D Hex 1066c4762a1bSJed Brown test: 1067c4762a1bSJed Brown suffix: hex_0 1068c4762a1bSJed Brown args: -dim 3 -cell_simplex 0 1069c4762a1bSJed Brown test: 1070c4762a1bSJed Brown suffix: hex_1 1071c4762a1bSJed Brown nsize: 2 1072c4762a1bSJed Brown args: -dim 3 -cell_simplex 0 1073c4762a1bSJed Brown test: 1074c4762a1bSJed Brown suffix: hex_t1_0 1075c4762a1bSJed Brown args: -dim 3 -cell_simplex 0 -test_num 1 1076c4762a1bSJed Brown test: 1077c4762a1bSJed Brown suffix: hex_t2_0 1078c4762a1bSJed Brown args: -dim 3 -cell_simplex 0 -test_num 2 1079c4762a1bSJed Brown 1080c4762a1bSJed Brown TEST*/ 1081