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