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