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