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