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 = DMPlexSetDimension(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, NULL, 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__ "DMPlexCreateSquareMesh" 440 /*@ 441 DMPlexCreateSquareMesh - Creates a 2D mesh for a square lattice. 442 443 Collective on MPI_Comm 444 445 Input Parameters: 446 + comm - The communicator for the DM object 447 . lower - The lower left corner coordinates 448 . upper - The upper right corner coordinates 449 . edges - The number of cells in each direction 450 . bdX - The boundary type for the X direction 451 - bdY - The boundary type for the Y direction 452 453 Output Parameter: 454 . dm - The DM object 455 456 Note: Here is the numbering returned for 2 cells in each direction: 457 $ 22--8-23--9--24 458 $ | | | 459 $ 13 2 14 3 15 460 $ | | | 461 $ 19--6-20--7--21 462 $ | | | 463 $ 10 0 11 1 12 464 $ | | | 465 $ 16--4-17--5--18 466 467 Level: beginner 468 469 .keywords: DM, create 470 .seealso: DMPlexCreateBoxMesh(), DMPlexCreateSquareBoundary(), DMPlexCreateCubeBoundary(), DMSetType(), DMCreate() 471 @*/ 472 PetscErrorCode DMPlexCreateSquareMesh(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[], DMBoundaryType bdX, DMBoundaryType bdY) 473 { 474 PetscInt markerTop = 1; 475 PetscInt markerBottom = 1; 476 PetscInt markerRight = 1; 477 PetscInt markerLeft = 1; 478 PetscBool markerSeparate = PETSC_FALSE; 479 PetscMPIInt rank; 480 PetscErrorCode ierr; 481 482 PetscFunctionBegin; 483 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 484 ierr = PetscOptionsGetBool(((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);CHKERRQ(ierr); 485 if (markerSeparate) { 486 markerTop = 3; 487 markerBottom = 1; 488 markerRight = 2; 489 markerLeft = 4; 490 } 491 { 492 const PetscInt numXEdges = !rank ? edges[0] : 0; 493 const PetscInt numYEdges = !rank ? edges[1] : 0; 494 const PetscInt numXVertices = !rank ? (bdX == DM_BOUNDARY_PERIODIC || bdX == DM_BOUNDARY_TWIST ? edges[0] : edges[0]+1) : 0; 495 const PetscInt numYVertices = !rank ? (bdY == DM_BOUNDARY_PERIODIC || bdY == DM_BOUNDARY_TWIST ? edges[1] : edges[1]+1) : 0; 496 const PetscInt numTotXEdges = numXEdges*numYVertices; 497 const PetscInt numTotYEdges = numYEdges*numXVertices; 498 const PetscInt numVertices = numXVertices*numYVertices; 499 const PetscInt numEdges = numTotXEdges + numTotYEdges; 500 const PetscInt numFaces = numXEdges*numYEdges; 501 const PetscInt firstVertex = numFaces; 502 const PetscInt firstXEdge = numFaces + numVertices; 503 const PetscInt firstYEdge = numFaces + numVertices + numTotXEdges; 504 Vec coordinates; 505 PetscSection coordSection; 506 PetscScalar *coords; 507 PetscInt coordSize; 508 PetscInt v, vx, vy; 509 PetscInt f, fx, fy, e, ex, ey; 510 511 ierr = DMPlexSetChart(dm, 0, numFaces+numEdges+numVertices);CHKERRQ(ierr); 512 for (f = 0; f < numFaces; ++f) { 513 ierr = DMPlexSetConeSize(dm, f, 4);CHKERRQ(ierr); 514 } 515 for (e = firstXEdge; e < firstXEdge+numEdges; ++e) { 516 ierr = DMPlexSetConeSize(dm, e, 2);CHKERRQ(ierr); 517 } 518 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 519 /* Build faces */ 520 for (fy = 0; fy < numYEdges; ++fy) { 521 for (fx = 0; fx < numXEdges; ++fx) { 522 PetscInt face = fy*numXEdges + fx; 523 PetscInt edgeL = firstYEdge + fx *numYEdges + fy; 524 PetscInt edgeR = firstYEdge + ((fx+1)%numXVertices)*numYEdges + fy; 525 PetscInt edgeB = firstXEdge + fy *numXEdges + fx; 526 PetscInt edgeT = firstXEdge + ((fy+1)%numYVertices)*numXEdges + fx; 527 PetscInt ornt[4] = {0, 0, -2, -2}; 528 PetscInt cone[4]; 529 530 if (bdX == DM_BOUNDARY_TWIST && fx == numXEdges-1) {edgeR += numYEdges-1-2*fy; ornt[1] = -2;} 531 if (bdY == DM_BOUNDARY_TWIST && fy == numYEdges-1) {edgeT += numXEdges-1-2*fx; ornt[2] = 0;} 532 cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL; 533 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr); 534 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr); 535 } 536 } 537 /* Build Y edges*/ 538 for (vx = 0; vx < numXVertices; vx++) { 539 for (ey = 0; ey < numYEdges; ey++) { 540 const PetscInt nextv = (bdY == DM_BOUNDARY_TWIST && ey == numYEdges-1) ? (numXVertices-vx-1) : ((ey+1)%numYVertices)*numXVertices + vx; 541 const PetscInt edge = firstYEdge + vx*numYEdges + ey; 542 const PetscInt vertexB = firstVertex + ey*numXVertices + vx; 543 const PetscInt vertexT = firstVertex + nextv; 544 PetscInt cone[2]; 545 546 cone[0] = vertexB; cone[1] = vertexT; 547 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 548 if ((bdX != DM_BOUNDARY_PERIODIC) && (bdX != DM_BOUNDARY_TWIST)) { 549 if (vx == numXVertices-1) { 550 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 551 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerRight);CHKERRQ(ierr); 552 if (ey == numYEdges-1) { 553 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerRight);CHKERRQ(ierr); 554 } 555 } else if (vx == 0) { 556 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 557 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerLeft);CHKERRQ(ierr); 558 if (ey == numYEdges-1) { 559 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerLeft);CHKERRQ(ierr); 560 } 561 } 562 } 563 } 564 } 565 /* Build X edges*/ 566 for (vy = 0; vy < numYVertices; vy++) { 567 for (ex = 0; ex < numXEdges; ex++) { 568 const PetscInt nextv = (bdX == DM_BOUNDARY_TWIST && ex == numXEdges-1) ? (numYVertices-vy-1)*numXVertices : vy*numXVertices + (ex+1)%numXVertices; 569 const PetscInt edge = firstXEdge + vy*numXEdges + ex; 570 const PetscInt vertexL = firstVertex + vy*numXVertices + ex; 571 const PetscInt vertexR = firstVertex + nextv; 572 PetscInt cone[2]; 573 574 cone[0] = vertexL; cone[1] = vertexR; 575 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 576 if ((bdY != DM_BOUNDARY_PERIODIC) && (bdY != DM_BOUNDARY_TWIST)) { 577 if (vy == numYVertices-1) { 578 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 579 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerTop);CHKERRQ(ierr); 580 if (ex == numXEdges-1) { 581 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerTop);CHKERRQ(ierr); 582 } 583 } else if (vy == 0) { 584 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 585 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerBottom);CHKERRQ(ierr); 586 if (ex == numXEdges-1) { 587 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerBottom);CHKERRQ(ierr); 588 } 589 } 590 } 591 } 592 } 593 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 594 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 595 /* Build coordinates */ 596 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 597 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVertices);CHKERRQ(ierr); 598 for (v = firstVertex; v < firstVertex+numVertices; ++v) { 599 ierr = PetscSectionSetDof(coordSection, v, 2);CHKERRQ(ierr); 600 } 601 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 602 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 603 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 604 ierr = VecSetBlockSize(coordinates, 2);CHKERRQ(ierr); 605 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 606 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 607 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 608 for (vy = 0; vy < numYVertices; ++vy) { 609 for (vx = 0; vx < numXVertices; ++vx) { 610 coords[(vy*numXVertices+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/numXEdges)*vx; 611 coords[(vy*numXVertices+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/numYEdges)*vy; 612 } 613 } 614 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 615 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 616 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 617 } 618 PetscFunctionReturn(0); 619 } 620 621 #undef __FUNCT__ 622 #define __FUNCT__ "DMPlexCreateBoxMesh" 623 /*@ 624 DMPlexCreateBoxMesh - Creates a mesh on the tensor product of unit intervals (box) using simplices. 625 626 Collective on MPI_Comm 627 628 Input Parameters: 629 + comm - The communicator for the DM object 630 . dim - The spatial dimension 631 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 632 633 Output Parameter: 634 . dm - The DM object 635 636 Level: beginner 637 638 .keywords: DM, create 639 .seealso: DMPlexCreateHexBoxMesh(), DMSetType(), DMCreate() 640 @*/ 641 PetscErrorCode DMPlexCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscBool interpolate, DM *dm) 642 { 643 DM boundary; 644 PetscErrorCode ierr; 645 646 PetscFunctionBegin; 647 PetscValidPointer(dm, 4); 648 ierr = DMCreate(comm, &boundary);CHKERRQ(ierr); 649 PetscValidLogicalCollectiveInt(boundary,dim,2); 650 ierr = DMSetType(boundary, DMPLEX);CHKERRQ(ierr); 651 ierr = DMPlexSetDimension(boundary, dim-1);CHKERRQ(ierr); 652 switch (dim) { 653 case 2: 654 { 655 PetscReal lower[2] = {0.0, 0.0}; 656 PetscReal upper[2] = {1.0, 1.0}; 657 PetscInt edges[2] = {2, 2}; 658 659 ierr = DMPlexCreateSquareBoundary(boundary, lower, upper, edges);CHKERRQ(ierr); 660 break; 661 } 662 case 3: 663 { 664 PetscReal lower[3] = {0.0, 0.0, 0.0}; 665 PetscReal upper[3] = {1.0, 1.0, 1.0}; 666 PetscInt faces[3] = {1, 1, 1}; 667 668 ierr = DMPlexCreateCubeBoundary(boundary, lower, upper, faces);CHKERRQ(ierr); 669 break; 670 } 671 default: 672 SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim); 673 } 674 ierr = DMPlexGenerate(boundary, NULL, interpolate, dm);CHKERRQ(ierr); 675 ierr = DMDestroy(&boundary);CHKERRQ(ierr); 676 PetscFunctionReturn(0); 677 } 678 679 #undef __FUNCT__ 680 #define __FUNCT__ "DMPlexCreateHexBoxMesh" 681 /*@ 682 DMPlexCreateHexBoxMesh - Creates a mesh on the tensor product of unit intervals (box) using hexahedra. 683 684 Collective on MPI_Comm 685 686 Input Parameters: 687 + comm - The communicator for the DM object 688 . dim - The spatial dimension 689 . periodicX - The boundary type for the X direction 690 . periodicY - The boundary type for the Y direction 691 . periodicZ - The boundary type for the Z direction 692 - cells - The number of cells in each direction 693 694 Output Parameter: 695 . dm - The DM object 696 697 Level: beginner 698 699 .keywords: DM, create 700 .seealso: DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 701 @*/ 702 PetscErrorCode DMPlexCreateHexBoxMesh(MPI_Comm comm, PetscInt dim, const PetscInt cells[], DMBoundaryType periodicX, DMBoundaryType periodicY, DMBoundaryType periodicZ, DM *dm) 703 { 704 PetscErrorCode ierr; 705 706 PetscFunctionBegin; 707 PetscValidPointer(dm, 4); 708 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 709 PetscValidLogicalCollectiveInt(*dm,dim,2); 710 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 711 ierr = DMPlexSetDimension(*dm, dim);CHKERRQ(ierr); 712 switch (dim) { 713 case 2: 714 { 715 PetscReal lower[2] = {0.0, 0.0}; 716 PetscReal upper[2] = {1.0, 1.0}; 717 718 ierr = DMPlexCreateSquareMesh(*dm, lower, upper, cells, periodicX, periodicY);CHKERRQ(ierr); 719 break; 720 } 721 #if 0 722 case 3: 723 { 724 PetscReal lower[3] = {0.0, 0.0, 0.0}; 725 PetscReal upper[3] = {1.0, 1.0, 1.0}; 726 727 ierr = DMPlexCreateCubeMesh(boundary, lower, upper, cells, periodicX, periodicY, periodicZ);CHKERRQ(ierr); 728 break; 729 } 730 #endif 731 default: 732 SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim); 733 } 734 PetscFunctionReturn(0); 735 } 736 737 /* External function declarations here */ 738 extern PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling); 739 extern PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, VecScatter *ctx); 740 extern PetscErrorCode DMCreateDefaultSection_Plex(DM dm); 741 extern PetscErrorCode DMCreateMatrix_Plex(DM dm, Mat *J); 742 extern PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm); 743 extern PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *dmRefined); 744 extern PetscErrorCode DMCoarsen_Plex(DM dm, MPI_Comm comm, DM *dmCoarsened); 745 extern PetscErrorCode DMRefineHierarchy_Plex(DM dm, PetscInt nlevels, DM dmRefined[]); 746 extern PetscErrorCode DMClone_Plex(DM dm, DM *newdm); 747 extern PetscErrorCode DMSetUp_Plex(DM dm); 748 extern PetscErrorCode DMDestroy_Plex(DM dm); 749 extern PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer); 750 extern PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer); 751 extern PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm); 752 extern PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, IS *cellIS); 753 754 #undef __FUNCT__ 755 #define __FUNCT__ "DMPlexReplace_Static" 756 /* Replace dm with the contents of dmNew 757 - Share the DM_Plex structure 758 - Share the coordinates 759 */ 760 static PetscErrorCode DMPlexReplace_Static(DM dm, DM dmNew) 761 { 762 PetscSection coordSection; 763 Vec coords; 764 PetscErrorCode ierr; 765 766 PetscFunctionBegin; 767 ierr = DMGetCoordinateSection(dmNew, &coordSection);CHKERRQ(ierr); 768 ierr = DMGetCoordinatesLocal(dmNew, &coords);CHKERRQ(ierr); 769 ierr = DMSetCoordinateSection(dm, coordSection);CHKERRQ(ierr); 770 ierr = DMSetCoordinatesLocal(dm, coords);CHKERRQ(ierr); 771 ierr = DMDestroy_Plex(dm);CHKERRQ(ierr); 772 dm->data = dmNew->data; 773 ((DM_Plex *) dmNew->data)->refct++; 774 PetscFunctionReturn(0); 775 } 776 777 #undef __FUNCT__ 778 #define __FUNCT__ "DMPlexSwap_Static" 779 /* Swap dm with the contents of dmNew 780 - Swap the DM_Plex structure 781 - Swap the coordinates 782 */ 783 static PetscErrorCode DMPlexSwap_Static(DM dmA, DM dmB) 784 { 785 DM coordDMA, coordDMB; 786 Vec coordsA, coordsB; 787 void *tmp; 788 PetscErrorCode ierr; 789 790 PetscFunctionBegin; 791 ierr = DMGetCoordinateDM(dmA, &coordDMA);CHKERRQ(ierr); 792 ierr = DMGetCoordinateDM(dmB, &coordDMB);CHKERRQ(ierr); 793 ierr = PetscObjectReference((PetscObject) coordDMA);CHKERRQ(ierr); 794 ierr = DMSetCoordinateDM(dmA, coordDMB);CHKERRQ(ierr); 795 ierr = DMSetCoordinateDM(dmB, coordDMA);CHKERRQ(ierr); 796 ierr = PetscObjectDereference((PetscObject) coordDMA);CHKERRQ(ierr); 797 798 ierr = DMGetCoordinatesLocal(dmA, &coordsA);CHKERRQ(ierr); 799 ierr = DMGetCoordinatesLocal(dmB, &coordsB);CHKERRQ(ierr); 800 ierr = PetscObjectReference((PetscObject) coordsA);CHKERRQ(ierr); 801 ierr = DMSetCoordinatesLocal(dmA, coordsB);CHKERRQ(ierr); 802 ierr = DMSetCoordinatesLocal(dmB, coordsA);CHKERRQ(ierr); 803 ierr = PetscObjectDereference((PetscObject) coordsA);CHKERRQ(ierr); 804 tmp = dmA->data; 805 dmA->data = dmB->data; 806 dmB->data = tmp; 807 PetscFunctionReturn(0); 808 } 809 810 #undef __FUNCT__ 811 #define __FUNCT__ "DMSetFromOptions_NonRefinement_Plex" 812 PetscErrorCode DMSetFromOptions_NonRefinement_Plex(DM dm) 813 { 814 DM_Plex *mesh = (DM_Plex*) dm->data; 815 DMBoundary b; 816 PetscErrorCode ierr; 817 818 PetscFunctionBegin; 819 /* Handle boundary conditions */ 820 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject) dm), NULL, "Boundary condition options", "");CHKERRQ(ierr); 821 for (b = mesh->boundary; b; b = b->next) { 822 char optname[1024]; 823 PetscInt ids[1024], len = 1024, i; 824 PetscBool flg; 825 826 ierr = PetscSNPrintf(optname, sizeof(optname), "-bc_%s", b->name);CHKERRQ(ierr); 827 ierr = PetscMemzero(ids, sizeof(ids));CHKERRQ(ierr); 828 ierr = PetscOptionsIntArray(optname, "List of boundary IDs", "", ids, &len, &flg);CHKERRQ(ierr); 829 if (flg) { 830 DMLabel label; 831 832 ierr = DMPlexGetLabel(dm, b->name, &label);CHKERRQ(ierr); 833 for (i = 0; i < len; ++i) { 834 PetscBool has; 835 836 ierr = DMLabelHasValue(label, ids[i], &has);CHKERRQ(ierr); 837 if (!has) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Boundary id %D is not present in the label %s", ids[i], b->name); 838 } 839 b->numids = len; 840 ierr = PetscFree(b->ids);CHKERRQ(ierr); 841 ierr = PetscMalloc1(len, &b->ids);CHKERRQ(ierr); 842 ierr = PetscMemcpy(b->ids, ids, len*sizeof(PetscInt));CHKERRQ(ierr); 843 } 844 } 845 ierr = PetscOptionsEnd();CHKERRQ(ierr); 846 /* Handle viewing */ 847 ierr = PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMView", PETSC_FALSE, &mesh->printSetValues, NULL);CHKERRQ(ierr); 848 ierr = PetscOptionsInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMView", 0, &mesh->printFEM, NULL);CHKERRQ(ierr); 849 ierr = PetscOptionsReal("-dm_plex_print_tol", "Tolerance for FEM output", "DMView", PETSC_FALSE, &mesh->printTol, NULL);CHKERRQ(ierr); 850 PetscFunctionReturn(0); 851 } 852 853 #undef __FUNCT__ 854 #define __FUNCT__ "DMSetFromOptions_Plex" 855 PetscErrorCode DMSetFromOptions_Plex(DM dm) 856 { 857 PetscInt refine = 0, r; 858 PetscBool isHierarchy; 859 PetscErrorCode ierr; 860 861 PetscFunctionBegin; 862 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 863 ierr = PetscOptionsHead("DMPlex Options");CHKERRQ(ierr); 864 /* Handle DMPlex refinement */ 865 ierr = PetscOptionsInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL);CHKERRQ(ierr); 866 ierr = PetscOptionsInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy);CHKERRQ(ierr); 867 ierr = DMPlexSetRefinementUniform(dm, refine ? PETSC_TRUE : PETSC_FALSE);CHKERRQ(ierr); 868 if (refine && isHierarchy) { 869 DM *dms; 870 871 ierr = PetscMalloc1(refine,&dms);CHKERRQ(ierr); 872 ierr = DMRefineHierarchy(dm, refine, dms);CHKERRQ(ierr); 873 /* Total hack since we do not pass in a pointer */ 874 ierr = DMPlexSwap_Static(dm, dms[refine-1]);CHKERRQ(ierr); 875 if (refine == 1) { 876 ierr = DMPlexSetCoarseDM(dm, dms[0]);CHKERRQ(ierr); 877 } else { 878 ierr = DMPlexSetCoarseDM(dm, dms[refine-2]);CHKERRQ(ierr); 879 ierr = DMPlexSetCoarseDM(dms[0], dms[refine-1]);CHKERRQ(ierr); 880 } 881 /* Free DMs */ 882 for (r = 0; r < refine; ++r) { 883 ierr = DMSetFromOptions_NonRefinement_Plex(dm);CHKERRQ(ierr); 884 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 885 } 886 ierr = PetscFree(dms);CHKERRQ(ierr); 887 } else { 888 for (r = 0; r < refine; ++r) { 889 DM refinedMesh; 890 891 ierr = DMSetFromOptions_NonRefinement_Plex(dm);CHKERRQ(ierr); 892 ierr = DMRefine(dm, PetscObjectComm((PetscObject) dm), &refinedMesh);CHKERRQ(ierr); 893 /* Total hack since we do not pass in a pointer */ 894 ierr = DMPlexReplace_Static(dm, refinedMesh);CHKERRQ(ierr); 895 ierr = DMDestroy(&refinedMesh);CHKERRQ(ierr); 896 } 897 } 898 ierr = DMSetFromOptions_NonRefinement_Plex(dm);CHKERRQ(ierr); 899 ierr = PetscOptionsTail();CHKERRQ(ierr); 900 PetscFunctionReturn(0); 901 } 902 903 #undef __FUNCT__ 904 #define __FUNCT__ "DMCreateGlobalVector_Plex" 905 static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec) 906 { 907 PetscErrorCode ierr; 908 909 PetscFunctionBegin; 910 ierr = DMCreateGlobalVector_Section_Private(dm,vec);CHKERRQ(ierr); 911 /* ierr = VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM);CHKERRQ(ierr); */ 912 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex);CHKERRQ(ierr); 913 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex);CHKERRQ(ierr); 914 PetscFunctionReturn(0); 915 } 916 917 #undef __FUNCT__ 918 #define __FUNCT__ "DMCreateLocalVector_Plex" 919 static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec) 920 { 921 PetscErrorCode ierr; 922 923 PetscFunctionBegin; 924 ierr = DMCreateLocalVector_Section_Private(dm,vec);CHKERRQ(ierr); 925 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex_Local);CHKERRQ(ierr); 926 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex_Local);CHKERRQ(ierr); 927 PetscFunctionReturn(0); 928 } 929 930 #undef __FUNCT__ 931 #define __FUNCT__ "DMInitialize_Plex" 932 PetscErrorCode DMInitialize_Plex(DM dm) 933 { 934 PetscFunctionBegin; 935 dm->ops->view = DMView_Plex; 936 dm->ops->load = DMLoad_Plex; 937 dm->ops->setfromoptions = DMSetFromOptions_Plex; 938 dm->ops->clone = DMClone_Plex; 939 dm->ops->setup = DMSetUp_Plex; 940 dm->ops->createdefaultsection = DMCreateDefaultSection_Plex; 941 dm->ops->createglobalvector = DMCreateGlobalVector_Plex; 942 dm->ops->createlocalvector = DMCreateLocalVector_Plex; 943 dm->ops->getlocaltoglobalmapping = NULL; 944 dm->ops->createfieldis = NULL; 945 dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex; 946 dm->ops->getcoloring = NULL; 947 dm->ops->creatematrix = DMCreateMatrix_Plex; 948 dm->ops->createinterpolation = DMCreateInterpolation_Plex; 949 dm->ops->getaggregates = NULL; 950 dm->ops->getinjection = DMCreateInjection_Plex; 951 dm->ops->refine = DMRefine_Plex; 952 dm->ops->coarsen = DMCoarsen_Plex; 953 dm->ops->refinehierarchy = DMRefineHierarchy_Plex; 954 dm->ops->coarsenhierarchy = NULL; 955 dm->ops->globaltolocalbegin = NULL; 956 dm->ops->globaltolocalend = NULL; 957 dm->ops->localtoglobalbegin = NULL; 958 dm->ops->localtoglobalend = NULL; 959 dm->ops->destroy = DMDestroy_Plex; 960 dm->ops->createsubdm = DMCreateSubDM_Plex; 961 dm->ops->locatepoints = DMLocatePoints_Plex; 962 PetscFunctionReturn(0); 963 } 964 965 #undef __FUNCT__ 966 #define __FUNCT__ "DMClone_Plex" 967 PetscErrorCode DMClone_Plex(DM dm, DM *newdm) 968 { 969 DM_Plex *mesh = (DM_Plex *) dm->data; 970 PetscErrorCode ierr; 971 972 PetscFunctionBegin; 973 mesh->refct++; 974 (*newdm)->data = mesh; 975 ierr = PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);CHKERRQ(ierr); 976 ierr = DMInitialize_Plex(*newdm);CHKERRQ(ierr); 977 PetscFunctionReturn(0); 978 } 979 980 /*MC 981 DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram. 982 In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is 983 specified by a PetscSection object. Ownership in the global representation is determined by 984 ownership of the underlying DMPlex points. This is specified by another PetscSection object. 985 986 Level: intermediate 987 988 .seealso: DMType, DMPlexCreate(), DMCreate(), DMSetType() 989 M*/ 990 991 #undef __FUNCT__ 992 #define __FUNCT__ "DMCreate_Plex" 993 PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm) 994 { 995 DM_Plex *mesh; 996 PetscInt unit, d; 997 PetscErrorCode ierr; 998 999 PetscFunctionBegin; 1000 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1001 ierr = PetscNewLog(dm,&mesh);CHKERRQ(ierr); 1002 dm->data = mesh; 1003 1004 mesh->refct = 1; 1005 mesh->dim = 0; 1006 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);CHKERRQ(ierr); 1007 mesh->maxConeSize = 0; 1008 mesh->cones = NULL; 1009 mesh->coneOrientations = NULL; 1010 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);CHKERRQ(ierr); 1011 mesh->maxSupportSize = 0; 1012 mesh->supports = NULL; 1013 mesh->refinementUniform = PETSC_TRUE; 1014 mesh->refinementLimit = -1.0; 1015 1016 mesh->facesTmp = NULL; 1017 1018 mesh->subpointMap = NULL; 1019 1020 for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0; 1021 1022 mesh->labels = NULL; 1023 mesh->depthLabel = NULL; 1024 mesh->globalVertexNumbers = NULL; 1025 mesh->globalCellNumbers = NULL; 1026 for (d = 0; d < 8; ++d) mesh->hybridPointMax[d] = PETSC_DETERMINE; 1027 mesh->vtkCellHeight = 0; 1028 mesh->useCone = PETSC_FALSE; 1029 mesh->useClosure = PETSC_TRUE; 1030 1031 mesh->printSetValues = PETSC_FALSE; 1032 mesh->printFEM = 0; 1033 mesh->printTol = 1.0e-10; 1034 1035 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 1036 PetscFunctionReturn(0); 1037 } 1038 1039 #undef __FUNCT__ 1040 #define __FUNCT__ "DMPlexCreate" 1041 /*@ 1042 DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram. 1043 1044 Collective on MPI_Comm 1045 1046 Input Parameter: 1047 . comm - The communicator for the DMPlex object 1048 1049 Output Parameter: 1050 . mesh - The DMPlex object 1051 1052 Level: beginner 1053 1054 .keywords: DMPlex, create 1055 @*/ 1056 PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh) 1057 { 1058 PetscErrorCode ierr; 1059 1060 PetscFunctionBegin; 1061 PetscValidPointer(mesh,2); 1062 ierr = DMCreate(comm, mesh);CHKERRQ(ierr); 1063 ierr = DMSetType(*mesh, DMPLEX);CHKERRQ(ierr); 1064 PetscFunctionReturn(0); 1065 } 1066 1067 #undef __FUNCT__ 1068 #define __FUNCT__ "DMPlexBuildFromCellList_Private" 1069 /* 1070 This takes as input the common mesh generator output, a list of the vertices for each cell 1071 */ 1072 PetscErrorCode DMPlexBuildFromCellList_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[]) 1073 { 1074 PetscInt *cone, c, p; 1075 PetscErrorCode ierr; 1076 1077 PetscFunctionBegin; 1078 ierr = DMPlexSetChart(dm, 0, numCells+numVertices);CHKERRQ(ierr); 1079 for (c = 0; c < numCells; ++c) { 1080 ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr); 1081 } 1082 ierr = DMSetUp(dm);CHKERRQ(ierr); 1083 ierr = DMGetWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 1084 for (c = 0; c < numCells; ++c) { 1085 for (p = 0; p < numCorners; ++p) { 1086 cone[p] = cells[c*numCorners+p]+numCells; 1087 } 1088 ierr = DMPlexSetCone(dm, c, cone);CHKERRQ(ierr); 1089 } 1090 ierr = DMRestoreWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 1091 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 1092 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 1093 PetscFunctionReturn(0); 1094 } 1095 1096 #undef __FUNCT__ 1097 #define __FUNCT__ "DMPlexBuildCoordinates_Private" 1098 /* 1099 This takes as input the coordinates for each vertex 1100 */ 1101 PetscErrorCode DMPlexBuildCoordinates_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, const double vertexCoords[]) 1102 { 1103 PetscSection coordSection; 1104 Vec coordinates; 1105 PetscScalar *coords; 1106 PetscInt coordSize, v, d; 1107 PetscErrorCode ierr; 1108 1109 PetscFunctionBegin; 1110 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1111 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1112 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 1113 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVertices);CHKERRQ(ierr); 1114 for (v = numCells; v < numCells+numVertices; ++v) { 1115 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 1116 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 1117 } 1118 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1119 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1120 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 1121 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 1122 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1123 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1124 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1125 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1126 for (v = 0; v < numVertices; ++v) { 1127 for (d = 0; d < spaceDim; ++d) { 1128 coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d]; 1129 } 1130 } 1131 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1132 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 1133 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1134 PetscFunctionReturn(0); 1135 } 1136 1137 #undef __FUNCT__ 1138 #define __FUNCT__ "DMPlexCreateFromCellList" 1139 /*@C 1140 DMPlexCreateFromCellList - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM 1141 1142 Input Parameters: 1143 + comm - The communicator 1144 . dim - The topological dimension of the mesh 1145 . numCells - The number of cells 1146 . numVertices - The number of vertices 1147 . numCorners - The number of vertices for each cell 1148 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 1149 . cells - An array of numCells*numCorners numbers, the vertices for each cell 1150 . spaceDim - The spatial dimension used for coordinates 1151 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 1152 1153 Output Parameter: 1154 . dm - The DM 1155 1156 Note: Two triangles sharing a face 1157 $ 1158 $ 2 1159 $ / | \ 1160 $ / | \ 1161 $ / | \ 1162 $ 0 0 | 1 3 1163 $ \ | / 1164 $ \ | / 1165 $ \ | / 1166 $ 1 1167 would have input 1168 $ numCells = 2, numVertices = 4 1169 $ cells = [0 1 2 1 3 2] 1170 $ 1171 which would result in the DMPlex 1172 $ 1173 $ 4 1174 $ / | \ 1175 $ / | \ 1176 $ / | \ 1177 $ 2 0 | 1 5 1178 $ \ | / 1179 $ \ | / 1180 $ \ | / 1181 $ 3 1182 1183 Level: beginner 1184 1185 .seealso: DMPlexCreateFromDAG(), DMPlexCreate() 1186 @*/ 1187 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) 1188 { 1189 PetscErrorCode ierr; 1190 1191 PetscFunctionBegin; 1192 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1193 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1194 ierr = DMPlexSetDimension(*dm, dim);CHKERRQ(ierr); 1195 ierr = DMPlexBuildFromCellList_Private(*dm, numCells, numVertices, numCorners, cells);CHKERRQ(ierr); 1196 if (interpolate) { 1197 DM idm = NULL; 1198 1199 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1200 ierr = DMDestroy(dm);CHKERRQ(ierr); 1201 *dm = idm; 1202 } 1203 ierr = DMPlexBuildCoordinates_Private(*dm, spaceDim, numCells, numVertices, vertexCoords);CHKERRQ(ierr); 1204 PetscFunctionReturn(0); 1205 } 1206 1207 #undef __FUNCT__ 1208 #define __FUNCT__ "DMPlexCreateFromDAG" 1209 /*@ 1210 DMPlexCreateFromDAG - This takes as input the adjacency-list representation of the Directed Acyclic Graph (Hasse Diagram) encoding a mesh, and produces a DM 1211 1212 Input Parameters: 1213 + dm - The empty DM object, usually from DMCreate() and DMPlexSetDimension() 1214 . depth - The depth of the DAG 1215 . numPoints - The number of points at each depth 1216 . coneSize - The cone size of each point 1217 . cones - The concatenation of the cone points for each point, the cone list must be oriented correctly for each point 1218 . coneOrientations - The orientation of each cone point 1219 - vertexCoords - An array of numVertices*dim numbers, the coordinates of each vertex 1220 1221 Output Parameter: 1222 . dm - The DM 1223 1224 Note: Two triangles sharing a face would have input 1225 $ depth = 1, numPoints = [4 2], coneSize = [3 3 0 0 0 0] 1226 $ cones = [2 3 4 3 5 4], coneOrientations = [0 0 0 0 0 0] 1227 $ vertexCoords = [-1.0 0.0 0.0 -1.0 0.0 1.0 1.0 0.0] 1228 $ 1229 which would result in the DMPlex 1230 $ 1231 $ 4 1232 $ / | \ 1233 $ / | \ 1234 $ / | \ 1235 $ 2 0 | 1 5 1236 $ \ | / 1237 $ \ | / 1238 $ \ | / 1239 $ 3 1240 $ 1241 $ Notice that all points are numbered consecutively, unlikely DMPlexCreateFromCellList() 1242 1243 Level: advanced 1244 1245 .seealso: DMPlexCreateFromCellList(), DMPlexCreate() 1246 @*/ 1247 PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[]) 1248 { 1249 Vec coordinates; 1250 PetscSection coordSection; 1251 PetscScalar *coords; 1252 PetscInt coordSize, firstVertex = numPoints[depth], pStart = 0, pEnd = 0, p, v, dim, d, off; 1253 PetscErrorCode ierr; 1254 1255 PetscFunctionBegin; 1256 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 1257 for (d = 0; d <= depth; ++d) pEnd += numPoints[d]; 1258 ierr = DMPlexSetChart(dm, pStart, pEnd);CHKERRQ(ierr); 1259 for (p = pStart; p < pEnd; ++p) { 1260 ierr = DMPlexSetConeSize(dm, p, coneSize[p-pStart]);CHKERRQ(ierr); 1261 } 1262 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 1263 for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) { 1264 ierr = DMPlexSetCone(dm, p, &cones[off]);CHKERRQ(ierr); 1265 ierr = DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);CHKERRQ(ierr); 1266 } 1267 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 1268 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 1269 /* Build coordinates */ 1270 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1271 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1272 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 1273 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0]);CHKERRQ(ierr); 1274 for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) { 1275 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 1276 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 1277 } 1278 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1279 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1280 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 1281 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 1282 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1283 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1284 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1285 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1286 for (v = 0; v < numPoints[0]; ++v) { 1287 PetscInt off; 1288 1289 ierr = PetscSectionGetOffset(coordSection, v+firstVertex, &off);CHKERRQ(ierr); 1290 for (d = 0; d < dim; ++d) { 1291 coords[off+d] = vertexCoords[v*dim+d]; 1292 } 1293 } 1294 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1295 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 1296 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1297 PetscFunctionReturn(0); 1298 } 1299