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