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