1 #define PETSCDM_DLL 2 #include <petsc-private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 3 #include <petscdmda.h> 4 #include <petscsf.h> 5 6 #undef __FUNCT__ 7 #define __FUNCT__ "DMSetFromOptions_Plex" 8 PetscErrorCode DMSetFromOptions_Plex(DM dm) 9 { 10 DM_Plex *mesh = (DM_Plex*) dm->data; 11 PetscErrorCode ierr; 12 13 PetscFunctionBegin; 14 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 15 ierr = PetscOptionsHead("DMPlex Options");CHKERRQ(ierr); 16 /* Handle DMPlex refinement */ 17 /* Handle associated vectors */ 18 /* Handle viewing */ 19 ierr = PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMView", PETSC_FALSE, &mesh->printSetValues, NULL);CHKERRQ(ierr); 20 ierr = PetscOptionsInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMView", 0, &mesh->printFEM, NULL);CHKERRQ(ierr); 21 ierr = PetscOptionsReal("-dm_plex_print_tol", "Tolerance for FEM output", "DMView", PETSC_FALSE, &mesh->printTol, NULL);CHKERRQ(ierr); 22 ierr = PetscOptionsTail();CHKERRQ(ierr); 23 PetscFunctionReturn(0); 24 } 25 26 #undef __FUNCT__ 27 #define __FUNCT__ "DMPlexCreateDoublet" 28 /*@ 29 DMPlexCreateDoublet - Creates a mesh of two cells of the specified type, optionally with later refinement. 30 31 Collective on MPI_Comm 32 33 Input Parameters: 34 + comm - The communicator for the DM object 35 . dim - The spatial dimension 36 . simplex - Flag for simplicial cells, otherwise they are tensor product cells 37 . interpolate - Flag to create intermediate mesh pieces (edges, faces) 38 . refinementUniform - Flag for uniform parallel refinement 39 - refinementLimit - A nonzero number indicates the largest admissible volume for a refined cell 40 41 Output Parameter: 42 . dm - The DM object 43 44 Level: beginner 45 46 .keywords: DM, create 47 .seealso: DMSetType(), DMCreate() 48 @*/ 49 PetscErrorCode DMPlexCreateDoublet(MPI_Comm comm, PetscInt dim, PetscBool simplex, PetscBool interpolate, PetscBool refinementUniform, PetscReal refinementLimit, DM *newdm) 50 { 51 DM dm; 52 PetscInt p; 53 PetscMPIInt rank; 54 PetscErrorCode ierr; 55 56 PetscFunctionBegin; 57 ierr = DMCreate(comm, &dm);CHKERRQ(ierr); 58 ierr = DMSetType(dm, DMPLEX);CHKERRQ(ierr); 59 ierr = DMPlexSetDimension(dm, dim);CHKERRQ(ierr); 60 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 61 switch (dim) { 62 case 2: 63 if (simplex) {ierr = PetscObjectSetName((PetscObject) dm, "triangular");CHKERRQ(ierr);} 64 else {ierr = PetscObjectSetName((PetscObject) dm, "quadrilateral");CHKERRQ(ierr);} 65 break; 66 case 3: 67 if (simplex) {ierr = PetscObjectSetName((PetscObject) dm, "tetrahedral");CHKERRQ(ierr);} 68 else {ierr = PetscObjectSetName((PetscObject) dm, "hexahedral");CHKERRQ(ierr);} 69 break; 70 default: 71 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %d", dim); 72 } 73 if (rank) { 74 PetscInt numPoints[2] = {0, 0}; 75 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 76 } else { 77 switch (dim) { 78 case 2: 79 if (simplex) { 80 PetscInt numPoints[2] = {4, 2}; 81 PetscInt coneSize[6] = {3, 3, 0, 0, 0, 0}; 82 PetscInt cones[6] = {2, 3, 4, 5, 4, 3}; 83 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 84 PetscScalar vertexCoords[8] = {-0.5, 0.5, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5}; 85 PetscInt markerPoints[8] = {2, 1, 3, 1, 4, 1, 5, 1}; 86 87 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 88 for (p = 0; p < 4; ++p) {ierr = DMPlexSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 89 } else { 90 PetscInt numPoints[2] = {6, 2}; 91 PetscInt coneSize[8] = {4, 4, 0, 0, 0, 0, 0, 0}; 92 PetscInt cones[8] = {2, 3, 4, 5, 3, 6, 7, 4}; 93 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 94 PetscScalar vertexCoords[12] = {-1.0, -0.5, 0.0, -0.5, 0.0, 0.5, -1.0, 0.5, 1.0, -0.5, 1.0, 0.5}; 95 96 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 97 } 98 break; 99 case 3: 100 if (simplex) { 101 PetscInt numPoints[2] = {5, 2}; 102 PetscInt coneSize[7] = {4, 4, 0, 0, 0, 0, 0}; 103 PetscInt cones[8] = {4, 3, 5, 2, 5, 3, 4, 6}; 104 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 105 PetscScalar vertexCoords[15] = {-1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0}; 106 PetscInt markerPoints[10] = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1}; 107 108 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 109 for (p = 0; p < 5; ++p) {ierr = DMPlexSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 110 } else { 111 PetscInt numPoints[2] = {12, 2}; 112 PetscInt coneSize[14] = {8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 113 PetscInt cones[16] = {2, 3, 4, 5, 6, 7, 8, 9, 5, 4, 10, 11, 7, 12, 13, 8}; 114 PetscInt coneOrientations[16] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 115 PetscScalar vertexCoords[36] = {-1.0, -0.5, -0.5, -1.0, 0.5, -0.5, 0.0, 0.5, -0.5, 0.0, -0.5, -0.5, 116 -1.0, -0.5, 0.5, 0.0, -0.5, 0.5, 0.0, 0.5, 0.5, -1.0, 0.5, 0.5, 117 1.0, 0.5, -0.5, 1.0, -0.5, -0.5, 1.0, -0.5, 0.5, 1.0, 0.5, 0.5}; 118 119 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 120 } 121 break; 122 default: 123 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %d", dim); 124 } 125 } 126 *newdm = dm; 127 if (refinementLimit > 0.0) { 128 DM rdm; 129 const char *name; 130 131 ierr = DMPlexSetRefinementUniform(*newdm, PETSC_FALSE);CHKERRQ(ierr); 132 ierr = DMPlexSetRefinementLimit(*newdm, refinementLimit);CHKERRQ(ierr); 133 ierr = DMRefine(*newdm, comm, &rdm);CHKERRQ(ierr); 134 ierr = PetscObjectGetName((PetscObject) *newdm, &name);CHKERRQ(ierr); 135 ierr = PetscObjectSetName((PetscObject) rdm, name);CHKERRQ(ierr); 136 ierr = DMDestroy(newdm);CHKERRQ(ierr); 137 *newdm = rdm; 138 } 139 if (interpolate) { 140 DM idm; 141 const char *name; 142 143 ierr = DMPlexInterpolate(*newdm, &idm);CHKERRQ(ierr); 144 ierr = PetscObjectGetName((PetscObject) *newdm, &name);CHKERRQ(ierr); 145 ierr = PetscObjectSetName((PetscObject) idm, name);CHKERRQ(ierr); 146 ierr = DMPlexCopyCoordinates(*newdm, idm);CHKERRQ(ierr); 147 ierr = DMPlexCopyLabels(*newdm, idm);CHKERRQ(ierr); 148 ierr = DMDestroy(newdm);CHKERRQ(ierr); 149 *newdm = idm; 150 } 151 { 152 DM refinedMesh = NULL; 153 DM distributedMesh = NULL; 154 155 /* Distribute mesh over processes */ 156 ierr = DMPlexDistribute(*newdm, NULL, 0, NULL, &distributedMesh);CHKERRQ(ierr); 157 if (distributedMesh) { 158 ierr = DMDestroy(newdm);CHKERRQ(ierr); 159 *newdm = distributedMesh; 160 } 161 if (refinementUniform) { 162 ierr = DMPlexSetRefinementUniform(*newdm, refinementUniform);CHKERRQ(ierr); 163 ierr = DMRefine(*newdm, comm, &refinedMesh);CHKERRQ(ierr); 164 if (refinedMesh) { 165 ierr = DMDestroy(newdm);CHKERRQ(ierr); 166 *newdm = refinedMesh; 167 } 168 } 169 } 170 PetscFunctionReturn(0); 171 } 172 173 #undef __FUNCT__ 174 #define __FUNCT__ "DMPlexCreateSquareBoundary" 175 /* 176 Simple square boundary: 177 178 18--5-17--4--16 179 | | | 180 6 10 3 181 | | | 182 19-11-20--9--15 183 | | | 184 7 8 2 185 | | | 186 12--0-13--1--14 187 */ 188 PetscErrorCode DMPlexCreateSquareBoundary(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[]) 189 { 190 PetscInt numVertices = (edges[0]+1)*(edges[1]+1); 191 PetscInt numEdges = edges[0]*(edges[1]+1) + (edges[0]+1)*edges[1]; 192 PetscInt markerTop = 1; 193 PetscInt markerBottom = 1; 194 PetscInt markerRight = 1; 195 PetscInt markerLeft = 1; 196 PetscBool markerSeparate = PETSC_FALSE; 197 Vec coordinates; 198 PetscSection coordSection; 199 PetscScalar *coords; 200 PetscInt coordSize; 201 PetscMPIInt rank; 202 PetscInt v, vx, vy; 203 PetscErrorCode ierr; 204 205 PetscFunctionBegin; 206 ierr = PetscOptionsGetBool(((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);CHKERRQ(ierr); 207 if (markerSeparate) { 208 markerTop = 1; 209 markerBottom = 0; 210 markerRight = 0; 211 markerLeft = 0; 212 } 213 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 214 if (!rank) { 215 PetscInt e, ex, ey; 216 217 ierr = DMPlexSetChart(dm, 0, numEdges+numVertices);CHKERRQ(ierr); 218 for (e = 0; e < numEdges; ++e) { 219 ierr = DMPlexSetConeSize(dm, e, 2);CHKERRQ(ierr); 220 } 221 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 222 for (vx = 0; vx <= edges[0]; vx++) { 223 for (ey = 0; ey < edges[1]; ey++) { 224 PetscInt edge = vx*edges[1] + ey + edges[0]*(edges[1]+1); 225 PetscInt vertex = ey*(edges[0]+1) + vx + numEdges; 226 PetscInt cone[2]; 227 228 cone[0] = vertex; cone[1] = vertex+edges[0]+1; 229 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 230 if (vx == edges[0]) { 231 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 232 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerRight);CHKERRQ(ierr); 233 if (ey == edges[1]-1) { 234 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerRight);CHKERRQ(ierr); 235 } 236 } else if (vx == 0) { 237 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 238 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerLeft);CHKERRQ(ierr); 239 if (ey == edges[1]-1) { 240 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerLeft);CHKERRQ(ierr); 241 } 242 } 243 } 244 } 245 for (vy = 0; vy <= edges[1]; vy++) { 246 for (ex = 0; ex < edges[0]; ex++) { 247 PetscInt edge = vy*edges[0] + ex; 248 PetscInt vertex = vy*(edges[0]+1) + ex + numEdges; 249 PetscInt cone[2]; 250 251 cone[0] = vertex; cone[1] = vertex+1; 252 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 253 if (vy == edges[1]) { 254 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 255 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerTop);CHKERRQ(ierr); 256 if (ex == edges[0]-1) { 257 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerTop);CHKERRQ(ierr); 258 } 259 } else if (vy == 0) { 260 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 261 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerBottom);CHKERRQ(ierr); 262 if (ex == edges[0]-1) { 263 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerBottom);CHKERRQ(ierr); 264 } 265 } 266 } 267 } 268 } 269 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 270 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 271 /* Build coordinates */ 272 ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 273 ierr = PetscSectionSetChart(coordSection, numEdges, numEdges + numVertices);CHKERRQ(ierr); 274 for (v = numEdges; v < numEdges+numVertices; ++v) { 275 ierr = PetscSectionSetDof(coordSection, v, 2);CHKERRQ(ierr); 276 } 277 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 278 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 279 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 280 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 281 ierr = VecSetFromOptions(coordinates);CHKERRQ(ierr); 282 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 283 for (vy = 0; vy <= edges[1]; ++vy) { 284 for (vx = 0; vx <= edges[0]; ++vx) { 285 coords[(vy*(edges[0]+1)+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/edges[0])*vx; 286 coords[(vy*(edges[0]+1)+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/edges[1])*vy; 287 } 288 } 289 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 290 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 291 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 292 PetscFunctionReturn(0); 293 } 294 295 #undef __FUNCT__ 296 #define __FUNCT__ "DMPlexCreateCubeBoundary" 297 /* 298 Simple cubic boundary: 299 300 2-------3 301 /| /| 302 6-------7 | 303 | | | | 304 | 0-----|-1 305 |/ |/ 306 4-------5 307 */ 308 PetscErrorCode DMPlexCreateCubeBoundary(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt faces[]) 309 { 310 PetscInt numVertices = (faces[0]+1)*(faces[1]+1)*(faces[2]+1); 311 PetscInt numFaces = 6; 312 Vec coordinates; 313 PetscSection coordSection; 314 PetscScalar *coords; 315 PetscInt coordSize; 316 PetscMPIInt rank; 317 PetscInt v, vx, vy, vz; 318 PetscErrorCode ierr; 319 320 PetscFunctionBegin; 321 if ((faces[0] < 1) || (faces[1] < 1) || (faces[2] < 1)) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Must have at least 1 face per side"); 322 if ((faces[0] > 1) || (faces[1] > 1) || (faces[2] > 1)) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Currently can't handle more than 1 face per side"); 323 ierr = PetscMalloc(numVertices*2 * sizeof(PetscReal), &coords);CHKERRQ(ierr); 324 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 325 if (!rank) { 326 PetscInt f; 327 328 ierr = DMPlexSetChart(dm, 0, numFaces+numVertices);CHKERRQ(ierr); 329 for (f = 0; f < numFaces; ++f) { 330 ierr = DMPlexSetConeSize(dm, f, 4);CHKERRQ(ierr); 331 } 332 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 333 for (v = 0; v < numFaces+numVertices; ++v) { 334 ierr = DMPlexSetLabelValue(dm, "marker", v, 1);CHKERRQ(ierr); 335 } 336 { /* Side 0 (Front) */ 337 PetscInt cone[4]; 338 cone[0] = numFaces+4; cone[1] = numFaces+5; cone[2] = numFaces+7; cone[3] = numFaces+6; 339 ierr = DMPlexSetCone(dm, 0, cone);CHKERRQ(ierr); 340 } 341 { /* Side 1 (Back) */ 342 PetscInt cone[4]; 343 cone[0] = numFaces+1; cone[1] = numFaces+0; cone[2] = numFaces+2; cone[3] = numFaces+3; 344 ierr = DMPlexSetCone(dm, 1, cone);CHKERRQ(ierr); 345 } 346 { /* Side 2 (Bottom) */ 347 PetscInt cone[4]; 348 cone[0] = numFaces+0; cone[1] = numFaces+1; cone[2] = numFaces+5; cone[3] = numFaces+4; 349 ierr = DMPlexSetCone(dm, 2, cone);CHKERRQ(ierr); 350 } 351 { /* Side 3 (Top) */ 352 PetscInt cone[4]; 353 cone[0] = numFaces+6; cone[1] = numFaces+7; cone[2] = numFaces+3; cone[3] = numFaces+2; 354 ierr = DMPlexSetCone(dm, 3, cone);CHKERRQ(ierr); 355 } 356 { /* Side 4 (Left) */ 357 PetscInt cone[4]; 358 cone[0] = numFaces+0; cone[1] = numFaces+4; cone[2] = numFaces+6; cone[3] = numFaces+2; 359 ierr = DMPlexSetCone(dm, 4, cone);CHKERRQ(ierr); 360 } 361 { /* Side 5 (Right) */ 362 PetscInt cone[4]; 363 cone[0] = numFaces+5; cone[1] = numFaces+1; cone[2] = numFaces+3; cone[3] = numFaces+7; 364 ierr = DMPlexSetCone(dm, 5, cone);CHKERRQ(ierr); 365 } 366 } 367 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 368 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 369 /* Build coordinates */ 370 ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 371 ierr = PetscSectionSetChart(coordSection, numFaces, numFaces + numVertices);CHKERRQ(ierr); 372 for (v = numFaces; v < numFaces+numVertices; ++v) { 373 ierr = PetscSectionSetDof(coordSection, v, 3);CHKERRQ(ierr); 374 } 375 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 376 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 377 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 378 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 379 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 380 ierr = VecSetFromOptions(coordinates);CHKERRQ(ierr); 381 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 382 for (vz = 0; vz <= faces[2]; ++vz) { 383 for (vy = 0; vy <= faces[1]; ++vy) { 384 for (vx = 0; vx <= faces[0]; ++vx) { 385 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+0] = lower[0] + ((upper[0] - lower[0])/faces[0])*vx; 386 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+1] = lower[1] + ((upper[1] - lower[1])/faces[1])*vy; 387 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+2] = lower[2] + ((upper[2] - lower[2])/faces[2])*vz; 388 } 389 } 390 } 391 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 392 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 393 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 394 PetscFunctionReturn(0); 395 } 396 397 #undef __FUNCT__ 398 #define __FUNCT__ "DMPlexCreateSquareMesh" 399 /* 400 Simple square mesh: 401 402 22--8-23--9--24 403 | | | 404 13 2 14 3 15 405 | | | 406 19--6-20--7--21 407 | | | 408 10 0 11 1 12 409 | | | 410 16--4-17--5--18 411 */ 412 PetscErrorCode DMPlexCreateSquareMesh(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[]) 413 { 414 PetscInt markerTop = 1; 415 PetscInt markerBottom = 1; 416 PetscInt markerRight = 1; 417 PetscInt markerLeft = 1; 418 PetscBool markerSeparate = PETSC_FALSE; 419 PetscMPIInt rank; 420 PetscErrorCode ierr; 421 422 PetscFunctionBegin; 423 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 424 ierr = PetscOptionsGetBool(((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);CHKERRQ(ierr); 425 if (markerSeparate) { 426 markerTop = 3; 427 markerBottom = 1; 428 markerRight = 2; 429 markerLeft = 4; 430 } 431 { 432 const PetscInt numXEdges = !rank ? edges[0] : 0; 433 const PetscInt numYEdges = !rank ? edges[1] : 0; 434 const PetscInt numXVertices = !rank ? edges[0]+1 : 0; 435 const PetscInt numYVertices = !rank ? edges[1]+1 : 0; 436 const PetscInt numTotXEdges = numXEdges*numYVertices; 437 const PetscInt numTotYEdges = numYEdges*numXVertices; 438 const PetscInt numVertices = numXVertices*numYVertices; 439 const PetscInt numEdges = numTotXEdges + numTotYEdges; 440 const PetscInt numFaces = numXEdges*numYEdges; 441 const PetscInt firstVertex = numFaces; 442 const PetscInt firstXEdge = numFaces + numVertices; 443 const PetscInt firstYEdge = numFaces + numVertices + numTotXEdges; 444 Vec coordinates; 445 PetscSection coordSection; 446 PetscScalar *coords; 447 PetscInt coordSize; 448 PetscInt v, vx, vy; 449 PetscInt f, fx, fy, e, ex, ey; 450 451 ierr = DMPlexSetChart(dm, 0, numFaces+numEdges+numVertices);CHKERRQ(ierr); 452 for (f = 0; f < numFaces; ++f) { 453 ierr = DMPlexSetConeSize(dm, f, 4);CHKERRQ(ierr); 454 } 455 for (e = firstXEdge; e < firstXEdge+numEdges; ++e) { 456 ierr = DMPlexSetConeSize(dm, e, 2);CHKERRQ(ierr); 457 } 458 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 459 /* Build faces */ 460 for (fy = 0; fy < numYEdges; fy++) { 461 for (fx = 0; fx < numXEdges; fx++) { 462 const PetscInt face = fy*numXEdges + fx; 463 const PetscInt edgeL = firstYEdge + fx*numYEdges + fy; 464 const PetscInt edgeB = firstXEdge + fy*numXEdges + fx; 465 const PetscInt ornt[4] = {0, 0, -2, -2}; 466 PetscInt cone[4]; 467 468 cone[0] = edgeB; cone[1] = edgeL+numYEdges; cone[2] = edgeB+numXEdges; cone[3] = edgeL; 469 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr); 470 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr); 471 } 472 } 473 /* Build Y edges*/ 474 for (vx = 0; vx < numXVertices; vx++) { 475 for (ey = 0; ey < numYEdges; ey++) { 476 const PetscInt edge = firstYEdge + vx*numYEdges + ey; 477 const PetscInt vertex = firstVertex + ey*numXVertices + vx; 478 PetscInt cone[2]; 479 480 cone[0] = vertex; cone[1] = vertex+numXVertices; 481 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 482 if (vx == numXVertices-1) { 483 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 484 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerRight);CHKERRQ(ierr); 485 if (ey == numYEdges-1) { 486 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerRight);CHKERRQ(ierr); 487 } 488 } else if (vx == 0) { 489 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 490 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerLeft);CHKERRQ(ierr); 491 if (ey == numYEdges-1) { 492 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerLeft);CHKERRQ(ierr); 493 } 494 } 495 } 496 } 497 /* Build X edges*/ 498 for (vy = 0; vy < numYVertices; vy++) { 499 for (ex = 0; ex < numXEdges; ex++) { 500 const PetscInt edge = firstXEdge + vy*numXEdges + ex; 501 const PetscInt vertex = firstVertex + vy*numXVertices + ex; 502 PetscInt cone[2]; 503 504 cone[0] = vertex; cone[1] = vertex+1; 505 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 506 if (vy == numYVertices-1) { 507 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 508 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerTop);CHKERRQ(ierr); 509 if (ex == numXEdges-1) { 510 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerTop);CHKERRQ(ierr); 511 } 512 } else if (vy == 0) { 513 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 514 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerBottom);CHKERRQ(ierr); 515 if (ex == numXEdges-1) { 516 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerBottom);CHKERRQ(ierr); 517 } 518 } 519 } 520 } 521 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 522 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 523 /* Build coordinates */ 524 ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 525 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVertices);CHKERRQ(ierr); 526 for (v = firstVertex; v < firstVertex+numVertices; ++v) { 527 ierr = PetscSectionSetDof(coordSection, v, 2);CHKERRQ(ierr); 528 } 529 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 530 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 531 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 532 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 533 ierr = VecSetFromOptions(coordinates);CHKERRQ(ierr); 534 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 535 for (vy = 0; vy < numYVertices; ++vy) { 536 for (vx = 0; vx < numXVertices; ++vx) { 537 coords[(vy*numXVertices+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/numXEdges)*vx; 538 coords[(vy*numXVertices+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/numYEdges)*vy; 539 } 540 } 541 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 542 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 543 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 544 } 545 PetscFunctionReturn(0); 546 } 547 548 #undef __FUNCT__ 549 #define __FUNCT__ "DMPlexCreateBoxMesh" 550 /*@ 551 DMPlexCreateBoxMesh - Creates a mesh on the tensor product of unit intervals (box). 552 553 Collective on MPI_Comm 554 555 Input Parameters: 556 + comm - The communicator for the DM object 557 . dim - The spatial dimension 558 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 559 560 Output Parameter: 561 . dm - The DM object 562 563 Level: beginner 564 565 .keywords: DM, create 566 .seealso: DMSetType(), DMCreate() 567 @*/ 568 PetscErrorCode DMPlexCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscBool interpolate, DM *dm) 569 { 570 DM boundary; 571 PetscErrorCode ierr; 572 573 PetscFunctionBegin; 574 PetscValidPointer(dm, 4); 575 ierr = DMCreate(comm, &boundary);CHKERRQ(ierr); 576 PetscValidLogicalCollectiveInt(boundary,dim,2); 577 ierr = DMSetType(boundary, DMPLEX);CHKERRQ(ierr); 578 ierr = DMPlexSetDimension(boundary, dim-1);CHKERRQ(ierr); 579 switch (dim) { 580 case 2: 581 { 582 PetscReal lower[2] = {0.0, 0.0}; 583 PetscReal upper[2] = {1.0, 1.0}; 584 PetscInt edges[2] = {2, 2}; 585 586 ierr = DMPlexCreateSquareBoundary(boundary, lower, upper, edges);CHKERRQ(ierr); 587 break; 588 } 589 case 3: 590 { 591 PetscReal lower[3] = {0.0, 0.0, 0.0}; 592 PetscReal upper[3] = {1.0, 1.0, 1.0}; 593 PetscInt faces[3] = {1, 1, 1}; 594 595 ierr = DMPlexCreateCubeBoundary(boundary, lower, upper, faces);CHKERRQ(ierr); 596 break; 597 } 598 default: 599 SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim); 600 } 601 ierr = DMPlexGenerate(boundary, NULL, interpolate, dm);CHKERRQ(ierr); 602 ierr = DMDestroy(&boundary);CHKERRQ(ierr); 603 PetscFunctionReturn(0); 604 } 605 606 #undef __FUNCT__ 607 #define __FUNCT__ "DMPlexCreateHexBoxMesh" 608 PetscErrorCode DMPlexCreateHexBoxMesh(MPI_Comm comm, PetscInt dim, const PetscInt cells[], DM *dm) 609 { 610 PetscErrorCode ierr; 611 612 PetscFunctionBegin; 613 PetscValidPointer(dm, 4); 614 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 615 PetscValidLogicalCollectiveInt(*dm,dim,2); 616 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 617 ierr = DMPlexSetDimension(*dm, dim);CHKERRQ(ierr); 618 switch (dim) { 619 case 2: 620 { 621 PetscReal lower[2] = {0.0, 0.0}; 622 PetscReal upper[2] = {1.0, 1.0}; 623 624 ierr = DMPlexCreateSquareMesh(*dm, lower, upper, cells);CHKERRQ(ierr); 625 break; 626 } 627 #if 0 628 case 3: 629 { 630 PetscReal lower[3] = {0.0, 0.0, 0.0}; 631 PetscReal upper[3] = {1.0, 1.0, 1.0}; 632 633 ierr = DMPlexCreateCubeMesh(boundary, lower, upper, cells);CHKERRQ(ierr); 634 break; 635 } 636 #endif 637 default: 638 SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim); 639 } 640 PetscFunctionReturn(0); 641 } 642 643 /* External function declarations here */ 644 extern PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling); 645 extern PetscErrorCode DMCreateMatrix_Plex(DM dm, MatType mtype, Mat *J); 646 extern PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm); 647 extern PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *dmRefined); 648 extern PetscErrorCode DMClone_Plex(DM dm, DM *newdm); 649 extern PetscErrorCode DMSetUp_Plex(DM dm); 650 extern PetscErrorCode DMDestroy_Plex(DM dm); 651 extern PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer); 652 extern PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm); 653 extern PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, IS *cellIS); 654 655 #undef __FUNCT__ 656 #define __FUNCT__ "DMCreateGlobalVector_Plex" 657 static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec) 658 { 659 PetscErrorCode ierr; 660 661 PetscFunctionBegin; 662 ierr = DMCreateGlobalVector_Section_Private(dm,vec);CHKERRQ(ierr); 663 /* ierr = VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM);CHKERRQ(ierr); */ 664 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void))VecView_Plex);CHKERRQ(ierr); 665 PetscFunctionReturn(0); 666 } 667 668 #undef __FUNCT__ 669 #define __FUNCT__ "DMCreateLocalVector_Plex" 670 static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec) 671 { 672 PetscErrorCode ierr; 673 674 PetscFunctionBegin; 675 ierr = DMCreateLocalVector_Section_Private(dm,vec);CHKERRQ(ierr); 676 ierr = VecSetOperation(*vec, VECOP_VIEW, (void(*)(void)) VecView_Plex_Local);CHKERRQ(ierr); 677 PetscFunctionReturn(0); 678 } 679 680 #undef __FUNCT__ 681 #define __FUNCT__ "DMInitialize_Plex" 682 PetscErrorCode DMInitialize_Plex(DM dm) 683 { 684 PetscErrorCode ierr; 685 686 PetscFunctionBegin; 687 ierr = PetscStrallocpy(VECSTANDARD, (char**)&dm->vectype);CHKERRQ(ierr); 688 689 dm->ops->view = DMView_Plex; 690 dm->ops->setfromoptions = DMSetFromOptions_Plex; 691 dm->ops->clone = DMClone_Plex; 692 dm->ops->setup = DMSetUp_Plex; 693 dm->ops->createglobalvector = DMCreateGlobalVector_Plex; 694 dm->ops->createlocalvector = DMCreateLocalVector_Plex; 695 dm->ops->getlocaltoglobalmapping = NULL; 696 dm->ops->getlocaltoglobalmappingblock = NULL; 697 dm->ops->createfieldis = NULL; 698 dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex; 699 dm->ops->getcoloring = 0; 700 dm->ops->creatematrix = DMCreateMatrix_Plex; 701 dm->ops->createinterpolation = 0; 702 dm->ops->getaggregates = 0; 703 dm->ops->getinjection = 0; 704 dm->ops->refine = DMRefine_Plex; 705 dm->ops->coarsen = 0; 706 dm->ops->refinehierarchy = 0; 707 dm->ops->coarsenhierarchy = 0; 708 dm->ops->globaltolocalbegin = NULL; 709 dm->ops->globaltolocalend = NULL; 710 dm->ops->localtoglobalbegin = NULL; 711 dm->ops->localtoglobalend = NULL; 712 dm->ops->destroy = DMDestroy_Plex; 713 dm->ops->createsubdm = DMCreateSubDM_Plex; 714 dm->ops->locatepoints = DMLocatePoints_Plex; 715 PetscFunctionReturn(0); 716 } 717 718 #undef __FUNCT__ 719 #define __FUNCT__ "DMClone_Plex" 720 PetscErrorCode DMClone_Plex(DM dm, DM *newdm) 721 { 722 DM_Plex *mesh = (DM_Plex *) dm->data; 723 PetscErrorCode ierr; 724 725 PetscFunctionBegin; 726 mesh->refct++; 727 (*newdm)->data = mesh; 728 ierr = PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);CHKERRQ(ierr); 729 ierr = DMInitialize_Plex(*newdm);CHKERRQ(ierr); 730 PetscFunctionReturn(0); 731 } 732 733 /*MC 734 DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram. 735 In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is 736 specified by a PetscSection object. Ownership in the global representation is determined by 737 ownership of the underlying DMPlex points. This is specified by another PetscSection object. 738 739 Level: intermediate 740 741 .seealso: DMType, DMPlexCreate(), DMCreate(), DMSetType() 742 M*/ 743 744 #undef __FUNCT__ 745 #define __FUNCT__ "DMCreate_Plex" 746 PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm) 747 { 748 DM_Plex *mesh; 749 PetscInt unit, d; 750 PetscErrorCode ierr; 751 752 PetscFunctionBegin; 753 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 754 ierr = PetscNewLog(dm, DM_Plex, &mesh);CHKERRQ(ierr); 755 dm->data = mesh; 756 757 mesh->refct = 1; 758 mesh->dim = 0; 759 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);CHKERRQ(ierr); 760 mesh->maxConeSize = 0; 761 mesh->cones = NULL; 762 mesh->coneOrientations = NULL; 763 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);CHKERRQ(ierr); 764 mesh->maxSupportSize = 0; 765 mesh->supports = NULL; 766 mesh->refinementUniform = PETSC_TRUE; 767 mesh->refinementLimit = -1.0; 768 769 mesh->facesTmp = NULL; 770 771 mesh->subpointMap = NULL; 772 773 for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0; 774 775 mesh->labels = NULL; 776 mesh->depthLabel = NULL; 777 mesh->globalVertexNumbers = NULL; 778 mesh->globalCellNumbers = NULL; 779 for (d = 0; d < 8; ++d) mesh->hybridPointMax[d] = PETSC_DETERMINE; 780 mesh->vtkCellHeight = 0; 781 mesh->preallocCenterDim = -1; 782 783 mesh->printSetValues = PETSC_FALSE; 784 mesh->printFEM = 0; 785 mesh->printTol = 1.0e-10; 786 787 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 788 PetscFunctionReturn(0); 789 } 790 791 #undef __FUNCT__ 792 #define __FUNCT__ "DMPlexCreate" 793 /*@ 794 DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram. 795 796 Collective on MPI_Comm 797 798 Input Parameter: 799 . comm - The communicator for the DMPlex object 800 801 Output Parameter: 802 . mesh - The DMPlex object 803 804 Level: beginner 805 806 .keywords: DMPlex, create 807 @*/ 808 PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh) 809 { 810 PetscErrorCode ierr; 811 812 PetscFunctionBegin; 813 PetscValidPointer(mesh,2); 814 ierr = DMCreate(comm, mesh);CHKERRQ(ierr); 815 ierr = DMSetType(*mesh, DMPLEX);CHKERRQ(ierr); 816 PetscFunctionReturn(0); 817 } 818 819 #undef __FUNCT__ 820 #define __FUNCT__ "DMPlexBuildFromCellList_Private" 821 /* 822 This takes as input the common mesh generator output, a list of the vertices for each cell 823 */ 824 PetscErrorCode DMPlexBuildFromCellList_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[]) 825 { 826 PetscInt *cone, c, p; 827 PetscErrorCode ierr; 828 829 PetscFunctionBegin; 830 ierr = DMPlexSetChart(dm, 0, numCells+numVertices);CHKERRQ(ierr); 831 for (c = 0; c < numCells; ++c) { 832 ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr); 833 } 834 ierr = DMSetUp(dm);CHKERRQ(ierr); 835 ierr = DMGetWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 836 for (c = 0; c < numCells; ++c) { 837 for (p = 0; p < numCorners; ++p) { 838 cone[p] = cells[c*numCorners+p]+numCells; 839 } 840 ierr = DMPlexSetCone(dm, c, cone);CHKERRQ(ierr); 841 } 842 ierr = DMRestoreWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 843 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 844 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 845 PetscFunctionReturn(0); 846 } 847 848 #undef __FUNCT__ 849 #define __FUNCT__ "DMPlexBuildCoordinates_Private" 850 /* 851 This takes as input the coordinates for each vertex 852 */ 853 PetscErrorCode DMPlexBuildCoordinates_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, const double vertexCoords[]) 854 { 855 PetscSection coordSection; 856 Vec coordinates; 857 PetscScalar *coords; 858 PetscInt coordSize, v, d; 859 PetscErrorCode ierr; 860 861 PetscFunctionBegin; 862 ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 863 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 864 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 865 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVertices);CHKERRQ(ierr); 866 for (v = numCells; v < numCells+numVertices; ++v) { 867 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 868 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 869 } 870 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 871 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 872 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 873 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 874 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 875 ierr = VecSetFromOptions(coordinates);CHKERRQ(ierr); 876 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 877 for (v = 0; v < numVertices; ++v) { 878 for (d = 0; d < spaceDim; ++d) { 879 coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d]; 880 } 881 } 882 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 883 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 884 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 885 PetscFunctionReturn(0); 886 } 887 888 #undef __FUNCT__ 889 #define __FUNCT__ "DMPlexCreateFromCellList" 890 /*@C 891 DMPlexCreateFromCellList - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM 892 893 Input Parameters: 894 + comm - The communicator 895 . dim - The topological dimension of the mesh 896 . numCells - The number of cells 897 . numVertices - The number of vertices 898 . numCorners - The number of vertices for each cell 899 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 900 . cells - An array of numCells*numCorners numbers, the vertices for each cell 901 . spaceDim - The spatial dimension used for coordinates 902 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 903 904 Output Parameter: 905 . dm - The DM 906 907 Note: Two triangles sharing a face 908 $ 909 $ 2 910 $ / | \ 911 $ / | \ 912 $ / | \ 913 $ 0 0 | 1 3 914 $ \ | / 915 $ \ | / 916 $ \ | / 917 $ 1 918 would have input 919 $ numCells = 2, numVertices = 4 920 $ cells = [0 1 2 1 3 2] 921 $ 922 which would result in the DMPlex 923 $ 924 $ 4 925 $ / | \ 926 $ / | \ 927 $ / | \ 928 $ 2 0 | 1 5 929 $ \ | / 930 $ \ | / 931 $ \ | / 932 $ 3 933 934 Level: beginner 935 936 .seealso: DMPlexCreate() 937 @*/ 938 PetscErrorCode DMPlexCreateFromCellList(MPI_Comm comm, PetscInt dim, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, PetscBool interpolate, const int cells[], PetscInt spaceDim, const double vertexCoords[], DM *dm) 939 { 940 PetscErrorCode ierr; 941 942 PetscFunctionBegin; 943 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 944 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 945 ierr = DMPlexSetDimension(*dm, dim);CHKERRQ(ierr); 946 ierr = DMPlexBuildFromCellList_Private(*dm, numCells, numVertices, numCorners, cells);CHKERRQ(ierr); 947 if (interpolate) { 948 DM idm; 949 950 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 951 ierr = DMDestroy(dm);CHKERRQ(ierr); 952 *dm = idm; 953 } 954 ierr = DMPlexBuildCoordinates_Private(*dm, spaceDim, numCells, numVertices, vertexCoords);CHKERRQ(ierr); 955 PetscFunctionReturn(0); 956 } 957 958 #undef __FUNCT__ 959 #define __FUNCT__ "DMPlexCreateFromDAG" 960 /* 961 This takes as input the raw Hasse Diagram data 962 */ 963 PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[]) 964 { 965 Vec coordinates; 966 PetscSection coordSection; 967 PetscScalar *coords; 968 PetscInt coordSize, firstVertex = numPoints[depth], pStart = 0, pEnd = 0, p, v, dim, d, off; 969 PetscErrorCode ierr; 970 971 PetscFunctionBegin; 972 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 973 for (d = 0; d <= depth; ++d) pEnd += numPoints[d]; 974 ierr = DMPlexSetChart(dm, pStart, pEnd);CHKERRQ(ierr); 975 for (p = pStart; p < pEnd; ++p) { 976 ierr = DMPlexSetConeSize(dm, p, coneSize[p-pStart]);CHKERRQ(ierr); 977 } 978 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 979 for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) { 980 ierr = DMPlexSetCone(dm, p, &cones[off]);CHKERRQ(ierr); 981 ierr = DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);CHKERRQ(ierr); 982 } 983 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 984 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 985 /* Build coordinates */ 986 ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 987 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 988 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 989 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0]);CHKERRQ(ierr); 990 for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) { 991 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 992 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 993 } 994 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 995 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 996 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 997 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 998 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 999 ierr = VecSetFromOptions(coordinates);CHKERRQ(ierr); 1000 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1001 for (v = 0; v < numPoints[0]; ++v) { 1002 PetscInt off; 1003 1004 ierr = PetscSectionGetOffset(coordSection, v+firstVertex, &off);CHKERRQ(ierr); 1005 for (d = 0; d < dim; ++d) { 1006 coords[off+d] = vertexCoords[v*dim+d]; 1007 } 1008 } 1009 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1010 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 1011 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1012 PetscFunctionReturn(0); 1013 } 1014