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,PETSC_TRUE,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,PETSC_TRUE,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, PETSC_TRUE, 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 /*@ 1332 DMPlexCreateWedgeCylinderMesh - Creates a mesh on the tensor product of the unit interval with the circle (cylinder) using wedges. 1333 1334 Collective on MPI_Comm 1335 1336 Input Parameters: 1337 + comm - The communicator for the DM object 1338 . n - The number of wedges around the origin 1339 - interpolate - Create edges and faces 1340 1341 Output Parameter: 1342 . dm - The DM object 1343 1344 Level: beginner 1345 1346 .keywords: DM, create 1347 .seealso: DMPlexCreateHexCylinderMesh(), DMPlexCreateHexBoxMesh(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 1348 @*/ 1349 PetscErrorCode DMPlexCreateWedgeCylinderMesh(MPI_Comm comm, PetscInt n, PetscBool interpolate, DM *dm) 1350 { 1351 const PetscInt dim = 3; 1352 PetscInt numCells, numVertices; 1353 PetscErrorCode ierr; 1354 1355 PetscFunctionBegin; 1356 PetscValidPointer(dm, 3); 1357 if (n < 0) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Number of wedges %D cannot be negative", n); 1358 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1359 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1360 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1361 /* Create topology */ 1362 { 1363 PetscInt cone[6], c; 1364 1365 numCells = n; 1366 numVertices = 2*(n+1); 1367 ierr = DMPlexSetChart(*dm, 0, numCells+numVertices);CHKERRQ(ierr); 1368 for (c = 0; c < numCells; c++) {ierr = DMPlexSetConeSize(*dm, c, 6);CHKERRQ(ierr);} 1369 ierr = DMSetUp(*dm);CHKERRQ(ierr); 1370 for (c = 0; c < numCells; c++) { 1371 cone[0] = c+n*1; cone[1] = (c+1)%n+n*1; cone[2] = 0+3*n; 1372 cone[3] = c+n*2; cone[4] = (c+1)%n+n*2; cone[5] = 1+3*n; 1373 ierr = DMPlexSetCone(*dm, c, cone);CHKERRQ(ierr); 1374 } 1375 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1376 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1377 } 1378 /* Interpolate */ 1379 if (interpolate) { 1380 DM idm = NULL; 1381 1382 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1383 ierr = DMDestroy(dm);CHKERRQ(ierr); 1384 *dm = idm; 1385 } 1386 /* Create cylinder geometry */ 1387 { 1388 Vec coordinates; 1389 PetscSection coordSection; 1390 PetscScalar *coords; 1391 PetscInt coordSize, v, c; 1392 1393 /* Build coordinates */ 1394 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 1395 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1396 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 1397 ierr = PetscSectionSetChart(coordSection, numCells, numCells+numVertices);CHKERRQ(ierr); 1398 for (v = numCells; v < numCells+numVertices; ++v) { 1399 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 1400 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 1401 } 1402 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1403 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1404 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 1405 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1406 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1407 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 1408 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1409 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1410 for (c = 0; c < numCells; c++) { 1411 coords[(c+0*n)*dim+0] = PetscCosReal(2.0*c*PETSC_PI/n); coords[(c+0*n)*dim+1] = PetscSinReal(2.0*c*PETSC_PI/n); coords[(c+0*n)*dim+2] = 1.0; 1412 coords[(c+1*n)*dim+0] = PetscCosReal(2.0*c*PETSC_PI/n); coords[(c+1*n)*dim+1] = PetscSinReal(2.0*c*PETSC_PI/n); coords[(c+1*n)*dim+2] = 0.0; 1413 } 1414 coords[(2*n+0)*dim+0] = 0.0; coords[(2*n+0)*dim+1] = 0.0; coords[(2*n+0)*dim+2] = 1.0; 1415 coords[(2*n+1)*dim+0] = 0.0; coords[(2*n+1)*dim+1] = 0.0; coords[(2*n+1)*dim+2] = 0.0; 1416 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1417 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 1418 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1419 } 1420 PetscFunctionReturn(0); 1421 } 1422 1423 PETSC_STATIC_INLINE PetscReal DiffNormReal(PetscInt dim, const PetscReal x[], const PetscReal y[]) 1424 { 1425 PetscReal prod = 0.0; 1426 PetscInt i; 1427 for (i = 0; i < dim; ++i) prod += PetscSqr(x[i] - y[i]); 1428 return PetscSqrtReal(prod); 1429 } 1430 PETSC_STATIC_INLINE PetscReal DotReal(PetscInt dim, const PetscReal x[], const PetscReal y[]) 1431 { 1432 PetscReal prod = 0.0; 1433 PetscInt i; 1434 for (i = 0; i < dim; ++i) prod += x[i]*y[i]; 1435 return prod; 1436 } 1437 1438 #undef __FUNCT__ 1439 #define __FUNCT__ "DMPlexCreateSphereMesh" 1440 /*@ 1441 DMPlexCreateSphereMesh - Creates a mesh on the d-dimensional sphere, S^d. 1442 1443 Collective on MPI_Comm 1444 1445 Input Parameters: 1446 . comm - The communicator for the DM object 1447 . dim - The dimension 1448 - simplex - Use simplices, or tensor product cells 1449 1450 Output Parameter: 1451 . dm - The DM object 1452 1453 Level: beginner 1454 1455 .keywords: DM, create 1456 .seealso: DMPlexCreateBoxMesh(), DMPlexCreateHexBoxMesh(), DMSetType(), DMCreate() 1457 @*/ 1458 PetscErrorCode DMPlexCreateSphereMesh(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *dm) 1459 { 1460 const PetscInt embedDim = dim+1; 1461 PetscSection coordSection; 1462 Vec coordinates; 1463 PetscScalar *coords; 1464 PetscReal *coordsIn; 1465 PetscInt numCells, numEdges, numVerts, firstVertex, v, firstEdge, coordSize, d, c, e; 1466 PetscMPIInt rank; 1467 PetscErrorCode ierr; 1468 1469 PetscFunctionBegin; 1470 PetscValidPointer(dm, 4); 1471 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1472 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1473 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1474 ierr = DMSetCoordinateDim(*dm, dim+1);CHKERRQ(ierr); 1475 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) *dm), &rank);CHKERRQ(ierr); 1476 switch (dim) { 1477 case 2: 1478 if (simplex) { 1479 DM idm = NULL; 1480 const PetscReal edgeLen = 2.0/(1.0 + PETSC_PHI); 1481 const PetscReal vertex[3] = {0.0, 1.0/(1.0 + PETSC_PHI), PETSC_PHI/(1.0 + PETSC_PHI)}; 1482 const PetscInt degree = 5; 1483 PetscInt s[3] = {1, 1, 1}; 1484 PetscInt cone[3]; 1485 PetscInt *graph, p, i, j, k; 1486 1487 numCells = !rank ? 20 : 0; 1488 numEdges = !rank ? 30 : 0; 1489 numVerts = !rank ? 12 : 0; 1490 firstVertex = numCells; 1491 firstEdge = numCells + numVerts; 1492 /* Use icosahedron, which for a unit sphere has coordinates which are all cyclic permutations of 1493 1494 (0, \pm 1/\phi+1, \pm \phi/\phi+1) 1495 1496 where \phi^2 - \phi - 1 = 0, meaning \phi is the golden ratio \frac{1 + \sqrt{5}}{2}. The edge 1497 length is then given by 2/\phi = 2 * 2.73606 = 5.47214. 1498 */ 1499 /* Construct vertices */ 1500 ierr = PetscCalloc1(numVerts * embedDim, &coordsIn);CHKERRQ(ierr); 1501 for (p = 0, i = 0; p < embedDim; ++p) { 1502 for (s[1] = -1; s[1] < 2; s[1] += 2) { 1503 for (s[2] = -1; s[2] < 2; s[2] += 2) { 1504 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[(d+p)%embedDim]*vertex[(d+p)%embedDim]; 1505 ++i; 1506 } 1507 } 1508 } 1509 /* Construct graph */ 1510 ierr = PetscCalloc1(numVerts * numVerts, &graph);CHKERRQ(ierr); 1511 for (i = 0; i < numVerts; ++i) { 1512 for (j = 0, k = 0; j < numVerts; ++j) { 1513 if (PetscAbsReal(DiffNormReal(embedDim, &coordsIn[i*embedDim], &coordsIn[j*embedDim]) - edgeLen) < PETSC_SMALL) {graph[i*numVerts+j] = 1; ++k;} 1514 } 1515 if (k != degree) SETERRQ3(comm, PETSC_ERR_PLIB, "Invalid icosahedron, vertex %D degree %D != %D", i, k, degree); 1516 } 1517 /* Build Topology */ 1518 ierr = DMPlexSetChart(*dm, 0, numCells+numVerts);CHKERRQ(ierr); 1519 for (c = 0; c < numCells; c++) { 1520 ierr = DMPlexSetConeSize(*dm, c, embedDim);CHKERRQ(ierr); 1521 } 1522 ierr = DMSetUp(*dm);CHKERRQ(ierr); /* Allocate space for cones */ 1523 /* Cells */ 1524 for (i = 0, c = 0; i < numVerts; ++i) { 1525 for (j = 0; j < i; ++j) { 1526 for (k = 0; k < j; ++k) { 1527 if (graph[i*numVerts+j] && graph[j*numVerts+k] && graph[k*numVerts+i]) { 1528 cone[0] = firstVertex+i; cone[1] = firstVertex+j; cone[2] = firstVertex+k; 1529 /* Check orientation */ 1530 { 1531 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}}}; 1532 PetscReal normal[3]; 1533 PetscInt e, f; 1534 1535 for (d = 0; d < embedDim; ++d) { 1536 normal[d] = 0.0; 1537 for (e = 0; e < embedDim; ++e) { 1538 for (f = 0; f < embedDim; ++f) { 1539 normal[d] += epsilon[d][e][f]*(coordsIn[j*embedDim+e] - coordsIn[i*embedDim+e])*(coordsIn[k*embedDim+f] - coordsIn[i*embedDim+f]); 1540 } 1541 } 1542 } 1543 if (DotReal(embedDim, normal, &coordsIn[i*embedDim]) < 0) {PetscInt tmp = cone[1]; cone[1] = cone[2]; cone[2] = tmp;} 1544 } 1545 ierr = DMPlexSetCone(*dm, c++, cone);CHKERRQ(ierr); 1546 } 1547 } 1548 } 1549 } 1550 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1551 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1552 ierr = PetscFree(graph);CHKERRQ(ierr); 1553 /* Interpolate mesh */ 1554 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1555 ierr = DMDestroy(dm);CHKERRQ(ierr); 1556 *dm = idm; 1557 } else { 1558 /* 1559 12-21--13 1560 | | 1561 25 4 24 1562 | | 1563 12-25--9-16--8-24--13 1564 | | | | 1565 23 5 17 0 15 3 22 1566 | | | | 1567 10-20--6-14--7-19--11 1568 | | 1569 20 1 19 1570 | | 1571 10-18--11 1572 | | 1573 23 2 22 1574 | | 1575 12-21--13 1576 */ 1577 const PetscReal dist = 1.0/PetscSqrtReal(3.0); 1578 PetscInt cone[4], ornt[4]; 1579 1580 numCells = !rank ? 6 : 0; 1581 numEdges = !rank ? 12 : 0; 1582 numVerts = !rank ? 8 : 0; 1583 firstVertex = numCells; 1584 firstEdge = numCells + numVerts; 1585 /* Build Topology */ 1586 ierr = DMPlexSetChart(*dm, 0, numCells+numEdges+numVerts);CHKERRQ(ierr); 1587 for (c = 0; c < numCells; c++) { 1588 ierr = DMPlexSetConeSize(*dm, c, 4);CHKERRQ(ierr); 1589 } 1590 for (e = firstEdge; e < firstEdge+numEdges; ++e) { 1591 ierr = DMPlexSetConeSize(*dm, e, 2);CHKERRQ(ierr); 1592 } 1593 ierr = DMSetUp(*dm);CHKERRQ(ierr); /* Allocate space for cones */ 1594 /* Cell 0 */ 1595 cone[0] = 14; cone[1] = 15; cone[2] = 16; cone[3] = 17; 1596 ierr = DMPlexSetCone(*dm, 0, cone);CHKERRQ(ierr); 1597 ornt[0] = 0; ornt[1] = 0; ornt[2] = 0; ornt[3] = 0; 1598 ierr = DMPlexSetConeOrientation(*dm, 0, ornt);CHKERRQ(ierr); 1599 /* Cell 1 */ 1600 cone[0] = 18; cone[1] = 19; cone[2] = 14; cone[3] = 20; 1601 ierr = DMPlexSetCone(*dm, 1, cone);CHKERRQ(ierr); 1602 ornt[0] = 0; ornt[1] = 0; ornt[2] = -2; ornt[3] = 0; 1603 ierr = DMPlexSetConeOrientation(*dm, 1, ornt);CHKERRQ(ierr); 1604 /* Cell 2 */ 1605 cone[0] = 21; cone[1] = 22; cone[2] = 18; cone[3] = 23; 1606 ierr = DMPlexSetCone(*dm, 2, cone);CHKERRQ(ierr); 1607 ornt[0] = 0; ornt[1] = 0; ornt[2] = -2; ornt[3] = 0; 1608 ierr = DMPlexSetConeOrientation(*dm, 2, ornt);CHKERRQ(ierr); 1609 /* Cell 3 */ 1610 cone[0] = 19; cone[1] = 22; cone[2] = 24; cone[3] = 15; 1611 ierr = DMPlexSetCone(*dm, 3, cone);CHKERRQ(ierr); 1612 ornt[0] = -2; ornt[1] = -2; ornt[2] = 0; ornt[3] = -2; 1613 ierr = DMPlexSetConeOrientation(*dm, 3, ornt);CHKERRQ(ierr); 1614 /* Cell 4 */ 1615 cone[0] = 16; cone[1] = 24; cone[2] = 21; cone[3] = 25; 1616 ierr = DMPlexSetCone(*dm, 4, cone);CHKERRQ(ierr); 1617 ornt[0] = -2; ornt[1] = -2; ornt[2] = -2; ornt[3] = 0; 1618 ierr = DMPlexSetConeOrientation(*dm, 4, ornt);CHKERRQ(ierr); 1619 /* Cell 5 */ 1620 cone[0] = 20; cone[1] = 17; cone[2] = 25; cone[3] = 23; 1621 ierr = DMPlexSetCone(*dm, 5, cone);CHKERRQ(ierr); 1622 ornt[0] = -2; ornt[1] = -2; ornt[2] = -2; ornt[3] = -2; 1623 ierr = DMPlexSetConeOrientation(*dm, 5, ornt);CHKERRQ(ierr); 1624 /* Edges */ 1625 cone[0] = 6; cone[1] = 7; 1626 ierr = DMPlexSetCone(*dm, 14, cone);CHKERRQ(ierr); 1627 cone[0] = 7; cone[1] = 8; 1628 ierr = DMPlexSetCone(*dm, 15, cone);CHKERRQ(ierr); 1629 cone[0] = 8; cone[1] = 9; 1630 ierr = DMPlexSetCone(*dm, 16, cone);CHKERRQ(ierr); 1631 cone[0] = 9; cone[1] = 6; 1632 ierr = DMPlexSetCone(*dm, 17, cone);CHKERRQ(ierr); 1633 cone[0] = 10; cone[1] = 11; 1634 ierr = DMPlexSetCone(*dm, 18, cone);CHKERRQ(ierr); 1635 cone[0] = 11; cone[1] = 7; 1636 ierr = DMPlexSetCone(*dm, 19, cone);CHKERRQ(ierr); 1637 cone[0] = 6; cone[1] = 10; 1638 ierr = DMPlexSetCone(*dm, 20, cone);CHKERRQ(ierr); 1639 cone[0] = 12; cone[1] = 13; 1640 ierr = DMPlexSetCone(*dm, 21, cone);CHKERRQ(ierr); 1641 cone[0] = 13; cone[1] = 11; 1642 ierr = DMPlexSetCone(*dm, 22, cone);CHKERRQ(ierr); 1643 cone[0] = 10; cone[1] = 12; 1644 ierr = DMPlexSetCone(*dm, 23, cone);CHKERRQ(ierr); 1645 cone[0] = 13; cone[1] = 8; 1646 ierr = DMPlexSetCone(*dm, 24, cone);CHKERRQ(ierr); 1647 cone[0] = 12; cone[1] = 9; 1648 ierr = DMPlexSetCone(*dm, 25, cone);CHKERRQ(ierr); 1649 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1650 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1651 /* Build coordinates */ 1652 ierr = PetscCalloc1(numVerts * embedDim, &coordsIn);CHKERRQ(ierr); 1653 coordsIn[0*embedDim+0] = -dist; coordsIn[0*embedDim+1] = dist; coordsIn[0*embedDim+2] = -dist; 1654 coordsIn[1*embedDim+0] = dist; coordsIn[1*embedDim+1] = dist; coordsIn[1*embedDim+2] = -dist; 1655 coordsIn[2*embedDim+0] = dist; coordsIn[2*embedDim+1] = -dist; coordsIn[2*embedDim+2] = -dist; 1656 coordsIn[3*embedDim+0] = -dist; coordsIn[3*embedDim+1] = -dist; coordsIn[3*embedDim+2] = -dist; 1657 coordsIn[4*embedDim+0] = -dist; coordsIn[4*embedDim+1] = dist; coordsIn[4*embedDim+2] = dist; 1658 coordsIn[5*embedDim+0] = dist; coordsIn[5*embedDim+1] = dist; coordsIn[5*embedDim+2] = dist; 1659 coordsIn[6*embedDim+0] = -dist; coordsIn[6*embedDim+1] = -dist; coordsIn[6*embedDim+2] = dist; 1660 coordsIn[7*embedDim+0] = dist; coordsIn[7*embedDim+1] = -dist; coordsIn[7*embedDim+2] = dist; 1661 } 1662 break; 1663 case 3: 1664 if (simplex) { 1665 DM idm = NULL; 1666 const PetscReal edgeLen = 1.0/PETSC_PHI; 1667 const PetscReal vertexA[4] = {0.5, 0.5, 0.5, 0.5}; 1668 const PetscReal vertexB[4] = {1.0, 0.0, 0.0, 0.0}; 1669 const PetscReal vertexC[4] = {0.5, 0.5*PETSC_PHI, 0.5/PETSC_PHI, 0.0}; 1670 const PetscInt degree = 12; 1671 PetscInt s[4] = {1, 1, 1}; 1672 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}, 1673 {2, 0, 1, 3}, {2, 1, 3, 0}, {2, 3, 0, 1}, {3, 0, 2, 1}, {3, 1, 0, 2}, {3, 2, 1, 0}}; 1674 PetscInt cone[4]; 1675 PetscInt *graph, p, i, j, k, l; 1676 1677 numCells = !rank ? 600 : 0; 1678 numVerts = !rank ? 120 : 0; 1679 firstVertex = numCells; 1680 /* Use the 600-cell, which for a unit sphere has coordinates which are 1681 1682 1/2 (\pm 1, \pm 1, \pm 1, \pm 1) 16 1683 (\pm 1, 0, 0, 0) all cyclic permutations 8 1684 1/2 (\pm 1, \pm phi, \pm 1/phi, 0) all even permutations 96 1685 1686 where \phi^2 - \phi - 1 = 0, meaning \phi is the golden ratio \frac{1 + \sqrt{5}}{2}. The edge 1687 length is then given by 1/\phi = 2.73606. 1688 1689 http://buzzard.pugetsound.edu/sage-practice/ch03s03.html 1690 http://mathworld.wolfram.com/600-Cell.html 1691 */ 1692 /* Construct vertices */ 1693 ierr = PetscCalloc1(numVerts * embedDim, &coordsIn);CHKERRQ(ierr); 1694 i = 0; 1695 for (s[0] = -1; s[0] < 2; s[0] += 2) { 1696 for (s[1] = -1; s[1] < 2; s[1] += 2) { 1697 for (s[2] = -1; s[2] < 2; s[2] += 2) { 1698 for (s[3] = -1; s[3] < 2; s[3] += 2) { 1699 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[d]*vertexA[d]; 1700 ++i; 1701 } 1702 } 1703 } 1704 } 1705 for (p = 0; p < embedDim; ++p) { 1706 s[1] = s[2] = s[3] = 1; 1707 for (s[0] = -1; s[0] < 2; s[0] += 2) { 1708 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[(d+p)%embedDim]*vertexB[(d+p)%embedDim]; 1709 ++i; 1710 } 1711 } 1712 for (p = 0; p < 12; ++p) { 1713 s[3] = 1; 1714 for (s[0] = -1; s[0] < 2; s[0] += 2) { 1715 for (s[1] = -1; s[1] < 2; s[1] += 2) { 1716 for (s[2] = -1; s[2] < 2; s[2] += 2) { 1717 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[evenPerm[p][d]]*vertexC[evenPerm[p][d]]; 1718 ++i; 1719 } 1720 } 1721 } 1722 } 1723 if (i != numVerts) SETERRQ2(comm, PETSC_ERR_PLIB, "Invalid 600-cell, vertices %D != %D", i, numVerts); 1724 /* Construct graph */ 1725 ierr = PetscCalloc1(numVerts * numVerts, &graph);CHKERRQ(ierr); 1726 for (i = 0; i < numVerts; ++i) { 1727 for (j = 0, k = 0; j < numVerts; ++j) { 1728 if (PetscAbsReal(DiffNormReal(embedDim, &coordsIn[i*embedDim], &coordsIn[j*embedDim]) - edgeLen) < PETSC_SMALL) {graph[i*numVerts+j] = 1; ++k;} 1729 } 1730 if (k != degree) SETERRQ3(comm, PETSC_ERR_PLIB, "Invalid 600-cell, vertex %D degree %D != %D", i, k, degree); 1731 } 1732 /* Build Topology */ 1733 ierr = DMPlexSetChart(*dm, 0, numCells+numVerts);CHKERRQ(ierr); 1734 for (c = 0; c < numCells; c++) { 1735 ierr = DMPlexSetConeSize(*dm, c, embedDim);CHKERRQ(ierr); 1736 } 1737 ierr = DMSetUp(*dm);CHKERRQ(ierr); /* Allocate space for cones */ 1738 /* Cells */ 1739 for (i = 0, c = 0; i < numVerts; ++i) { 1740 for (j = 0; j < i; ++j) { 1741 for (k = 0; k < j; ++k) { 1742 for (l = 0; l < k; ++l) { 1743 if (graph[i*numVerts+j] && graph[j*numVerts+k] && graph[k*numVerts+i] && 1744 graph[l*numVerts+i] && graph[l*numVerts+j] && graph[l*numVerts+k]) { 1745 cone[0] = firstVertex+i; cone[1] = firstVertex+j; cone[2] = firstVertex+k; cone[3] = firstVertex+l; 1746 /* Check orientation: https://ef.gy/linear-algebra:normal-vectors-in-higher-dimensional-spaces */ 1747 { 1748 const PetscInt epsilon[4][4][4][4] = {{{{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 1749 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 1}, { 0, 0, -1, 0}}, 1750 {{0, 0, 0, 0}, { 0, 0, 0, -1}, { 0, 0, 0, 0}, { 0, 1, 0, 0}}, 1751 {{0, 0, 0, 0}, { 0, 0, 1, 0}, { 0, -1, 0, 0}, { 0, 0, 0, 0}}}, 1752 1753 {{{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, -1}, { 0, 0, 1, 0}}, 1754 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 1755 {{0, 0, 0, 1}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, {-1, 0, 0, 0}}, 1756 {{0, 0, -1, 0}, { 0, 0, 0, 0}, { 1, 0, 0, 0}, { 0, 0, 0, 0}}}, 1757 1758 {{{0, 0, 0, 0}, { 0, 0, 0, 1}, { 0, 0, 0, 0}, { 0, -1, 0, 0}}, 1759 {{0, 0, 0, -1}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 1, 0, 0, 0}}, 1760 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 1761 {{0, 1, 0, 0}, {-1, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}}, 1762 1763 {{{0, 0, 0, 0}, { 0, 0, -1, 0}, { 0, 1, 0, 0}, { 0, 0, 0, 0}}, 1764 {{0, 0, 1, 0}, { 0, 0, 0, 0}, {-1, 0, 0, 0}, { 0, 0, 0, 0}}, 1765 {{0, -1, 0, 0}, { 1, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 1766 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}}}; 1767 PetscReal normal[4]; 1768 PetscInt e, f, g; 1769 1770 for (d = 0; d < embedDim; ++d) { 1771 normal[d] = 0.0; 1772 for (e = 0; e < embedDim; ++e) { 1773 for (f = 0; f < embedDim; ++f) { 1774 for (g = 0; g < embedDim; ++g) { 1775 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]); 1776 } 1777 } 1778 } 1779 } 1780 if (DotReal(embedDim, normal, &coordsIn[i*embedDim]) < 0) {PetscInt tmp = cone[1]; cone[1] = cone[2]; cone[2] = tmp;} 1781 } 1782 ierr = DMPlexSetCone(*dm, c++, cone);CHKERRQ(ierr); 1783 } 1784 } 1785 } 1786 } 1787 } 1788 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1789 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1790 ierr = PetscFree(graph);CHKERRQ(ierr); 1791 /* Interpolate mesh */ 1792 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1793 ierr = DMDestroy(dm);CHKERRQ(ierr); 1794 *dm = idm; 1795 break; 1796 } 1797 default: SETERRQ1(comm, PETSC_ERR_SUP, "Unsupported dimension for sphere: %D", dim); 1798 } 1799 /* Create coordinates */ 1800 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 1801 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1802 ierr = PetscSectionSetFieldComponents(coordSection, 0, embedDim);CHKERRQ(ierr); 1803 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVerts);CHKERRQ(ierr); 1804 for (v = firstVertex; v < firstVertex+numVerts; ++v) { 1805 ierr = PetscSectionSetDof(coordSection, v, embedDim);CHKERRQ(ierr); 1806 ierr = PetscSectionSetFieldDof(coordSection, v, 0, embedDim);CHKERRQ(ierr); 1807 } 1808 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1809 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1810 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 1811 ierr = VecSetBlockSize(coordinates, embedDim);CHKERRQ(ierr); 1812 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1813 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1814 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1815 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1816 for (v = 0; v < numVerts; ++v) for (d = 0; d < embedDim; ++d) {coords[v*embedDim+d] = coordsIn[v*embedDim+d];} 1817 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1818 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 1819 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1820 ierr = PetscFree(coordsIn);CHKERRQ(ierr); 1821 PetscFunctionReturn(0); 1822 } 1823 1824 /* External function declarations here */ 1825 extern PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling); 1826 extern PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, Mat *mat); 1827 extern PetscErrorCode DMCreateDefaultSection_Plex(DM dm); 1828 extern PetscErrorCode DMCreateDefaultConstraints_Plex(DM dm); 1829 extern PetscErrorCode DMCreateMatrix_Plex(DM dm, Mat *J); 1830 extern PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm); 1831 PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm); 1832 extern PetscErrorCode DMSetUp_Plex(DM dm); 1833 extern PetscErrorCode DMDestroy_Plex(DM dm); 1834 extern PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer); 1835 extern PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer); 1836 extern PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm); 1837 static PetscErrorCode DMInitialize_Plex(DM dm); 1838 1839 /* Replace dm with the contents of dmNew 1840 - Share the DM_Plex structure 1841 - Share the coordinates 1842 - Share the SF 1843 */ 1844 static PetscErrorCode DMPlexReplace_Static(DM dm, DM dmNew) 1845 { 1846 PetscSF sf; 1847 DM coordDM, coarseDM; 1848 Vec coords; 1849 PetscBool isper; 1850 const PetscReal *maxCell, *L; 1851 const DMBoundaryType *bd; 1852 PetscErrorCode ierr; 1853 1854 PetscFunctionBegin; 1855 ierr = DMGetPointSF(dmNew, &sf);CHKERRQ(ierr); 1856 ierr = DMSetPointSF(dm, sf);CHKERRQ(ierr); 1857 ierr = DMGetCoordinateDM(dmNew, &coordDM);CHKERRQ(ierr); 1858 ierr = DMGetCoordinatesLocal(dmNew, &coords);CHKERRQ(ierr); 1859 ierr = DMSetCoordinateDM(dm, coordDM);CHKERRQ(ierr); 1860 ierr = DMSetCoordinatesLocal(dm, coords);CHKERRQ(ierr); 1861 ierr = DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);CHKERRQ(ierr); 1862 ierr = DMSetPeriodicity(dmNew, isper, maxCell, L, bd);CHKERRQ(ierr); 1863 ierr = DMDestroy_Plex(dm);CHKERRQ(ierr); 1864 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 1865 dm->data = dmNew->data; 1866 ((DM_Plex *) dmNew->data)->refct++; 1867 dmNew->labels->refct++; 1868 if (!--(dm->labels->refct)) { 1869 DMLabelLink next = dm->labels->next; 1870 1871 /* destroy the labels */ 1872 while (next) { 1873 DMLabelLink tmp = next->next; 1874 1875 ierr = DMLabelDestroy(&next->label);CHKERRQ(ierr); 1876 ierr = PetscFree(next);CHKERRQ(ierr); 1877 next = tmp; 1878 } 1879 ierr = PetscFree(dm->labels);CHKERRQ(ierr); 1880 } 1881 dm->labels = dmNew->labels; 1882 dm->depthLabel = dmNew->depthLabel; 1883 ierr = DMGetCoarseDM(dmNew,&coarseDM);CHKERRQ(ierr); 1884 ierr = DMSetCoarseDM(dm,coarseDM);CHKERRQ(ierr); 1885 PetscFunctionReturn(0); 1886 } 1887 1888 /* Swap dm with the contents of dmNew 1889 - Swap the DM_Plex structure 1890 - Swap the coordinates 1891 - Swap the point PetscSF 1892 */ 1893 static PetscErrorCode DMPlexSwap_Static(DM dmA, DM dmB) 1894 { 1895 DM coordDMA, coordDMB; 1896 Vec coordsA, coordsB; 1897 PetscSF sfA, sfB; 1898 void *tmp; 1899 DMLabelLinkList listTmp; 1900 DMLabel depthTmp; 1901 PetscErrorCode ierr; 1902 1903 PetscFunctionBegin; 1904 ierr = DMGetPointSF(dmA, &sfA);CHKERRQ(ierr); 1905 ierr = DMGetPointSF(dmB, &sfB);CHKERRQ(ierr); 1906 ierr = PetscObjectReference((PetscObject) sfA);CHKERRQ(ierr); 1907 ierr = DMSetPointSF(dmA, sfB);CHKERRQ(ierr); 1908 ierr = DMSetPointSF(dmB, sfA);CHKERRQ(ierr); 1909 ierr = PetscObjectDereference((PetscObject) sfA);CHKERRQ(ierr); 1910 1911 ierr = DMGetCoordinateDM(dmA, &coordDMA);CHKERRQ(ierr); 1912 ierr = DMGetCoordinateDM(dmB, &coordDMB);CHKERRQ(ierr); 1913 ierr = PetscObjectReference((PetscObject) coordDMA);CHKERRQ(ierr); 1914 ierr = DMSetCoordinateDM(dmA, coordDMB);CHKERRQ(ierr); 1915 ierr = DMSetCoordinateDM(dmB, coordDMA);CHKERRQ(ierr); 1916 ierr = PetscObjectDereference((PetscObject) coordDMA);CHKERRQ(ierr); 1917 1918 ierr = DMGetCoordinatesLocal(dmA, &coordsA);CHKERRQ(ierr); 1919 ierr = DMGetCoordinatesLocal(dmB, &coordsB);CHKERRQ(ierr); 1920 ierr = PetscObjectReference((PetscObject) coordsA);CHKERRQ(ierr); 1921 ierr = DMSetCoordinatesLocal(dmA, coordsB);CHKERRQ(ierr); 1922 ierr = DMSetCoordinatesLocal(dmB, coordsA);CHKERRQ(ierr); 1923 ierr = PetscObjectDereference((PetscObject) coordsA);CHKERRQ(ierr); 1924 1925 tmp = dmA->data; 1926 dmA->data = dmB->data; 1927 dmB->data = tmp; 1928 listTmp = dmA->labels; 1929 dmA->labels = dmB->labels; 1930 dmB->labels = listTmp; 1931 depthTmp = dmA->depthLabel; 1932 dmA->depthLabel = dmB->depthLabel; 1933 dmB->depthLabel = depthTmp; 1934 PetscFunctionReturn(0); 1935 } 1936 1937 PetscErrorCode DMSetFromOptions_NonRefinement_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 1938 { 1939 DM_Plex *mesh = (DM_Plex*) dm->data; 1940 PetscErrorCode ierr; 1941 1942 PetscFunctionBegin; 1943 /* Handle viewing */ 1944 ierr = PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMView", PETSC_FALSE, &mesh->printSetValues, NULL);CHKERRQ(ierr); 1945 ierr = PetscOptionsInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMView", 0, &mesh->printFEM, NULL);CHKERRQ(ierr); 1946 ierr = PetscOptionsReal("-dm_plex_print_tol", "Tolerance for FEM output", "DMView", mesh->printTol, &mesh->printTol, NULL);CHKERRQ(ierr); 1947 /* Point Location */ 1948 ierr = PetscOptionsBool("-dm_plex_hash_location", "Use grid hashing for point location", "DMView", PETSC_FALSE, &mesh->useHashLocation, NULL);CHKERRQ(ierr); 1949 /* Generation and remeshing */ 1950 ierr = PetscOptionsBool("-dm_plex_remesh_bd", "Allow changes to the boundary on remeshing", "DMView", PETSC_FALSE, &mesh->remeshBd, NULL);CHKERRQ(ierr); 1951 /* Projection behavior */ 1952 ierr = PetscOptionsInt("-dm_plex_max_projection_height", "Maxmimum mesh point height used to project locally", "DMPlexSetMaxProjectionHeight", 0, &mesh->maxProjectionHeight, NULL);CHKERRQ(ierr); 1953 ierr = PetscOptionsBool("-dm_plex_regular_refinement", "Use special nested projection algorithm for regular refinement", "DMPlexSetRegularRefinement", mesh->regularRefinement, &mesh->regularRefinement, NULL);CHKERRQ(ierr); 1954 1955 ierr = PetscPartitionerSetFromOptions(mesh->partitioner);CHKERRQ(ierr); 1956 PetscFunctionReturn(0); 1957 } 1958 1959 static PetscErrorCode DMSetFromOptions_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 1960 { 1961 PetscInt refine = 0, coarsen = 0, r; 1962 PetscBool isHierarchy; 1963 PetscErrorCode ierr; 1964 1965 PetscFunctionBegin; 1966 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1967 ierr = PetscOptionsHead(PetscOptionsObject,"DMPlex Options");CHKERRQ(ierr); 1968 /* Handle DMPlex refinement */ 1969 ierr = PetscOptionsInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL);CHKERRQ(ierr); 1970 ierr = PetscOptionsInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy);CHKERRQ(ierr); 1971 if (refine) {ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE);CHKERRQ(ierr);} 1972 if (refine && isHierarchy) { 1973 DM *dms, coarseDM; 1974 1975 ierr = DMGetCoarseDM(dm, &coarseDM);CHKERRQ(ierr); 1976 ierr = PetscObjectReference((PetscObject)coarseDM);CHKERRQ(ierr); 1977 ierr = PetscMalloc1(refine,&dms);CHKERRQ(ierr); 1978 ierr = DMRefineHierarchy(dm, refine, dms);CHKERRQ(ierr); 1979 /* Total hack since we do not pass in a pointer */ 1980 ierr = DMPlexSwap_Static(dm, dms[refine-1]);CHKERRQ(ierr); 1981 if (refine == 1) { 1982 ierr = DMSetCoarseDM(dm, dms[0]);CHKERRQ(ierr); 1983 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr); 1984 } else { 1985 ierr = DMSetCoarseDM(dm, dms[refine-2]);CHKERRQ(ierr); 1986 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr); 1987 ierr = DMSetCoarseDM(dms[0], dms[refine-1]);CHKERRQ(ierr); 1988 ierr = DMPlexSetRegularRefinement(dms[0], PETSC_TRUE);CHKERRQ(ierr); 1989 } 1990 ierr = DMSetCoarseDM(dms[refine-1], coarseDM);CHKERRQ(ierr); 1991 ierr = PetscObjectDereference((PetscObject)coarseDM);CHKERRQ(ierr); 1992 /* Free DMs */ 1993 for (r = 0; r < refine; ++r) { 1994 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr); 1995 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 1996 } 1997 ierr = PetscFree(dms);CHKERRQ(ierr); 1998 } else { 1999 for (r = 0; r < refine; ++r) { 2000 DM refinedMesh; 2001 2002 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2003 ierr = DMRefine(dm, PetscObjectComm((PetscObject) dm), &refinedMesh);CHKERRQ(ierr); 2004 /* Total hack since we do not pass in a pointer */ 2005 ierr = DMPlexReplace_Static(dm, refinedMesh);CHKERRQ(ierr); 2006 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2007 ierr = DMDestroy(&refinedMesh);CHKERRQ(ierr); 2008 } 2009 } 2010 /* Handle DMPlex coarsening */ 2011 ierr = PetscOptionsInt("-dm_coarsen", "Coarsen the mesh", "DMCreate", coarsen, &coarsen, NULL);CHKERRQ(ierr); 2012 ierr = PetscOptionsInt("-dm_coarsen_hierarchy", "The number of coarsenings", "DMCreate", coarsen, &coarsen, &isHierarchy);CHKERRQ(ierr); 2013 if (coarsen && isHierarchy) { 2014 DM *dms; 2015 2016 ierr = PetscMalloc1(coarsen, &dms);CHKERRQ(ierr); 2017 ierr = DMCoarsenHierarchy(dm, coarsen, dms);CHKERRQ(ierr); 2018 /* Free DMs */ 2019 for (r = 0; r < coarsen; ++r) { 2020 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr); 2021 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 2022 } 2023 ierr = PetscFree(dms);CHKERRQ(ierr); 2024 } else { 2025 for (r = 0; r < coarsen; ++r) { 2026 DM coarseMesh; 2027 2028 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2029 ierr = DMCoarsen(dm, PetscObjectComm((PetscObject) dm), &coarseMesh);CHKERRQ(ierr); 2030 /* Total hack since we do not pass in a pointer */ 2031 ierr = DMPlexReplace_Static(dm, coarseMesh);CHKERRQ(ierr); 2032 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2033 ierr = DMDestroy(&coarseMesh);CHKERRQ(ierr); 2034 } 2035 } 2036 /* Handle */ 2037 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2038 ierr = PetscOptionsTail();CHKERRQ(ierr); 2039 PetscFunctionReturn(0); 2040 } 2041 2042 static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec) 2043 { 2044 PetscErrorCode ierr; 2045 2046 PetscFunctionBegin; 2047 ierr = DMCreateGlobalVector_Section_Private(dm,vec);CHKERRQ(ierr); 2048 /* ierr = VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM);CHKERRQ(ierr); */ 2049 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex);CHKERRQ(ierr); 2050 ierr = VecSetOperation(*vec, VECOP_VIEWNATIVE, (void (*)(void)) VecView_Plex_Native);CHKERRQ(ierr); 2051 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex);CHKERRQ(ierr); 2052 ierr = VecSetOperation(*vec, VECOP_LOADNATIVE, (void (*)(void)) VecLoad_Plex_Native);CHKERRQ(ierr); 2053 PetscFunctionReturn(0); 2054 } 2055 2056 static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec) 2057 { 2058 PetscErrorCode ierr; 2059 2060 PetscFunctionBegin; 2061 ierr = DMCreateLocalVector_Section_Private(dm,vec);CHKERRQ(ierr); 2062 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex_Local);CHKERRQ(ierr); 2063 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex_Local);CHKERRQ(ierr); 2064 PetscFunctionReturn(0); 2065 } 2066 2067 static PetscErrorCode DMGetDimPoints_Plex(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd) 2068 { 2069 PetscInt depth, d; 2070 PetscErrorCode ierr; 2071 2072 PetscFunctionBegin; 2073 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 2074 if (depth == 1) { 2075 ierr = DMGetDimension(dm, &d);CHKERRQ(ierr); 2076 if (dim == 0) {ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr);} 2077 else if (dim == d) {ierr = DMPlexGetDepthStratum(dm, 1, pStart, pEnd);CHKERRQ(ierr);} 2078 else {*pStart = 0; *pEnd = 0;} 2079 } else { 2080 ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr); 2081 } 2082 PetscFunctionReturn(0); 2083 } 2084 2085 static PetscErrorCode DMGetNeighors_Plex(DM dm, PetscInt *nranks, const PetscMPIInt *ranks[]) 2086 { 2087 PetscSF sf; 2088 PetscErrorCode ierr; 2089 2090 PetscFunctionBegin; 2091 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 2092 ierr = PetscSFGetRanks(sf, nranks, ranks, NULL, NULL, NULL);CHKERRQ(ierr); 2093 PetscFunctionReturn(0); 2094 } 2095 2096 static PetscErrorCode DMInitialize_Plex(DM dm) 2097 { 2098 PetscErrorCode ierr; 2099 2100 PetscFunctionBegin; 2101 dm->ops->view = DMView_Plex; 2102 dm->ops->load = DMLoad_Plex; 2103 dm->ops->setfromoptions = DMSetFromOptions_Plex; 2104 dm->ops->clone = DMClone_Plex; 2105 dm->ops->setup = DMSetUp_Plex; 2106 dm->ops->createdefaultsection = DMCreateDefaultSection_Plex; 2107 dm->ops->createdefaultconstraints = DMCreateDefaultConstraints_Plex; 2108 dm->ops->createglobalvector = DMCreateGlobalVector_Plex; 2109 dm->ops->createlocalvector = DMCreateLocalVector_Plex; 2110 dm->ops->getlocaltoglobalmapping = NULL; 2111 dm->ops->createfieldis = NULL; 2112 dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex; 2113 dm->ops->getcoloring = NULL; 2114 dm->ops->creatematrix = DMCreateMatrix_Plex; 2115 dm->ops->createinterpolation = DMCreateInterpolation_Plex; 2116 dm->ops->getaggregates = NULL; 2117 dm->ops->getinjection = DMCreateInjection_Plex; 2118 dm->ops->refine = DMRefine_Plex; 2119 dm->ops->coarsen = DMCoarsen_Plex; 2120 dm->ops->refinehierarchy = DMRefineHierarchy_Plex; 2121 dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Plex; 2122 dm->ops->adaptlabel = DMAdaptLabel_Plex; 2123 dm->ops->adaptmetric = DMAdaptMetric_Plex; 2124 dm->ops->globaltolocalbegin = NULL; 2125 dm->ops->globaltolocalend = NULL; 2126 dm->ops->localtoglobalbegin = NULL; 2127 dm->ops->localtoglobalend = NULL; 2128 dm->ops->destroy = DMDestroy_Plex; 2129 dm->ops->createsubdm = DMCreateSubDM_Plex; 2130 dm->ops->getdimpoints = DMGetDimPoints_Plex; 2131 dm->ops->locatepoints = DMLocatePoints_Plex; 2132 dm->ops->projectfunctionlocal = DMProjectFunctionLocal_Plex; 2133 dm->ops->projectfunctionlabellocal = DMProjectFunctionLabelLocal_Plex; 2134 dm->ops->projectfieldlocal = DMProjectFieldLocal_Plex; 2135 dm->ops->projectfieldlabellocal = DMProjectFieldLabelLocal_Plex; 2136 dm->ops->computel2diff = DMComputeL2Diff_Plex; 2137 dm->ops->computel2gradientdiff = DMComputeL2GradientDiff_Plex; 2138 dm->ops->computel2fielddiff = DMComputeL2FieldDiff_Plex; 2139 dm->ops->getneighbors = DMGetNeighors_Plex; 2140 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMPlexInsertBoundaryValues_C",DMPlexInsertBoundaryValues_Plex);CHKERRQ(ierr); 2141 PetscFunctionReturn(0); 2142 } 2143 2144 PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm) 2145 { 2146 DM_Plex *mesh = (DM_Plex *) dm->data; 2147 PetscErrorCode ierr; 2148 2149 PetscFunctionBegin; 2150 mesh->refct++; 2151 (*newdm)->data = mesh; 2152 ierr = PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);CHKERRQ(ierr); 2153 ierr = DMInitialize_Plex(*newdm);CHKERRQ(ierr); 2154 PetscFunctionReturn(0); 2155 } 2156 2157 /*MC 2158 DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram. 2159 In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is 2160 specified by a PetscSection object. Ownership in the global representation is determined by 2161 ownership of the underlying DMPlex points. This is specified by another PetscSection object. 2162 2163 Level: intermediate 2164 2165 .seealso: DMType, DMPlexCreate(), DMCreate(), DMSetType() 2166 M*/ 2167 2168 PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm) 2169 { 2170 DM_Plex *mesh; 2171 PetscInt unit, d; 2172 PetscErrorCode ierr; 2173 2174 PetscFunctionBegin; 2175 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2176 ierr = PetscNewLog(dm,&mesh);CHKERRQ(ierr); 2177 dm->dim = 0; 2178 dm->data = mesh; 2179 2180 mesh->refct = 1; 2181 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);CHKERRQ(ierr); 2182 mesh->maxConeSize = 0; 2183 mesh->cones = NULL; 2184 mesh->coneOrientations = NULL; 2185 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);CHKERRQ(ierr); 2186 mesh->maxSupportSize = 0; 2187 mesh->supports = NULL; 2188 mesh->refinementUniform = PETSC_TRUE; 2189 mesh->refinementLimit = -1.0; 2190 2191 mesh->facesTmp = NULL; 2192 2193 mesh->tetgenOpts = NULL; 2194 mesh->triangleOpts = NULL; 2195 ierr = PetscPartitionerCreate(PetscObjectComm((PetscObject)dm), &mesh->partitioner);CHKERRQ(ierr); 2196 mesh->remeshBd = PETSC_FALSE; 2197 2198 mesh->subpointMap = NULL; 2199 2200 for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0; 2201 2202 mesh->regularRefinement = PETSC_FALSE; 2203 mesh->depthState = -1; 2204 mesh->globalVertexNumbers = NULL; 2205 mesh->globalCellNumbers = NULL; 2206 mesh->anchorSection = NULL; 2207 mesh->anchorIS = NULL; 2208 mesh->createanchors = NULL; 2209 mesh->computeanchormatrix = NULL; 2210 mesh->parentSection = NULL; 2211 mesh->parents = NULL; 2212 mesh->childIDs = NULL; 2213 mesh->childSection = NULL; 2214 mesh->children = NULL; 2215 mesh->referenceTree = NULL; 2216 mesh->getchildsymmetry = NULL; 2217 for (d = 0; d < 8; ++d) mesh->hybridPointMax[d] = PETSC_DETERMINE; 2218 mesh->vtkCellHeight = 0; 2219 mesh->useCone = PETSC_FALSE; 2220 mesh->useClosure = PETSC_TRUE; 2221 mesh->useAnchors = PETSC_FALSE; 2222 2223 mesh->maxProjectionHeight = 0; 2224 2225 mesh->printSetValues = PETSC_FALSE; 2226 mesh->printFEM = 0; 2227 mesh->printTol = 1.0e-10; 2228 2229 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 2230 PetscFunctionReturn(0); 2231 } 2232 2233 /*@ 2234 DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram. 2235 2236 Collective on MPI_Comm 2237 2238 Input Parameter: 2239 . comm - The communicator for the DMPlex object 2240 2241 Output Parameter: 2242 . mesh - The DMPlex object 2243 2244 Level: beginner 2245 2246 .keywords: DMPlex, create 2247 @*/ 2248 PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh) 2249 { 2250 PetscErrorCode ierr; 2251 2252 PetscFunctionBegin; 2253 PetscValidPointer(mesh,2); 2254 ierr = DMCreate(comm, mesh);CHKERRQ(ierr); 2255 ierr = DMSetType(*mesh, DMPLEX);CHKERRQ(ierr); 2256 PetscFunctionReturn(0); 2257 } 2258 2259 /* 2260 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 2261 */ 2262 static PetscErrorCode DMPlexBuildFromCellList_Parallel_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[], PetscSF *sfVert) 2263 { 2264 PetscSF sfPoint; 2265 PetscLayout vLayout; 2266 PetscHashI vhash; 2267 PetscSFNode *remoteVerticesAdj, *vertexLocal, *vertexOwner, *remoteVertex; 2268 const PetscInt *vrange; 2269 PetscInt numVerticesAdj, off, *verticesAdj, numVerticesGhost = 0, *localVertex, *cone, c, p, v, g; 2270 PETSC_UNUSED PetscHashIIter ret, iter; 2271 PetscMPIInt rank, size; 2272 PetscErrorCode ierr; 2273 2274 PetscFunctionBegin; 2275 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr); 2276 ierr = MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);CHKERRQ(ierr); 2277 /* Partition vertices */ 2278 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) dm), &vLayout);CHKERRQ(ierr); 2279 ierr = PetscLayoutSetLocalSize(vLayout, numVertices);CHKERRQ(ierr); 2280 ierr = PetscLayoutSetBlockSize(vLayout, 1);CHKERRQ(ierr); 2281 ierr = PetscLayoutSetUp(vLayout);CHKERRQ(ierr); 2282 ierr = PetscLayoutGetRanges(vLayout, &vrange);CHKERRQ(ierr); 2283 /* Count vertices and map them to procs */ 2284 PetscHashICreate(vhash); 2285 for (c = 0; c < numCells; ++c) { 2286 for (p = 0; p < numCorners; ++p) { 2287 PetscHashIPut(vhash, cells[c*numCorners+p], ret, iter); 2288 } 2289 } 2290 PetscHashISize(vhash, numVerticesAdj); 2291 ierr = PetscMalloc1(numVerticesAdj, &verticesAdj);CHKERRQ(ierr); 2292 off = 0; ierr = PetscHashIGetKeys(vhash, &off, verticesAdj);CHKERRQ(ierr); 2293 if (off != numVerticesAdj) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid number of local vertices %D should be %D", off, numVerticesAdj); 2294 ierr = PetscSortInt(numVerticesAdj, verticesAdj);CHKERRQ(ierr); 2295 ierr = PetscMalloc1(numVerticesAdj, &remoteVerticesAdj);CHKERRQ(ierr); 2296 for (v = 0; v < numVerticesAdj; ++v) { 2297 const PetscInt gv = verticesAdj[v]; 2298 PetscInt vrank; 2299 2300 ierr = PetscFindInt(gv, size+1, vrange, &vrank);CHKERRQ(ierr); 2301 vrank = vrank < 0 ? -(vrank+2) : vrank; 2302 remoteVerticesAdj[v].index = gv - vrange[vrank]; 2303 remoteVerticesAdj[v].rank = vrank; 2304 } 2305 /* Create cones */ 2306 ierr = DMPlexSetChart(dm, 0, numCells+numVerticesAdj);CHKERRQ(ierr); 2307 for (c = 0; c < numCells; ++c) {ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr);} 2308 ierr = DMSetUp(dm);CHKERRQ(ierr); 2309 ierr = DMGetWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 2310 for (c = 0; c < numCells; ++c) { 2311 for (p = 0; p < numCorners; ++p) { 2312 const PetscInt gv = cells[c*numCorners+p]; 2313 PetscInt lv; 2314 2315 ierr = PetscFindInt(gv, numVerticesAdj, verticesAdj, &lv);CHKERRQ(ierr); 2316 if (lv < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not find global vertex %D in local connectivity", gv); 2317 cone[p] = lv+numCells; 2318 } 2319 ierr = DMPlexSetCone(dm, c, cone);CHKERRQ(ierr); 2320 } 2321 ierr = DMRestoreWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 2322 /* Create SF for vertices */ 2323 ierr = PetscSFCreate(PetscObjectComm((PetscObject)dm), sfVert);CHKERRQ(ierr); 2324 ierr = PetscObjectSetName((PetscObject) *sfVert, "Vertex Ownership SF");CHKERRQ(ierr); 2325 ierr = PetscSFSetFromOptions(*sfVert);CHKERRQ(ierr); 2326 ierr = PetscSFSetGraph(*sfVert, numVertices, numVerticesAdj, NULL, PETSC_OWN_POINTER, remoteVerticesAdj, PETSC_OWN_POINTER);CHKERRQ(ierr); 2327 ierr = PetscFree(verticesAdj);CHKERRQ(ierr); 2328 /* Build pointSF */ 2329 ierr = PetscMalloc2(numVerticesAdj, &vertexLocal, numVertices, &vertexOwner);CHKERRQ(ierr); 2330 for (v = 0; v < numVerticesAdj; ++v) {vertexLocal[v].index = v+numCells; vertexLocal[v].rank = rank;} 2331 for (v = 0; v < numVertices; ++v) {vertexOwner[v].index = -1; vertexOwner[v].rank = -1;} 2332 ierr = PetscSFReduceBegin(*sfVert, MPIU_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);CHKERRQ(ierr); 2333 ierr = PetscSFReduceEnd(*sfVert, MPIU_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);CHKERRQ(ierr); 2334 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); 2335 ierr = PetscSFBcastBegin(*sfVert, MPIU_2INT, vertexOwner, vertexLocal);CHKERRQ(ierr); 2336 ierr = PetscSFBcastEnd(*sfVert, MPIU_2INT, vertexOwner, vertexLocal);CHKERRQ(ierr); 2337 for (v = 0; v < numVerticesAdj; ++v) if (vertexLocal[v].rank != rank) ++numVerticesGhost; 2338 ierr = PetscMalloc1(numVerticesGhost, &localVertex);CHKERRQ(ierr); 2339 ierr = PetscMalloc1(numVerticesGhost, &remoteVertex);CHKERRQ(ierr); 2340 for (v = 0, g = 0; v < numVerticesAdj; ++v) { 2341 if (vertexLocal[v].rank != rank) { 2342 localVertex[g] = v+numCells; 2343 remoteVertex[g].index = vertexLocal[v].index; 2344 remoteVertex[g].rank = vertexLocal[v].rank; 2345 ++g; 2346 } 2347 } 2348 ierr = PetscFree2(vertexLocal, vertexOwner);CHKERRQ(ierr); 2349 if (g != numVerticesGhost) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of vertex ghosts %D should be %D", g, numVerticesGhost); 2350 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 2351 ierr = PetscObjectSetName((PetscObject) sfPoint, "point SF");CHKERRQ(ierr); 2352 ierr = PetscSFSetGraph(sfPoint, numCells+numVerticesAdj, numVerticesGhost, localVertex, PETSC_OWN_POINTER, remoteVertex, PETSC_OWN_POINTER);CHKERRQ(ierr); 2353 ierr = PetscLayoutDestroy(&vLayout);CHKERRQ(ierr); 2354 PetscHashIDestroy(vhash); 2355 /* Fill in the rest of the topology structure */ 2356 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 2357 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 2358 PetscFunctionReturn(0); 2359 } 2360 2361 /* 2362 This takes as input the coordinates for each owned vertex 2363 */ 2364 static PetscErrorCode DMPlexBuildCoordinates_Parallel_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numV, PetscSF sfVert, const PetscReal vertexCoords[]) 2365 { 2366 PetscSection coordSection; 2367 Vec coordinates; 2368 PetscScalar *coords; 2369 PetscInt numVertices, numVerticesAdj, coordSize, v; 2370 PetscErrorCode ierr; 2371 2372 PetscFunctionBegin; 2373 ierr = PetscSFGetGraph(sfVert, &numVertices, &numVerticesAdj, NULL, NULL);CHKERRQ(ierr); 2374 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2375 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 2376 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 2377 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVerticesAdj);CHKERRQ(ierr); 2378 for (v = numCells; v < numCells+numVerticesAdj; ++v) { 2379 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 2380 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 2381 } 2382 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 2383 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 2384 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 2385 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 2386 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 2387 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 2388 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 2389 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2390 { 2391 MPI_Datatype coordtype; 2392 2393 /* Need a temp buffer for coords if we have complex/single */ 2394 ierr = MPI_Type_contiguous(spaceDim, MPIU_SCALAR, &coordtype);CHKERRQ(ierr); 2395 ierr = MPI_Type_commit(&coordtype);CHKERRQ(ierr); 2396 #if defined(PETSC_USE_COMPLEX) 2397 { 2398 PetscScalar *svertexCoords; 2399 PetscInt i; 2400 ierr = PetscMalloc1(numV*spaceDim,&svertexCoords);CHKERRQ(ierr); 2401 for (i=0; i<numV*spaceDim; i++) svertexCoords[i] = vertexCoords[i]; 2402 ierr = PetscSFBcastBegin(sfVert, coordtype, svertexCoords, coords);CHKERRQ(ierr); 2403 ierr = PetscSFBcastEnd(sfVert, coordtype, svertexCoords, coords);CHKERRQ(ierr); 2404 ierr = PetscFree(svertexCoords);CHKERRQ(ierr); 2405 } 2406 #else 2407 ierr = PetscSFBcastBegin(sfVert, coordtype, vertexCoords, coords);CHKERRQ(ierr); 2408 ierr = PetscSFBcastEnd(sfVert, coordtype, vertexCoords, coords);CHKERRQ(ierr); 2409 #endif 2410 ierr = MPI_Type_free(&coordtype);CHKERRQ(ierr); 2411 } 2412 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2413 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 2414 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 2415 PetscFunctionReturn(0); 2416 } 2417 2418 /*@C 2419 DMPlexCreateFromCellListParallel - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM 2420 2421 Input Parameters: 2422 + comm - The communicator 2423 . dim - The topological dimension of the mesh 2424 . numCells - The number of cells owned by this process 2425 . numVertices - The number of vertices owned by this process 2426 . numCorners - The number of vertices for each cell 2427 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 2428 . cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 2429 . spaceDim - The spatial dimension used for coordinates 2430 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 2431 2432 Output Parameter: 2433 + dm - The DM 2434 - vertexSF - Optional, SF describing complete vertex ownership 2435 2436 Note: Two triangles sharing a face 2437 $ 2438 $ 2 2439 $ / | \ 2440 $ / | \ 2441 $ / | \ 2442 $ 0 0 | 1 3 2443 $ \ | / 2444 $ \ | / 2445 $ \ | / 2446 $ 1 2447 would have input 2448 $ numCells = 2, numVertices = 4 2449 $ cells = [0 1 2 1 3 2] 2450 $ 2451 which would result in the DMPlex 2452 $ 2453 $ 4 2454 $ / | \ 2455 $ / | \ 2456 $ / | \ 2457 $ 2 0 | 1 5 2458 $ \ | / 2459 $ \ | / 2460 $ \ | / 2461 $ 3 2462 2463 Level: beginner 2464 2465 .seealso: DMPlexCreateFromCellList(), DMPlexCreateFromDAG(), DMPlexCreate() 2466 @*/ 2467 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) 2468 { 2469 PetscSF sfVert; 2470 PetscErrorCode ierr; 2471 2472 PetscFunctionBegin; 2473 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 2474 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 2475 PetscValidLogicalCollectiveInt(*dm, dim, 2); 2476 PetscValidLogicalCollectiveInt(*dm, spaceDim, 8); 2477 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 2478 ierr = DMPlexBuildFromCellList_Parallel_Private(*dm, numCells, numVertices, numCorners, cells, &sfVert);CHKERRQ(ierr); 2479 if (interpolate) { 2480 DM idm = NULL; 2481 2482 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 2483 ierr = DMDestroy(dm);CHKERRQ(ierr); 2484 *dm = idm; 2485 } 2486 ierr = DMPlexBuildCoordinates_Parallel_Private(*dm, spaceDim, numCells, numVertices,sfVert, vertexCoords);CHKERRQ(ierr); 2487 if (vertexSF) *vertexSF = sfVert; 2488 else {ierr = PetscSFDestroy(&sfVert);CHKERRQ(ierr);} 2489 PetscFunctionReturn(0); 2490 } 2491 2492 /* 2493 This takes as input the common mesh generator output, a list of the vertices for each cell 2494 */ 2495 static PetscErrorCode DMPlexBuildFromCellList_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[]) 2496 { 2497 PetscInt *cone, c, p; 2498 PetscErrorCode ierr; 2499 2500 PetscFunctionBegin; 2501 ierr = DMPlexSetChart(dm, 0, numCells+numVertices);CHKERRQ(ierr); 2502 for (c = 0; c < numCells; ++c) { 2503 ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr); 2504 } 2505 ierr = DMSetUp(dm);CHKERRQ(ierr); 2506 ierr = DMGetWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 2507 for (c = 0; c < numCells; ++c) { 2508 for (p = 0; p < numCorners; ++p) { 2509 cone[p] = cells[c*numCorners+p]+numCells; 2510 } 2511 ierr = DMPlexSetCone(dm, c, cone);CHKERRQ(ierr); 2512 } 2513 ierr = DMRestoreWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 2514 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 2515 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 2516 PetscFunctionReturn(0); 2517 } 2518 2519 /* 2520 This takes as input the coordinates for each vertex 2521 */ 2522 static PetscErrorCode DMPlexBuildCoordinates_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, const double vertexCoords[]) 2523 { 2524 PetscSection coordSection; 2525 Vec coordinates; 2526 DM cdm; 2527 PetscScalar *coords; 2528 PetscInt v, d; 2529 PetscErrorCode ierr; 2530 2531 PetscFunctionBegin; 2532 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2533 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 2534 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 2535 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVertices);CHKERRQ(ierr); 2536 for (v = numCells; v < numCells+numVertices; ++v) { 2537 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 2538 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 2539 } 2540 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 2541 2542 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 2543 ierr = DMCreateLocalVector(cdm, &coordinates);CHKERRQ(ierr); 2544 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 2545 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 2546 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2547 for (v = 0; v < numVertices; ++v) { 2548 for (d = 0; d < spaceDim; ++d) { 2549 coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d]; 2550 } 2551 } 2552 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2553 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 2554 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 2555 PetscFunctionReturn(0); 2556 } 2557 2558 /*@C 2559 DMPlexCreateFromCellList - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM 2560 2561 Input Parameters: 2562 + comm - The communicator 2563 . dim - The topological dimension of the mesh 2564 . numCells - The number of cells 2565 . numVertices - The number of vertices 2566 . numCorners - The number of vertices for each cell 2567 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 2568 . cells - An array of numCells*numCorners numbers, the vertices for each cell 2569 . spaceDim - The spatial dimension used for coordinates 2570 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 2571 2572 Output Parameter: 2573 . dm - The DM 2574 2575 Note: Two triangles sharing a face 2576 $ 2577 $ 2 2578 $ / | \ 2579 $ / | \ 2580 $ / | \ 2581 $ 0 0 | 1 3 2582 $ \ | / 2583 $ \ | / 2584 $ \ | / 2585 $ 1 2586 would have input 2587 $ numCells = 2, numVertices = 4 2588 $ cells = [0 1 2 1 3 2] 2589 $ 2590 which would result in the DMPlex 2591 $ 2592 $ 4 2593 $ / | \ 2594 $ / | \ 2595 $ / | \ 2596 $ 2 0 | 1 5 2597 $ \ | / 2598 $ \ | / 2599 $ \ | / 2600 $ 3 2601 2602 Level: beginner 2603 2604 .seealso: DMPlexCreateFromDAG(), DMPlexCreate() 2605 @*/ 2606 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) 2607 { 2608 PetscErrorCode ierr; 2609 2610 PetscFunctionBegin; 2611 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 2612 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 2613 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 2614 ierr = DMPlexBuildFromCellList_Private(*dm, numCells, numVertices, numCorners, cells);CHKERRQ(ierr); 2615 if (interpolate) { 2616 DM idm = NULL; 2617 2618 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 2619 ierr = DMDestroy(dm);CHKERRQ(ierr); 2620 *dm = idm; 2621 } 2622 ierr = DMPlexBuildCoordinates_Private(*dm, spaceDim, numCells, numVertices, vertexCoords);CHKERRQ(ierr); 2623 PetscFunctionReturn(0); 2624 } 2625 2626 /*@ 2627 DMPlexCreateFromDAG - This takes as input the adjacency-list representation of the Directed Acyclic Graph (Hasse Diagram) encoding a mesh, and produces a DM 2628 2629 Input Parameters: 2630 + dm - The empty DM object, usually from DMCreate() and DMSetDimension() 2631 . depth - The depth of the DAG 2632 . numPoints - The number of points at each depth 2633 . coneSize - The cone size of each point 2634 . cones - The concatenation of the cone points for each point, the cone list must be oriented correctly for each point 2635 . coneOrientations - The orientation of each cone point 2636 - vertexCoords - An array of numVertices*dim numbers, the coordinates of each vertex 2637 2638 Output Parameter: 2639 . dm - The DM 2640 2641 Note: Two triangles sharing a face would have input 2642 $ depth = 1, numPoints = [4 2], coneSize = [3 3 0 0 0 0] 2643 $ cones = [2 3 4 3 5 4], coneOrientations = [0 0 0 0 0 0] 2644 $ vertexCoords = [-1.0 0.0 0.0 -1.0 0.0 1.0 1.0 0.0] 2645 $ 2646 which would result in the DMPlex 2647 $ 2648 $ 4 2649 $ / | \ 2650 $ / | \ 2651 $ / | \ 2652 $ 2 0 | 1 5 2653 $ \ | / 2654 $ \ | / 2655 $ \ | / 2656 $ 3 2657 $ 2658 $ Notice that all points are numbered consecutively, unlikely DMPlexCreateFromCellList() 2659 2660 Level: advanced 2661 2662 .seealso: DMPlexCreateFromCellList(), DMPlexCreate() 2663 @*/ 2664 PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[]) 2665 { 2666 Vec coordinates; 2667 PetscSection coordSection; 2668 PetscScalar *coords; 2669 PetscInt coordSize, firstVertex = -1, pStart = 0, pEnd = 0, p, v, dim, dimEmbed, d, off; 2670 PetscErrorCode ierr; 2671 2672 PetscFunctionBegin; 2673 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2674 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 2675 if (dimEmbed < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Embedding dimension %d cannot be less than intrinsic dimension %d",dimEmbed,dim); 2676 for (d = 0; d <= depth; ++d) pEnd += numPoints[d]; 2677 ierr = DMPlexSetChart(dm, pStart, pEnd);CHKERRQ(ierr); 2678 for (p = pStart; p < pEnd; ++p) { 2679 ierr = DMPlexSetConeSize(dm, p, coneSize[p-pStart]);CHKERRQ(ierr); 2680 if (firstVertex < 0 && !coneSize[p - pStart]) { 2681 firstVertex = p - pStart; 2682 } 2683 } 2684 if (firstVertex < 0 && numPoints[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected %d vertices but could not find any", numPoints[0]); 2685 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 2686 for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) { 2687 ierr = DMPlexSetCone(dm, p, &cones[off]);CHKERRQ(ierr); 2688 ierr = DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);CHKERRQ(ierr); 2689 } 2690 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 2691 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 2692 /* Build coordinates */ 2693 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2694 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 2695 ierr = PetscSectionSetFieldComponents(coordSection, 0, dimEmbed);CHKERRQ(ierr); 2696 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0]);CHKERRQ(ierr); 2697 for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) { 2698 ierr = PetscSectionSetDof(coordSection, v, dimEmbed);CHKERRQ(ierr); 2699 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed);CHKERRQ(ierr); 2700 } 2701 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 2702 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 2703 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 2704 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 2705 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 2706 ierr = VecSetBlockSize(coordinates, dimEmbed);CHKERRQ(ierr); 2707 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 2708 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2709 for (v = 0; v < numPoints[0]; ++v) { 2710 PetscInt off; 2711 2712 ierr = PetscSectionGetOffset(coordSection, v+firstVertex, &off);CHKERRQ(ierr); 2713 for (d = 0; d < dimEmbed; ++d) { 2714 coords[off+d] = vertexCoords[v*dimEmbed+d]; 2715 } 2716 } 2717 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2718 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 2719 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 2720 PetscFunctionReturn(0); 2721 } 2722 2723 /*@C 2724 DMPlexCreateFromFile - This takes a filename and produces a DM 2725 2726 Input Parameters: 2727 + comm - The communicator 2728 . filename - A file name 2729 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 2730 2731 Output Parameter: 2732 . dm - The DM 2733 2734 Level: beginner 2735 2736 .seealso: DMPlexCreateFromDAG(), DMPlexCreateFromCellList(), DMPlexCreate() 2737 @*/ 2738 PetscErrorCode DMPlexCreateFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 2739 { 2740 const char *extGmsh = ".msh"; 2741 const char *extCGNS = ".cgns"; 2742 const char *extExodus = ".exo"; 2743 const char *extFluent = ".cas"; 2744 const char *extHDF5 = ".h5"; 2745 const char *extMed = ".med"; 2746 const char *extPLY = ".ply"; 2747 size_t len; 2748 PetscBool isGmsh, isCGNS, isExodus, isFluent, isHDF5, isMed, isPLY; 2749 PetscMPIInt rank; 2750 PetscErrorCode ierr; 2751 2752 PetscFunctionBegin; 2753 PetscValidPointer(filename, 2); 2754 PetscValidPointer(dm, 4); 2755 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 2756 ierr = PetscStrlen(filename, &len);CHKERRQ(ierr); 2757 if (!len) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Filename must be a valid path"); 2758 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGmsh, 4, &isGmsh);CHKERRQ(ierr); 2759 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extCGNS, 5, &isCGNS);CHKERRQ(ierr); 2760 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extExodus, 4, &isExodus);CHKERRQ(ierr); 2761 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extFluent, 4, &isFluent);CHKERRQ(ierr); 2762 ierr = PetscStrncmp(&filename[PetscMax(0,len-3)], extHDF5, 3, &isHDF5);CHKERRQ(ierr); 2763 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extMed, 4, &isMed);CHKERRQ(ierr); 2764 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extPLY, 4, &isPLY);CHKERRQ(ierr); 2765 if (isGmsh) { 2766 ierr = DMPlexCreateGmshFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2767 } else if (isCGNS) { 2768 ierr = DMPlexCreateCGNSFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2769 } else if (isExodus) { 2770 ierr = DMPlexCreateExodusFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2771 } else if (isFluent) { 2772 ierr = DMPlexCreateFluentFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2773 } else if (isHDF5) { 2774 PetscViewer viewer; 2775 2776 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 2777 ierr = PetscViewerSetType(viewer, PETSCVIEWERHDF5);CHKERRQ(ierr); 2778 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 2779 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 2780 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 2781 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 2782 ierr = DMLoad(*dm, viewer);CHKERRQ(ierr); 2783 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 2784 } else if (isMed) { 2785 ierr = DMPlexCreateMedFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2786 } else if (isPLY) { 2787 ierr = DMPlexCreatePLYFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2788 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot load file %s: unrecognized extension", filename); 2789 PetscFunctionReturn(0); 2790 } 2791 2792 /*@ 2793 DMPlexCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell 2794 2795 Collective on comm 2796 2797 Input Parameters: 2798 + comm - The communicator 2799 . dim - The spatial dimension 2800 - simplex - Flag for simplex, otherwise use a tensor-product cell 2801 2802 Output Parameter: 2803 . refdm - The reference cell 2804 2805 Level: intermediate 2806 2807 .keywords: reference cell 2808 .seealso: 2809 @*/ 2810 PetscErrorCode DMPlexCreateReferenceCell(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *refdm) 2811 { 2812 DM rdm; 2813 Vec coords; 2814 PetscErrorCode ierr; 2815 2816 PetscFunctionBeginUser; 2817 ierr = DMCreate(comm, &rdm);CHKERRQ(ierr); 2818 ierr = DMSetType(rdm, DMPLEX);CHKERRQ(ierr); 2819 ierr = DMSetDimension(rdm, dim);CHKERRQ(ierr); 2820 switch (dim) { 2821 case 0: 2822 { 2823 PetscInt numPoints[1] = {1}; 2824 PetscInt coneSize[1] = {0}; 2825 PetscInt cones[1] = {0}; 2826 PetscInt coneOrientations[1] = {0}; 2827 PetscScalar vertexCoords[1] = {0.0}; 2828 2829 ierr = DMPlexCreateFromDAG(rdm, 0, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2830 } 2831 break; 2832 case 1: 2833 { 2834 PetscInt numPoints[2] = {2, 1}; 2835 PetscInt coneSize[3] = {2, 0, 0}; 2836 PetscInt cones[2] = {1, 2}; 2837 PetscInt coneOrientations[2] = {0, 0}; 2838 PetscScalar vertexCoords[2] = {-1.0, 1.0}; 2839 2840 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2841 } 2842 break; 2843 case 2: 2844 if (simplex) { 2845 PetscInt numPoints[2] = {3, 1}; 2846 PetscInt coneSize[4] = {3, 0, 0, 0}; 2847 PetscInt cones[3] = {1, 2, 3}; 2848 PetscInt coneOrientations[3] = {0, 0, 0}; 2849 PetscScalar vertexCoords[6] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0}; 2850 2851 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2852 } else { 2853 PetscInt numPoints[2] = {4, 1}; 2854 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 2855 PetscInt cones[4] = {1, 2, 3, 4}; 2856 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 2857 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0}; 2858 2859 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2860 } 2861 break; 2862 case 3: 2863 if (simplex) { 2864 PetscInt numPoints[2] = {4, 1}; 2865 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 2866 PetscInt cones[4] = {1, 3, 2, 4}; 2867 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 2868 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}; 2869 2870 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2871 } else { 2872 PetscInt numPoints[2] = {8, 1}; 2873 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0}; 2874 PetscInt cones[8] = {1, 4, 3, 2, 5, 6, 7, 8}; 2875 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 2876 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, 2877 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0}; 2878 2879 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2880 } 2881 break; 2882 default: 2883 SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Cannot create reference cell for dimension %d", dim); 2884 } 2885 *refdm = NULL; 2886 ierr = DMPlexInterpolate(rdm, refdm);CHKERRQ(ierr); 2887 if (rdm->coordinateDM) { 2888 DM ncdm; 2889 PetscSection cs; 2890 PetscInt pEnd = -1; 2891 2892 ierr = DMGetDefaultSection(rdm->coordinateDM, &cs);CHKERRQ(ierr); 2893 if (cs) {ierr = PetscSectionGetChart(cs, NULL, &pEnd);CHKERRQ(ierr);} 2894 if (pEnd >= 0) { 2895 ierr = DMClone(rdm->coordinateDM, &ncdm);CHKERRQ(ierr); 2896 ierr = DMSetDefaultSection(ncdm, cs);CHKERRQ(ierr); 2897 ierr = DMSetCoordinateDM(*refdm, ncdm);CHKERRQ(ierr); 2898 ierr = DMDestroy(&ncdm);CHKERRQ(ierr); 2899 } 2900 } 2901 ierr = DMGetCoordinatesLocal(rdm, &coords);CHKERRQ(ierr); 2902 if (coords) { 2903 ierr = DMSetCoordinatesLocal(*refdm, coords);CHKERRQ(ierr); 2904 } else { 2905 ierr = DMGetCoordinates(rdm, &coords);CHKERRQ(ierr); 2906 if (coords) {ierr = DMSetCoordinates(*refdm, coords);CHKERRQ(ierr);} 2907 } 2908 ierr = DMDestroy(&rdm);CHKERRQ(ierr); 2909 PetscFunctionReturn(0); 2910 } 2911