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