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__ "DMPlexCreateDoublet" 8 /*@ 9 DMPlexCreateDoublet - Creates a mesh of two cells of the specified type, optionally with later refinement. 10 11 Collective on MPI_Comm 12 13 Input Parameters: 14 + comm - The communicator for the DM object 15 . dim - The spatial dimension 16 . simplex - Flag for simplicial cells, otherwise they are tensor product cells 17 . interpolate - Flag to create intermediate mesh pieces (edges, faces) 18 . refinementUniform - Flag for uniform parallel refinement 19 - refinementLimit - A nonzero number indicates the largest admissible volume for a refined cell 20 21 Output Parameter: 22 . dm - The DM object 23 24 Level: beginner 25 26 .keywords: DM, create 27 .seealso: DMSetType(), DMCreate() 28 @*/ 29 PetscErrorCode DMPlexCreateDoublet(MPI_Comm comm, PetscInt dim, PetscBool simplex, PetscBool interpolate, PetscBool refinementUniform, PetscReal refinementLimit, DM *newdm) 30 { 31 DM dm; 32 PetscInt p; 33 PetscMPIInt rank; 34 PetscErrorCode ierr; 35 36 PetscFunctionBegin; 37 ierr = DMCreate(comm, &dm);CHKERRQ(ierr); 38 ierr = DMSetType(dm, DMPLEX);CHKERRQ(ierr); 39 ierr = DMSetDimension(dm, dim);CHKERRQ(ierr); 40 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 41 switch (dim) { 42 case 2: 43 if (simplex) {ierr = PetscObjectSetName((PetscObject) dm, "triangular");CHKERRQ(ierr);} 44 else {ierr = PetscObjectSetName((PetscObject) dm, "quadrilateral");CHKERRQ(ierr);} 45 break; 46 case 3: 47 if (simplex) {ierr = PetscObjectSetName((PetscObject) dm, "tetrahedral");CHKERRQ(ierr);} 48 else {ierr = PetscObjectSetName((PetscObject) dm, "hexahedral");CHKERRQ(ierr);} 49 break; 50 default: 51 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %d", dim); 52 } 53 if (rank) { 54 PetscInt numPoints[2] = {0, 0}; 55 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 56 } else { 57 switch (dim) { 58 case 2: 59 if (simplex) { 60 PetscInt numPoints[2] = {4, 2}; 61 PetscInt coneSize[6] = {3, 3, 0, 0, 0, 0}; 62 PetscInt cones[6] = {2, 3, 4, 5, 4, 3}; 63 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 64 PetscScalar vertexCoords[8] = {-0.5, 0.5, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5}; 65 PetscInt markerPoints[8] = {2, 1, 3, 1, 4, 1, 5, 1}; 66 67 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 68 for (p = 0; p < 4; ++p) {ierr = DMPlexSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 69 } else { 70 PetscInt numPoints[2] = {6, 2}; 71 PetscInt coneSize[8] = {4, 4, 0, 0, 0, 0, 0, 0}; 72 PetscInt cones[8] = {2, 3, 4, 5, 3, 6, 7, 4}; 73 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 74 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}; 75 76 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 77 } 78 break; 79 case 3: 80 if (simplex) { 81 PetscInt numPoints[2] = {5, 2}; 82 PetscInt coneSize[7] = {4, 4, 0, 0, 0, 0, 0}; 83 PetscInt cones[8] = {4, 3, 5, 2, 5, 3, 4, 6}; 84 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 85 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}; 86 PetscInt markerPoints[10] = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1}; 87 88 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 89 for (p = 0; p < 5; ++p) {ierr = DMPlexSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 90 } else { 91 PetscInt numPoints[2] = {12, 2}; 92 PetscInt coneSize[14] = {8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 93 PetscInt cones[16] = {2, 3, 4, 5, 6, 7, 8, 9, 5, 4, 10, 11, 7, 12, 13, 8}; 94 PetscInt coneOrientations[16] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 95 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, 96 -1.0, -0.5, 0.5, 0.0, -0.5, 0.5, 0.0, 0.5, 0.5, -1.0, 0.5, 0.5, 97 1.0, 0.5, -0.5, 1.0, -0.5, -0.5, 1.0, -0.5, 0.5, 1.0, 0.5, 0.5}; 98 99 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 100 } 101 break; 102 default: 103 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %d", dim); 104 } 105 } 106 *newdm = dm; 107 if (refinementLimit > 0.0) { 108 DM rdm; 109 const char *name; 110 111 ierr = DMPlexSetRefinementUniform(*newdm, PETSC_FALSE);CHKERRQ(ierr); 112 ierr = DMPlexSetRefinementLimit(*newdm, refinementLimit);CHKERRQ(ierr); 113 ierr = DMRefine(*newdm, comm, &rdm);CHKERRQ(ierr); 114 ierr = PetscObjectGetName((PetscObject) *newdm, &name);CHKERRQ(ierr); 115 ierr = PetscObjectSetName((PetscObject) rdm, name);CHKERRQ(ierr); 116 ierr = DMDestroy(newdm);CHKERRQ(ierr); 117 *newdm = rdm; 118 } 119 if (interpolate) { 120 DM idm = NULL; 121 const char *name; 122 123 ierr = DMPlexInterpolate(*newdm, &idm);CHKERRQ(ierr); 124 ierr = PetscObjectGetName((PetscObject) *newdm, &name);CHKERRQ(ierr); 125 ierr = PetscObjectSetName((PetscObject) idm, name);CHKERRQ(ierr); 126 ierr = DMPlexCopyCoordinates(*newdm, idm);CHKERRQ(ierr); 127 ierr = DMPlexCopyLabels(*newdm, idm);CHKERRQ(ierr); 128 ierr = DMDestroy(newdm);CHKERRQ(ierr); 129 *newdm = idm; 130 } 131 { 132 DM refinedMesh = NULL; 133 DM distributedMesh = NULL; 134 135 /* Distribute mesh over processes */ 136 ierr = DMPlexDistribute(*newdm, 0, NULL, &distributedMesh);CHKERRQ(ierr); 137 if (distributedMesh) { 138 ierr = DMDestroy(newdm);CHKERRQ(ierr); 139 *newdm = distributedMesh; 140 } 141 if (refinementUniform) { 142 ierr = DMPlexSetRefinementUniform(*newdm, refinementUniform);CHKERRQ(ierr); 143 ierr = DMRefine(*newdm, comm, &refinedMesh);CHKERRQ(ierr); 144 if (refinedMesh) { 145 ierr = DMDestroy(newdm);CHKERRQ(ierr); 146 *newdm = refinedMesh; 147 } 148 } 149 } 150 PetscFunctionReturn(0); 151 } 152 153 #undef __FUNCT__ 154 #define __FUNCT__ "DMPlexCreateSquareBoundary" 155 /*@ 156 DMPlexCreateSquareBoundary - Creates a 1D mesh the is the boundary of a square lattice. 157 158 Collective on MPI_Comm 159 160 Input Parameters: 161 + comm - The communicator for the DM object 162 . lower - The lower left corner coordinates 163 . upper - The upper right corner coordinates 164 - edges - The number of cells in each direction 165 166 Output Parameter: 167 . dm - The DM object 168 169 Note: Here is the numbering returned for 2 cells in each direction: 170 $ 18--5-17--4--16 171 $ | | | 172 $ 6 10 3 173 $ | | | 174 $ 19-11-20--9--15 175 $ | | | 176 $ 7 8 2 177 $ | | | 178 $ 12--0-13--1--14 179 180 Level: beginner 181 182 .keywords: DM, create 183 .seealso: DMPlexCreateBoxMesh(), DMPlexCreateCubeBoundary(), DMSetType(), DMCreate() 184 @*/ 185 PetscErrorCode DMPlexCreateSquareBoundary(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[]) 186 { 187 PetscInt numVertices = (edges[0]+1)*(edges[1]+1); 188 PetscInt numEdges = edges[0]*(edges[1]+1) + (edges[0]+1)*edges[1]; 189 PetscInt markerTop = 1; 190 PetscInt markerBottom = 1; 191 PetscInt markerRight = 1; 192 PetscInt markerLeft = 1; 193 PetscBool markerSeparate = PETSC_FALSE; 194 Vec coordinates; 195 PetscSection coordSection; 196 PetscScalar *coords; 197 PetscInt coordSize; 198 PetscMPIInt rank; 199 PetscInt v, vx, vy; 200 PetscErrorCode ierr; 201 202 PetscFunctionBegin; 203 ierr = PetscOptionsGetBool(((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);CHKERRQ(ierr); 204 if (markerSeparate) { 205 markerTop = 1; 206 markerBottom = 0; 207 markerRight = 0; 208 markerLeft = 0; 209 } 210 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 211 if (!rank) { 212 PetscInt e, ex, ey; 213 214 ierr = DMPlexSetChart(dm, 0, numEdges+numVertices);CHKERRQ(ierr); 215 for (e = 0; e < numEdges; ++e) { 216 ierr = DMPlexSetConeSize(dm, e, 2);CHKERRQ(ierr); 217 } 218 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 219 for (vx = 0; vx <= edges[0]; vx++) { 220 for (ey = 0; ey < edges[1]; ey++) { 221 PetscInt edge = vx*edges[1] + ey + edges[0]*(edges[1]+1); 222 PetscInt vertex = ey*(edges[0]+1) + vx + numEdges; 223 PetscInt cone[2]; 224 225 cone[0] = vertex; cone[1] = vertex+edges[0]+1; 226 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 227 if (vx == edges[0]) { 228 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 229 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerRight);CHKERRQ(ierr); 230 if (ey == edges[1]-1) { 231 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerRight);CHKERRQ(ierr); 232 } 233 } else if (vx == 0) { 234 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 235 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerLeft);CHKERRQ(ierr); 236 if (ey == edges[1]-1) { 237 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerLeft);CHKERRQ(ierr); 238 } 239 } 240 } 241 } 242 for (vy = 0; vy <= edges[1]; vy++) { 243 for (ex = 0; ex < edges[0]; ex++) { 244 PetscInt edge = vy*edges[0] + ex; 245 PetscInt vertex = vy*(edges[0]+1) + ex + numEdges; 246 PetscInt cone[2]; 247 248 cone[0] = vertex; cone[1] = vertex+1; 249 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 250 if (vy == edges[1]) { 251 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 252 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerTop);CHKERRQ(ierr); 253 if (ex == edges[0]-1) { 254 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerTop);CHKERRQ(ierr); 255 } 256 } else if (vy == 0) { 257 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 258 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerBottom);CHKERRQ(ierr); 259 if (ex == edges[0]-1) { 260 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerBottom);CHKERRQ(ierr); 261 } 262 } 263 } 264 } 265 } 266 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 267 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 268 /* Build coordinates */ 269 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 270 ierr = PetscSectionSetChart(coordSection, numEdges, numEdges + numVertices);CHKERRQ(ierr); 271 for (v = numEdges; v < numEdges+numVertices; ++v) { 272 ierr = PetscSectionSetDof(coordSection, v, 2);CHKERRQ(ierr); 273 } 274 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 275 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 276 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 277 ierr = VecSetBlockSize(coordinates, 2);CHKERRQ(ierr); 278 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 279 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 280 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 281 for (vy = 0; vy <= edges[1]; ++vy) { 282 for (vx = 0; vx <= edges[0]; ++vx) { 283 coords[(vy*(edges[0]+1)+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/edges[0])*vx; 284 coords[(vy*(edges[0]+1)+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/edges[1])*vy; 285 } 286 } 287 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 288 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 289 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 290 PetscFunctionReturn(0); 291 } 292 293 #undef __FUNCT__ 294 #define __FUNCT__ "DMPlexCreateCubeBoundary" 295 /*@ 296 DMPlexCreateCubeBoundary - Creates a 2D mesh the is the boundary of a cubic lattice. 297 298 Collective on MPI_Comm 299 300 Input Parameters: 301 + comm - The communicator for the DM object 302 . lower - The lower left front corner coordinates 303 . upper - The upper right back corner coordinates 304 - edges - The number of cells in each direction 305 306 Output Parameter: 307 . dm - The DM object 308 309 Level: beginner 310 311 .keywords: DM, create 312 .seealso: DMPlexCreateBoxMesh(), DMPlexCreateSquareBoundary(), DMSetType(), DMCreate() 313 @*/ 314 PetscErrorCode DMPlexCreateCubeBoundary(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt faces[]) 315 { 316 PetscInt vertices[3], numVertices; 317 PetscInt numFaces = 2*faces[0]*faces[1] + 2*faces[1]*faces[2] + 2*faces[0]*faces[2]; 318 Vec coordinates; 319 PetscSection coordSection; 320 PetscScalar *coords; 321 PetscInt coordSize; 322 PetscMPIInt rank; 323 PetscInt v, vx, vy, vz; 324 PetscInt voffset, iface=0, cone[4]; 325 PetscErrorCode ierr; 326 327 PetscFunctionBegin; 328 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"); 329 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 330 vertices[0] = faces[0]+1; vertices[1] = faces[1]+1; vertices[2] = faces[2]+1; 331 numVertices = vertices[0]*vertices[1]*vertices[2]; 332 if (!rank) { 333 PetscInt f; 334 335 ierr = DMPlexSetChart(dm, 0, numFaces+numVertices);CHKERRQ(ierr); 336 for (f = 0; f < numFaces; ++f) { 337 ierr = DMPlexSetConeSize(dm, f, 4);CHKERRQ(ierr); 338 } 339 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 340 for (v = 0; v < numFaces+numVertices; ++v) { 341 ierr = DMPlexSetLabelValue(dm, "marker", v, 1);CHKERRQ(ierr); 342 } 343 344 /* Side 0 (Top) */ 345 for (vy = 0; vy < faces[1]; vy++) { 346 for (vx = 0; vx < faces[0]; vx++) { 347 voffset = numFaces + vertices[0]*vertices[1]*(vertices[2]-1) + vy*vertices[0] + vx; 348 cone[0] = voffset; cone[1] = voffset+1; cone[2] = voffset+vertices[0]+1; cone[3] = voffset+vertices[0]; 349 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 350 iface++; 351 } 352 } 353 354 /* Side 1 (Bottom) */ 355 for (vy = 0; vy < faces[1]; vy++) { 356 for (vx = 0; vx < faces[0]; vx++) { 357 voffset = numFaces + vy*(faces[0]+1) + vx; 358 cone[0] = voffset+1; cone[1] = voffset; cone[2] = voffset+vertices[0]; cone[3] = voffset+vertices[0]+1; 359 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 360 iface++; 361 } 362 } 363 364 /* Side 2 (Front) */ 365 for (vz = 0; vz < faces[2]; vz++) { 366 for (vx = 0; vx < faces[0]; vx++) { 367 voffset = numFaces + vz*vertices[0]*vertices[1] + vx; 368 cone[0] = voffset; cone[1] = voffset+1; cone[2] = voffset+vertices[0]*vertices[1]+1; cone[3] = voffset+vertices[0]*vertices[1]; 369 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 370 iface++; 371 } 372 } 373 374 /* Side 3 (Back) */ 375 for (vz = 0; vz < faces[2]; vz++) { 376 for (vx = 0; vx < faces[0]; vx++) { 377 voffset = numFaces + vz*vertices[0]*vertices[1] + vertices[0]*(vertices[1]-1) + vx; 378 cone[0] = voffset+vertices[0]*vertices[1]; cone[1] = voffset+vertices[0]*vertices[1]+1; 379 cone[2] = voffset+1; cone[3] = voffset; 380 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 381 iface++; 382 } 383 } 384 385 /* Side 4 (Left) */ 386 for (vz = 0; vz < faces[2]; vz++) { 387 for (vy = 0; vy < faces[1]; vy++) { 388 voffset = numFaces + vz*vertices[0]*vertices[1] + vy*vertices[0]; 389 cone[0] = voffset; cone[1] = voffset+vertices[0]*vertices[1]; 390 cone[2] = voffset+vertices[0]*vertices[1]+vertices[0]; cone[3] = voffset+vertices[0]; 391 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 392 iface++; 393 } 394 } 395 396 /* Side 5 (Right) */ 397 for (vz = 0; vz < faces[2]; vz++) { 398 for (vy = 0; vy < faces[1]; vy++) { 399 voffset = numFaces + vz*vertices[0]*vertices[1] + vy*vertices[0] + vx; 400 cone[0] = voffset+vertices[0]*vertices[1]; cone[1] = voffset; 401 cone[2] = voffset+vertices[0]; cone[3] = voffset+vertices[0]*vertices[1]+vertices[0]; 402 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 403 iface++; 404 } 405 } 406 } 407 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 408 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 409 /* Build coordinates */ 410 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 411 ierr = PetscSectionSetChart(coordSection, numFaces, numFaces + numVertices);CHKERRQ(ierr); 412 for (v = numFaces; v < numFaces+numVertices; ++v) { 413 ierr = PetscSectionSetDof(coordSection, v, 3);CHKERRQ(ierr); 414 } 415 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 416 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 417 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 418 ierr = VecSetBlockSize(coordinates, 3);CHKERRQ(ierr); 419 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 420 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 421 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 422 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 423 for (vz = 0; vz <= faces[2]; ++vz) { 424 for (vy = 0; vy <= faces[1]; ++vy) { 425 for (vx = 0; vx <= faces[0]; ++vx) { 426 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+0] = lower[0] + ((upper[0] - lower[0])/faces[0])*vx; 427 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+1] = lower[1] + ((upper[1] - lower[1])/faces[1])*vy; 428 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+2] = lower[2] + ((upper[2] - lower[2])/faces[2])*vz; 429 } 430 } 431 } 432 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 433 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 434 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 435 PetscFunctionReturn(0); 436 } 437 438 #undef __FUNCT__ 439 #define __FUNCT__ "DMPlexCreateCubeMesh_Internal" 440 static PetscErrorCode DMPlexCreateCubeMesh_Internal(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[], DMBoundaryType bdX, DMBoundaryType bdY, DMBoundaryType bdZ) 441 { 442 PetscInt markerTop = 1; 443 PetscInt markerBottom = 1; 444 PetscInt markerFront = 1; 445 PetscInt markerBack = 1; 446 PetscInt markerRight = 1; 447 PetscInt markerLeft = 1; 448 PetscInt dim; 449 PetscBool markerSeparate = PETSC_FALSE; 450 PetscMPIInt rank; 451 PetscErrorCode ierr; 452 453 PetscFunctionBegin; 454 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 455 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 456 ierr = PetscOptionsGetBool(((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);CHKERRQ(ierr); 457 if (markerSeparate) { 458 switch (dim) { 459 case 2: 460 markerTop = 3; 461 markerBottom = 1; 462 markerRight = 2; 463 markerLeft = 4; 464 break; 465 case 3: 466 markerBottom = 1; 467 markerTop = 2; 468 markerFront = 3; 469 markerBack = 4; 470 markerRight = 5; 471 markerLeft = 6; 472 break; 473 default: 474 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Dimension %d not supported",dim); 475 break; 476 } 477 } 478 { 479 const PetscInt numXEdges = !rank ? edges[0] : 0; 480 const PetscInt numYEdges = !rank ? edges[1] : 0; 481 const PetscInt numZEdges = !rank ? edges[2] : 0; 482 const PetscInt numXVertices = !rank ? (bdX == DM_BOUNDARY_PERIODIC || bdX == DM_BOUNDARY_TWIST ? edges[0] : edges[0]+1) : 0; 483 const PetscInt numYVertices = !rank ? (bdY == DM_BOUNDARY_PERIODIC || bdY == DM_BOUNDARY_TWIST ? edges[1] : edges[1]+1) : 0; 484 const PetscInt numZVertices = !rank ? (bdZ == DM_BOUNDARY_PERIODIC || bdY == DM_BOUNDARY_TWIST ? edges[2] : edges[2]+1) : 0; 485 const PetscInt numCells = numXEdges*numYEdges*numZEdges; 486 const PetscInt numXFaces = numYEdges*numZEdges; 487 const PetscInt numYFaces = numXEdges*numZEdges; 488 const PetscInt numZFaces = numXEdges*numYEdges; 489 const PetscInt numTotXFaces = numXVertices*numXFaces; 490 const PetscInt numTotYFaces = numYVertices*numYFaces; 491 const PetscInt numTotZFaces = numZVertices*numZFaces; 492 const PetscInt numFaces = numTotXFaces + numTotYFaces + numTotZFaces; 493 const PetscInt numTotXEdges = numXEdges*numYVertices*numZVertices; 494 const PetscInt numTotYEdges = numYEdges*numXVertices*numZVertices; 495 const PetscInt numTotZEdges = numZEdges*numXVertices*numYVertices; 496 const PetscInt numVertices = numXVertices*numYVertices*numZVertices; 497 const PetscInt numEdges = numTotXEdges + numTotYEdges + numTotZEdges; 498 const PetscInt firstVertex = (dim == 2) ? numFaces : numCells; 499 const PetscInt firstXFace = (dim == 2) ? 0 : numCells + numVertices; 500 const PetscInt firstYFace = firstXFace + numTotXFaces; 501 const PetscInt firstZFace = firstYFace + numTotYFaces; 502 const PetscInt firstXEdge = numCells + numFaces + numVertices; 503 const PetscInt firstYEdge = firstXEdge + numTotXEdges; 504 const PetscInt firstZEdge = firstYEdge + numTotYEdges; 505 Vec coordinates; 506 PetscSection coordSection; 507 PetscScalar *coords; 508 PetscInt coordSize; 509 PetscInt v, vx, vy, vz; 510 PetscInt c, f, fx, fy, fz, e, ex, ey, ez; 511 512 ierr = DMPlexSetChart(dm, 0, numCells+numFaces+numEdges+numVertices);CHKERRQ(ierr); 513 for (c = 0; c < numCells; c++) { 514 ierr = DMPlexSetConeSize(dm, c, 6);CHKERRQ(ierr); 515 } 516 for (f = firstXFace; f < firstXFace+numFaces; ++f) { 517 ierr = DMPlexSetConeSize(dm, f, 4);CHKERRQ(ierr); 518 } 519 for (e = firstXEdge; e < firstXEdge+numEdges; ++e) { 520 ierr = DMPlexSetConeSize(dm, e, 2);CHKERRQ(ierr); 521 } 522 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 523 /* Build cells */ 524 for (fz = 0; fz < numZEdges; ++fz) { 525 for (fy = 0; fy < numYEdges; ++fy) { 526 for (fx = 0; fx < numXEdges; ++fx) { 527 PetscInt cell = (fz*numYEdges + fy)*numXEdges + fx; 528 PetscInt faceB = firstZFace + (fy*numXEdges+fx)*numZVertices + fz; 529 PetscInt faceT = firstZFace + (fy*numXEdges+fx)*numZVertices + ((fz+1)%numZVertices); 530 PetscInt faceF = firstYFace + (fz*numXEdges+fx)*numYVertices + fy; 531 PetscInt faceK = firstYFace + (fz*numXEdges+fx)*numYVertices + ((fy+1)%numYVertices); 532 PetscInt faceL = firstXFace + (fz*numYEdges+fy)*numXVertices + fx; 533 PetscInt faceR = firstXFace + (fz*numYEdges+fy)*numXVertices + ((fx+1)%numXVertices); 534 /* B, T, F, K, R, L */ 535 PetscInt ornt[8] = {-4, 0, 0, -1, 0, -4}; /* ??? */ 536 PetscInt cone[8]; 537 538 /* no boundary twisting in 3D */ 539 cone[0] = faceB; cone[1] = faceT; cone[2] = faceF; cone[3] = faceK; cone[4] = faceR; cone[5] = faceL; 540 ierr = DMPlexSetCone(dm, cell, cone);CHKERRQ(ierr); 541 ierr = DMPlexSetConeOrientation(dm, cell, ornt);CHKERRQ(ierr); 542 } 543 } 544 } 545 /* Build x faces */ 546 for (fz = 0; fz < numZEdges; ++fz) { 547 for (fy = 0; fy < numYEdges; ++fy) { 548 for (fx = 0; fx < numXVertices; ++fx) { 549 PetscInt face = firstXFace + (fz*numYEdges+fy)*numXVertices + fx; 550 PetscInt edgeL = firstZEdge + ( fy* numXVertices+fx)*numZEdges + fz; 551 PetscInt edgeR = firstZEdge + (((fy+1)%numYVertices)*numXVertices+fx)*numZEdges + fz; 552 PetscInt edgeB = firstYEdge + ( fz* numXVertices+fx)*numYEdges + fy; 553 PetscInt edgeT = firstYEdge + (((fz+1)%numZVertices)*numXVertices+fx)*numYEdges + fy; 554 PetscInt ornt[4] = {0, 0, -2, -2}; 555 PetscInt cone[4]; 556 557 if (dim == 3) { 558 /* markers */ 559 if (bdX != DM_BOUNDARY_PERIODIC) { 560 if (fx == numXVertices-1) { 561 ierr = DMPlexSetLabelValue(dm, "marker", face, markerRight);CHKERRQ(ierr); 562 } 563 else if (fx == 0) { 564 ierr = DMPlexSetLabelValue(dm, "marker", face, markerLeft);CHKERRQ(ierr); 565 } 566 } 567 } 568 cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL; 569 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr); 570 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr); 571 } 572 } 573 } 574 /* Build y faces */ 575 for (fz = 0; fz < numZEdges; ++fz) { 576 for (fx = 0; fx < numYEdges; ++fx) { 577 for (fy = 0; fy < numYVertices; ++fy) { 578 PetscInt face = firstYFace + (fz*numXEdges+fx)*numYVertices + fy; 579 PetscInt edgeL = firstZEdge + (fy*numXVertices+ fx )*numZEdges + fz; 580 PetscInt edgeR = firstZEdge + (fy*numXVertices+((fx+1)%numXVertices))*numZEdges + fz; 581 PetscInt edgeB = firstXEdge + ( fz *numYVertices+fy)*numXEdges + fx; 582 PetscInt edgeT = firstXEdge + (((fz+1)%numZVertices)*numYVertices+fy)*numXEdges + fx; 583 PetscInt ornt[4] = {0, 0, -2, -2}; 584 PetscInt cone[4]; 585 586 if (dim == 3) { 587 /* markers */ 588 if (bdY != DM_BOUNDARY_PERIODIC) { 589 if (fy == numYVertices-1) { 590 ierr = DMPlexSetLabelValue(dm, "marker", face, markerBack);CHKERRQ(ierr); 591 } 592 else if (fy == 0) { 593 ierr = DMPlexSetLabelValue(dm, "marker", face, markerFront);CHKERRQ(ierr); 594 } 595 } 596 } 597 cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL; 598 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr); 599 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr); 600 } 601 } 602 } 603 /* Build z faces */ 604 for (fy = 0; fy < numYEdges; ++fy) { 605 for (fx = 0; fx < numXEdges; ++fx) { 606 for (fz = 0; fz < numZVertices; fz++) { 607 PetscInt face = firstZFace + (fy*numXEdges+fx)*numZVertices + fz; 608 PetscInt edgeL = firstYEdge + (fz*numXVertices+ fx )*numYEdges + fy; 609 PetscInt edgeR = firstYEdge + (fz*numXVertices+((fx+1)%numXVertices))*numYEdges + fy; 610 PetscInt edgeB = firstXEdge + (fz*numYVertices+ fy )*numXEdges + fx; 611 PetscInt edgeT = firstXEdge + (fz*numYVertices+((fy+1)%numYVertices))*numXEdges + fx; 612 PetscInt ornt[4] = {0, 0, -2, -2}; 613 PetscInt cone[4]; 614 615 if (dim == 2) { 616 if (bdX == DM_BOUNDARY_TWIST && fx == numXEdges-1) {edgeR += numYEdges-1-2*fy; ornt[1] = -2;} 617 if (bdY == DM_BOUNDARY_TWIST && fy == numYEdges-1) {edgeT += numXEdges-1-2*fx; ornt[2] = 0;} 618 } 619 else { 620 /* markers */ 621 if (bdZ != DM_BOUNDARY_PERIODIC) { 622 if (fz == numZVertices-1) { 623 ierr = DMPlexSetLabelValue(dm, "marker", face, markerTop);CHKERRQ(ierr); 624 } 625 else if (fz == 0) { 626 ierr = DMPlexSetLabelValue(dm, "marker", face, markerBottom);CHKERRQ(ierr); 627 } 628 } 629 } 630 cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL; 631 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr); 632 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr); 633 } 634 } 635 } 636 /* Build Z edges*/ 637 for (vy = 0; vy < numYVertices; vy++) { 638 for (vx = 0; vx < numXVertices; vx++) { 639 for (ez = 0; ez < numZEdges; ez++) { 640 const PetscInt edge = firstZEdge + (vy*numXVertices+vx)*numZEdges + ez; 641 const PetscInt vertexB = firstVertex + ( ez *numYVertices+vy)*numXVertices + vx; 642 const PetscInt vertexT = firstVertex + (((ez+1)%numZVertices)*numYVertices+vy)*numXVertices + vx; 643 PetscInt cone[2]; 644 645 if (dim == 3) { 646 if (bdX != DM_BOUNDARY_PERIODIC) { 647 if (vx == numXVertices-1) { 648 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 649 } 650 else if (vx == 0) { 651 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 652 } 653 } 654 if (bdY != DM_BOUNDARY_PERIODIC) { 655 if (vy == numYVertices-1) { 656 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerBack);CHKERRQ(ierr); 657 } 658 else if (vy == 0) { 659 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerFront);CHKERRQ(ierr); 660 } 661 } 662 } 663 cone[0] = vertexB; cone[1] = vertexT; 664 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 665 } 666 } 667 } 668 /* Build Y edges*/ 669 for (vz = 0; vz < numZVertices; vz++) { 670 for (vx = 0; vx < numXVertices; vx++) { 671 for (ey = 0; ey < numYEdges; ey++) { 672 const PetscInt nextv = (dim == 2 && bdY == DM_BOUNDARY_TWIST && ey == numYEdges-1) ? (numXVertices-vx-1) : (vz*numYVertices+((ey+1)%numYVertices))*numXVertices + vx; 673 const PetscInt edge = firstYEdge + (vz*numXVertices+vx)*numYEdges + ey; 674 const PetscInt vertexF = firstVertex + (vz*numYVertices+ey)*numXVertices + vx; 675 const PetscInt vertexK = firstVertex + nextv; 676 PetscInt cone[2]; 677 678 cone[0] = vertexF; cone[1] = vertexK; 679 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 680 if (dim == 2) { 681 if ((bdX != DM_BOUNDARY_PERIODIC) && (bdX != DM_BOUNDARY_TWIST)) { 682 if (vx == numXVertices-1) { 683 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 684 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerRight);CHKERRQ(ierr); 685 if (ey == numYEdges-1) { 686 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerRight);CHKERRQ(ierr); 687 } 688 } 689 else if (vx == 0) { 690 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 691 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerLeft);CHKERRQ(ierr); 692 if (ey == numYEdges-1) { 693 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerLeft);CHKERRQ(ierr); 694 } 695 } 696 } 697 } 698 else { 699 if (bdX != DM_BOUNDARY_PERIODIC) { 700 if (vx == numXVertices-1) { 701 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 702 } 703 else if (vx == 0) { 704 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 705 } 706 } 707 if (bdZ != DM_BOUNDARY_PERIODIC) { 708 if (vz == numZVertices-1) { 709 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 710 } 711 else if (vz == 0) { 712 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 713 } 714 } 715 } 716 } 717 } 718 } 719 /* Build X edges*/ 720 for (vz = 0; vz < numZVertices; vz++) { 721 for (vy = 0; vy < numYVertices; vy++) { 722 for (ex = 0; ex < numXEdges; ex++) { 723 const PetscInt nextv = (dim == 2 && bdX == DM_BOUNDARY_TWIST && ex == numXEdges-1) ? (numYVertices-vy-1)*numXVertices : (vz*numYVertices+vy)*numXVertices + (ex+1)%numXVertices; 724 const PetscInt edge = firstXEdge + (vz*numYVertices+vy)*numXEdges + ex; 725 const PetscInt vertexL = firstVertex + (vz*numYVertices+vy)*numXVertices + ex; 726 const PetscInt vertexR = firstVertex + nextv; 727 PetscInt cone[2]; 728 729 cone[0] = vertexL; cone[1] = vertexR; 730 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 731 if (dim == 2) { 732 if ((bdY != DM_BOUNDARY_PERIODIC) && (bdY != DM_BOUNDARY_TWIST)) { 733 if (vy == numYVertices-1) { 734 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 735 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerTop);CHKERRQ(ierr); 736 if (ex == numXEdges-1) { 737 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerTop);CHKERRQ(ierr); 738 } 739 } 740 else if (vy == 0) { 741 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 742 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerBottom);CHKERRQ(ierr); 743 if (ex == numXEdges-1) { 744 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerBottom);CHKERRQ(ierr); 745 } 746 } 747 } 748 } 749 else { 750 if (bdY != DM_BOUNDARY_PERIODIC) { 751 if (vy == numYVertices-1) { 752 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerBack);CHKERRQ(ierr); 753 } 754 else if (vy == 0) { 755 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerFront);CHKERRQ(ierr); 756 } 757 } 758 if (bdZ != DM_BOUNDARY_PERIODIC) { 759 if (vz == numZVertices-1) { 760 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 761 } 762 else if (vz == 0) { 763 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 764 } 765 } 766 } 767 } 768 } 769 } 770 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 771 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 772 /* Build coordinates */ 773 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 774 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 775 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 776 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVertices);CHKERRQ(ierr); 777 for (v = firstVertex; v < firstVertex+numVertices; ++v) { 778 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 779 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 780 } 781 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 782 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 783 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 784 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 785 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 786 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 787 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 788 for (vz = 0; vz < numZVertices; ++vz) { 789 for (vy = 0; vy < numYVertices; ++vy) { 790 for (vx = 0; vx < numXVertices; ++vx) { 791 coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+0] = lower[0] + ((upper[0] - lower[0])/numXEdges)*vx; 792 coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+1] = lower[1] + ((upper[1] - lower[1])/numYEdges)*vy; 793 if (dim == 3) { 794 coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+2] = lower[2] + ((upper[2] - lower[2])/numZEdges)*vz; 795 } 796 } 797 } 798 } 799 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 800 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 801 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 802 } 803 PetscFunctionReturn(0); 804 } 805 806 #undef __FUNCT__ 807 #define __FUNCT__ "DMPlexCreateSquareMesh" 808 /*@ 809 DMPlexCreateSquareMesh - Creates a 2D mesh for a square lattice. 810 811 Collective on MPI_Comm 812 813 Input Parameters: 814 + comm - The communicator for the DM object 815 . lower - The lower left corner coordinates 816 . upper - The upper right corner coordinates 817 . edges - The number of cells in each direction 818 . bdX - The boundary type for the X direction 819 - bdY - The boundary type for the Y direction 820 821 Output Parameter: 822 . dm - The DM object 823 824 Note: Here is the numbering returned for 2 cells in each direction: 825 $ 22--8-23--9--24 826 $ | | | 827 $ 13 2 14 3 15 828 $ | | | 829 $ 19--6-20--7--21 830 $ | | | 831 $ 10 0 11 1 12 832 $ | | | 833 $ 16--4-17--5--18 834 835 Level: beginner 836 837 .keywords: DM, create 838 .seealso: DMPlexCreateBoxMesh(), DMPlexCreateSquareBoundary(), DMPlexCreateCubeBoundary(), DMSetType(), DMCreate() 839 @*/ 840 PetscErrorCode DMPlexCreateSquareMesh(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[], DMBoundaryType bdX, DMBoundaryType bdY) 841 { 842 PetscReal lower3[3], upper3[3]; 843 PetscInt edges3[3]; 844 PetscErrorCode ierr; 845 846 PetscFunctionBegin; 847 lower3[0] = lower[0]; lower3[1] = lower[1]; lower3[2] = 0.; 848 upper3[0] = upper[0]; upper3[1] = upper[1]; upper3[2] = 0.; 849 edges3[0] = edges[0]; edges3[1] = edges[1]; edges3[2] = 0; 850 ierr = DMPlexCreateCubeMesh_Internal(dm, lower3, upper3, edges3, bdX, bdY, DM_BOUNDARY_NONE);CHKERRQ(ierr); 851 PetscFunctionReturn(0); 852 } 853 854 #undef __FUNCT__ 855 #define __FUNCT__ "DMPlexCreateBoxMesh" 856 /*@ 857 DMPlexCreateBoxMesh - Creates a mesh on the tensor product of unit intervals (box) using simplices. 858 859 Collective on MPI_Comm 860 861 Input Parameters: 862 + comm - The communicator for the DM object 863 . dim - The spatial dimension 864 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 865 866 Output Parameter: 867 . dm - The DM object 868 869 Level: beginner 870 871 .keywords: DM, create 872 .seealso: DMPlexCreateHexBoxMesh(), DMSetType(), DMCreate() 873 @*/ 874 PetscErrorCode DMPlexCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscBool interpolate, DM *dm) 875 { 876 DM boundary; 877 PetscErrorCode ierr; 878 879 PetscFunctionBegin; 880 PetscValidPointer(dm, 4); 881 ierr = DMCreate(comm, &boundary);CHKERRQ(ierr); 882 PetscValidLogicalCollectiveInt(boundary,dim,2); 883 ierr = DMSetType(boundary, DMPLEX);CHKERRQ(ierr); 884 ierr = DMSetDimension(boundary, dim-1);CHKERRQ(ierr); 885 switch (dim) { 886 case 2: 887 { 888 PetscReal lower[2] = {0.0, 0.0}; 889 PetscReal upper[2] = {1.0, 1.0}; 890 PetscInt edges[2] = {2, 2}; 891 892 ierr = DMPlexCreateSquareBoundary(boundary, lower, upper, edges);CHKERRQ(ierr); 893 break; 894 } 895 case 3: 896 { 897 PetscReal lower[3] = {0.0, 0.0, 0.0}; 898 PetscReal upper[3] = {1.0, 1.0, 1.0}; 899 PetscInt faces[3] = {1, 1, 1}; 900 901 ierr = DMPlexCreateCubeBoundary(boundary, lower, upper, faces);CHKERRQ(ierr); 902 break; 903 } 904 default: 905 SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim); 906 } 907 ierr = DMPlexGenerate(boundary, NULL, interpolate, dm);CHKERRQ(ierr); 908 ierr = DMDestroy(&boundary);CHKERRQ(ierr); 909 PetscFunctionReturn(0); 910 } 911 912 #undef __FUNCT__ 913 #define __FUNCT__ "DMPlexCreateHexBoxMesh" 914 /*@ 915 DMPlexCreateHexBoxMesh - Creates a mesh on the tensor product of unit intervals (box) using hexahedra. 916 917 Collective on MPI_Comm 918 919 Input Parameters: 920 + comm - The communicator for the DM object 921 . dim - The spatial dimension 922 . periodicX - The boundary type for the X direction 923 . periodicY - The boundary type for the Y direction 924 . periodicZ - The boundary type for the Z direction 925 - cells - The number of cells in each direction 926 927 Output Parameter: 928 . dm - The DM object 929 930 Level: beginner 931 932 .keywords: DM, create 933 .seealso: DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 934 @*/ 935 PetscErrorCode DMPlexCreateHexBoxMesh(MPI_Comm comm, PetscInt dim, const PetscInt cells[], DMBoundaryType periodicX, DMBoundaryType periodicY, DMBoundaryType periodicZ, DM *dm) 936 { 937 PetscErrorCode ierr; 938 939 PetscFunctionBegin; 940 PetscValidPointer(dm, 4); 941 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 942 PetscValidLogicalCollectiveInt(*dm,dim,2); 943 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 944 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 945 switch (dim) { 946 case 2: 947 { 948 PetscReal lower[2] = {0.0, 0.0}; 949 PetscReal upper[2] = {1.0, 1.0}; 950 951 ierr = DMPlexCreateSquareMesh(*dm, lower, upper, cells, periodicX, periodicY);CHKERRQ(ierr); 952 break; 953 } 954 case 3: 955 { 956 PetscReal lower[3] = {0.0, 0.0, 0.0}; 957 PetscReal upper[3] = {1.0, 1.0, 1.0}; 958 959 ierr = DMPlexCreateCubeMesh_Internal(*dm, lower, upper, cells, periodicX, periodicY, periodicZ);CHKERRQ(ierr); 960 break; 961 } 962 default: 963 SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim); 964 } 965 PetscFunctionReturn(0); 966 } 967 968 /* External function declarations here */ 969 extern PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling); 970 extern PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, VecScatter *ctx); 971 extern PetscErrorCode DMCreateDefaultSection_Plex(DM dm); 972 extern PetscErrorCode DMCreateDefaultConstraints_Plex(DM dm); 973 extern PetscErrorCode DMCreateMatrix_Plex(DM dm, Mat *J); 974 extern PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm); 975 extern PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *dmRefined); 976 extern PetscErrorCode DMCoarsen_Plex(DM dm, MPI_Comm comm, DM *dmCoarsened); 977 extern PetscErrorCode DMRefineHierarchy_Plex(DM dm, PetscInt nlevels, DM dmRefined[]); 978 extern PetscErrorCode DMClone_Plex(DM dm, DM *newdm); 979 extern PetscErrorCode DMSetUp_Plex(DM dm); 980 extern PetscErrorCode DMDestroy_Plex(DM dm); 981 extern PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer); 982 extern PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer); 983 extern PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm); 984 extern PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, IS *cellIS); 985 986 #undef __FUNCT__ 987 #define __FUNCT__ "DMPlexReplace_Static" 988 /* Replace dm with the contents of dmNew 989 - Share the DM_Plex structure 990 - Share the coordinates 991 - Share the SF 992 */ 993 static PetscErrorCode DMPlexReplace_Static(DM dm, DM dmNew) 994 { 995 PetscSF sf; 996 DM coordDM; 997 Vec coords; 998 const PetscReal *maxCell, *L; 999 PetscErrorCode ierr; 1000 1001 PetscFunctionBegin; 1002 ierr = DMGetPointSF(dmNew, &sf);CHKERRQ(ierr); 1003 ierr = DMSetPointSF(dm, sf);CHKERRQ(ierr); 1004 ierr = DMGetCoordinateDM(dmNew, &coordDM);CHKERRQ(ierr); 1005 ierr = DMGetCoordinatesLocal(dmNew, &coords);CHKERRQ(ierr); 1006 ierr = DMSetCoordinateDM(dm, coordDM);CHKERRQ(ierr); 1007 ierr = DMSetCoordinatesLocal(dm, coords);CHKERRQ(ierr); 1008 ierr = DMGetPeriodicity(dm, &maxCell, &L);CHKERRQ(ierr); 1009 if (L) {ierr = DMSetPeriodicity(dmNew, maxCell, L);CHKERRQ(ierr);} 1010 ierr = DMDestroy_Plex(dm);CHKERRQ(ierr); 1011 dm->data = dmNew->data; 1012 ((DM_Plex *) dmNew->data)->refct++; 1013 PetscFunctionReturn(0); 1014 } 1015 1016 #undef __FUNCT__ 1017 #define __FUNCT__ "DMPlexSwap_Static" 1018 /* Swap dm with the contents of dmNew 1019 - Swap the DM_Plex structure 1020 - Swap the coordinates 1021 */ 1022 static PetscErrorCode DMPlexSwap_Static(DM dmA, DM dmB) 1023 { 1024 DM coordDMA, coordDMB; 1025 Vec coordsA, coordsB; 1026 void *tmp; 1027 PetscErrorCode ierr; 1028 1029 PetscFunctionBegin; 1030 ierr = DMGetCoordinateDM(dmA, &coordDMA);CHKERRQ(ierr); 1031 ierr = DMGetCoordinateDM(dmB, &coordDMB);CHKERRQ(ierr); 1032 ierr = PetscObjectReference((PetscObject) coordDMA);CHKERRQ(ierr); 1033 ierr = DMSetCoordinateDM(dmA, coordDMB);CHKERRQ(ierr); 1034 ierr = DMSetCoordinateDM(dmB, coordDMA);CHKERRQ(ierr); 1035 ierr = PetscObjectDereference((PetscObject) coordDMA);CHKERRQ(ierr); 1036 1037 ierr = DMGetCoordinatesLocal(dmA, &coordsA);CHKERRQ(ierr); 1038 ierr = DMGetCoordinatesLocal(dmB, &coordsB);CHKERRQ(ierr); 1039 ierr = PetscObjectReference((PetscObject) coordsA);CHKERRQ(ierr); 1040 ierr = DMSetCoordinatesLocal(dmA, coordsB);CHKERRQ(ierr); 1041 ierr = DMSetCoordinatesLocal(dmB, coordsA);CHKERRQ(ierr); 1042 ierr = PetscObjectDereference((PetscObject) coordsA);CHKERRQ(ierr); 1043 tmp = dmA->data; 1044 dmA->data = dmB->data; 1045 dmB->data = tmp; 1046 PetscFunctionReturn(0); 1047 } 1048 1049 #undef __FUNCT__ 1050 #define __FUNCT__ "DMSetFromOptions_NonRefinement_Plex" 1051 PetscErrorCode DMSetFromOptions_NonRefinement_Plex(DM dm) 1052 { 1053 DM_Plex *mesh = (DM_Plex*) dm->data; 1054 DMBoundary b; 1055 PetscErrorCode ierr; 1056 1057 PetscFunctionBegin; 1058 /* Handle boundary conditions */ 1059 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject) dm), NULL, "Boundary condition options", "");CHKERRQ(ierr); 1060 for (b = mesh->boundary; b; b = b->next) { 1061 char optname[1024]; 1062 PetscInt ids[1024], len = 1024, i; 1063 PetscBool flg; 1064 1065 ierr = PetscSNPrintf(optname, sizeof(optname), "-bc_%s", b->name);CHKERRQ(ierr); 1066 ierr = PetscMemzero(ids, sizeof(ids));CHKERRQ(ierr); 1067 ierr = PetscOptionsIntArray(optname, "List of boundary IDs", "", ids, &len, &flg);CHKERRQ(ierr); 1068 if (flg) { 1069 DMLabel label; 1070 1071 ierr = DMPlexGetLabel(dm, b->labelname, &label);CHKERRQ(ierr); 1072 for (i = 0; i < len; ++i) { 1073 PetscBool has; 1074 1075 ierr = DMLabelHasValue(label, ids[i], &has);CHKERRQ(ierr); 1076 if (!has) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Boundary id %D is not present in the label %s", ids[i], b->name); 1077 } 1078 b->numids = len; 1079 ierr = PetscFree(b->ids);CHKERRQ(ierr); 1080 ierr = PetscMalloc1(len, &b->ids);CHKERRQ(ierr); 1081 ierr = PetscMemcpy(b->ids, ids, len*sizeof(PetscInt));CHKERRQ(ierr); 1082 } 1083 } 1084 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1085 /* Handle viewing */ 1086 ierr = PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMView", PETSC_FALSE, &mesh->printSetValues, NULL);CHKERRQ(ierr); 1087 ierr = PetscOptionsInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMView", 0, &mesh->printFEM, NULL);CHKERRQ(ierr); 1088 ierr = PetscOptionsReal("-dm_plex_print_tol", "Tolerance for FEM output", "DMView", PETSC_FALSE, &mesh->printTol, NULL);CHKERRQ(ierr); 1089 PetscFunctionReturn(0); 1090 } 1091 1092 #undef __FUNCT__ 1093 #define __FUNCT__ "DMSetFromOptions_Plex" 1094 PetscErrorCode DMSetFromOptions_Plex(DM dm) 1095 { 1096 PetscInt refine = 0, r; 1097 PetscBool isHierarchy; 1098 PetscErrorCode ierr; 1099 1100 PetscFunctionBegin; 1101 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1102 ierr = PetscOptionsHead("DMPlex Options");CHKERRQ(ierr); 1103 /* Handle DMPlex refinement */ 1104 ierr = PetscOptionsInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL);CHKERRQ(ierr); 1105 ierr = PetscOptionsInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy);CHKERRQ(ierr); 1106 if (refine) {ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE);CHKERRQ(ierr);} 1107 if (refine && isHierarchy) { 1108 DM *dms; 1109 1110 ierr = PetscMalloc1(refine,&dms);CHKERRQ(ierr); 1111 ierr = DMRefineHierarchy(dm, refine, dms);CHKERRQ(ierr); 1112 /* Total hack since we do not pass in a pointer */ 1113 ierr = DMPlexSwap_Static(dm, dms[refine-1]);CHKERRQ(ierr); 1114 if (refine == 1) { 1115 ierr = DMPlexSetCoarseDM(dm, dms[0]);CHKERRQ(ierr); 1116 } else { 1117 ierr = DMPlexSetCoarseDM(dm, dms[refine-2]);CHKERRQ(ierr); 1118 ierr = DMPlexSetCoarseDM(dms[0], dms[refine-1]);CHKERRQ(ierr); 1119 } 1120 /* Free DMs */ 1121 for (r = 0; r < refine; ++r) { 1122 ierr = DMSetFromOptions_NonRefinement_Plex(dm);CHKERRQ(ierr); 1123 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 1124 } 1125 ierr = PetscFree(dms);CHKERRQ(ierr); 1126 } else { 1127 for (r = 0; r < refine; ++r) { 1128 DM refinedMesh; 1129 1130 ierr = DMSetFromOptions_NonRefinement_Plex(dm);CHKERRQ(ierr); 1131 ierr = DMRefine(dm, PetscObjectComm((PetscObject) dm), &refinedMesh);CHKERRQ(ierr); 1132 /* Total hack since we do not pass in a pointer */ 1133 ierr = DMPlexReplace_Static(dm, refinedMesh);CHKERRQ(ierr); 1134 ierr = DMDestroy(&refinedMesh);CHKERRQ(ierr); 1135 } 1136 } 1137 ierr = DMSetFromOptions_NonRefinement_Plex(dm);CHKERRQ(ierr); 1138 ierr = PetscOptionsTail();CHKERRQ(ierr); 1139 PetscFunctionReturn(0); 1140 } 1141 1142 #undef __FUNCT__ 1143 #define __FUNCT__ "DMCreateGlobalVector_Plex" 1144 static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec) 1145 { 1146 PetscErrorCode ierr; 1147 1148 PetscFunctionBegin; 1149 ierr = DMCreateGlobalVector_Section_Private(dm,vec);CHKERRQ(ierr); 1150 /* ierr = VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM);CHKERRQ(ierr); */ 1151 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex);CHKERRQ(ierr); 1152 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex);CHKERRQ(ierr); 1153 PetscFunctionReturn(0); 1154 } 1155 1156 #undef __FUNCT__ 1157 #define __FUNCT__ "DMCreateLocalVector_Plex" 1158 static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec) 1159 { 1160 PetscErrorCode ierr; 1161 1162 PetscFunctionBegin; 1163 ierr = DMCreateLocalVector_Section_Private(dm,vec);CHKERRQ(ierr); 1164 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex_Local);CHKERRQ(ierr); 1165 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex_Local);CHKERRQ(ierr); 1166 PetscFunctionReturn(0); 1167 } 1168 1169 #undef __FUNCT__ 1170 #define __FUNCT__ "DMGetDimPoints_Plex" 1171 static PetscErrorCode DMGetDimPoints_Plex(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd) 1172 { 1173 PetscInt depth, d; 1174 PetscErrorCode ierr; 1175 1176 PetscFunctionBegin; 1177 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1178 if (depth == 1) { 1179 ierr = DMGetDimension(dm, &d);CHKERRQ(ierr); 1180 if (dim == 0) {ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr);} 1181 else if (dim == d) {ierr = DMPlexGetDepthStratum(dm, 1, pStart, pEnd);CHKERRQ(ierr);} 1182 else {*pStart = 0; *pEnd = 0;} 1183 } else { 1184 ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr); 1185 } 1186 PetscFunctionReturn(0); 1187 } 1188 1189 #undef __FUNCT__ 1190 #define __FUNCT__ "DMInitialize_Plex" 1191 PetscErrorCode DMInitialize_Plex(DM dm) 1192 { 1193 PetscFunctionBegin; 1194 dm->ops->view = DMView_Plex; 1195 dm->ops->load = DMLoad_Plex; 1196 dm->ops->setfromoptions = DMSetFromOptions_Plex; 1197 dm->ops->clone = DMClone_Plex; 1198 dm->ops->setup = DMSetUp_Plex; 1199 dm->ops->createdefaultsection = DMCreateDefaultSection_Plex; 1200 dm->ops->createdefaultconstraints = DMCreateDefaultConstraints_Plex; 1201 dm->ops->createglobalvector = DMCreateGlobalVector_Plex; 1202 dm->ops->createlocalvector = DMCreateLocalVector_Plex; 1203 dm->ops->getlocaltoglobalmapping = NULL; 1204 dm->ops->createfieldis = NULL; 1205 dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex; 1206 dm->ops->getcoloring = NULL; 1207 dm->ops->creatematrix = DMCreateMatrix_Plex; 1208 dm->ops->createinterpolation = DMCreateInterpolation_Plex; 1209 dm->ops->getaggregates = NULL; 1210 dm->ops->getinjection = DMCreateInjection_Plex; 1211 dm->ops->refine = DMRefine_Plex; 1212 dm->ops->coarsen = DMCoarsen_Plex; 1213 dm->ops->refinehierarchy = DMRefineHierarchy_Plex; 1214 dm->ops->coarsenhierarchy = NULL; 1215 dm->ops->globaltolocalbegin = NULL; 1216 dm->ops->globaltolocalend = NULL; 1217 dm->ops->localtoglobalbegin = NULL; 1218 dm->ops->localtoglobalend = NULL; 1219 dm->ops->destroy = DMDestroy_Plex; 1220 dm->ops->createsubdm = DMCreateSubDM_Plex; 1221 dm->ops->getdimpoints = DMGetDimPoints_Plex; 1222 dm->ops->locatepoints = DMLocatePoints_Plex; 1223 PetscFunctionReturn(0); 1224 } 1225 1226 #undef __FUNCT__ 1227 #define __FUNCT__ "DMClone_Plex" 1228 PetscErrorCode DMClone_Plex(DM dm, DM *newdm) 1229 { 1230 DM_Plex *mesh = (DM_Plex *) dm->data; 1231 PetscErrorCode ierr; 1232 1233 PetscFunctionBegin; 1234 mesh->refct++; 1235 (*newdm)->data = mesh; 1236 ierr = PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);CHKERRQ(ierr); 1237 ierr = DMInitialize_Plex(*newdm);CHKERRQ(ierr); 1238 PetscFunctionReturn(0); 1239 } 1240 1241 /*MC 1242 DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram. 1243 In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is 1244 specified by a PetscSection object. Ownership in the global representation is determined by 1245 ownership of the underlying DMPlex points. This is specified by another PetscSection object. 1246 1247 Level: intermediate 1248 1249 .seealso: DMType, DMPlexCreate(), DMCreate(), DMSetType() 1250 M*/ 1251 1252 #undef __FUNCT__ 1253 #define __FUNCT__ "DMCreate_Plex" 1254 PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm) 1255 { 1256 DM_Plex *mesh; 1257 PetscInt unit, d; 1258 PetscErrorCode ierr; 1259 1260 PetscFunctionBegin; 1261 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1262 ierr = PetscNewLog(dm,&mesh);CHKERRQ(ierr); 1263 dm->dim = 0; 1264 dm->data = mesh; 1265 1266 mesh->refct = 1; 1267 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);CHKERRQ(ierr); 1268 mesh->maxConeSize = 0; 1269 mesh->cones = NULL; 1270 mesh->coneOrientations = NULL; 1271 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);CHKERRQ(ierr); 1272 mesh->maxSupportSize = 0; 1273 mesh->supports = NULL; 1274 mesh->refinementUniform = PETSC_TRUE; 1275 mesh->refinementLimit = -1.0; 1276 1277 mesh->facesTmp = NULL; 1278 1279 mesh->tetgenOpts = NULL; 1280 mesh->triangleOpts = NULL; 1281 ierr = PetscPartitionerCreate(PetscObjectComm((PetscObject)dm), &mesh->partitioner);CHKERRQ(ierr); 1282 ierr = PetscPartitionerSetTypeFromOptions_Internal(mesh->partitioner);CHKERRQ(ierr); 1283 1284 mesh->subpointMap = NULL; 1285 1286 for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0; 1287 1288 mesh->labels = NULL; 1289 mesh->depthLabel = NULL; 1290 mesh->globalVertexNumbers = NULL; 1291 mesh->globalCellNumbers = NULL; 1292 mesh->anchorSection = NULL; 1293 mesh->anchorIS = NULL; 1294 mesh->createanchors = NULL; 1295 mesh->computeanchormatrix = NULL; 1296 mesh->parentSection = NULL; 1297 mesh->parents = NULL; 1298 mesh->childIDs = NULL; 1299 mesh->childSection = NULL; 1300 mesh->children = NULL; 1301 mesh->referenceTree = NULL; 1302 mesh->getchildsymmetry = NULL; 1303 for (d = 0; d < 8; ++d) mesh->hybridPointMax[d] = PETSC_DETERMINE; 1304 mesh->vtkCellHeight = 0; 1305 mesh->useCone = PETSC_FALSE; 1306 mesh->useClosure = PETSC_TRUE; 1307 mesh->useAnchors = PETSC_FALSE; 1308 1309 mesh->printSetValues = PETSC_FALSE; 1310 mesh->printFEM = 0; 1311 mesh->printTol = 1.0e-10; 1312 1313 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 1314 PetscFunctionReturn(0); 1315 } 1316 1317 #undef __FUNCT__ 1318 #define __FUNCT__ "DMPlexCreate" 1319 /*@ 1320 DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram. 1321 1322 Collective on MPI_Comm 1323 1324 Input Parameter: 1325 . comm - The communicator for the DMPlex object 1326 1327 Output Parameter: 1328 . mesh - The DMPlex object 1329 1330 Level: beginner 1331 1332 .keywords: DMPlex, create 1333 @*/ 1334 PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh) 1335 { 1336 PetscErrorCode ierr; 1337 1338 PetscFunctionBegin; 1339 PetscValidPointer(mesh,2); 1340 ierr = DMCreate(comm, mesh);CHKERRQ(ierr); 1341 ierr = DMSetType(*mesh, DMPLEX);CHKERRQ(ierr); 1342 PetscFunctionReturn(0); 1343 } 1344 1345 #undef __FUNCT__ 1346 #define __FUNCT__ "DMPlexBuildFromCellList_Private" 1347 /* 1348 This takes as input the common mesh generator output, a list of the vertices for each cell 1349 */ 1350 PetscErrorCode DMPlexBuildFromCellList_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[]) 1351 { 1352 PetscInt *cone, c, p; 1353 PetscErrorCode ierr; 1354 1355 PetscFunctionBegin; 1356 ierr = DMPlexSetChart(dm, 0, numCells+numVertices);CHKERRQ(ierr); 1357 for (c = 0; c < numCells; ++c) { 1358 ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr); 1359 } 1360 ierr = DMSetUp(dm);CHKERRQ(ierr); 1361 ierr = DMGetWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 1362 for (c = 0; c < numCells; ++c) { 1363 for (p = 0; p < numCorners; ++p) { 1364 cone[p] = cells[c*numCorners+p]+numCells; 1365 } 1366 ierr = DMPlexSetCone(dm, c, cone);CHKERRQ(ierr); 1367 } 1368 ierr = DMRestoreWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 1369 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 1370 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 1371 PetscFunctionReturn(0); 1372 } 1373 1374 #undef __FUNCT__ 1375 #define __FUNCT__ "DMPlexBuildCoordinates_Private" 1376 /* 1377 This takes as input the coordinates for each vertex 1378 */ 1379 PetscErrorCode DMPlexBuildCoordinates_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, const double vertexCoords[]) 1380 { 1381 PetscSection coordSection; 1382 Vec coordinates; 1383 PetscScalar *coords; 1384 PetscInt coordSize, v, d; 1385 PetscErrorCode ierr; 1386 1387 PetscFunctionBegin; 1388 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1389 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1390 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 1391 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVertices);CHKERRQ(ierr); 1392 for (v = numCells; v < numCells+numVertices; ++v) { 1393 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 1394 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 1395 } 1396 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1397 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1398 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 1399 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 1400 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1401 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1402 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1403 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1404 for (v = 0; v < numVertices; ++v) { 1405 for (d = 0; d < spaceDim; ++d) { 1406 coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d]; 1407 } 1408 } 1409 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1410 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 1411 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1412 PetscFunctionReturn(0); 1413 } 1414 1415 #undef __FUNCT__ 1416 #define __FUNCT__ "DMPlexCreateFromCellList" 1417 /*@C 1418 DMPlexCreateFromCellList - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM 1419 1420 Input Parameters: 1421 + comm - The communicator 1422 . dim - The topological dimension of the mesh 1423 . numCells - The number of cells 1424 . numVertices - The number of vertices 1425 . numCorners - The number of vertices for each cell 1426 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 1427 . cells - An array of numCells*numCorners numbers, the vertices for each cell 1428 . spaceDim - The spatial dimension used for coordinates 1429 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 1430 1431 Output Parameter: 1432 . dm - The DM 1433 1434 Note: Two triangles sharing a face 1435 $ 1436 $ 2 1437 $ / | \ 1438 $ / | \ 1439 $ / | \ 1440 $ 0 0 | 1 3 1441 $ \ | / 1442 $ \ | / 1443 $ \ | / 1444 $ 1 1445 would have input 1446 $ numCells = 2, numVertices = 4 1447 $ cells = [0 1 2 1 3 2] 1448 $ 1449 which would result in the DMPlex 1450 $ 1451 $ 4 1452 $ / | \ 1453 $ / | \ 1454 $ / | \ 1455 $ 2 0 | 1 5 1456 $ \ | / 1457 $ \ | / 1458 $ \ | / 1459 $ 3 1460 1461 Level: beginner 1462 1463 .seealso: DMPlexCreateFromDAG(), DMPlexCreate() 1464 @*/ 1465 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) 1466 { 1467 PetscErrorCode ierr; 1468 1469 PetscFunctionBegin; 1470 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1471 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1472 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1473 ierr = DMPlexBuildFromCellList_Private(*dm, numCells, numVertices, numCorners, cells);CHKERRQ(ierr); 1474 if (interpolate) { 1475 DM idm = NULL; 1476 1477 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1478 ierr = DMDestroy(dm);CHKERRQ(ierr); 1479 *dm = idm; 1480 } 1481 ierr = DMPlexBuildCoordinates_Private(*dm, spaceDim, numCells, numVertices, vertexCoords);CHKERRQ(ierr); 1482 PetscFunctionReturn(0); 1483 } 1484 1485 #undef __FUNCT__ 1486 #define __FUNCT__ "DMPlexCreateFromDAG" 1487 /*@ 1488 DMPlexCreateFromDAG - This takes as input the adjacency-list representation of the Directed Acyclic Graph (Hasse Diagram) encoding a mesh, and produces a DM 1489 1490 Input Parameters: 1491 + dm - The empty DM object, usually from DMCreate() and DMSetDimension() 1492 . depth - The depth of the DAG 1493 . numPoints - The number of points at each depth 1494 . coneSize - The cone size of each point 1495 . cones - The concatenation of the cone points for each point, the cone list must be oriented correctly for each point 1496 . coneOrientations - The orientation of each cone point 1497 - vertexCoords - An array of numVertices*dim numbers, the coordinates of each vertex 1498 1499 Output Parameter: 1500 . dm - The DM 1501 1502 Note: Two triangles sharing a face would have input 1503 $ depth = 1, numPoints = [4 2], coneSize = [3 3 0 0 0 0] 1504 $ cones = [2 3 4 3 5 4], coneOrientations = [0 0 0 0 0 0] 1505 $ vertexCoords = [-1.0 0.0 0.0 -1.0 0.0 1.0 1.0 0.0] 1506 $ 1507 which would result in the DMPlex 1508 $ 1509 $ 4 1510 $ / | \ 1511 $ / | \ 1512 $ / | \ 1513 $ 2 0 | 1 5 1514 $ \ | / 1515 $ \ | / 1516 $ \ | / 1517 $ 3 1518 $ 1519 $ Notice that all points are numbered consecutively, unlikely DMPlexCreateFromCellList() 1520 1521 Level: advanced 1522 1523 .seealso: DMPlexCreateFromCellList(), DMPlexCreate() 1524 @*/ 1525 PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[]) 1526 { 1527 Vec coordinates; 1528 PetscSection coordSection; 1529 PetscScalar *coords; 1530 PetscInt coordSize, firstVertex = -1, pStart = 0, pEnd = 0, p, v, dim, dimEmbed, d, off; 1531 PetscErrorCode ierr; 1532 1533 PetscFunctionBegin; 1534 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1535 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 1536 if (dimEmbed < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Embedding dimension %d cannot be less than intrinsic dimension %d",dimEmbed,dim); 1537 for (d = 0; d <= depth; ++d) pEnd += numPoints[d]; 1538 ierr = DMPlexSetChart(dm, pStart, pEnd);CHKERRQ(ierr); 1539 for (p = pStart; p < pEnd; ++p) { 1540 ierr = DMPlexSetConeSize(dm, p, coneSize[p-pStart]);CHKERRQ(ierr); 1541 if (firstVertex < 0 && !coneSize[p - pStart]) { 1542 firstVertex = p - pStart; 1543 } 1544 } 1545 if (firstVertex < 0 && numPoints[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected %d vertices but could not find any", numPoints[0]); 1546 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 1547 for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) { 1548 ierr = DMPlexSetCone(dm, p, &cones[off]);CHKERRQ(ierr); 1549 ierr = DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);CHKERRQ(ierr); 1550 } 1551 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 1552 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 1553 /* Build coordinates */ 1554 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1555 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1556 ierr = PetscSectionSetFieldComponents(coordSection, 0, dimEmbed);CHKERRQ(ierr); 1557 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0]);CHKERRQ(ierr); 1558 for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) { 1559 ierr = PetscSectionSetDof(coordSection, v, dimEmbed);CHKERRQ(ierr); 1560 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed);CHKERRQ(ierr); 1561 } 1562 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1563 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1564 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 1565 ierr = VecSetBlockSize(coordinates, dimEmbed);CHKERRQ(ierr); 1566 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1567 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1568 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1569 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1570 for (v = 0; v < numPoints[0]; ++v) { 1571 PetscInt off; 1572 1573 ierr = PetscSectionGetOffset(coordSection, v+firstVertex, &off);CHKERRQ(ierr); 1574 for (d = 0; d < dimEmbed; ++d) { 1575 coords[off+d] = vertexCoords[v*dimEmbed+d]; 1576 } 1577 } 1578 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1579 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 1580 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1581 PetscFunctionReturn(0); 1582 } 1583 1584 #undef __FUNCT__ 1585 #define __FUNCT__ "DMPlexCreateFromFile" 1586 /*@C 1587 DMPlexCreateFromFile - This takes a filename and produces a DM 1588 1589 Input Parameters: 1590 + comm - The communicator 1591 . filename - A file name 1592 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 1593 1594 Output Parameter: 1595 . dm - The DM 1596 1597 Level: beginner 1598 1599 .seealso: DMPlexCreateFromDAG(), DMPlexCreateFromCellList(), DMPlexCreate() 1600 @*/ 1601 PetscErrorCode DMPlexCreateFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 1602 { 1603 const char *extGmsh = ".msh"; 1604 const char *extCGNS = ".cgns"; 1605 const char *extExodus = ".exo"; 1606 size_t len; 1607 PetscBool isGmsh, isCGNS, isExodus; 1608 PetscMPIInt rank; 1609 PetscErrorCode ierr; 1610 1611 PetscFunctionBegin; 1612 PetscValidPointer(filename, 2); 1613 PetscValidPointer(dm, 4); 1614 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 1615 ierr = PetscStrlen(filename, &len);CHKERRQ(ierr); 1616 if (!len) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Filename must be a valid path"); 1617 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGmsh, 4, &isGmsh);CHKERRQ(ierr); 1618 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extCGNS, 5, &isCGNS);CHKERRQ(ierr); 1619 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extExodus, 4, &isExodus);CHKERRQ(ierr); 1620 if (isGmsh) { 1621 PetscViewer viewer; 1622 1623 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 1624 ierr = PetscViewerSetType(viewer, PETSCVIEWERASCII);CHKERRQ(ierr); 1625 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 1626 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 1627 ierr = DMPlexCreateGmsh(comm, viewer, interpolate, dm);CHKERRQ(ierr); 1628 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 1629 } else if (isCGNS) { 1630 ierr = DMPlexCreateCGNSFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 1631 } else if (isExodus) { 1632 ierr = DMPlexCreateExodusFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 1633 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot load file %s: unrecognized extension", filename); 1634 PetscFunctionReturn(0); 1635 } 1636 1637 #undef __FUNCT__ 1638 #define __FUNCT__ "DMPlexCreateReferenceCell" 1639 /*@ 1640 DMPlexCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell 1641 1642 Collective on comm 1643 1644 Input Parameters: 1645 + comm - The communicator 1646 . dim - The spatial dimension 1647 - simplex - Flag for simplex, otherwise use a tensor-product cell 1648 1649 Output Parameter: 1650 . refdm - The reference cell 1651 1652 Level: intermediate 1653 1654 .keywords: reference cell 1655 .seealso: 1656 @*/ 1657 PetscErrorCode DMPlexCreateReferenceCell(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *refdm) 1658 { 1659 DM rdm; 1660 PetscErrorCode ierr; 1661 1662 PetscFunctionBeginUser; 1663 ierr = DMCreate(comm, &rdm);CHKERRQ(ierr); 1664 ierr = DMSetType(rdm, DMPLEX);CHKERRQ(ierr); 1665 ierr = DMSetDimension(rdm, dim);CHKERRQ(ierr); 1666 switch (dim) { 1667 case 0: 1668 { 1669 PetscInt numPoints[1] = {1}; 1670 PetscInt coneSize[1] = {0}; 1671 PetscInt cones[1] = {0}; 1672 PetscInt coneOrientations[1] = {0}; 1673 PetscScalar vertexCoords[1] = {0.0}; 1674 1675 ierr = DMPlexCreateFromDAG(rdm, 0, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 1676 } 1677 break; 1678 case 1: 1679 { 1680 PetscInt numPoints[2] = {2, 1}; 1681 PetscInt coneSize[3] = {2, 0, 0}; 1682 PetscInt cones[2] = {1, 2}; 1683 PetscInt coneOrientations[2] = {0, 0}; 1684 PetscScalar vertexCoords[2] = {-1.0, 1.0}; 1685 1686 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 1687 } 1688 break; 1689 case 2: 1690 if (simplex) { 1691 PetscInt numPoints[2] = {3, 1}; 1692 PetscInt coneSize[4] = {3, 0, 0, 0}; 1693 PetscInt cones[3] = {1, 2, 3}; 1694 PetscInt coneOrientations[3] = {0, 0, 0}; 1695 PetscScalar vertexCoords[6] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0}; 1696 1697 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 1698 } else { 1699 PetscInt numPoints[2] = {4, 1}; 1700 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 1701 PetscInt cones[4] = {1, 2, 3, 4}; 1702 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 1703 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0}; 1704 1705 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 1706 } 1707 break; 1708 case 3: 1709 if (simplex) { 1710 PetscInt numPoints[2] = {4, 1}; 1711 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 1712 PetscInt cones[4] = {1, 3, 2, 4}; 1713 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 1714 PetscScalar vertexCoords[12] = {-1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0}; 1715 1716 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 1717 } else { 1718 PetscInt numPoints[2] = {8, 1}; 1719 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0}; 1720 PetscInt cones[8] = {1, 4, 3, 2, 5, 6, 7, 8}; 1721 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 1722 PetscScalar vertexCoords[24] = {-1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1723 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0}; 1724 1725 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 1726 } 1727 break; 1728 default: 1729 SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Cannot create reference cell for dimension %d", dim); 1730 } 1731 *refdm = NULL; 1732 ierr = DMPlexInterpolate(rdm, refdm);CHKERRQ(ierr); 1733 ierr = DMPlexCopyCoordinates(rdm, *refdm);CHKERRQ(ierr); 1734 ierr = DMDestroy(&rdm);CHKERRQ(ierr); 1735 PetscFunctionReturn(0); 1736 } 1737