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 = DMGetCoordinateSection(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 = VecSetType(coordinates,VECSTANDARD);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 = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 324 if (!rank) { 325 PetscInt f; 326 327 ierr = DMPlexSetChart(dm, 0, numFaces+numVertices);CHKERRQ(ierr); 328 for (f = 0; f < numFaces; ++f) { 329 ierr = DMPlexSetConeSize(dm, f, 4);CHKERRQ(ierr); 330 } 331 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 332 for (v = 0; v < numFaces+numVertices; ++v) { 333 ierr = DMPlexSetLabelValue(dm, "marker", v, 1);CHKERRQ(ierr); 334 } 335 { /* Side 0 (Front) */ 336 PetscInt cone[4]; 337 cone[0] = numFaces+4; cone[1] = numFaces+5; cone[2] = numFaces+7; cone[3] = numFaces+6; 338 ierr = DMPlexSetCone(dm, 0, cone);CHKERRQ(ierr); 339 } 340 { /* Side 1 (Back) */ 341 PetscInt cone[4]; 342 cone[0] = numFaces+1; cone[1] = numFaces+0; cone[2] = numFaces+2; cone[3] = numFaces+3; 343 ierr = DMPlexSetCone(dm, 1, cone);CHKERRQ(ierr); 344 } 345 { /* Side 2 (Bottom) */ 346 PetscInt cone[4]; 347 cone[0] = numFaces+0; cone[1] = numFaces+1; cone[2] = numFaces+5; cone[3] = numFaces+4; 348 ierr = DMPlexSetCone(dm, 2, cone);CHKERRQ(ierr); 349 } 350 { /* Side 3 (Top) */ 351 PetscInt cone[4]; 352 cone[0] = numFaces+6; cone[1] = numFaces+7; cone[2] = numFaces+3; cone[3] = numFaces+2; 353 ierr = DMPlexSetCone(dm, 3, cone);CHKERRQ(ierr); 354 } 355 { /* Side 4 (Left) */ 356 PetscInt cone[4]; 357 cone[0] = numFaces+0; cone[1] = numFaces+4; cone[2] = numFaces+6; cone[3] = numFaces+2; 358 ierr = DMPlexSetCone(dm, 4, cone);CHKERRQ(ierr); 359 } 360 { /* Side 5 (Right) */ 361 PetscInt cone[4]; 362 cone[0] = numFaces+5; cone[1] = numFaces+1; cone[2] = numFaces+3; cone[3] = numFaces+7; 363 ierr = DMPlexSetCone(dm, 5, cone);CHKERRQ(ierr); 364 } 365 } 366 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 367 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 368 /* Build coordinates */ 369 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 370 ierr = PetscSectionSetChart(coordSection, numFaces, numFaces + numVertices);CHKERRQ(ierr); 371 for (v = numFaces; v < numFaces+numVertices; ++v) { 372 ierr = PetscSectionSetDof(coordSection, v, 3);CHKERRQ(ierr); 373 } 374 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 375 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 376 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 377 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 378 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 379 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 380 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 381 for (vz = 0; vz <= faces[2]; ++vz) { 382 for (vy = 0; vy <= faces[1]; ++vy) { 383 for (vx = 0; vx <= faces[0]; ++vx) { 384 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+0] = lower[0] + ((upper[0] - lower[0])/faces[0])*vx; 385 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+1] = lower[1] + ((upper[1] - lower[1])/faces[1])*vy; 386 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+2] = lower[2] + ((upper[2] - lower[2])/faces[2])*vz; 387 } 388 } 389 } 390 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 391 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 392 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 393 PetscFunctionReturn(0); 394 } 395 396 #undef __FUNCT__ 397 #define __FUNCT__ "DMPlexCreateSquareMesh" 398 /* 399 Simple square mesh: 400 401 22--8-23--9--24 402 | | | 403 13 2 14 3 15 404 | | | 405 19--6-20--7--21 406 | | | 407 10 0 11 1 12 408 | | | 409 16--4-17--5--18 410 */ 411 PetscErrorCode DMPlexCreateSquareMesh(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[]) 412 { 413 PetscInt markerTop = 1; 414 PetscInt markerBottom = 1; 415 PetscInt markerRight = 1; 416 PetscInt markerLeft = 1; 417 PetscBool markerSeparate = PETSC_FALSE; 418 PetscMPIInt rank; 419 PetscErrorCode ierr; 420 421 PetscFunctionBegin; 422 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 423 ierr = PetscOptionsGetBool(((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);CHKERRQ(ierr); 424 if (markerSeparate) { 425 markerTop = 3; 426 markerBottom = 1; 427 markerRight = 2; 428 markerLeft = 4; 429 } 430 { 431 const PetscInt numXEdges = !rank ? edges[0] : 0; 432 const PetscInt numYEdges = !rank ? edges[1] : 0; 433 const PetscInt numXVertices = !rank ? edges[0]+1 : 0; 434 const PetscInt numYVertices = !rank ? edges[1]+1 : 0; 435 const PetscInt numTotXEdges = numXEdges*numYVertices; 436 const PetscInt numTotYEdges = numYEdges*numXVertices; 437 const PetscInt numVertices = numXVertices*numYVertices; 438 const PetscInt numEdges = numTotXEdges + numTotYEdges; 439 const PetscInt numFaces = numXEdges*numYEdges; 440 const PetscInt firstVertex = numFaces; 441 const PetscInt firstXEdge = numFaces + numVertices; 442 const PetscInt firstYEdge = numFaces + numVertices + numTotXEdges; 443 Vec coordinates; 444 PetscSection coordSection; 445 PetscScalar *coords; 446 PetscInt coordSize; 447 PetscInt v, vx, vy; 448 PetscInt f, fx, fy, e, ex, ey; 449 450 ierr = DMPlexSetChart(dm, 0, numFaces+numEdges+numVertices);CHKERRQ(ierr); 451 for (f = 0; f < numFaces; ++f) { 452 ierr = DMPlexSetConeSize(dm, f, 4);CHKERRQ(ierr); 453 } 454 for (e = firstXEdge; e < firstXEdge+numEdges; ++e) { 455 ierr = DMPlexSetConeSize(dm, e, 2);CHKERRQ(ierr); 456 } 457 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 458 /* Build faces */ 459 for (fy = 0; fy < numYEdges; fy++) { 460 for (fx = 0; fx < numXEdges; fx++) { 461 const PetscInt face = fy*numXEdges + fx; 462 const PetscInt edgeL = firstYEdge + fx*numYEdges + fy; 463 const PetscInt edgeB = firstXEdge + fy*numXEdges + fx; 464 const PetscInt ornt[4] = {0, 0, -2, -2}; 465 PetscInt cone[4]; 466 467 cone[0] = edgeB; cone[1] = edgeL+numYEdges; cone[2] = edgeB+numXEdges; cone[3] = edgeL; 468 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr); 469 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr); 470 } 471 } 472 /* Build Y edges*/ 473 for (vx = 0; vx < numXVertices; vx++) { 474 for (ey = 0; ey < numYEdges; ey++) { 475 const PetscInt edge = firstYEdge + vx*numYEdges + ey; 476 const PetscInt vertex = firstVertex + ey*numXVertices + vx; 477 PetscInt cone[2]; 478 479 cone[0] = vertex; cone[1] = vertex+numXVertices; 480 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 481 if (vx == numXVertices-1) { 482 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 483 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerRight);CHKERRQ(ierr); 484 if (ey == numYEdges-1) { 485 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerRight);CHKERRQ(ierr); 486 } 487 } else if (vx == 0) { 488 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 489 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerLeft);CHKERRQ(ierr); 490 if (ey == numYEdges-1) { 491 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerLeft);CHKERRQ(ierr); 492 } 493 } 494 } 495 } 496 /* Build X edges*/ 497 for (vy = 0; vy < numYVertices; vy++) { 498 for (ex = 0; ex < numXEdges; ex++) { 499 const PetscInt edge = firstXEdge + vy*numXEdges + ex; 500 const PetscInt vertex = firstVertex + vy*numXVertices + ex; 501 PetscInt cone[2]; 502 503 cone[0] = vertex; cone[1] = vertex+1; 504 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 505 if (vy == numYVertices-1) { 506 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 507 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerTop);CHKERRQ(ierr); 508 if (ex == numXEdges-1) { 509 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerTop);CHKERRQ(ierr); 510 } 511 } else if (vy == 0) { 512 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 513 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerBottom);CHKERRQ(ierr); 514 if (ex == numXEdges-1) { 515 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerBottom);CHKERRQ(ierr); 516 } 517 } 518 } 519 } 520 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 521 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 522 /* Build coordinates */ 523 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 524 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVertices);CHKERRQ(ierr); 525 for (v = firstVertex; v < firstVertex+numVertices; ++v) { 526 ierr = PetscSectionSetDof(coordSection, v, 2);CHKERRQ(ierr); 527 } 528 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 529 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 530 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 531 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 532 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 533 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 534 for (vy = 0; vy < numYVertices; ++vy) { 535 for (vx = 0; vx < numXVertices; ++vx) { 536 coords[(vy*numXVertices+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/numXEdges)*vx; 537 coords[(vy*numXVertices+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/numYEdges)*vy; 538 } 539 } 540 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 541 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 542 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 543 } 544 PetscFunctionReturn(0); 545 } 546 547 #undef __FUNCT__ 548 #define __FUNCT__ "DMPlexCreateBoxMesh" 549 /*@ 550 DMPlexCreateBoxMesh - Creates a mesh on the tensor product of unit intervals (box). 551 552 Collective on MPI_Comm 553 554 Input Parameters: 555 + comm - The communicator for the DM object 556 . dim - The spatial dimension 557 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 558 559 Output Parameter: 560 . dm - The DM object 561 562 Level: beginner 563 564 .keywords: DM, create 565 .seealso: DMSetType(), DMCreate() 566 @*/ 567 PetscErrorCode DMPlexCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscBool interpolate, DM *dm) 568 { 569 DM boundary; 570 PetscErrorCode ierr; 571 572 PetscFunctionBegin; 573 PetscValidPointer(dm, 4); 574 ierr = DMCreate(comm, &boundary);CHKERRQ(ierr); 575 PetscValidLogicalCollectiveInt(boundary,dim,2); 576 ierr = DMSetType(boundary, DMPLEX);CHKERRQ(ierr); 577 ierr = DMPlexSetDimension(boundary, dim-1);CHKERRQ(ierr); 578 switch (dim) { 579 case 2: 580 { 581 PetscReal lower[2] = {0.0, 0.0}; 582 PetscReal upper[2] = {1.0, 1.0}; 583 PetscInt edges[2] = {2, 2}; 584 585 ierr = DMPlexCreateSquareBoundary(boundary, lower, upper, edges);CHKERRQ(ierr); 586 break; 587 } 588 case 3: 589 { 590 PetscReal lower[3] = {0.0, 0.0, 0.0}; 591 PetscReal upper[3] = {1.0, 1.0, 1.0}; 592 PetscInt faces[3] = {1, 1, 1}; 593 594 ierr = DMPlexCreateCubeBoundary(boundary, lower, upper, faces);CHKERRQ(ierr); 595 break; 596 } 597 default: 598 SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim); 599 } 600 ierr = DMPlexGenerate(boundary, NULL, interpolate, dm);CHKERRQ(ierr); 601 ierr = DMDestroy(&boundary);CHKERRQ(ierr); 602 PetscFunctionReturn(0); 603 } 604 605 #undef __FUNCT__ 606 #define __FUNCT__ "DMPlexCreateHexBoxMesh" 607 PetscErrorCode DMPlexCreateHexBoxMesh(MPI_Comm comm, PetscInt dim, const PetscInt cells[], DM *dm) 608 { 609 PetscErrorCode ierr; 610 611 PetscFunctionBegin; 612 PetscValidPointer(dm, 4); 613 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 614 PetscValidLogicalCollectiveInt(*dm,dim,2); 615 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 616 ierr = DMPlexSetDimension(*dm, dim);CHKERRQ(ierr); 617 switch (dim) { 618 case 2: 619 { 620 PetscReal lower[2] = {0.0, 0.0}; 621 PetscReal upper[2] = {1.0, 1.0}; 622 623 ierr = DMPlexCreateSquareMesh(*dm, lower, upper, cells);CHKERRQ(ierr); 624 break; 625 } 626 #if 0 627 case 3: 628 { 629 PetscReal lower[3] = {0.0, 0.0, 0.0}; 630 PetscReal upper[3] = {1.0, 1.0, 1.0}; 631 632 ierr = DMPlexCreateCubeMesh(boundary, lower, upper, cells);CHKERRQ(ierr); 633 break; 634 } 635 #endif 636 default: 637 SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim); 638 } 639 PetscFunctionReturn(0); 640 } 641 642 /* External function declarations here */ 643 extern PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling); 644 extern PetscErrorCode DMCreateMatrix_Plex(DM dm, Mat *J); 645 extern PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm); 646 extern PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *dmRefined); 647 extern PetscErrorCode DMClone_Plex(DM dm, DM *newdm); 648 extern PetscErrorCode DMSetUp_Plex(DM dm); 649 extern PetscErrorCode DMDestroy_Plex(DM dm); 650 extern PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer); 651 extern PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm); 652 extern PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, IS *cellIS); 653 654 #undef __FUNCT__ 655 #define __FUNCT__ "DMCreateGlobalVector_Plex" 656 static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec) 657 { 658 PetscErrorCode ierr; 659 660 PetscFunctionBegin; 661 ierr = DMCreateGlobalVector_Section_Private(dm,vec);CHKERRQ(ierr); 662 /* ierr = VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM);CHKERRQ(ierr); */ 663 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void))VecView_Plex);CHKERRQ(ierr); 664 PetscFunctionReturn(0); 665 } 666 667 #undef __FUNCT__ 668 #define __FUNCT__ "DMCreateLocalVector_Plex" 669 static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec) 670 { 671 PetscErrorCode ierr; 672 673 PetscFunctionBegin; 674 ierr = DMCreateLocalVector_Section_Private(dm,vec);CHKERRQ(ierr); 675 ierr = VecSetOperation(*vec, VECOP_VIEW, (void(*)(void)) VecView_Plex_Local);CHKERRQ(ierr); 676 PetscFunctionReturn(0); 677 } 678 679 #undef __FUNCT__ 680 #define __FUNCT__ "DMInitialize_Plex" 681 PetscErrorCode DMInitialize_Plex(DM dm) 682 { 683 PetscFunctionBegin; 684 dm->ops->view = DMView_Plex; 685 dm->ops->setfromoptions = DMSetFromOptions_Plex; 686 dm->ops->clone = DMClone_Plex; 687 dm->ops->setup = DMSetUp_Plex; 688 dm->ops->createglobalvector = DMCreateGlobalVector_Plex; 689 dm->ops->createlocalvector = DMCreateLocalVector_Plex; 690 dm->ops->getlocaltoglobalmapping = NULL; 691 dm->ops->getlocaltoglobalmappingblock = NULL; 692 dm->ops->createfieldis = NULL; 693 dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex; 694 dm->ops->getcoloring = 0; 695 dm->ops->creatematrix = DMCreateMatrix_Plex; 696 dm->ops->createinterpolation = 0; 697 dm->ops->getaggregates = 0; 698 dm->ops->getinjection = 0; 699 dm->ops->refine = DMRefine_Plex; 700 dm->ops->coarsen = 0; 701 dm->ops->refinehierarchy = 0; 702 dm->ops->coarsenhierarchy = 0; 703 dm->ops->globaltolocalbegin = NULL; 704 dm->ops->globaltolocalend = NULL; 705 dm->ops->localtoglobalbegin = NULL; 706 dm->ops->localtoglobalend = NULL; 707 dm->ops->destroy = DMDestroy_Plex; 708 dm->ops->createsubdm = DMCreateSubDM_Plex; 709 dm->ops->locatepoints = DMLocatePoints_Plex; 710 PetscFunctionReturn(0); 711 } 712 713 #undef __FUNCT__ 714 #define __FUNCT__ "DMClone_Plex" 715 PetscErrorCode DMClone_Plex(DM dm, DM *newdm) 716 { 717 DM_Plex *mesh = (DM_Plex *) dm->data; 718 PetscErrorCode ierr; 719 720 PetscFunctionBegin; 721 mesh->refct++; 722 (*newdm)->data = mesh; 723 ierr = PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);CHKERRQ(ierr); 724 ierr = DMInitialize_Plex(*newdm);CHKERRQ(ierr); 725 PetscFunctionReturn(0); 726 } 727 728 /*MC 729 DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram. 730 In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is 731 specified by a PetscSection object. Ownership in the global representation is determined by 732 ownership of the underlying DMPlex points. This is specified by another PetscSection object. 733 734 Level: intermediate 735 736 .seealso: DMType, DMPlexCreate(), DMCreate(), DMSetType() 737 M*/ 738 739 #undef __FUNCT__ 740 #define __FUNCT__ "DMCreate_Plex" 741 PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm) 742 { 743 DM_Plex *mesh; 744 PetscInt unit, d; 745 PetscErrorCode ierr; 746 747 PetscFunctionBegin; 748 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 749 ierr = PetscNewLog(dm,&mesh);CHKERRQ(ierr); 750 dm->data = mesh; 751 752 mesh->refct = 1; 753 mesh->dim = 0; 754 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);CHKERRQ(ierr); 755 mesh->maxConeSize = 0; 756 mesh->cones = NULL; 757 mesh->coneOrientations = NULL; 758 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);CHKERRQ(ierr); 759 mesh->maxSupportSize = 0; 760 mesh->supports = NULL; 761 mesh->refinementUniform = PETSC_TRUE; 762 mesh->refinementLimit = -1.0; 763 764 mesh->facesTmp = NULL; 765 766 mesh->subpointMap = NULL; 767 768 for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0; 769 770 mesh->labels = NULL; 771 mesh->depthLabel = NULL; 772 mesh->globalVertexNumbers = NULL; 773 mesh->globalCellNumbers = NULL; 774 for (d = 0; d < 8; ++d) mesh->hybridPointMax[d] = PETSC_DETERMINE; 775 mesh->vtkCellHeight = 0; 776 mesh->preallocCenterDim = -1; 777 778 mesh->printSetValues = PETSC_FALSE; 779 mesh->printFEM = 0; 780 mesh->printTol = 1.0e-10; 781 782 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 783 PetscFunctionReturn(0); 784 } 785 786 #undef __FUNCT__ 787 #define __FUNCT__ "DMPlexCreate" 788 /*@ 789 DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram. 790 791 Collective on MPI_Comm 792 793 Input Parameter: 794 . comm - The communicator for the DMPlex object 795 796 Output Parameter: 797 . mesh - The DMPlex object 798 799 Level: beginner 800 801 .keywords: DMPlex, create 802 @*/ 803 PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh) 804 { 805 PetscErrorCode ierr; 806 807 PetscFunctionBegin; 808 PetscValidPointer(mesh,2); 809 ierr = DMCreate(comm, mesh);CHKERRQ(ierr); 810 ierr = DMSetType(*mesh, DMPLEX);CHKERRQ(ierr); 811 PetscFunctionReturn(0); 812 } 813 814 #undef __FUNCT__ 815 #define __FUNCT__ "DMPlexBuildFromCellList_Private" 816 /* 817 This takes as input the common mesh generator output, a list of the vertices for each cell 818 */ 819 PetscErrorCode DMPlexBuildFromCellList_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[]) 820 { 821 PetscInt *cone, c, p; 822 PetscErrorCode ierr; 823 824 PetscFunctionBegin; 825 ierr = DMPlexSetChart(dm, 0, numCells+numVertices);CHKERRQ(ierr); 826 for (c = 0; c < numCells; ++c) { 827 ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr); 828 } 829 ierr = DMSetUp(dm);CHKERRQ(ierr); 830 ierr = DMGetWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 831 for (c = 0; c < numCells; ++c) { 832 for (p = 0; p < numCorners; ++p) { 833 cone[p] = cells[c*numCorners+p]+numCells; 834 } 835 ierr = DMPlexSetCone(dm, c, cone);CHKERRQ(ierr); 836 } 837 ierr = DMRestoreWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 838 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 839 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 840 PetscFunctionReturn(0); 841 } 842 843 #undef __FUNCT__ 844 #define __FUNCT__ "DMPlexBuildCoordinates_Private" 845 /* 846 This takes as input the coordinates for each vertex 847 */ 848 PetscErrorCode DMPlexBuildCoordinates_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, const double vertexCoords[]) 849 { 850 PetscSection coordSection; 851 Vec coordinates; 852 PetscScalar *coords; 853 PetscInt coordSize, v, d; 854 PetscErrorCode ierr; 855 856 PetscFunctionBegin; 857 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 858 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 859 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 860 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVertices);CHKERRQ(ierr); 861 for (v = numCells; v < numCells+numVertices; ++v) { 862 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 863 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 864 } 865 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 866 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 867 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 868 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 869 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 870 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 871 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 872 for (v = 0; v < numVertices; ++v) { 873 for (d = 0; d < spaceDim; ++d) { 874 coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d]; 875 } 876 } 877 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 878 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 879 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 880 PetscFunctionReturn(0); 881 } 882 883 #undef __FUNCT__ 884 #define __FUNCT__ "DMPlexCreateFromCellList" 885 /*@C 886 DMPlexCreateFromCellList - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM 887 888 Input Parameters: 889 + comm - The communicator 890 . dim - The topological dimension of the mesh 891 . numCells - The number of cells 892 . numVertices - The number of vertices 893 . numCorners - The number of vertices for each cell 894 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 895 . cells - An array of numCells*numCorners numbers, the vertices for each cell 896 . spaceDim - The spatial dimension used for coordinates 897 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 898 899 Output Parameter: 900 . dm - The DM 901 902 Note: Two triangles sharing a face 903 $ 904 $ 2 905 $ / | \ 906 $ / | \ 907 $ / | \ 908 $ 0 0 | 1 3 909 $ \ | / 910 $ \ | / 911 $ \ | / 912 $ 1 913 would have input 914 $ numCells = 2, numVertices = 4 915 $ cells = [0 1 2 1 3 2] 916 $ 917 which would result in the DMPlex 918 $ 919 $ 4 920 $ / | \ 921 $ / | \ 922 $ / | \ 923 $ 2 0 | 1 5 924 $ \ | / 925 $ \ | / 926 $ \ | / 927 $ 3 928 929 Level: beginner 930 931 .seealso: DMPlexCreateFromDAG(), DMPlexCreate() 932 @*/ 933 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) 934 { 935 PetscErrorCode ierr; 936 937 PetscFunctionBegin; 938 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 939 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 940 ierr = DMPlexSetDimension(*dm, dim);CHKERRQ(ierr); 941 ierr = DMPlexBuildFromCellList_Private(*dm, numCells, numVertices, numCorners, cells);CHKERRQ(ierr); 942 if (interpolate) { 943 DM idm; 944 945 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 946 ierr = DMDestroy(dm);CHKERRQ(ierr); 947 *dm = idm; 948 } 949 ierr = DMPlexBuildCoordinates_Private(*dm, spaceDim, numCells, numVertices, vertexCoords);CHKERRQ(ierr); 950 PetscFunctionReturn(0); 951 } 952 953 #undef __FUNCT__ 954 #define __FUNCT__ "DMPlexCreateFromDAG" 955 /*@ 956 DMPlexCreateFromDAG - This takes as input the adjacency-list representation of the Directed Acyclic Graph (Hasse Diagram) encoding a mesh, and produces a DM 957 958 Input Parameters: 959 + dm - The empty DM object, usually from DMCreate() and DMPlexSetDimension() 960 . depth - The depth of the DAG 961 . numPoints - The number of points at each depth 962 . coneSize - The cone size of each point 963 . cones - The concatenation of the cone points for each point, the cone list must be oriented correctly for each point 964 . coneOrientations - The orientation of each cone point 965 - vertexCoords - An array of numVertices*dim numbers, the coordinates of each vertex 966 967 Output Parameter: 968 . dm - The DM 969 970 Note: Two triangles sharing a face would have input 971 $ depth = 1, numPoints = [4 2], coneSize = [3 3 0 0 0 0] 972 $ cones = [2 3 4 3 5 4], coneOrientations = [0 0 0 0 0 0] 973 $ vertexCoords = [-1.0 0.0 0.0 -1.0 0.0 1.0 1.0 0.0] 974 $ 975 which would result in the DMPlex 976 $ 977 $ 4 978 $ / | \ 979 $ / | \ 980 $ / | \ 981 $ 2 0 | 1 5 982 $ \ | / 983 $ \ | / 984 $ \ | / 985 $ 3 986 $ 987 $ Notice that all points are numbered consecutively, unlikely DMPlexCreateFromCellList() 988 989 Level: advanced 990 991 .seealso: DMPlexCreateFromCellList(), DMPlexCreate() 992 @*/ 993 PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[]) 994 { 995 Vec coordinates; 996 PetscSection coordSection; 997 PetscScalar *coords; 998 PetscInt coordSize, firstVertex = numPoints[depth], pStart = 0, pEnd = 0, p, v, dim, d, off; 999 PetscErrorCode ierr; 1000 1001 PetscFunctionBegin; 1002 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 1003 for (d = 0; d <= depth; ++d) pEnd += numPoints[d]; 1004 ierr = DMPlexSetChart(dm, pStart, pEnd);CHKERRQ(ierr); 1005 for (p = pStart; p < pEnd; ++p) { 1006 ierr = DMPlexSetConeSize(dm, p, coneSize[p-pStart]);CHKERRQ(ierr); 1007 } 1008 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 1009 for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) { 1010 ierr = DMPlexSetCone(dm, p, &cones[off]);CHKERRQ(ierr); 1011 ierr = DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);CHKERRQ(ierr); 1012 } 1013 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 1014 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 1015 /* Build coordinates */ 1016 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1017 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1018 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 1019 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0]);CHKERRQ(ierr); 1020 for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) { 1021 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 1022 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 1023 } 1024 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1025 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1026 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 1027 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1028 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1029 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1030 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1031 for (v = 0; v < numPoints[0]; ++v) { 1032 PetscInt off; 1033 1034 ierr = PetscSectionGetOffset(coordSection, v+firstVertex, &off);CHKERRQ(ierr); 1035 for (d = 0; d < dim; ++d) { 1036 coords[off+d] = vertexCoords[v*dim+d]; 1037 } 1038 } 1039 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1040 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 1041 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1042 PetscFunctionReturn(0); 1043 } 1044