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