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