1 #define PETSCDM_DLL 2 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 3 #include <petscdmda.h> 4 #include <petscsf.h> 5 6 /*@ 7 DMPlexCreateDoublet - Creates a mesh of two cells of the specified type, optionally with later refinement. 8 9 Collective on MPI_Comm 10 11 Input Parameters: 12 + comm - The communicator for the DM object 13 . dim - The spatial dimension 14 . simplex - Flag for simplicial cells, otherwise they are tensor product cells 15 . interpolate - Flag to create intermediate mesh pieces (edges, faces) 16 . refinementUniform - Flag for uniform parallel refinement 17 - refinementLimit - A nonzero number indicates the largest admissible volume for a refined cell 18 19 Output Parameter: 20 . dm - The DM object 21 22 Level: beginner 23 24 .keywords: DM, create 25 .seealso: DMSetType(), DMCreate() 26 @*/ 27 PetscErrorCode DMPlexCreateDoublet(MPI_Comm comm, PetscInt dim, PetscBool simplex, PetscBool interpolate, PetscBool refinementUniform, PetscReal refinementLimit, DM *newdm) 28 { 29 DM dm; 30 PetscInt p; 31 PetscMPIInt rank; 32 PetscErrorCode ierr; 33 34 PetscFunctionBegin; 35 ierr = DMCreate(comm, &dm);CHKERRQ(ierr); 36 ierr = DMSetType(dm, DMPLEX);CHKERRQ(ierr); 37 ierr = DMSetDimension(dm, dim);CHKERRQ(ierr); 38 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 39 switch (dim) { 40 case 2: 41 if (simplex) {ierr = PetscObjectSetName((PetscObject) dm, "triangular");CHKERRQ(ierr);} 42 else {ierr = PetscObjectSetName((PetscObject) dm, "quadrilateral");CHKERRQ(ierr);} 43 break; 44 case 3: 45 if (simplex) {ierr = PetscObjectSetName((PetscObject) dm, "tetrahedral");CHKERRQ(ierr);} 46 else {ierr = PetscObjectSetName((PetscObject) dm, "hexahedral");CHKERRQ(ierr);} 47 break; 48 default: 49 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %d", dim); 50 } 51 if (rank) { 52 PetscInt numPoints[2] = {0, 0}; 53 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 54 } else { 55 switch (dim) { 56 case 2: 57 if (simplex) { 58 PetscInt numPoints[2] = {4, 2}; 59 PetscInt coneSize[6] = {3, 3, 0, 0, 0, 0}; 60 PetscInt cones[6] = {2, 3, 4, 5, 4, 3}; 61 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 62 PetscScalar vertexCoords[8] = {-0.5, 0.5, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5}; 63 PetscInt markerPoints[8] = {2, 1, 3, 1, 4, 1, 5, 1}; 64 65 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 66 for (p = 0; p < 4; ++p) {ierr = DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 67 } else { 68 PetscInt numPoints[2] = {6, 2}; 69 PetscInt coneSize[8] = {4, 4, 0, 0, 0, 0, 0, 0}; 70 PetscInt cones[8] = {2, 3, 4, 5, 3, 6, 7, 4}; 71 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 72 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}; 73 74 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 75 } 76 break; 77 case 3: 78 if (simplex) { 79 PetscInt numPoints[2] = {5, 2}; 80 PetscInt coneSize[7] = {4, 4, 0, 0, 0, 0, 0}; 81 PetscInt cones[8] = {4, 3, 5, 2, 5, 3, 4, 6}; 82 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 83 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}; 84 PetscInt markerPoints[10] = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1}; 85 86 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 87 for (p = 0; p < 5; ++p) {ierr = DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 88 } else { 89 PetscInt numPoints[2] = {12, 2}; 90 PetscInt coneSize[14] = {8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 91 PetscInt cones[16] = {2, 3, 4, 5, 6, 7, 8, 9, 5, 4, 10, 11, 7, 12, 13, 8}; 92 PetscInt coneOrientations[16] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 93 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, 94 -1.0, -0.5, 0.5, 0.0, -0.5, 0.5, 0.0, 0.5, 0.5, -1.0, 0.5, 0.5, 95 1.0, 0.5, -0.5, 1.0, -0.5, -0.5, 1.0, -0.5, 0.5, 1.0, 0.5, 0.5}; 96 97 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 98 } 99 break; 100 default: 101 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %d", dim); 102 } 103 } 104 *newdm = dm; 105 if (refinementLimit > 0.0) { 106 DM rdm; 107 const char *name; 108 109 ierr = DMPlexSetRefinementUniform(*newdm, PETSC_FALSE);CHKERRQ(ierr); 110 ierr = DMPlexSetRefinementLimit(*newdm, refinementLimit);CHKERRQ(ierr); 111 ierr = DMRefine(*newdm, comm, &rdm);CHKERRQ(ierr); 112 ierr = PetscObjectGetName((PetscObject) *newdm, &name);CHKERRQ(ierr); 113 ierr = PetscObjectSetName((PetscObject) rdm, name);CHKERRQ(ierr); 114 ierr = DMDestroy(newdm);CHKERRQ(ierr); 115 *newdm = rdm; 116 } 117 if (interpolate) { 118 DM idm = NULL; 119 const char *name; 120 121 ierr = DMPlexInterpolate(*newdm, &idm);CHKERRQ(ierr); 122 ierr = PetscObjectGetName((PetscObject) *newdm, &name);CHKERRQ(ierr); 123 ierr = PetscObjectSetName((PetscObject) idm, name);CHKERRQ(ierr); 124 ierr = DMPlexCopyCoordinates(*newdm, idm);CHKERRQ(ierr); 125 ierr = DMCopyLabels(*newdm, idm);CHKERRQ(ierr); 126 ierr = DMDestroy(newdm);CHKERRQ(ierr); 127 *newdm = idm; 128 } 129 { 130 DM refinedMesh = NULL; 131 DM distributedMesh = NULL; 132 133 /* Distribute mesh over processes */ 134 ierr = DMPlexDistribute(*newdm, 0, NULL, &distributedMesh);CHKERRQ(ierr); 135 if (distributedMesh) { 136 ierr = DMDestroy(newdm);CHKERRQ(ierr); 137 *newdm = distributedMesh; 138 } 139 if (refinementUniform) { 140 ierr = DMPlexSetRefinementUniform(*newdm, refinementUniform);CHKERRQ(ierr); 141 ierr = DMRefine(*newdm, comm, &refinedMesh);CHKERRQ(ierr); 142 if (refinedMesh) { 143 ierr = DMDestroy(newdm);CHKERRQ(ierr); 144 *newdm = refinedMesh; 145 } 146 } 147 } 148 PetscFunctionReturn(0); 149 } 150 151 /*@ 152 DMPlexCreateSquareBoundary - Creates a 1D mesh the is the boundary of a square lattice. 153 154 Collective on MPI_Comm 155 156 Input Parameters: 157 + comm - The communicator for the DM object 158 . lower - The lower left corner coordinates 159 . upper - The upper right corner coordinates 160 - edges - The number of cells in each direction 161 162 Output Parameter: 163 . dm - The DM object 164 165 Note: Here is the numbering returned for 2 cells in each direction: 166 $ 18--5-17--4--16 167 $ | | | 168 $ 6 10 3 169 $ | | | 170 $ 19-11-20--9--15 171 $ | | | 172 $ 7 8 2 173 $ | | | 174 $ 12--0-13--1--14 175 176 Level: beginner 177 178 .keywords: DM, create 179 .seealso: DMPlexCreateBoxMesh(), DMPlexCreateCubeBoundary(), DMSetType(), DMCreate() 180 @*/ 181 PetscErrorCode DMPlexCreateSquareBoundary(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[]) 182 { 183 const PetscInt numVertices = (edges[0]+1)*(edges[1]+1); 184 const PetscInt numEdges = edges[0]*(edges[1]+1) + (edges[0]+1)*edges[1]; 185 PetscInt markerTop = 1; 186 PetscInt markerBottom = 1; 187 PetscInt markerRight = 1; 188 PetscInt markerLeft = 1; 189 PetscBool markerSeparate = PETSC_FALSE; 190 Vec coordinates; 191 PetscSection coordSection; 192 PetscScalar *coords; 193 PetscInt coordSize; 194 PetscMPIInt rank; 195 PetscInt v, vx, vy; 196 PetscErrorCode ierr; 197 198 PetscFunctionBegin; 199 ierr = PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);CHKERRQ(ierr); 200 if (markerSeparate) { 201 markerTop = 3; 202 markerBottom = 1; 203 markerRight = 2; 204 markerLeft = 4; 205 } 206 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 207 if (!rank) { 208 PetscInt e, ex, ey; 209 210 ierr = DMPlexSetChart(dm, 0, numEdges+numVertices);CHKERRQ(ierr); 211 for (e = 0; e < numEdges; ++e) { 212 ierr = DMPlexSetConeSize(dm, e, 2);CHKERRQ(ierr); 213 } 214 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 215 for (vx = 0; vx <= edges[0]; vx++) { 216 for (ey = 0; ey < edges[1]; ey++) { 217 PetscInt edge = vx*edges[1] + ey + edges[0]*(edges[1]+1); 218 PetscInt vertex = ey*(edges[0]+1) + vx + numEdges; 219 PetscInt cone[2]; 220 221 cone[0] = vertex; cone[1] = vertex+edges[0]+1; 222 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 223 if (vx == edges[0]) { 224 ierr = DMSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 225 ierr = DMSetLabelValue(dm, "marker", cone[0], markerRight);CHKERRQ(ierr); 226 if (ey == edges[1]-1) { 227 ierr = DMSetLabelValue(dm, "marker", cone[1], markerRight);CHKERRQ(ierr); 228 } 229 } else if (vx == 0) { 230 ierr = DMSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 231 ierr = DMSetLabelValue(dm, "marker", cone[0], markerLeft);CHKERRQ(ierr); 232 if (ey == edges[1]-1) { 233 ierr = DMSetLabelValue(dm, "marker", cone[1], markerLeft);CHKERRQ(ierr); 234 } 235 } 236 } 237 } 238 for (vy = 0; vy <= edges[1]; vy++) { 239 for (ex = 0; ex < edges[0]; ex++) { 240 PetscInt edge = vy*edges[0] + ex; 241 PetscInt vertex = vy*(edges[0]+1) + ex + numEdges; 242 PetscInt cone[2]; 243 244 cone[0] = vertex; cone[1] = vertex+1; 245 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 246 if (vy == edges[1]) { 247 ierr = DMSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 248 ierr = DMSetLabelValue(dm, "marker", cone[0], markerTop);CHKERRQ(ierr); 249 if (ex == edges[0]-1) { 250 ierr = DMSetLabelValue(dm, "marker", cone[1], markerTop);CHKERRQ(ierr); 251 } 252 } else if (vy == 0) { 253 ierr = DMSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 254 ierr = DMSetLabelValue(dm, "marker", cone[0], markerBottom);CHKERRQ(ierr); 255 if (ex == edges[0]-1) { 256 ierr = DMSetLabelValue(dm, "marker", cone[1], markerBottom);CHKERRQ(ierr); 257 } 258 } 259 } 260 } 261 } 262 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 263 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 264 /* Build coordinates */ 265 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 266 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 267 ierr = PetscSectionSetChart(coordSection, numEdges, numEdges + numVertices);CHKERRQ(ierr); 268 ierr = PetscSectionSetFieldComponents(coordSection, 0, 2);CHKERRQ(ierr); 269 for (v = numEdges; v < numEdges+numVertices; ++v) { 270 ierr = PetscSectionSetDof(coordSection, v, 2);CHKERRQ(ierr); 271 ierr = PetscSectionSetFieldDof(coordSection, v, 0, 2);CHKERRQ(ierr); 272 } 273 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 274 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 275 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 276 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 277 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 278 ierr = VecSetBlockSize(coordinates, 2);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 /*@ 294 DMPlexCreateCubeBoundary - Creates a 2D mesh the is the boundary of a cubic lattice. 295 296 Collective on MPI_Comm 297 298 Input Parameters: 299 + comm - The communicator for the DM object 300 . lower - The lower left front corner coordinates 301 . upper - The upper right back corner coordinates 302 - edges - The number of cells in each direction 303 304 Output Parameter: 305 . dm - The DM object 306 307 Level: beginner 308 309 .keywords: DM, create 310 .seealso: DMPlexCreateBoxMesh(), DMPlexCreateSquareBoundary(), DMSetType(), DMCreate() 311 @*/ 312 PetscErrorCode DMPlexCreateCubeBoundary(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt faces[]) 313 { 314 PetscInt vertices[3], numVertices; 315 PetscInt numFaces = 2*faces[0]*faces[1] + 2*faces[1]*faces[2] + 2*faces[0]*faces[2]; 316 Vec coordinates; 317 PetscSection coordSection; 318 PetscScalar *coords; 319 PetscInt coordSize; 320 PetscMPIInt rank; 321 PetscInt v, vx, vy, vz; 322 PetscInt voffset, iface=0, cone[4]; 323 PetscErrorCode ierr; 324 325 PetscFunctionBegin; 326 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"); 327 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 328 vertices[0] = faces[0]+1; vertices[1] = faces[1]+1; vertices[2] = faces[2]+1; 329 numVertices = vertices[0]*vertices[1]*vertices[2]; 330 if (!rank) { 331 PetscInt f; 332 333 ierr = DMPlexSetChart(dm, 0, numFaces+numVertices);CHKERRQ(ierr); 334 for (f = 0; f < numFaces; ++f) { 335 ierr = DMPlexSetConeSize(dm, f, 4);CHKERRQ(ierr); 336 } 337 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 338 for (v = 0; v < numFaces+numVertices; ++v) { 339 ierr = DMSetLabelValue(dm, "marker", v, 1);CHKERRQ(ierr); 340 } 341 342 /* Side 0 (Top) */ 343 for (vy = 0; vy < faces[1]; vy++) { 344 for (vx = 0; vx < faces[0]; vx++) { 345 voffset = numFaces + vertices[0]*vertices[1]*(vertices[2]-1) + vy*vertices[0] + vx; 346 cone[0] = voffset; cone[1] = voffset+1; cone[2] = voffset+vertices[0]+1; cone[3] = voffset+vertices[0]; 347 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 348 iface++; 349 } 350 } 351 352 /* Side 1 (Bottom) */ 353 for (vy = 0; vy < faces[1]; vy++) { 354 for (vx = 0; vx < faces[0]; vx++) { 355 voffset = numFaces + vy*(faces[0]+1) + vx; 356 cone[0] = voffset+1; cone[1] = voffset; cone[2] = voffset+vertices[0]; cone[3] = voffset+vertices[0]+1; 357 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 358 iface++; 359 } 360 } 361 362 /* Side 2 (Front) */ 363 for (vz = 0; vz < faces[2]; vz++) { 364 for (vx = 0; vx < faces[0]; vx++) { 365 voffset = numFaces + vz*vertices[0]*vertices[1] + vx; 366 cone[0] = voffset; cone[1] = voffset+1; cone[2] = voffset+vertices[0]*vertices[1]+1; cone[3] = voffset+vertices[0]*vertices[1]; 367 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 368 iface++; 369 } 370 } 371 372 /* Side 3 (Back) */ 373 for (vz = 0; vz < faces[2]; vz++) { 374 for (vx = 0; vx < faces[0]; vx++) { 375 voffset = numFaces + vz*vertices[0]*vertices[1] + vertices[0]*(vertices[1]-1) + vx; 376 cone[0] = voffset+vertices[0]*vertices[1]; cone[1] = voffset+vertices[0]*vertices[1]+1; 377 cone[2] = voffset+1; cone[3] = voffset; 378 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 379 iface++; 380 } 381 } 382 383 /* Side 4 (Left) */ 384 for (vz = 0; vz < faces[2]; vz++) { 385 for (vy = 0; vy < faces[1]; vy++) { 386 voffset = numFaces + vz*vertices[0]*vertices[1] + vy*vertices[0]; 387 cone[0] = voffset; cone[1] = voffset+vertices[0]*vertices[1]; 388 cone[2] = voffset+vertices[0]*vertices[1]+vertices[0]; cone[3] = voffset+vertices[0]; 389 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 390 iface++; 391 } 392 } 393 394 /* Side 5 (Right) */ 395 for (vz = 0; vz < faces[2]; vz++) { 396 for (vy = 0; vy < faces[1]; vy++) { 397 voffset = numFaces + vz*vertices[0]*vertices[1] + vy*vertices[0] + vx; 398 cone[0] = voffset+vertices[0]*vertices[1]; cone[1] = voffset; 399 cone[2] = voffset+vertices[0]; cone[3] = voffset+vertices[0]*vertices[1]+vertices[0]; 400 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 401 iface++; 402 } 403 } 404 } 405 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 406 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 407 /* Build coordinates */ 408 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 409 ierr = PetscSectionSetChart(coordSection, numFaces, numFaces + numVertices);CHKERRQ(ierr); 410 for (v = numFaces; v < numFaces+numVertices; ++v) { 411 ierr = PetscSectionSetDof(coordSection, v, 3);CHKERRQ(ierr); 412 } 413 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 414 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 415 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 416 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 417 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 418 ierr = VecSetBlockSize(coordinates, 3);CHKERRQ(ierr); 419 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 420 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 421 for (vz = 0; vz <= faces[2]; ++vz) { 422 for (vy = 0; vy <= faces[1]; ++vy) { 423 for (vx = 0; vx <= faces[0]; ++vx) { 424 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+0] = lower[0] + ((upper[0] - lower[0])/faces[0])*vx; 425 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+1] = lower[1] + ((upper[1] - lower[1])/faces[1])*vy; 426 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+2] = lower[2] + ((upper[2] - lower[2])/faces[2])*vz; 427 } 428 } 429 } 430 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 431 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 432 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 433 PetscFunctionReturn(0); 434 } 435 436 /*@ 437 DMPlexCreateBoxMesh - Creates a mesh on the tensor product of unit intervals (box) using simplices. 438 439 Collective on MPI_Comm 440 441 Input Parameters: 442 + comm - The communicator for the DM object 443 . dim - The spatial dimension 444 . numFaces - Number of faces per dimension 445 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 446 447 Output Parameter: 448 . dm - The DM object 449 450 Level: beginner 451 452 .keywords: DM, create 453 .seealso: DMPlexCreateHexBoxMesh(), DMSetType(), DMCreate() 454 @*/ 455 PetscErrorCode DMPlexCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscInt numFaces, PetscBool interpolate, DM *dm) 456 { 457 DM boundary; 458 PetscErrorCode ierr; 459 460 PetscFunctionBegin; 461 PetscValidPointer(dm, 4); 462 ierr = DMCreate(comm, &boundary);CHKERRQ(ierr); 463 PetscValidLogicalCollectiveInt(boundary,dim,2); 464 ierr = DMSetType(boundary, DMPLEX);CHKERRQ(ierr); 465 ierr = DMSetDimension(boundary, dim-1);CHKERRQ(ierr); 466 ierr = DMSetCoordinateDim(boundary, dim);CHKERRQ(ierr); 467 switch (dim) { 468 case 2: 469 { 470 PetscReal lower[2] = {0.0, 0.0}; 471 PetscReal upper[2] = {1.0, 1.0}; 472 PetscInt edges[2]; 473 474 edges[0] = numFaces; edges[1] = numFaces; 475 ierr = DMPlexCreateSquareBoundary(boundary, lower, upper, edges);CHKERRQ(ierr); 476 break; 477 } 478 case 3: 479 { 480 PetscReal lower[3] = {0.0, 0.0, 0.0}; 481 PetscReal upper[3] = {1.0, 1.0, 1.0}; 482 PetscInt faces[3]; 483 484 faces[0] = numFaces; faces[1] = numFaces; faces[2] = numFaces; 485 ierr = DMPlexCreateCubeBoundary(boundary, lower, upper, faces);CHKERRQ(ierr); 486 break; 487 } 488 default: 489 SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim); 490 } 491 ierr = DMPlexGenerate(boundary, NULL, interpolate, dm);CHKERRQ(ierr); 492 ierr = DMDestroy(&boundary);CHKERRQ(ierr); 493 PetscFunctionReturn(0); 494 } 495 496 static PetscErrorCode DMPlexCreateCubeMesh_Internal(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[], DMBoundaryType bdX, DMBoundaryType bdY, DMBoundaryType bdZ) 497 { 498 PetscInt markerTop = 1, faceMarkerTop = 1; 499 PetscInt markerBottom = 1, faceMarkerBottom = 1; 500 PetscInt markerFront = 1, faceMarkerFront = 1; 501 PetscInt markerBack = 1, faceMarkerBack = 1; 502 PetscInt markerRight = 1, faceMarkerRight = 1; 503 PetscInt markerLeft = 1, faceMarkerLeft = 1; 504 PetscInt dim; 505 PetscBool markerSeparate = PETSC_FALSE, cutMarker = PETSC_FALSE; 506 PetscMPIInt rank; 507 PetscErrorCode ierr; 508 509 PetscFunctionBegin; 510 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 511 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 512 ierr = DMCreateLabel(dm,"marker");CHKERRQ(ierr); 513 ierr = DMCreateLabel(dm,"Face Sets");CHKERRQ(ierr); 514 if (bdX == DM_BOUNDARY_PERIODIC || bdX == DM_BOUNDARY_TWIST || 515 bdY == DM_BOUNDARY_PERIODIC || bdY == DM_BOUNDARY_TWIST || 516 bdZ == DM_BOUNDARY_PERIODIC || bdZ == DM_BOUNDARY_TWIST) { 517 ierr = PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_periodic_cut", &cutMarker, NULL);CHKERRQ(ierr); 518 if (cutMarker) {ierr = DMCreateLabel(dm,"periodic_cut");CHKERRQ(ierr);} 519 } 520 switch (dim) { 521 case 2: 522 faceMarkerTop = 3; 523 faceMarkerBottom = 1; 524 faceMarkerRight = 2; 525 faceMarkerLeft = 4; 526 break; 527 case 3: 528 faceMarkerBottom = 1; 529 faceMarkerTop = 2; 530 faceMarkerFront = 3; 531 faceMarkerBack = 4; 532 faceMarkerRight = 5; 533 faceMarkerLeft = 6; 534 break; 535 default: 536 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Dimension %d not supported",dim); 537 break; 538 } 539 ierr = PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);CHKERRQ(ierr); 540 if (markerSeparate) { 541 markerBottom = faceMarkerBottom; 542 markerTop = faceMarkerTop; 543 markerFront = faceMarkerFront; 544 markerBack = faceMarkerBack; 545 markerRight = faceMarkerRight; 546 markerLeft = faceMarkerLeft; 547 } 548 { 549 const PetscInt numXEdges = !rank ? edges[0] : 0; 550 const PetscInt numYEdges = !rank ? edges[1] : 0; 551 const PetscInt numZEdges = !rank ? edges[2] : 0; 552 const PetscInt numXVertices = !rank ? (bdX == DM_BOUNDARY_PERIODIC || bdX == DM_BOUNDARY_TWIST ? edges[0] : edges[0]+1) : 0; 553 const PetscInt numYVertices = !rank ? (bdY == DM_BOUNDARY_PERIODIC || bdY == DM_BOUNDARY_TWIST ? edges[1] : edges[1]+1) : 0; 554 const PetscInt numZVertices = !rank ? (bdZ == DM_BOUNDARY_PERIODIC || bdZ == DM_BOUNDARY_TWIST ? edges[2] : edges[2]+1) : 0; 555 const PetscInt numCells = numXEdges*numYEdges*numZEdges; 556 const PetscInt numXFaces = numYEdges*numZEdges; 557 const PetscInt numYFaces = numXEdges*numZEdges; 558 const PetscInt numZFaces = numXEdges*numYEdges; 559 const PetscInt numTotXFaces = numXVertices*numXFaces; 560 const PetscInt numTotYFaces = numYVertices*numYFaces; 561 const PetscInt numTotZFaces = numZVertices*numZFaces; 562 const PetscInt numFaces = numTotXFaces + numTotYFaces + numTotZFaces; 563 const PetscInt numTotXEdges = numXEdges*numYVertices*numZVertices; 564 const PetscInt numTotYEdges = numYEdges*numXVertices*numZVertices; 565 const PetscInt numTotZEdges = numZEdges*numXVertices*numYVertices; 566 const PetscInt numVertices = numXVertices*numYVertices*numZVertices; 567 const PetscInt numEdges = numTotXEdges + numTotYEdges + numTotZEdges; 568 const PetscInt firstVertex = (dim == 2) ? numFaces : numCells; 569 const PetscInt firstXFace = (dim == 2) ? 0 : numCells + numVertices; 570 const PetscInt firstYFace = firstXFace + numTotXFaces; 571 const PetscInt firstZFace = firstYFace + numTotYFaces; 572 const PetscInt firstXEdge = numCells + numFaces + numVertices; 573 const PetscInt firstYEdge = firstXEdge + numTotXEdges; 574 const PetscInt firstZEdge = firstYEdge + numTotYEdges; 575 Vec coordinates; 576 PetscSection coordSection; 577 PetscScalar *coords; 578 PetscInt coordSize; 579 PetscInt v, vx, vy, vz; 580 PetscInt c, f, fx, fy, fz, e, ex, ey, ez; 581 582 ierr = DMPlexSetChart(dm, 0, numCells+numFaces+numEdges+numVertices);CHKERRQ(ierr); 583 for (c = 0; c < numCells; c++) { 584 ierr = DMPlexSetConeSize(dm, c, 6);CHKERRQ(ierr); 585 } 586 for (f = firstXFace; f < firstXFace+numFaces; ++f) { 587 ierr = DMPlexSetConeSize(dm, f, 4);CHKERRQ(ierr); 588 } 589 for (e = firstXEdge; e < firstXEdge+numEdges; ++e) { 590 ierr = DMPlexSetConeSize(dm, e, 2);CHKERRQ(ierr); 591 } 592 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 593 /* Build cells */ 594 for (fz = 0; fz < numZEdges; ++fz) { 595 for (fy = 0; fy < numYEdges; ++fy) { 596 for (fx = 0; fx < numXEdges; ++fx) { 597 PetscInt cell = (fz*numYEdges + fy)*numXEdges + fx; 598 PetscInt faceB = firstZFace + (fy*numXEdges+fx)*numZVertices + fz; 599 PetscInt faceT = firstZFace + (fy*numXEdges+fx)*numZVertices + ((fz+1)%numZVertices); 600 PetscInt faceF = firstYFace + (fz*numXEdges+fx)*numYVertices + fy; 601 PetscInt faceK = firstYFace + (fz*numXEdges+fx)*numYVertices + ((fy+1)%numYVertices); 602 PetscInt faceL = firstXFace + (fz*numYEdges+fy)*numXVertices + fx; 603 PetscInt faceR = firstXFace + (fz*numYEdges+fy)*numXVertices + ((fx+1)%numXVertices); 604 /* B, T, F, K, R, L */ 605 PetscInt ornt[6] = {-4, 0, 0, -1, 0, -4}; /* ??? */ 606 PetscInt cone[6]; 607 608 /* no boundary twisting in 3D */ 609 cone[0] = faceB; cone[1] = faceT; cone[2] = faceF; cone[3] = faceK; cone[4] = faceR; cone[5] = faceL; 610 ierr = DMPlexSetCone(dm, cell, cone);CHKERRQ(ierr); 611 ierr = DMPlexSetConeOrientation(dm, cell, ornt);CHKERRQ(ierr); 612 } 613 } 614 } 615 /* Build x faces */ 616 for (fz = 0; fz < numZEdges; ++fz) { 617 for (fy = 0; fy < numYEdges; ++fy) { 618 for (fx = 0; fx < numXVertices; ++fx) { 619 PetscInt face = firstXFace + (fz*numYEdges+fy)*numXVertices + fx; 620 PetscInt edgeL = firstZEdge + ( fy* numXVertices+fx)*numZEdges + fz; 621 PetscInt edgeR = firstZEdge + (((fy+1)%numYVertices)*numXVertices+fx)*numZEdges + fz; 622 PetscInt edgeB = firstYEdge + ( fz* numXVertices+fx)*numYEdges + fy; 623 PetscInt edgeT = firstYEdge + (((fz+1)%numZVertices)*numXVertices+fx)*numYEdges + fy; 624 PetscInt ornt[4] = {0, 0, -2, -2}; 625 PetscInt cone[4]; 626 627 if (dim == 3) { 628 /* markers */ 629 if (bdX != DM_BOUNDARY_PERIODIC) { 630 if (fx == numXVertices-1) { 631 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerRight);CHKERRQ(ierr); 632 ierr = DMSetLabelValue(dm, "marker", face, markerRight);CHKERRQ(ierr); 633 } 634 else if (fx == 0) { 635 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerLeft);CHKERRQ(ierr); 636 ierr = DMSetLabelValue(dm, "marker", face, markerLeft);CHKERRQ(ierr); 637 } 638 } 639 } 640 cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL; 641 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr); 642 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr); 643 } 644 } 645 } 646 /* Build y faces */ 647 for (fz = 0; fz < numZEdges; ++fz) { 648 for (fx = 0; fx < numXEdges; ++fx) { 649 for (fy = 0; fy < numYVertices; ++fy) { 650 PetscInt face = firstYFace + (fz*numXEdges+fx)*numYVertices + fy; 651 PetscInt edgeL = firstZEdge + (fy*numXVertices+ fx )*numZEdges + fz; 652 PetscInt edgeR = firstZEdge + (fy*numXVertices+((fx+1)%numXVertices))*numZEdges + fz; 653 PetscInt edgeB = firstXEdge + ( fz *numYVertices+fy)*numXEdges + fx; 654 PetscInt edgeT = firstXEdge + (((fz+1)%numZVertices)*numYVertices+fy)*numXEdges + fx; 655 PetscInt ornt[4] = {0, 0, -2, -2}; 656 PetscInt cone[4]; 657 658 if (dim == 3) { 659 /* markers */ 660 if (bdY != DM_BOUNDARY_PERIODIC) { 661 if (fy == numYVertices-1) { 662 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerBack);CHKERRQ(ierr); 663 ierr = DMSetLabelValue(dm, "marker", face, markerBack);CHKERRQ(ierr); 664 } 665 else if (fy == 0) { 666 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerFront);CHKERRQ(ierr); 667 ierr = DMSetLabelValue(dm, "marker", face, markerFront);CHKERRQ(ierr); 668 } 669 } 670 } 671 cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL; 672 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr); 673 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr); 674 } 675 } 676 } 677 /* Build z faces */ 678 for (fy = 0; fy < numYEdges; ++fy) { 679 for (fx = 0; fx < numXEdges; ++fx) { 680 for (fz = 0; fz < numZVertices; fz++) { 681 PetscInt face = firstZFace + (fy*numXEdges+fx)*numZVertices + fz; 682 PetscInt edgeL = firstYEdge + (fz*numXVertices+ fx )*numYEdges + fy; 683 PetscInt edgeR = firstYEdge + (fz*numXVertices+((fx+1)%numXVertices))*numYEdges + fy; 684 PetscInt edgeB = firstXEdge + (fz*numYVertices+ fy )*numXEdges + fx; 685 PetscInt edgeT = firstXEdge + (fz*numYVertices+((fy+1)%numYVertices))*numXEdges + fx; 686 PetscInt ornt[4] = {0, 0, -2, -2}; 687 PetscInt cone[4]; 688 689 if (dim == 2) { 690 if (bdX == DM_BOUNDARY_TWIST && fx == numXEdges-1) {edgeR += numYEdges-1-2*fy; ornt[1] = -2;} 691 if (bdY == DM_BOUNDARY_TWIST && fy == numYEdges-1) {edgeT += numXEdges-1-2*fx; ornt[2] = 0;} 692 } 693 else { 694 /* markers */ 695 if (bdZ != DM_BOUNDARY_PERIODIC) { 696 if (fz == numZVertices-1) { 697 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerTop);CHKERRQ(ierr); 698 ierr = DMSetLabelValue(dm, "marker", face, markerTop);CHKERRQ(ierr); 699 } 700 else if (fz == 0) { 701 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerBottom);CHKERRQ(ierr); 702 ierr = DMSetLabelValue(dm, "marker", face, markerBottom);CHKERRQ(ierr); 703 } 704 } 705 } 706 cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL; 707 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr); 708 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr); 709 } 710 } 711 } 712 /* Build Z edges*/ 713 for (vy = 0; vy < numYVertices; vy++) { 714 for (vx = 0; vx < numXVertices; vx++) { 715 for (ez = 0; ez < numZEdges; ez++) { 716 const PetscInt edge = firstZEdge + (vy*numXVertices+vx)*numZEdges + ez; 717 const PetscInt vertexB = firstVertex + ( ez *numYVertices+vy)*numXVertices + vx; 718 const PetscInt vertexT = firstVertex + (((ez+1)%numZVertices)*numYVertices+vy)*numXVertices + vx; 719 PetscInt cone[2]; 720 721 if (dim == 3) { 722 if (bdX != DM_BOUNDARY_PERIODIC) { 723 if (vx == numXVertices-1) { 724 ierr = DMSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 725 } 726 else if (vx == 0) { 727 ierr = DMSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 728 } 729 } 730 if (bdY != DM_BOUNDARY_PERIODIC) { 731 if (vy == numYVertices-1) { 732 ierr = DMSetLabelValue(dm, "marker", edge, markerBack);CHKERRQ(ierr); 733 } 734 else if (vy == 0) { 735 ierr = DMSetLabelValue(dm, "marker", edge, markerFront);CHKERRQ(ierr); 736 } 737 } 738 } 739 cone[0] = vertexB; cone[1] = vertexT; 740 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 741 } 742 } 743 } 744 /* Build Y edges*/ 745 for (vz = 0; vz < numZVertices; vz++) { 746 for (vx = 0; vx < numXVertices; vx++) { 747 for (ey = 0; ey < numYEdges; ey++) { 748 const PetscInt nextv = (dim == 2 && bdY == DM_BOUNDARY_TWIST && ey == numYEdges-1) ? (numXVertices-vx-1) : (vz*numYVertices+((ey+1)%numYVertices))*numXVertices + vx; 749 const PetscInt edge = firstYEdge + (vz*numXVertices+vx)*numYEdges + ey; 750 const PetscInt vertexF = firstVertex + (vz*numYVertices+ey)*numXVertices + vx; 751 const PetscInt vertexK = firstVertex + nextv; 752 PetscInt cone[2]; 753 754 cone[0] = vertexF; cone[1] = vertexK; 755 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 756 if (dim == 2) { 757 if ((bdX != DM_BOUNDARY_PERIODIC) && (bdX != DM_BOUNDARY_TWIST)) { 758 if (vx == numXVertices-1) { 759 ierr = DMSetLabelValue(dm, "Face Sets", edge, faceMarkerRight);CHKERRQ(ierr); 760 ierr = DMSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 761 ierr = DMSetLabelValue(dm, "marker", cone[0], markerRight);CHKERRQ(ierr); 762 if (ey == numYEdges-1) { 763 ierr = DMSetLabelValue(dm, "marker", cone[1], markerRight);CHKERRQ(ierr); 764 } 765 } else if (vx == 0) { 766 ierr = DMSetLabelValue(dm, "Face Sets", edge, faceMarkerLeft);CHKERRQ(ierr); 767 ierr = DMSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 768 ierr = DMSetLabelValue(dm, "marker", cone[0], markerLeft);CHKERRQ(ierr); 769 if (ey == numYEdges-1) { 770 ierr = DMSetLabelValue(dm, "marker", cone[1], markerLeft);CHKERRQ(ierr); 771 } 772 } 773 } else { 774 if (vx == 0 && cutMarker) { 775 ierr = DMSetLabelValue(dm, "periodic_cut", edge, 1);CHKERRQ(ierr); 776 ierr = DMSetLabelValue(dm, "periodic_cut", cone[0], 1);CHKERRQ(ierr); 777 if (ey == numYEdges-1) { 778 ierr = DMSetLabelValue(dm, "periodic_cut", cone[1], 1);CHKERRQ(ierr); 779 } 780 } 781 } 782 } else { 783 if (bdX != DM_BOUNDARY_PERIODIC) { 784 if (vx == numXVertices-1) { 785 ierr = DMSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 786 } else if (vx == 0) { 787 ierr = DMSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 788 } 789 } 790 if (bdZ != DM_BOUNDARY_PERIODIC) { 791 if (vz == numZVertices-1) { 792 ierr = DMSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 793 } else if (vz == 0) { 794 ierr = DMSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 795 } 796 } 797 } 798 } 799 } 800 } 801 /* Build X edges*/ 802 for (vz = 0; vz < numZVertices; vz++) { 803 for (vy = 0; vy < numYVertices; vy++) { 804 for (ex = 0; ex < numXEdges; ex++) { 805 const PetscInt nextv = (dim == 2 && bdX == DM_BOUNDARY_TWIST && ex == numXEdges-1) ? (numYVertices-vy-1)*numXVertices : (vz*numYVertices+vy)*numXVertices + (ex+1)%numXVertices; 806 const PetscInt edge = firstXEdge + (vz*numYVertices+vy)*numXEdges + ex; 807 const PetscInt vertexL = firstVertex + (vz*numYVertices+vy)*numXVertices + ex; 808 const PetscInt vertexR = firstVertex + nextv; 809 PetscInt cone[2]; 810 811 cone[0] = vertexL; cone[1] = vertexR; 812 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 813 if (dim == 2) { 814 if ((bdY != DM_BOUNDARY_PERIODIC) && (bdY != DM_BOUNDARY_TWIST)) { 815 if (vy == numYVertices-1) { 816 ierr = DMSetLabelValue(dm, "Face Sets", edge, faceMarkerTop);CHKERRQ(ierr); 817 ierr = DMSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 818 ierr = DMSetLabelValue(dm, "marker", cone[0], markerTop);CHKERRQ(ierr); 819 if (ex == numXEdges-1) { 820 ierr = DMSetLabelValue(dm, "marker", cone[1], markerTop);CHKERRQ(ierr); 821 } 822 } else if (vy == 0) { 823 ierr = DMSetLabelValue(dm, "Face Sets", edge, faceMarkerBottom);CHKERRQ(ierr); 824 ierr = DMSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 825 ierr = DMSetLabelValue(dm, "marker", cone[0], markerBottom);CHKERRQ(ierr); 826 if (ex == numXEdges-1) { 827 ierr = DMSetLabelValue(dm, "marker", cone[1], markerBottom);CHKERRQ(ierr); 828 } 829 } 830 } else { 831 if (vy == 0 && cutMarker) { 832 ierr = DMSetLabelValue(dm, "periodic_cut", edge, 1);CHKERRQ(ierr); 833 ierr = DMSetLabelValue(dm, "periodic_cut", cone[0], 1);CHKERRQ(ierr); 834 if (ex == numXEdges-1) { 835 ierr = DMSetLabelValue(dm, "periodic_cut", cone[1], 1);CHKERRQ(ierr); 836 } 837 } 838 } 839 } else { 840 if (bdY != DM_BOUNDARY_PERIODIC) { 841 if (vy == numYVertices-1) { 842 ierr = DMSetLabelValue(dm, "marker", edge, markerBack);CHKERRQ(ierr); 843 } 844 else if (vy == 0) { 845 ierr = DMSetLabelValue(dm, "marker", edge, markerFront);CHKERRQ(ierr); 846 } 847 } 848 if (bdZ != DM_BOUNDARY_PERIODIC) { 849 if (vz == numZVertices-1) { 850 ierr = DMSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 851 } 852 else if (vz == 0) { 853 ierr = DMSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 854 } 855 } 856 } 857 } 858 } 859 } 860 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 861 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 862 /* Build coordinates */ 863 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 864 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 865 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 866 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVertices);CHKERRQ(ierr); 867 for (v = firstVertex; v < firstVertex+numVertices; ++v) { 868 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 869 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 870 } 871 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 872 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 873 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 874 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 875 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 876 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 877 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 878 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 879 for (vz = 0; vz < numZVertices; ++vz) { 880 for (vy = 0; vy < numYVertices; ++vy) { 881 for (vx = 0; vx < numXVertices; ++vx) { 882 coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+0] = lower[0] + ((upper[0] - lower[0])/numXEdges)*vx; 883 coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+1] = lower[1] + ((upper[1] - lower[1])/numYEdges)*vy; 884 if (dim == 3) { 885 coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+2] = lower[2] + ((upper[2] - lower[2])/numZEdges)*vz; 886 } 887 } 888 } 889 } 890 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 891 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 892 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 893 } 894 PetscFunctionReturn(0); 895 } 896 897 /*@ 898 DMPlexCreateHexBoxMesh - Creates a mesh on the tensor product of unit intervals (box) using hexahedra. 899 900 Collective on MPI_Comm 901 902 Input Parameters: 903 + comm - The communicator for the DM object 904 . dim - The spatial dimension 905 . periodicX - The boundary type for the X direction 906 . periodicY - The boundary type for the Y direction 907 . periodicZ - The boundary type for the Z direction 908 - cells - The number of cells in each direction 909 910 Output Parameter: 911 . dm - The DM object 912 913 Note: Here is the numbering returned for 2 cells in each direction: 914 $ 22--8-23--9--24 915 $ | | | 916 $ 13 2 14 3 15 917 $ | | | 918 $ 19--6-20--7--21 919 $ | | | 920 $ 10 0 11 1 12 921 $ | | | 922 $ 16--4-17--5--18 923 924 Level: beginner 925 926 .keywords: DM, create 927 .seealso: DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 928 @*/ 929 PetscErrorCode DMPlexCreateHexBoxMesh(MPI_Comm comm, PetscInt dim, const PetscInt cells[], DMBoundaryType periodicX, DMBoundaryType periodicY, DMBoundaryType periodicZ, DM *dm) 930 { 931 PetscInt i; 932 PetscErrorCode ierr; 933 934 PetscFunctionBegin; 935 PetscValidPointer(dm, 7); 936 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 937 PetscValidLogicalCollectiveInt(*dm,dim,2); 938 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 939 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 940 switch (dim) { 941 case 2: 942 { 943 PetscReal lower[3] = {0.0, 0.0, 0.0}; 944 PetscReal upper[3] = {1.0, 1.0, 0.0}; 945 PetscInt edges[3] = {cells[0], cells[1], 0}; 946 947 ierr = DMPlexCreateCubeMesh_Internal(*dm, lower, upper, edges, periodicX, periodicY, DM_BOUNDARY_NONE);CHKERRQ(ierr); 948 if (periodicX == DM_BOUNDARY_PERIODIC || periodicX == DM_BOUNDARY_TWIST || 949 periodicY == DM_BOUNDARY_PERIODIC || periodicY == DM_BOUNDARY_TWIST) { 950 PetscReal L[2]; 951 PetscReal maxCell[2]; 952 DMBoundaryType bdType[2]; 953 954 bdType[0] = periodicX; 955 bdType[1] = periodicY; 956 for (i = 0; i < dim; i++) { 957 L[i] = upper[i] - lower[i]; 958 maxCell[i] = 1.1 * (L[i] / cells[i]); 959 } 960 961 ierr = DMSetPeriodicity(*dm,maxCell,L,bdType);CHKERRQ(ierr); 962 } 963 break; 964 } 965 case 3: 966 { 967 PetscReal lower[3] = {0.0, 0.0, 0.0}; 968 PetscReal upper[3] = {1.0, 1.0, 1.0}; 969 970 ierr = DMPlexCreateCubeMesh_Internal(*dm, lower, upper, cells, periodicX, periodicY, periodicZ);CHKERRQ(ierr); 971 if (periodicX == DM_BOUNDARY_PERIODIC || periodicX == DM_BOUNDARY_TWIST || 972 periodicY == DM_BOUNDARY_PERIODIC || periodicY == DM_BOUNDARY_TWIST || 973 periodicZ == DM_BOUNDARY_PERIODIC || periodicZ == DM_BOUNDARY_TWIST) { 974 PetscReal L[3]; 975 PetscReal maxCell[3]; 976 DMBoundaryType bdType[3]; 977 978 bdType[0] = periodicX; 979 bdType[1] = periodicY; 980 bdType[2] = periodicZ; 981 for (i = 0; i < dim; i++) { 982 L[i] = upper[i] - lower[i]; 983 maxCell[i] = 1.1 * (L[i] / cells[i]); 984 } 985 986 ierr = DMSetPeriodicity(*dm,maxCell,L,bdType);CHKERRQ(ierr); 987 } 988 break; 989 } 990 default: 991 SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim); 992 } 993 PetscFunctionReturn(0); 994 } 995 996 /*@C 997 DMPlexSetOptionsPrefix - Sets the prefix used for searching for all DM options in the database. 998 999 Logically Collective on DM 1000 1001 Input Parameters: 1002 + dm - the DM context 1003 - prefix - the prefix to prepend to all option names 1004 1005 Notes: 1006 A hyphen (-) must NOT be given at the beginning of the prefix name. 1007 The first character of all runtime options is AUTOMATICALLY the hyphen. 1008 1009 Level: advanced 1010 1011 .seealso: SNESSetFromOptions() 1012 @*/ 1013 PetscErrorCode DMPlexSetOptionsPrefix(DM dm, const char prefix[]) 1014 { 1015 DM_Plex *mesh = (DM_Plex *) dm->data; 1016 PetscErrorCode ierr; 1017 1018 PetscFunctionBegin; 1019 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1020 ierr = PetscObjectSetOptionsPrefix((PetscObject) dm, prefix);CHKERRQ(ierr); 1021 ierr = PetscObjectSetOptionsPrefix((PetscObject) mesh->partitioner, prefix);CHKERRQ(ierr); 1022 PetscFunctionReturn(0); 1023 } 1024 1025 /*@ 1026 DMPlexCreateHexCylinderMesh - Creates a mesh on the tensor product of the unit interval with the circle (cylinder) using hexahedra. 1027 1028 Collective on MPI_Comm 1029 1030 Input Parameters: 1031 + comm - The communicator for the DM object 1032 . numRefine - The number of regular refinements to the basic 5 cell structure 1033 - periodicZ - The boundary type for the Z direction 1034 1035 Output Parameter: 1036 . dm - The DM object 1037 1038 Note: Here is the output numbering looking from the bottom of the cylinder: 1039 $ 17-----14 1040 $ | | 1041 $ | 2 | 1042 $ | | 1043 $ 17-----8-----7-----14 1044 $ | | | | 1045 $ | 3 | 0 | 1 | 1046 $ | | | | 1047 $ 19-----5-----6-----13 1048 $ | | 1049 $ | 4 | 1050 $ | | 1051 $ 19-----13 1052 $ 1053 $ and up through the top 1054 $ 1055 $ 18-----16 1056 $ | | 1057 $ | 2 | 1058 $ | | 1059 $ 18----10----11-----16 1060 $ | | | | 1061 $ | 3 | 0 | 1 | 1062 $ | | | | 1063 $ 20-----9----12-----15 1064 $ | | 1065 $ | 4 | 1066 $ | | 1067 $ 20-----15 1068 1069 Level: beginner 1070 1071 .keywords: DM, create 1072 .seealso: DMPlexCreateHexBoxMesh(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 1073 @*/ 1074 PetscErrorCode DMPlexCreateHexCylinderMesh(MPI_Comm comm, PetscInt numRefine, DMBoundaryType periodicZ, DM *dm) 1075 { 1076 const PetscInt dim = 3; 1077 PetscInt numCells, numVertices, r; 1078 PetscErrorCode ierr; 1079 1080 PetscFunctionBegin; 1081 PetscValidPointer(dm, 4); 1082 if (numRefine < 0) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Number of refinements %D cannot be negative", numRefine); 1083 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1084 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1085 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1086 /* Create topology */ 1087 { 1088 PetscInt cone[8], c; 1089 1090 numCells = 5; 1091 numVertices = 16; 1092 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1093 numCells *= 3; 1094 numVertices = 24; 1095 } 1096 ierr = DMPlexSetChart(*dm, 0, numCells+numVertices);CHKERRQ(ierr); 1097 for (c = 0; c < numCells; c++) {ierr = DMPlexSetConeSize(*dm, c, 8);CHKERRQ(ierr);} 1098 ierr = DMSetUp(*dm);CHKERRQ(ierr); 1099 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1100 cone[0] = 15; cone[1] = 18; cone[2] = 17; cone[3] = 16; 1101 cone[4] = 31; cone[5] = 32; cone[6] = 33; cone[7] = 34; 1102 ierr = DMPlexSetCone(*dm, 0, cone);CHKERRQ(ierr); 1103 cone[0] = 16; cone[1] = 17; cone[2] = 24; cone[3] = 23; 1104 cone[4] = 32; cone[5] = 36; cone[6] = 37; cone[7] = 33; /* 22 25 26 21 */ 1105 ierr = DMPlexSetCone(*dm, 1, cone);CHKERRQ(ierr); 1106 cone[0] = 18; cone[1] = 27; cone[2] = 24; cone[3] = 17; 1107 cone[4] = 34; cone[5] = 33; cone[6] = 37; cone[7] = 38; 1108 ierr = DMPlexSetCone(*dm, 2, cone);CHKERRQ(ierr); 1109 cone[0] = 29; cone[1] = 27; cone[2] = 18; cone[3] = 15; 1110 cone[4] = 35; cone[5] = 31; cone[6] = 34; cone[7] = 38; 1111 ierr = DMPlexSetCone(*dm, 3, cone);CHKERRQ(ierr); 1112 cone[0] = 29; cone[1] = 15; cone[2] = 16; cone[3] = 23; 1113 cone[4] = 35; cone[5] = 36; cone[6] = 32; cone[7] = 31; 1114 ierr = DMPlexSetCone(*dm, 4, cone);CHKERRQ(ierr); 1115 1116 cone[0] = 31; cone[1] = 34; cone[2] = 33; cone[3] = 32; 1117 cone[4] = 19; cone[5] = 22; cone[6] = 21; cone[7] = 20; 1118 ierr = DMPlexSetCone(*dm, 5, cone);CHKERRQ(ierr); 1119 cone[0] = 32; cone[1] = 33; cone[2] = 37; cone[3] = 36; 1120 cone[4] = 22; cone[5] = 25; cone[6] = 26; cone[7] = 21; 1121 ierr = DMPlexSetCone(*dm, 6, cone);CHKERRQ(ierr); 1122 cone[0] = 34; cone[1] = 38; cone[2] = 37; cone[3] = 33; 1123 cone[4] = 20; cone[5] = 21; cone[6] = 26; cone[7] = 28; 1124 ierr = DMPlexSetCone(*dm, 7, cone);CHKERRQ(ierr); 1125 cone[0] = 35; cone[1] = 38; cone[2] = 34; cone[3] = 31; 1126 cone[4] = 30; cone[5] = 19; cone[6] = 20; cone[7] = 28; 1127 ierr = DMPlexSetCone(*dm, 8, cone);CHKERRQ(ierr); 1128 cone[0] = 35; cone[1] = 31; cone[2] = 32; cone[3] = 36; 1129 cone[4] = 30; cone[5] = 25; cone[6] = 22; cone[7] = 19; 1130 ierr = DMPlexSetCone(*dm, 9, cone);CHKERRQ(ierr); 1131 1132 cone[0] = 19; cone[1] = 20; cone[2] = 21; cone[3] = 22; 1133 cone[4] = 15; cone[5] = 16; cone[6] = 17; cone[7] = 18; 1134 ierr = DMPlexSetCone(*dm, 10, cone);CHKERRQ(ierr); 1135 cone[0] = 22; cone[1] = 21; cone[2] = 26; cone[3] = 25; 1136 cone[4] = 16; cone[5] = 23; cone[6] = 24; cone[7] = 17; 1137 ierr = DMPlexSetCone(*dm, 11, cone);CHKERRQ(ierr); 1138 cone[0] = 20; cone[1] = 28; cone[2] = 26; cone[3] = 21; 1139 cone[4] = 18; cone[5] = 17; cone[6] = 24; cone[7] = 27; 1140 ierr = DMPlexSetCone(*dm, 12, cone);CHKERRQ(ierr); 1141 cone[0] = 30; cone[1] = 28; cone[2] = 20; cone[3] = 19; 1142 cone[4] = 29; cone[5] = 15; cone[6] = 18; cone[7] = 27; 1143 ierr = DMPlexSetCone(*dm, 13, cone);CHKERRQ(ierr); 1144 cone[0] = 30; cone[1] = 19; cone[2] = 22; cone[3] = 25; 1145 cone[4] = 29; cone[5] = 23; cone[6] = 16; cone[7] = 15; 1146 ierr = DMPlexSetCone(*dm, 14, cone);CHKERRQ(ierr); 1147 } else { 1148 cone[0] = 5; cone[1] = 8; cone[2] = 7; cone[3] = 6; 1149 cone[4] = 9; cone[5] = 12; cone[6] = 11; cone[7] = 10; 1150 ierr = DMPlexSetCone(*dm, 0, cone);CHKERRQ(ierr); 1151 cone[0] = 6; cone[1] = 7; cone[2] = 14; cone[3] = 13; 1152 cone[4] = 12; cone[5] = 15; cone[6] = 16; cone[7] = 11; 1153 ierr = DMPlexSetCone(*dm, 1, cone);CHKERRQ(ierr); 1154 cone[0] = 8; cone[1] = 17; cone[2] = 14; cone[3] = 7; 1155 cone[4] = 10; cone[5] = 11; cone[6] = 16; cone[7] = 18; 1156 ierr = DMPlexSetCone(*dm, 2, cone);CHKERRQ(ierr); 1157 cone[0] = 19; cone[1] = 17; cone[2] = 8; cone[3] = 5; 1158 cone[4] = 20; cone[5] = 9; cone[6] = 10; cone[7] = 18; 1159 ierr = DMPlexSetCone(*dm, 3, cone);CHKERRQ(ierr); 1160 cone[0] = 19; cone[1] = 5; cone[2] = 6; cone[3] = 13; 1161 cone[4] = 20; cone[5] = 15; cone[6] = 12; cone[7] = 9; 1162 ierr = DMPlexSetCone(*dm, 4, cone);CHKERRQ(ierr); 1163 } 1164 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1165 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1166 } 1167 /* Interpolate */ 1168 { 1169 DM idm = NULL; 1170 1171 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1172 ierr = DMDestroy(dm);CHKERRQ(ierr); 1173 *dm = idm; 1174 } 1175 /* Create cube geometry */ 1176 { 1177 Vec coordinates; 1178 PetscSection coordSection; 1179 PetscScalar *coords; 1180 PetscInt coordSize, v; 1181 const PetscReal dis = 1.0/PetscSqrtReal(2.0); 1182 const PetscReal ds2 = dis/2.0; 1183 1184 /* Build coordinates */ 1185 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 1186 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1187 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 1188 ierr = PetscSectionSetChart(coordSection, numCells, numCells+numVertices);CHKERRQ(ierr); 1189 for (v = numCells; v < numCells+numVertices; ++v) { 1190 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 1191 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 1192 } 1193 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1194 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1195 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 1196 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1197 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1198 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 1199 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1200 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1201 coords[0*dim+0] = -ds2; coords[0*dim+1] = -ds2; coords[0*dim+2] = 0.0; 1202 coords[1*dim+0] = ds2; coords[1*dim+1] = -ds2; coords[1*dim+2] = 0.0; 1203 coords[2*dim+0] = ds2; coords[2*dim+1] = ds2; coords[2*dim+2] = 0.0; 1204 coords[3*dim+0] = -ds2; coords[3*dim+1] = ds2; coords[3*dim+2] = 0.0; 1205 coords[4*dim+0] = -ds2; coords[4*dim+1] = -ds2; coords[4*dim+2] = 1.0; 1206 coords[5*dim+0] = -ds2; coords[5*dim+1] = ds2; coords[5*dim+2] = 1.0; 1207 coords[6*dim+0] = ds2; coords[6*dim+1] = ds2; coords[6*dim+2] = 1.0; 1208 coords[7*dim+0] = ds2; coords[7*dim+1] = -ds2; coords[7*dim+2] = 1.0; 1209 coords[ 8*dim+0] = dis; coords[ 8*dim+1] = -dis; coords[ 8*dim+2] = 0.0; 1210 coords[ 9*dim+0] = dis; coords[ 9*dim+1] = dis; coords[ 9*dim+2] = 0.0; 1211 coords[10*dim+0] = dis; coords[10*dim+1] = -dis; coords[10*dim+2] = 1.0; 1212 coords[11*dim+0] = dis; coords[11*dim+1] = dis; coords[11*dim+2] = 1.0; 1213 coords[12*dim+0] = -dis; coords[12*dim+1] = dis; coords[12*dim+2] = 0.0; 1214 coords[13*dim+0] = -dis; coords[13*dim+1] = dis; coords[13*dim+2] = 1.0; 1215 coords[14*dim+0] = -dis; coords[14*dim+1] = -dis; coords[14*dim+2] = 0.0; 1216 coords[15*dim+0] = -dis; coords[15*dim+1] = -dis; coords[15*dim+2] = 1.0; 1217 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1218 /* 15 31 19 */ coords[16*dim+0] = -ds2; coords[16*dim+1] = -ds2; coords[16*dim+2] = 0.5; 1219 /* 16 32 22 */ coords[17*dim+0] = ds2; coords[17*dim+1] = -ds2; coords[17*dim+2] = 0.5; 1220 /* 17 33 21 */ coords[18*dim+0] = ds2; coords[18*dim+1] = ds2; coords[18*dim+2] = 0.5; 1221 /* 18 34 20 */ coords[19*dim+0] = -ds2; coords[19*dim+1] = ds2; coords[19*dim+2] = 0.5; 1222 /* 29 35 30 */ coords[20*dim+0] = -dis; coords[20*dim+1] = -dis; coords[20*dim+2] = 0.5; 1223 /* 23 36 25 */ coords[21*dim+0] = dis; coords[21*dim+1] = -dis; coords[21*dim+2] = 0.5; 1224 /* 24 37 26 */ coords[22*dim+0] = dis; coords[22*dim+1] = dis; coords[22*dim+2] = 0.5; 1225 /* 27 38 28 */ coords[23*dim+0] = -dis; coords[23*dim+1] = dis; coords[23*dim+2] = 0.5; 1226 } 1227 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1228 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 1229 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1230 } 1231 /* Create periodicity */ 1232 if (periodicZ == DM_BOUNDARY_PERIODIC || periodicZ == DM_BOUNDARY_TWIST) { 1233 PetscReal L[3]; 1234 PetscReal maxCell[3]; 1235 DMBoundaryType bdType[3]; 1236 PetscReal lower[3] = {0.0, 0.0, 0.0}; 1237 PetscReal upper[3] = {1.0, 1.0, 1.5}; 1238 PetscInt i, numZCells = 3; 1239 1240 bdType[0] = DM_BOUNDARY_NONE; 1241 bdType[1] = DM_BOUNDARY_NONE; 1242 bdType[2] = periodicZ; 1243 for (i = 0; i < dim; i++) { 1244 L[i] = upper[i] - lower[i]; 1245 maxCell[i] = 1.1 * (L[i] / numZCells); 1246 } 1247 ierr = DMSetPeriodicity(*dm, maxCell, L, bdType);CHKERRQ(ierr); 1248 } 1249 /* Refine topology */ 1250 for (r = 0; r < numRefine; ++r) { 1251 DM rdm = NULL; 1252 1253 ierr = DMRefine(*dm, comm, &rdm);CHKERRQ(ierr); 1254 ierr = DMDestroy(dm);CHKERRQ(ierr); 1255 *dm = rdm; 1256 } 1257 /* Remap geometry to cylinder 1258 Interior square: Linear interpolation is correct 1259 The other cells all have vertices on rays from the origin. We want to uniformly expand the spacing 1260 such that the last vertex is on the unit circle. So the closest and farthest vertices are at distance 1261 1262 phi = arctan(y/x) 1263 d_close = sqrt(1/8 + 1/4 sin^2(phi)) 1264 d_far = sqrt(1/2 + sin^2(phi)) 1265 1266 so we remap them using 1267 1268 x_new = x_close + (x - x_close) (1 - d_close) / (d_far - d_close) 1269 y_new = y_close + (y - y_close) (1 - d_close) / (d_far - d_close) 1270 1271 If pi/4 < phi < 3pi/4 or -3pi/4 < phi < -pi/4, then we switch x and y. 1272 */ 1273 { 1274 Vec coordinates; 1275 PetscSection coordSection; 1276 PetscScalar *coords; 1277 PetscInt vStart, vEnd, v; 1278 const PetscReal dis = 1.0/PetscSqrtReal(2.0); 1279 const PetscReal ds2 = 0.5*dis; 1280 1281 ierr = DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1282 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 1283 ierr = DMGetCoordinatesLocal(*dm, &coordinates);CHKERRQ(ierr); 1284 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1285 for (v = vStart; v < vEnd; ++v) { 1286 PetscReal phi, sinp, cosp, dc, df, x, y, xc, yc; 1287 PetscInt off; 1288 1289 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr); 1290 if ((PetscAbsScalar(coords[off+0]) <= ds2) && (PetscAbsScalar(coords[off+1]) <= ds2)) continue; 1291 x = PetscRealPart(coords[off]); 1292 y = PetscRealPart(coords[off+1]); 1293 phi = PetscAtan2Real(y, x); 1294 sinp = PetscSinReal(phi); 1295 cosp = PetscCosReal(phi); 1296 if ((PetscAbsReal(phi) > PETSC_PI/4.0) && (PetscAbsReal(phi) < 3.0*PETSC_PI/4.0)) { 1297 dc = PetscAbsReal(ds2/sinp); 1298 df = PetscAbsReal(dis/sinp); 1299 xc = ds2*x/PetscAbsScalar(y); 1300 yc = ds2*PetscSignReal(y); 1301 } else { 1302 dc = PetscAbsReal(ds2/cosp); 1303 df = PetscAbsReal(dis/cosp); 1304 xc = ds2*PetscSignReal(x); 1305 yc = ds2*y/PetscAbsScalar(x); 1306 } 1307 coords[off+0] = xc + (coords[off+0] - xc)*(1.0 - dc)/(df - dc); 1308 coords[off+1] = yc + (coords[off+1] - yc)*(1.0 - dc)/(df - dc); 1309 } 1310 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1311 if (periodicZ == DM_BOUNDARY_PERIODIC || periodicZ == DM_BOUNDARY_TWIST) { 1312 ierr = DMLocalizeCoordinates(*dm);CHKERRQ(ierr); 1313 } 1314 } 1315 PetscFunctionReturn(0); 1316 } 1317 1318 /* External function declarations here */ 1319 extern PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling); 1320 extern PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, Mat *mat); 1321 extern PetscErrorCode DMCreateDefaultSection_Plex(DM dm); 1322 extern PetscErrorCode DMCreateDefaultConstraints_Plex(DM dm); 1323 extern PetscErrorCode DMCreateMatrix_Plex(DM dm, Mat *J); 1324 extern PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm); 1325 PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm); 1326 extern PetscErrorCode DMSetUp_Plex(DM dm); 1327 extern PetscErrorCode DMDestroy_Plex(DM dm); 1328 extern PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer); 1329 extern PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer); 1330 extern PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm); 1331 1332 /* Replace dm with the contents of dmNew 1333 - Share the DM_Plex structure 1334 - Share the coordinates 1335 - Share the SF 1336 */ 1337 static PetscErrorCode DMPlexReplace_Static(DM dm, DM dmNew) 1338 { 1339 PetscSF sf; 1340 DM coordDM, coarseDM; 1341 Vec coords; 1342 const PetscReal *maxCell, *L; 1343 const DMBoundaryType *bd; 1344 PetscErrorCode ierr; 1345 1346 PetscFunctionBegin; 1347 ierr = DMGetPointSF(dmNew, &sf);CHKERRQ(ierr); 1348 ierr = DMSetPointSF(dm, sf);CHKERRQ(ierr); 1349 ierr = DMGetCoordinateDM(dmNew, &coordDM);CHKERRQ(ierr); 1350 ierr = DMGetCoordinatesLocal(dmNew, &coords);CHKERRQ(ierr); 1351 ierr = DMSetCoordinateDM(dm, coordDM);CHKERRQ(ierr); 1352 ierr = DMSetCoordinatesLocal(dm, coords);CHKERRQ(ierr); 1353 ierr = DMGetPeriodicity(dm, &maxCell, &L, &bd);CHKERRQ(ierr); 1354 if (L) {ierr = DMSetPeriodicity(dmNew, maxCell, L, bd);CHKERRQ(ierr);} 1355 ierr = DMDestroy_Plex(dm);CHKERRQ(ierr); 1356 dm->data = dmNew->data; 1357 ((DM_Plex *) dmNew->data)->refct++; 1358 dmNew->labels->refct++; 1359 if (!--(dm->labels->refct)) { 1360 DMLabelLink next = dm->labels->next; 1361 1362 /* destroy the labels */ 1363 while (next) { 1364 DMLabelLink tmp = next->next; 1365 1366 ierr = DMLabelDestroy(&next->label);CHKERRQ(ierr); 1367 ierr = PetscFree(next);CHKERRQ(ierr); 1368 next = tmp; 1369 } 1370 ierr = PetscFree(dm->labels);CHKERRQ(ierr); 1371 } 1372 dm->labels = dmNew->labels; 1373 dm->depthLabel = dmNew->depthLabel; 1374 ierr = DMGetCoarseDM(dmNew,&coarseDM);CHKERRQ(ierr); 1375 ierr = DMSetCoarseDM(dm,coarseDM);CHKERRQ(ierr); 1376 PetscFunctionReturn(0); 1377 } 1378 1379 /* Swap dm with the contents of dmNew 1380 - Swap the DM_Plex structure 1381 - Swap the coordinates 1382 - Swap the point PetscSF 1383 */ 1384 static PetscErrorCode DMPlexSwap_Static(DM dmA, DM dmB) 1385 { 1386 DM coordDMA, coordDMB; 1387 Vec coordsA, coordsB; 1388 PetscSF sfA, sfB; 1389 void *tmp; 1390 DMLabelLinkList listTmp; 1391 DMLabel depthTmp; 1392 PetscErrorCode ierr; 1393 1394 PetscFunctionBegin; 1395 ierr = DMGetPointSF(dmA, &sfA);CHKERRQ(ierr); 1396 ierr = DMGetPointSF(dmB, &sfB);CHKERRQ(ierr); 1397 ierr = PetscObjectReference((PetscObject) sfA);CHKERRQ(ierr); 1398 ierr = DMSetPointSF(dmA, sfB);CHKERRQ(ierr); 1399 ierr = DMSetPointSF(dmB, sfA);CHKERRQ(ierr); 1400 ierr = PetscObjectDereference((PetscObject) sfA);CHKERRQ(ierr); 1401 1402 ierr = DMGetCoordinateDM(dmA, &coordDMA);CHKERRQ(ierr); 1403 ierr = DMGetCoordinateDM(dmB, &coordDMB);CHKERRQ(ierr); 1404 ierr = PetscObjectReference((PetscObject) coordDMA);CHKERRQ(ierr); 1405 ierr = DMSetCoordinateDM(dmA, coordDMB);CHKERRQ(ierr); 1406 ierr = DMSetCoordinateDM(dmB, coordDMA);CHKERRQ(ierr); 1407 ierr = PetscObjectDereference((PetscObject) coordDMA);CHKERRQ(ierr); 1408 1409 ierr = DMGetCoordinatesLocal(dmA, &coordsA);CHKERRQ(ierr); 1410 ierr = DMGetCoordinatesLocal(dmB, &coordsB);CHKERRQ(ierr); 1411 ierr = PetscObjectReference((PetscObject) coordsA);CHKERRQ(ierr); 1412 ierr = DMSetCoordinatesLocal(dmA, coordsB);CHKERRQ(ierr); 1413 ierr = DMSetCoordinatesLocal(dmB, coordsA);CHKERRQ(ierr); 1414 ierr = PetscObjectDereference((PetscObject) coordsA);CHKERRQ(ierr); 1415 1416 tmp = dmA->data; 1417 dmA->data = dmB->data; 1418 dmB->data = tmp; 1419 listTmp = dmA->labels; 1420 dmA->labels = dmB->labels; 1421 dmB->labels = listTmp; 1422 depthTmp = dmA->depthLabel; 1423 dmA->depthLabel = dmB->depthLabel; 1424 dmB->depthLabel = depthTmp; 1425 PetscFunctionReturn(0); 1426 } 1427 1428 PetscErrorCode DMSetFromOptions_NonRefinement_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 1429 { 1430 DM_Plex *mesh = (DM_Plex*) dm->data; 1431 PetscErrorCode ierr; 1432 1433 PetscFunctionBegin; 1434 /* Handle viewing */ 1435 ierr = PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMView", PETSC_FALSE, &mesh->printSetValues, NULL);CHKERRQ(ierr); 1436 ierr = PetscOptionsInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMView", 0, &mesh->printFEM, NULL);CHKERRQ(ierr); 1437 ierr = PetscOptionsReal("-dm_plex_print_tol", "Tolerance for FEM output", "DMView", mesh->printTol, &mesh->printTol, NULL);CHKERRQ(ierr); 1438 /* Point Location */ 1439 ierr = PetscOptionsBool("-dm_plex_hash_location", "Use grid hashing for point location", "DMView", PETSC_FALSE, &mesh->useHashLocation, NULL);CHKERRQ(ierr); 1440 /* Generation and remeshing */ 1441 ierr = PetscOptionsBool("-dm_plex_remesh_bd", "Allow changes to the boundary on remeshing", "DMView", PETSC_FALSE, &mesh->remeshBd, NULL);CHKERRQ(ierr); 1442 /* Projection behavior */ 1443 ierr = PetscOptionsInt("-dm_plex_max_projection_height", "Maxmimum mesh point height used to project locally", "DMPlexSetMaxProjectionHeight", 0, &mesh->maxProjectionHeight, NULL);CHKERRQ(ierr); 1444 ierr = PetscOptionsBool("-dm_plex_regular_refinement", "Use special nested projection algorithm for regular refinement", "DMPlexSetRegularRefinement", mesh->regularRefinement, &mesh->regularRefinement, NULL);CHKERRQ(ierr); 1445 1446 ierr = PetscPartitionerSetFromOptions(mesh->partitioner);CHKERRQ(ierr); 1447 PetscFunctionReturn(0); 1448 } 1449 1450 static PetscErrorCode DMSetFromOptions_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 1451 { 1452 PetscInt refine = 0, coarsen = 0, r; 1453 PetscBool isHierarchy; 1454 PetscErrorCode ierr; 1455 1456 PetscFunctionBegin; 1457 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1458 ierr = PetscOptionsHead(PetscOptionsObject,"DMPlex Options");CHKERRQ(ierr); 1459 /* Handle DMPlex refinement */ 1460 ierr = PetscOptionsInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL);CHKERRQ(ierr); 1461 ierr = PetscOptionsInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy);CHKERRQ(ierr); 1462 if (refine) {ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE);CHKERRQ(ierr);} 1463 if (refine && isHierarchy) { 1464 DM *dms, coarseDM; 1465 1466 ierr = DMGetCoarseDM(dm, &coarseDM);CHKERRQ(ierr); 1467 ierr = PetscObjectReference((PetscObject)coarseDM);CHKERRQ(ierr); 1468 ierr = PetscMalloc1(refine,&dms);CHKERRQ(ierr); 1469 ierr = DMRefineHierarchy(dm, refine, dms);CHKERRQ(ierr); 1470 /* Total hack since we do not pass in a pointer */ 1471 ierr = DMPlexSwap_Static(dm, dms[refine-1]);CHKERRQ(ierr); 1472 if (refine == 1) { 1473 ierr = DMSetCoarseDM(dm, dms[0]);CHKERRQ(ierr); 1474 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr); 1475 } else { 1476 ierr = DMSetCoarseDM(dm, dms[refine-2]);CHKERRQ(ierr); 1477 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr); 1478 ierr = DMSetCoarseDM(dms[0], dms[refine-1]);CHKERRQ(ierr); 1479 ierr = DMPlexSetRegularRefinement(dms[0], PETSC_TRUE);CHKERRQ(ierr); 1480 } 1481 ierr = DMSetCoarseDM(dms[refine-1], coarseDM);CHKERRQ(ierr); 1482 ierr = PetscObjectDereference((PetscObject)coarseDM);CHKERRQ(ierr); 1483 /* Free DMs */ 1484 for (r = 0; r < refine; ++r) { 1485 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr); 1486 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 1487 } 1488 ierr = PetscFree(dms);CHKERRQ(ierr); 1489 } else { 1490 for (r = 0; r < refine; ++r) { 1491 DM refinedMesh; 1492 1493 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 1494 ierr = DMRefine(dm, PetscObjectComm((PetscObject) dm), &refinedMesh);CHKERRQ(ierr); 1495 /* Total hack since we do not pass in a pointer */ 1496 ierr = DMPlexReplace_Static(dm, refinedMesh);CHKERRQ(ierr); 1497 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 1498 ierr = DMDestroy(&refinedMesh);CHKERRQ(ierr); 1499 } 1500 } 1501 /* Handle DMPlex coarsening */ 1502 ierr = PetscOptionsInt("-dm_coarsen", "Coarsen the mesh", "DMCreate", coarsen, &coarsen, NULL);CHKERRQ(ierr); 1503 ierr = PetscOptionsInt("-dm_coarsen_hierarchy", "The number of coarsenings", "DMCreate", coarsen, &coarsen, &isHierarchy);CHKERRQ(ierr); 1504 if (coarsen && isHierarchy) { 1505 DM *dms; 1506 1507 ierr = PetscMalloc1(coarsen, &dms);CHKERRQ(ierr); 1508 ierr = DMCoarsenHierarchy(dm, coarsen, dms);CHKERRQ(ierr); 1509 /* Free DMs */ 1510 for (r = 0; r < coarsen; ++r) { 1511 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr); 1512 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 1513 } 1514 ierr = PetscFree(dms);CHKERRQ(ierr); 1515 } else { 1516 for (r = 0; r < coarsen; ++r) { 1517 DM coarseMesh; 1518 1519 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 1520 ierr = DMCoarsen(dm, PetscObjectComm((PetscObject) dm), &coarseMesh);CHKERRQ(ierr); 1521 /* Total hack since we do not pass in a pointer */ 1522 ierr = DMPlexReplace_Static(dm, coarseMesh);CHKERRQ(ierr); 1523 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 1524 ierr = DMDestroy(&coarseMesh);CHKERRQ(ierr); 1525 } 1526 } 1527 /* Handle */ 1528 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 1529 ierr = PetscOptionsTail();CHKERRQ(ierr); 1530 PetscFunctionReturn(0); 1531 } 1532 1533 static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec) 1534 { 1535 PetscErrorCode ierr; 1536 1537 PetscFunctionBegin; 1538 ierr = DMCreateGlobalVector_Section_Private(dm,vec);CHKERRQ(ierr); 1539 /* ierr = VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM);CHKERRQ(ierr); */ 1540 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex);CHKERRQ(ierr); 1541 ierr = VecSetOperation(*vec, VECOP_VIEWNATIVE, (void (*)(void)) VecView_Plex_Native);CHKERRQ(ierr); 1542 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex);CHKERRQ(ierr); 1543 ierr = VecSetOperation(*vec, VECOP_LOADNATIVE, (void (*)(void)) VecLoad_Plex_Native);CHKERRQ(ierr); 1544 PetscFunctionReturn(0); 1545 } 1546 1547 static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec) 1548 { 1549 PetscErrorCode ierr; 1550 1551 PetscFunctionBegin; 1552 ierr = DMCreateLocalVector_Section_Private(dm,vec);CHKERRQ(ierr); 1553 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex_Local);CHKERRQ(ierr); 1554 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex_Local);CHKERRQ(ierr); 1555 PetscFunctionReturn(0); 1556 } 1557 1558 static PetscErrorCode DMGetDimPoints_Plex(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd) 1559 { 1560 PetscInt depth, d; 1561 PetscErrorCode ierr; 1562 1563 PetscFunctionBegin; 1564 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1565 if (depth == 1) { 1566 ierr = DMGetDimension(dm, &d);CHKERRQ(ierr); 1567 if (dim == 0) {ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr);} 1568 else if (dim == d) {ierr = DMPlexGetDepthStratum(dm, 1, pStart, pEnd);CHKERRQ(ierr);} 1569 else {*pStart = 0; *pEnd = 0;} 1570 } else { 1571 ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr); 1572 } 1573 PetscFunctionReturn(0); 1574 } 1575 1576 static PetscErrorCode DMGetNeighors_Plex(DM dm, PetscInt *nranks, const PetscMPIInt *ranks[]) 1577 { 1578 PetscSF sf; 1579 PetscErrorCode ierr; 1580 1581 PetscFunctionBegin; 1582 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 1583 ierr = PetscSFGetRanks(sf, nranks, ranks, NULL, NULL, NULL);CHKERRQ(ierr); 1584 PetscFunctionReturn(0); 1585 } 1586 1587 static PetscErrorCode DMInitialize_Plex(DM dm) 1588 { 1589 PetscErrorCode ierr; 1590 1591 PetscFunctionBegin; 1592 dm->ops->view = DMView_Plex; 1593 dm->ops->load = DMLoad_Plex; 1594 dm->ops->setfromoptions = DMSetFromOptions_Plex; 1595 dm->ops->clone = DMClone_Plex; 1596 dm->ops->setup = DMSetUp_Plex; 1597 dm->ops->createdefaultsection = DMCreateDefaultSection_Plex; 1598 dm->ops->createdefaultconstraints = DMCreateDefaultConstraints_Plex; 1599 dm->ops->createglobalvector = DMCreateGlobalVector_Plex; 1600 dm->ops->createlocalvector = DMCreateLocalVector_Plex; 1601 dm->ops->getlocaltoglobalmapping = NULL; 1602 dm->ops->createfieldis = NULL; 1603 dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex; 1604 dm->ops->getcoloring = NULL; 1605 dm->ops->creatematrix = DMCreateMatrix_Plex; 1606 dm->ops->createinterpolation = DMCreateInterpolation_Plex; 1607 dm->ops->getaggregates = NULL; 1608 dm->ops->getinjection = DMCreateInjection_Plex; 1609 dm->ops->refine = DMRefine_Plex; 1610 dm->ops->coarsen = DMCoarsen_Plex; 1611 dm->ops->refinehierarchy = DMRefineHierarchy_Plex; 1612 dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Plex; 1613 dm->ops->globaltolocalbegin = NULL; 1614 dm->ops->globaltolocalend = NULL; 1615 dm->ops->localtoglobalbegin = NULL; 1616 dm->ops->localtoglobalend = NULL; 1617 dm->ops->destroy = DMDestroy_Plex; 1618 dm->ops->createsubdm = DMCreateSubDM_Plex; 1619 dm->ops->getdimpoints = DMGetDimPoints_Plex; 1620 dm->ops->locatepoints = DMLocatePoints_Plex; 1621 dm->ops->projectfunctionlocal = DMProjectFunctionLocal_Plex; 1622 dm->ops->projectfunctionlabellocal = DMProjectFunctionLabelLocal_Plex; 1623 dm->ops->projectfieldlocal = DMProjectFieldLocal_Plex; 1624 dm->ops->projectfieldlabellocal = DMProjectFieldLabelLocal_Plex; 1625 dm->ops->computel2diff = DMComputeL2Diff_Plex; 1626 dm->ops->computel2gradientdiff = DMComputeL2GradientDiff_Plex; 1627 dm->ops->computel2fielddiff = DMComputeL2FieldDiff_Plex; 1628 dm->ops->getneighbors = DMGetNeighors_Plex; 1629 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMAdaptLabel_C",DMAdaptLabel_Plex);CHKERRQ(ierr); 1630 PetscFunctionReturn(0); 1631 } 1632 1633 PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm) 1634 { 1635 DM_Plex *mesh = (DM_Plex *) dm->data; 1636 PetscErrorCode ierr; 1637 1638 PetscFunctionBegin; 1639 mesh->refct++; 1640 (*newdm)->data = mesh; 1641 ierr = PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);CHKERRQ(ierr); 1642 ierr = DMInitialize_Plex(*newdm);CHKERRQ(ierr); 1643 PetscFunctionReturn(0); 1644 } 1645 1646 /*MC 1647 DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram. 1648 In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is 1649 specified by a PetscSection object. Ownership in the global representation is determined by 1650 ownership of the underlying DMPlex points. This is specified by another PetscSection object. 1651 1652 Level: intermediate 1653 1654 .seealso: DMType, DMPlexCreate(), DMCreate(), DMSetType() 1655 M*/ 1656 1657 PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm) 1658 { 1659 DM_Plex *mesh; 1660 PetscInt unit, d; 1661 PetscErrorCode ierr; 1662 1663 PetscFunctionBegin; 1664 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1665 ierr = PetscNewLog(dm,&mesh);CHKERRQ(ierr); 1666 dm->dim = 0; 1667 dm->data = mesh; 1668 1669 mesh->refct = 1; 1670 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);CHKERRQ(ierr); 1671 mesh->maxConeSize = 0; 1672 mesh->cones = NULL; 1673 mesh->coneOrientations = NULL; 1674 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);CHKERRQ(ierr); 1675 mesh->maxSupportSize = 0; 1676 mesh->supports = NULL; 1677 mesh->refinementUniform = PETSC_TRUE; 1678 mesh->refinementLimit = -1.0; 1679 1680 mesh->facesTmp = NULL; 1681 1682 mesh->tetgenOpts = NULL; 1683 mesh->triangleOpts = NULL; 1684 ierr = PetscPartitionerCreate(PetscObjectComm((PetscObject)dm), &mesh->partitioner);CHKERRQ(ierr); 1685 mesh->remeshBd = PETSC_FALSE; 1686 1687 mesh->subpointMap = NULL; 1688 1689 for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0; 1690 1691 mesh->regularRefinement = PETSC_FALSE; 1692 mesh->depthState = -1; 1693 mesh->globalVertexNumbers = NULL; 1694 mesh->globalCellNumbers = NULL; 1695 mesh->anchorSection = NULL; 1696 mesh->anchorIS = NULL; 1697 mesh->createanchors = NULL; 1698 mesh->computeanchormatrix = NULL; 1699 mesh->parentSection = NULL; 1700 mesh->parents = NULL; 1701 mesh->childIDs = NULL; 1702 mesh->childSection = NULL; 1703 mesh->children = NULL; 1704 mesh->referenceTree = NULL; 1705 mesh->getchildsymmetry = NULL; 1706 for (d = 0; d < 8; ++d) mesh->hybridPointMax[d] = PETSC_DETERMINE; 1707 mesh->vtkCellHeight = 0; 1708 mesh->useCone = PETSC_FALSE; 1709 mesh->useClosure = PETSC_TRUE; 1710 mesh->useAnchors = PETSC_FALSE; 1711 1712 mesh->maxProjectionHeight = 0; 1713 1714 mesh->printSetValues = PETSC_FALSE; 1715 mesh->printFEM = 0; 1716 mesh->printTol = 1.0e-10; 1717 1718 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 1719 PetscFunctionReturn(0); 1720 } 1721 1722 /*@ 1723 DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram. 1724 1725 Collective on MPI_Comm 1726 1727 Input Parameter: 1728 . comm - The communicator for the DMPlex object 1729 1730 Output Parameter: 1731 . mesh - The DMPlex object 1732 1733 Level: beginner 1734 1735 .keywords: DMPlex, create 1736 @*/ 1737 PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh) 1738 { 1739 PetscErrorCode ierr; 1740 1741 PetscFunctionBegin; 1742 PetscValidPointer(mesh,2); 1743 ierr = DMCreate(comm, mesh);CHKERRQ(ierr); 1744 ierr = DMSetType(*mesh, DMPLEX);CHKERRQ(ierr); 1745 PetscFunctionReturn(0); 1746 } 1747 1748 /* 1749 This takes as input the common mesh generator output, a list of the vertices for each cell, but vertex numbers are global and an SF is built for them 1750 */ 1751 static PetscErrorCode DMPlexBuildFromCellList_Parallel_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[], PetscSF *sfVert) 1752 { 1753 PetscSF sfPoint; 1754 PetscLayout vLayout; 1755 PetscHashI vhash; 1756 PetscSFNode *remoteVerticesAdj, *vertexLocal, *vertexOwner, *remoteVertex; 1757 const PetscInt *vrange; 1758 PetscInt numVerticesAdj, off, *verticesAdj, numVerticesGhost = 0, *localVertex, *cone, c, p, v, g; 1759 PETSC_UNUSED PetscHashIIter ret, iter; 1760 PetscMPIInt rank, size; 1761 PetscErrorCode ierr; 1762 1763 PetscFunctionBegin; 1764 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr); 1765 ierr = MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);CHKERRQ(ierr); 1766 /* Partition vertices */ 1767 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) dm), &vLayout);CHKERRQ(ierr); 1768 ierr = PetscLayoutSetLocalSize(vLayout, numVertices);CHKERRQ(ierr); 1769 ierr = PetscLayoutSetBlockSize(vLayout, 1);CHKERRQ(ierr); 1770 ierr = PetscLayoutSetUp(vLayout);CHKERRQ(ierr); 1771 ierr = PetscLayoutGetRanges(vLayout, &vrange);CHKERRQ(ierr); 1772 /* Count vertices and map them to procs */ 1773 PetscHashICreate(vhash); 1774 for (c = 0; c < numCells; ++c) { 1775 for (p = 0; p < numCorners; ++p) { 1776 PetscHashIPut(vhash, cells[c*numCorners+p], ret, iter); 1777 } 1778 } 1779 PetscHashISize(vhash, numVerticesAdj); 1780 ierr = PetscMalloc1(numVerticesAdj, &verticesAdj);CHKERRQ(ierr); 1781 off = 0; ierr = PetscHashIGetKeys(vhash, &off, verticesAdj);CHKERRQ(ierr); 1782 if (off != numVerticesAdj) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid number of local vertices %D should be %D", off, numVerticesAdj); 1783 ierr = PetscSortInt(numVerticesAdj, verticesAdj);CHKERRQ(ierr); 1784 ierr = PetscMalloc1(numVerticesAdj, &remoteVerticesAdj);CHKERRQ(ierr); 1785 for (v = 0; v < numVerticesAdj; ++v) { 1786 const PetscInt gv = verticesAdj[v]; 1787 PetscInt vrank; 1788 1789 ierr = PetscFindInt(gv, size+1, vrange, &vrank);CHKERRQ(ierr); 1790 vrank = vrank < 0 ? -(vrank+2) : vrank; 1791 remoteVerticesAdj[v].index = gv - vrange[vrank]; 1792 remoteVerticesAdj[v].rank = vrank; 1793 } 1794 /* Create cones */ 1795 ierr = DMPlexSetChart(dm, 0, numCells+numVerticesAdj);CHKERRQ(ierr); 1796 for (c = 0; c < numCells; ++c) {ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr);} 1797 ierr = DMSetUp(dm);CHKERRQ(ierr); 1798 ierr = DMGetWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 1799 for (c = 0; c < numCells; ++c) { 1800 for (p = 0; p < numCorners; ++p) { 1801 const PetscInt gv = cells[c*numCorners+p]; 1802 PetscInt lv; 1803 1804 ierr = PetscFindInt(gv, numVerticesAdj, verticesAdj, &lv);CHKERRQ(ierr); 1805 if (lv < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not find global vertex %D in local connectivity", gv); 1806 cone[p] = lv+numCells; 1807 } 1808 ierr = DMPlexSetCone(dm, c, cone);CHKERRQ(ierr); 1809 } 1810 ierr = DMRestoreWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 1811 /* Create SF for vertices */ 1812 ierr = PetscSFCreate(PetscObjectComm((PetscObject)dm), sfVert);CHKERRQ(ierr); 1813 ierr = PetscObjectSetName((PetscObject) *sfVert, "Vertex Ownership SF");CHKERRQ(ierr); 1814 ierr = PetscSFSetFromOptions(*sfVert);CHKERRQ(ierr); 1815 ierr = PetscSFSetGraph(*sfVert, numVertices, numVerticesAdj, NULL, PETSC_OWN_POINTER, remoteVerticesAdj, PETSC_OWN_POINTER);CHKERRQ(ierr); 1816 ierr = PetscFree(verticesAdj);CHKERRQ(ierr); 1817 /* Build pointSF */ 1818 ierr = PetscMalloc2(numVerticesAdj, &vertexLocal, numVertices, &vertexOwner);CHKERRQ(ierr); 1819 for (v = 0; v < numVerticesAdj; ++v) {vertexLocal[v].index = v+numCells; vertexLocal[v].rank = rank;} 1820 for (v = 0; v < numVertices; ++v) {vertexOwner[v].index = -1; vertexOwner[v].rank = -1;} 1821 ierr = PetscSFReduceBegin(*sfVert, MPI_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);CHKERRQ(ierr); 1822 ierr = PetscSFReduceEnd(*sfVert, MPI_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);CHKERRQ(ierr); 1823 for (v = 0; v < numVertices; ++v) if (vertexOwner[v].rank < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Global vertex %d on rank %d was unclaimed", v + vrange[rank], rank); 1824 ierr = PetscSFBcastBegin(*sfVert, MPI_2INT, vertexOwner, vertexLocal);CHKERRQ(ierr); 1825 ierr = PetscSFBcastEnd(*sfVert, MPI_2INT, vertexOwner, vertexLocal);CHKERRQ(ierr); 1826 for (v = 0; v < numVerticesAdj; ++v) if (vertexLocal[v].rank != rank) ++numVerticesGhost; 1827 ierr = PetscMalloc1(numVerticesGhost, &localVertex);CHKERRQ(ierr); 1828 ierr = PetscMalloc1(numVerticesGhost, &remoteVertex);CHKERRQ(ierr); 1829 for (v = 0, g = 0; v < numVerticesAdj; ++v) { 1830 if (vertexLocal[v].rank != rank) { 1831 localVertex[g] = v+numCells; 1832 remoteVertex[g].index = vertexLocal[v].index; 1833 remoteVertex[g].rank = vertexLocal[v].rank; 1834 ++g; 1835 } 1836 } 1837 ierr = PetscFree2(vertexLocal, vertexOwner);CHKERRQ(ierr); 1838 if (g != numVerticesGhost) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of vertex ghosts %D should be %D", g, numVerticesGhost); 1839 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 1840 ierr = PetscObjectSetName((PetscObject) sfPoint, "point SF");CHKERRQ(ierr); 1841 ierr = PetscSFSetGraph(sfPoint, numCells+numVerticesAdj, numVerticesGhost, localVertex, PETSC_OWN_POINTER, remoteVertex, PETSC_OWN_POINTER);CHKERRQ(ierr); 1842 ierr = PetscLayoutDestroy(&vLayout);CHKERRQ(ierr); 1843 PetscHashIDestroy(vhash); 1844 /* Fill in the rest of the topology structure */ 1845 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 1846 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 1847 PetscFunctionReturn(0); 1848 } 1849 1850 /* 1851 This takes as input the coordinates for each owned vertex 1852 */ 1853 static PetscErrorCode DMPlexBuildCoordinates_Parallel_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numV, PetscSF sfVert, const PetscReal vertexCoords[]) 1854 { 1855 PetscSection coordSection; 1856 Vec coordinates; 1857 PetscScalar *coords; 1858 PetscInt numVertices, numVerticesAdj, coordSize, v; 1859 PetscErrorCode ierr; 1860 1861 PetscFunctionBegin; 1862 ierr = PetscSFGetGraph(sfVert, &numVertices, &numVerticesAdj, NULL, NULL);CHKERRQ(ierr); 1863 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1864 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1865 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 1866 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVerticesAdj);CHKERRQ(ierr); 1867 for (v = numCells; v < numCells+numVerticesAdj; ++v) { 1868 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 1869 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 1870 } 1871 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1872 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1873 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 1874 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 1875 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1876 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1877 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1878 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1879 { 1880 MPI_Datatype coordtype; 1881 1882 /* Need a temp buffer for coords if we have complex/single */ 1883 ierr = MPI_Type_contiguous(spaceDim, MPIU_SCALAR, &coordtype);CHKERRQ(ierr); 1884 ierr = MPI_Type_commit(&coordtype);CHKERRQ(ierr); 1885 #if defined(PETSC_USE_COMPLEX) 1886 { 1887 PetscScalar *svertexCoords; 1888 PetscInt i; 1889 ierr = PetscMalloc1(numV*spaceDim,&svertexCoords);CHKERRQ(ierr); 1890 for (i=0; i<numV*spaceDim; i++) svertexCoords[i] = vertexCoords[i]; 1891 ierr = PetscSFBcastBegin(sfVert, coordtype, svertexCoords, coords);CHKERRQ(ierr); 1892 ierr = PetscSFBcastEnd(sfVert, coordtype, svertexCoords, coords);CHKERRQ(ierr); 1893 ierr = PetscFree(svertexCoords);CHKERRQ(ierr); 1894 } 1895 #else 1896 ierr = PetscSFBcastBegin(sfVert, coordtype, vertexCoords, coords);CHKERRQ(ierr); 1897 ierr = PetscSFBcastEnd(sfVert, coordtype, vertexCoords, coords);CHKERRQ(ierr); 1898 #endif 1899 ierr = MPI_Type_free(&coordtype);CHKERRQ(ierr); 1900 } 1901 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1902 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 1903 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1904 PetscFunctionReturn(0); 1905 } 1906 1907 /*@C 1908 DMPlexCreateFromCellListParallel - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM 1909 1910 Input Parameters: 1911 + comm - The communicator 1912 . dim - The topological dimension of the mesh 1913 . numCells - The number of cells owned by this process 1914 . numVertices - The number of vertices owned by this process 1915 . numCorners - The number of vertices for each cell 1916 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 1917 . cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 1918 . spaceDim - The spatial dimension used for coordinates 1919 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 1920 1921 Output Parameter: 1922 + dm - The DM 1923 - vertexSF - Optional, SF describing complete vertex ownership 1924 1925 Note: Two triangles sharing a face 1926 $ 1927 $ 2 1928 $ / | \ 1929 $ / | \ 1930 $ / | \ 1931 $ 0 0 | 1 3 1932 $ \ | / 1933 $ \ | / 1934 $ \ | / 1935 $ 1 1936 would have input 1937 $ numCells = 2, numVertices = 4 1938 $ cells = [0 1 2 1 3 2] 1939 $ 1940 which would result in the DMPlex 1941 $ 1942 $ 4 1943 $ / | \ 1944 $ / | \ 1945 $ / | \ 1946 $ 2 0 | 1 5 1947 $ \ | / 1948 $ \ | / 1949 $ \ | / 1950 $ 3 1951 1952 Level: beginner 1953 1954 .seealso: DMPlexCreateFromCellList(), DMPlexCreateFromDAG(), DMPlexCreate() 1955 @*/ 1956 PetscErrorCode DMPlexCreateFromCellListParallel(MPI_Comm comm, PetscInt dim, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, PetscBool interpolate, const int cells[], PetscInt spaceDim, const PetscReal vertexCoords[], PetscSF *vertexSF, DM *dm) 1957 { 1958 PetscSF sfVert; 1959 PetscErrorCode ierr; 1960 1961 PetscFunctionBegin; 1962 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1963 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1964 PetscValidLogicalCollectiveInt(*dm, dim, 2); 1965 PetscValidLogicalCollectiveInt(*dm, spaceDim, 8); 1966 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1967 ierr = DMPlexBuildFromCellList_Parallel_Private(*dm, numCells, numVertices, numCorners, cells, &sfVert);CHKERRQ(ierr); 1968 if (interpolate) { 1969 DM idm = NULL; 1970 1971 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1972 ierr = DMDestroy(dm);CHKERRQ(ierr); 1973 *dm = idm; 1974 } 1975 ierr = DMPlexBuildCoordinates_Parallel_Private(*dm, spaceDim, numCells, numVertices,sfVert, vertexCoords);CHKERRQ(ierr); 1976 if (vertexSF) *vertexSF = sfVert; 1977 else {ierr = PetscSFDestroy(&sfVert);CHKERRQ(ierr);} 1978 PetscFunctionReturn(0); 1979 } 1980 1981 /* 1982 This takes as input the common mesh generator output, a list of the vertices for each cell 1983 */ 1984 static PetscErrorCode DMPlexBuildFromCellList_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[]) 1985 { 1986 PetscInt *cone, c, p; 1987 PetscErrorCode ierr; 1988 1989 PetscFunctionBegin; 1990 ierr = DMPlexSetChart(dm, 0, numCells+numVertices);CHKERRQ(ierr); 1991 for (c = 0; c < numCells; ++c) { 1992 ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr); 1993 } 1994 ierr = DMSetUp(dm);CHKERRQ(ierr); 1995 ierr = DMGetWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 1996 for (c = 0; c < numCells; ++c) { 1997 for (p = 0; p < numCorners; ++p) { 1998 cone[p] = cells[c*numCorners+p]+numCells; 1999 } 2000 ierr = DMPlexSetCone(dm, c, cone);CHKERRQ(ierr); 2001 } 2002 ierr = DMRestoreWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 2003 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 2004 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 2005 PetscFunctionReturn(0); 2006 } 2007 2008 /* 2009 This takes as input the coordinates for each vertex 2010 */ 2011 static PetscErrorCode DMPlexBuildCoordinates_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, const double vertexCoords[]) 2012 { 2013 PetscSection coordSection; 2014 Vec coordinates; 2015 DM cdm; 2016 PetscScalar *coords; 2017 PetscInt v, d; 2018 PetscErrorCode ierr; 2019 2020 PetscFunctionBegin; 2021 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2022 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 2023 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 2024 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVertices);CHKERRQ(ierr); 2025 for (v = numCells; v < numCells+numVertices; ++v) { 2026 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 2027 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 2028 } 2029 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 2030 2031 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 2032 ierr = DMCreateLocalVector(cdm, &coordinates);CHKERRQ(ierr); 2033 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 2034 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 2035 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2036 for (v = 0; v < numVertices; ++v) { 2037 for (d = 0; d < spaceDim; ++d) { 2038 coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d]; 2039 } 2040 } 2041 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2042 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 2043 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 2044 PetscFunctionReturn(0); 2045 } 2046 2047 /*@C 2048 DMPlexCreateFromCellList - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM 2049 2050 Input Parameters: 2051 + comm - The communicator 2052 . dim - The topological dimension of the mesh 2053 . numCells - The number of cells 2054 . numVertices - The number of vertices 2055 . numCorners - The number of vertices for each cell 2056 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 2057 . cells - An array of numCells*numCorners numbers, the vertices for each cell 2058 . spaceDim - The spatial dimension used for coordinates 2059 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 2060 2061 Output Parameter: 2062 . dm - The DM 2063 2064 Note: Two triangles sharing a face 2065 $ 2066 $ 2 2067 $ / | \ 2068 $ / | \ 2069 $ / | \ 2070 $ 0 0 | 1 3 2071 $ \ | / 2072 $ \ | / 2073 $ \ | / 2074 $ 1 2075 would have input 2076 $ numCells = 2, numVertices = 4 2077 $ cells = [0 1 2 1 3 2] 2078 $ 2079 which would result in the DMPlex 2080 $ 2081 $ 4 2082 $ / | \ 2083 $ / | \ 2084 $ / | \ 2085 $ 2 0 | 1 5 2086 $ \ | / 2087 $ \ | / 2088 $ \ | / 2089 $ 3 2090 2091 Level: beginner 2092 2093 .seealso: DMPlexCreateFromDAG(), DMPlexCreate() 2094 @*/ 2095 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) 2096 { 2097 PetscErrorCode ierr; 2098 2099 PetscFunctionBegin; 2100 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 2101 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 2102 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 2103 ierr = DMPlexBuildFromCellList_Private(*dm, numCells, numVertices, numCorners, cells);CHKERRQ(ierr); 2104 if (interpolate) { 2105 DM idm = NULL; 2106 2107 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 2108 ierr = DMDestroy(dm);CHKERRQ(ierr); 2109 *dm = idm; 2110 } 2111 ierr = DMPlexBuildCoordinates_Private(*dm, spaceDim, numCells, numVertices, vertexCoords);CHKERRQ(ierr); 2112 PetscFunctionReturn(0); 2113 } 2114 2115 /*@ 2116 DMPlexCreateFromDAG - This takes as input the adjacency-list representation of the Directed Acyclic Graph (Hasse Diagram) encoding a mesh, and produces a DM 2117 2118 Input Parameters: 2119 + dm - The empty DM object, usually from DMCreate() and DMSetDimension() 2120 . depth - The depth of the DAG 2121 . numPoints - The number of points at each depth 2122 . coneSize - The cone size of each point 2123 . cones - The concatenation of the cone points for each point, the cone list must be oriented correctly for each point 2124 . coneOrientations - The orientation of each cone point 2125 - vertexCoords - An array of numVertices*dim numbers, the coordinates of each vertex 2126 2127 Output Parameter: 2128 . dm - The DM 2129 2130 Note: Two triangles sharing a face would have input 2131 $ depth = 1, numPoints = [4 2], coneSize = [3 3 0 0 0 0] 2132 $ cones = [2 3 4 3 5 4], coneOrientations = [0 0 0 0 0 0] 2133 $ vertexCoords = [-1.0 0.0 0.0 -1.0 0.0 1.0 1.0 0.0] 2134 $ 2135 which would result in the DMPlex 2136 $ 2137 $ 4 2138 $ / | \ 2139 $ / | \ 2140 $ / | \ 2141 $ 2 0 | 1 5 2142 $ \ | / 2143 $ \ | / 2144 $ \ | / 2145 $ 3 2146 $ 2147 $ Notice that all points are numbered consecutively, unlikely DMPlexCreateFromCellList() 2148 2149 Level: advanced 2150 2151 .seealso: DMPlexCreateFromCellList(), DMPlexCreate() 2152 @*/ 2153 PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[]) 2154 { 2155 Vec coordinates; 2156 PetscSection coordSection; 2157 PetscScalar *coords; 2158 PetscInt coordSize, firstVertex = -1, pStart = 0, pEnd = 0, p, v, dim, dimEmbed, d, off; 2159 PetscErrorCode ierr; 2160 2161 PetscFunctionBegin; 2162 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2163 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 2164 if (dimEmbed < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Embedding dimension %d cannot be less than intrinsic dimension %d",dimEmbed,dim); 2165 for (d = 0; d <= depth; ++d) pEnd += numPoints[d]; 2166 ierr = DMPlexSetChart(dm, pStart, pEnd);CHKERRQ(ierr); 2167 for (p = pStart; p < pEnd; ++p) { 2168 ierr = DMPlexSetConeSize(dm, p, coneSize[p-pStart]);CHKERRQ(ierr); 2169 if (firstVertex < 0 && !coneSize[p - pStart]) { 2170 firstVertex = p - pStart; 2171 } 2172 } 2173 if (firstVertex < 0 && numPoints[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected %d vertices but could not find any", numPoints[0]); 2174 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 2175 for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) { 2176 ierr = DMPlexSetCone(dm, p, &cones[off]);CHKERRQ(ierr); 2177 ierr = DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);CHKERRQ(ierr); 2178 } 2179 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 2180 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 2181 /* Build coordinates */ 2182 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2183 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 2184 ierr = PetscSectionSetFieldComponents(coordSection, 0, dimEmbed);CHKERRQ(ierr); 2185 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0]);CHKERRQ(ierr); 2186 for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) { 2187 ierr = PetscSectionSetDof(coordSection, v, dimEmbed);CHKERRQ(ierr); 2188 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed);CHKERRQ(ierr); 2189 } 2190 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 2191 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 2192 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 2193 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 2194 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 2195 ierr = VecSetBlockSize(coordinates, dimEmbed);CHKERRQ(ierr); 2196 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 2197 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2198 for (v = 0; v < numPoints[0]; ++v) { 2199 PetscInt off; 2200 2201 ierr = PetscSectionGetOffset(coordSection, v+firstVertex, &off);CHKERRQ(ierr); 2202 for (d = 0; d < dimEmbed; ++d) { 2203 coords[off+d] = vertexCoords[v*dimEmbed+d]; 2204 } 2205 } 2206 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2207 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 2208 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 2209 PetscFunctionReturn(0); 2210 } 2211 2212 /*@C 2213 DMPlexCreateFromFile - This takes a filename and produces a DM 2214 2215 Input Parameters: 2216 + comm - The communicator 2217 . filename - A file name 2218 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 2219 2220 Output Parameter: 2221 . dm - The DM 2222 2223 Level: beginner 2224 2225 .seealso: DMPlexCreateFromDAG(), DMPlexCreateFromCellList(), DMPlexCreate() 2226 @*/ 2227 PetscErrorCode DMPlexCreateFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 2228 { 2229 const char *extGmsh = ".msh"; 2230 const char *extCGNS = ".cgns"; 2231 const char *extExodus = ".exo"; 2232 const char *extFluent = ".cas"; 2233 const char *extHDF5 = ".h5"; 2234 const char *extMed = ".med"; 2235 size_t len; 2236 PetscBool isGmsh, isCGNS, isExodus, isFluent, isHDF5, isMed; 2237 PetscMPIInt rank; 2238 PetscErrorCode ierr; 2239 2240 PetscFunctionBegin; 2241 PetscValidPointer(filename, 2); 2242 PetscValidPointer(dm, 4); 2243 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 2244 ierr = PetscStrlen(filename, &len);CHKERRQ(ierr); 2245 if (!len) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Filename must be a valid path"); 2246 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGmsh, 4, &isGmsh);CHKERRQ(ierr); 2247 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extCGNS, 5, &isCGNS);CHKERRQ(ierr); 2248 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extExodus, 4, &isExodus);CHKERRQ(ierr); 2249 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extFluent, 4, &isFluent);CHKERRQ(ierr); 2250 ierr = PetscStrncmp(&filename[PetscMax(0,len-3)], extHDF5, 3, &isHDF5);CHKERRQ(ierr); 2251 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extMed, 4, &isMed);CHKERRQ(ierr); 2252 if (isGmsh) { 2253 ierr = DMPlexCreateGmshFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2254 } else if (isCGNS) { 2255 ierr = DMPlexCreateCGNSFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2256 } else if (isExodus) { 2257 ierr = DMPlexCreateExodusFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2258 } else if (isFluent) { 2259 ierr = DMPlexCreateFluentFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2260 } else if (isHDF5) { 2261 PetscViewer viewer; 2262 2263 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 2264 ierr = PetscViewerSetType(viewer, PETSCVIEWERHDF5);CHKERRQ(ierr); 2265 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 2266 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 2267 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 2268 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 2269 ierr = DMLoad(*dm, viewer);CHKERRQ(ierr); 2270 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 2271 } else if (isMed) { 2272 ierr = DMPlexCreateMedFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2273 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot load file %s: unrecognized extension", filename); 2274 PetscFunctionReturn(0); 2275 } 2276 2277 /*@ 2278 DMPlexCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell 2279 2280 Collective on comm 2281 2282 Input Parameters: 2283 + comm - The communicator 2284 . dim - The spatial dimension 2285 - simplex - Flag for simplex, otherwise use a tensor-product cell 2286 2287 Output Parameter: 2288 . refdm - The reference cell 2289 2290 Level: intermediate 2291 2292 .keywords: reference cell 2293 .seealso: 2294 @*/ 2295 PetscErrorCode DMPlexCreateReferenceCell(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *refdm) 2296 { 2297 DM rdm; 2298 Vec coords; 2299 PetscErrorCode ierr; 2300 2301 PetscFunctionBeginUser; 2302 ierr = DMCreate(comm, &rdm);CHKERRQ(ierr); 2303 ierr = DMSetType(rdm, DMPLEX);CHKERRQ(ierr); 2304 ierr = DMSetDimension(rdm, dim);CHKERRQ(ierr); 2305 switch (dim) { 2306 case 0: 2307 { 2308 PetscInt numPoints[1] = {1}; 2309 PetscInt coneSize[1] = {0}; 2310 PetscInt cones[1] = {0}; 2311 PetscInt coneOrientations[1] = {0}; 2312 PetscScalar vertexCoords[1] = {0.0}; 2313 2314 ierr = DMPlexCreateFromDAG(rdm, 0, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2315 } 2316 break; 2317 case 1: 2318 { 2319 PetscInt numPoints[2] = {2, 1}; 2320 PetscInt coneSize[3] = {2, 0, 0}; 2321 PetscInt cones[2] = {1, 2}; 2322 PetscInt coneOrientations[2] = {0, 0}; 2323 PetscScalar vertexCoords[2] = {-1.0, 1.0}; 2324 2325 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2326 } 2327 break; 2328 case 2: 2329 if (simplex) { 2330 PetscInt numPoints[2] = {3, 1}; 2331 PetscInt coneSize[4] = {3, 0, 0, 0}; 2332 PetscInt cones[3] = {1, 2, 3}; 2333 PetscInt coneOrientations[3] = {0, 0, 0}; 2334 PetscScalar vertexCoords[6] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0}; 2335 2336 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2337 } else { 2338 PetscInt numPoints[2] = {4, 1}; 2339 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 2340 PetscInt cones[4] = {1, 2, 3, 4}; 2341 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 2342 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0}; 2343 2344 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2345 } 2346 break; 2347 case 3: 2348 if (simplex) { 2349 PetscInt numPoints[2] = {4, 1}; 2350 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 2351 PetscInt cones[4] = {1, 3, 2, 4}; 2352 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 2353 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}; 2354 2355 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2356 } else { 2357 PetscInt numPoints[2] = {8, 1}; 2358 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0}; 2359 PetscInt cones[8] = {1, 4, 3, 2, 5, 6, 7, 8}; 2360 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 2361 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, 2362 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0}; 2363 2364 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2365 } 2366 break; 2367 default: 2368 SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Cannot create reference cell for dimension %d", dim); 2369 } 2370 *refdm = NULL; 2371 ierr = DMPlexInterpolate(rdm, refdm);CHKERRQ(ierr); 2372 if (rdm->coordinateDM) { 2373 DM ncdm; 2374 PetscSection cs; 2375 PetscInt pEnd = -1; 2376 2377 ierr = DMGetDefaultSection(rdm->coordinateDM, &cs);CHKERRQ(ierr); 2378 if (cs) {ierr = PetscSectionGetChart(cs, NULL, &pEnd);CHKERRQ(ierr);} 2379 if (pEnd >= 0) { 2380 ierr = DMClone(rdm->coordinateDM, &ncdm);CHKERRQ(ierr); 2381 ierr = DMSetDefaultSection(ncdm, cs);CHKERRQ(ierr); 2382 ierr = DMSetCoordinateDM(*refdm, ncdm);CHKERRQ(ierr); 2383 ierr = DMDestroy(&ncdm);CHKERRQ(ierr); 2384 } 2385 } 2386 ierr = DMGetCoordinatesLocal(rdm, &coords);CHKERRQ(ierr); 2387 if (coords) { 2388 ierr = DMSetCoordinatesLocal(*refdm, coords);CHKERRQ(ierr); 2389 } else { 2390 ierr = DMGetCoordinates(rdm, &coords);CHKERRQ(ierr); 2391 if (coords) {ierr = DMSetCoordinates(*refdm, coords);CHKERRQ(ierr);} 2392 } 2393 ierr = DMDestroy(&rdm);CHKERRQ(ierr); 2394 PetscFunctionReturn(0); 2395 } 2396