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 Options Database Keys: 979 + -dm_plex_box_lower <x,y,z> - Specify lower-left-bottom coordinates for the box 980 - -dm_plex_box_upper <x,y,z> - Specify upper-right-top coordinates for the box 981 982 Note: Here is the numbering returned for 2 faces in each direction for tensor cells: 983 $ 10---17---11---18----12 984 $ | | | 985 $ | | | 986 $ 20 2 22 3 24 987 $ | | | 988 $ | | | 989 $ 7---15----8---16----9 990 $ | | | 991 $ | | | 992 $ 19 0 21 1 23 993 $ | | | 994 $ | | | 995 $ 4---13----5---14----6 996 997 and for simplicial cells 998 999 $ 14----8---15----9----16 1000 $ |\ 5 |\ 7 | 1001 $ | \ | \ | 1002 $ 13 2 14 3 15 1003 $ | 4 \ | 6 \ | 1004 $ | \ | \ | 1005 $ 11----6---12----7----13 1006 $ |\ |\ | 1007 $ | \ 1 | \ 3 | 1008 $ 10 0 11 1 12 1009 $ | 0 \ | 2 \ | 1010 $ | \ | \ | 1011 $ 8----4----9----5----10 1012 1013 Level: beginner 1014 1015 .keywords: DM, create 1016 .seealso: DMPlexCreateFromFile(), DMPlexCreateHexCylinderMesh(), DMSetType(), DMCreate() 1017 @*/ 1018 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) 1019 { 1020 PetscInt fac[3] = {0, 0, 0}; 1021 PetscReal low[3] = {0, 0, 0}; 1022 PetscReal upp[3] = {1, 1, 1}; 1023 DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE}; 1024 PetscInt i, n; 1025 PetscBool flg; 1026 PetscErrorCode ierr; 1027 1028 PetscFunctionBegin; 1029 for (i = 0; i < dim; ++i) fac[i] = faces ? faces[i] : (dim == 1 ? 1 : 4-dim); 1030 if (lower) for (i = 0; i < dim; ++i) low[i] = lower[i]; 1031 if (upper) for (i = 0; i < dim; ++i) upp[i] = upper[i]; 1032 if (periodicity) for (i = 0; i < dim; ++i) bdt[i] = periodicity[i]; 1033 /* Allow bounds to be specified from the command line */ 1034 n = 3; 1035 ierr = PetscOptionsGetRealArray(NULL, NULL, "-dm_plex_box_lower", low, &n, &flg);CHKERRQ(ierr); 1036 if (flg && (n != dim)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Lower box point had %D values, should have been %D", n, dim); 1037 n = 3; 1038 ierr = PetscOptionsGetRealArray(NULL, NULL, "-dm_plex_box_upper", upp, &n, &flg);CHKERRQ(ierr); 1039 if (flg && (n != dim)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Upper box point had %D values, should have been %D", n, dim); 1040 1041 if (dim == 1) {ierr = DMPlexCreateLineMesh_Internal(comm, fac[0], low[0], upp[0], bdt[0], dm);CHKERRQ(ierr);} 1042 else if (simplex) {ierr = DMPlexCreateBoxMesh_Simplex_Internal(comm, dim, fac, low, upp, bdt, interpolate, dm);CHKERRQ(ierr);} 1043 else {ierr = DMPlexCreateBoxMesh_Tensor_Internal(comm, dim, fac, low, upp, bdt, interpolate, dm);CHKERRQ(ierr);} 1044 PetscFunctionReturn(0); 1045 } 1046 1047 /*@ 1048 DMPlexCreateWedgeBoxMesh - Creates a 3-D mesh tesselating the (x,y) plane and extruding in the third direction using wedge cells. 1049 1050 Collective on MPI_Comm 1051 1052 Input Parameters: 1053 + comm - The communicator for the DM object 1054 . faces - Number of faces per dimension, or NULL for (1, 1, 1) 1055 . lower - The lower left corner, or NULL for (0, 0, 0) 1056 . upper - The upper right corner, or NULL for (1, 1, 1) 1057 . periodicity - The boundary type for the X,Y,Z direction, or NULL for DM_BOUNDARY_NONE 1058 . ordExt - If PETSC_TRUE, orders the extruded cells in the height first. Otherwise, orders the cell on the layers first 1059 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 1060 1061 Output Parameter: 1062 . dm - The DM object 1063 1064 Level: beginner 1065 1066 .keywords: DM, create 1067 .seealso: DMPlexCreateHexCylinderMesh(), DMPlexCreateWedgeCylinderMesh(), DMPlexExtrude(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 1068 @*/ 1069 PetscErrorCode DMPlexCreateWedgeBoxMesh(MPI_Comm comm, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool ordExt, PetscBool interpolate, DM *dm) 1070 { 1071 DM bdm, botdm; 1072 PetscInt i; 1073 PetscInt fac[3] = {0, 0, 0}; 1074 PetscReal low[3] = {0, 0, 0}; 1075 PetscReal upp[3] = {1, 1, 1}; 1076 DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE}; 1077 PetscErrorCode ierr; 1078 1079 PetscFunctionBegin; 1080 for (i = 0; i < 3; ++i) fac[i] = faces ? (faces[i] > 0 ? faces[i] : 1) : 1; 1081 if (lower) for (i = 0; i < 3; ++i) low[i] = lower[i]; 1082 if (upper) for (i = 0; i < 3; ++i) upp[i] = upper[i]; 1083 if (periodicity) for (i = 0; i < 3; ++i) bdt[i] = periodicity[i]; 1084 for (i = 0; i < 3; ++i) if (bdt[i] != DM_BOUNDARY_NONE) SETERRQ(comm, PETSC_ERR_SUP, "Periodicity not yet supported"); 1085 1086 ierr = DMCreate(comm, &bdm);CHKERRQ(ierr); 1087 ierr = DMSetType(bdm, DMPLEX);CHKERRQ(ierr); 1088 ierr = DMSetDimension(bdm, 1);CHKERRQ(ierr); 1089 ierr = DMSetCoordinateDim(bdm, 2);CHKERRQ(ierr); 1090 ierr = DMPlexCreateSquareBoundary(bdm, low, upp, fac);CHKERRQ(ierr); 1091 ierr = DMPlexGenerate(bdm, NULL, PETSC_FALSE, &botdm);CHKERRQ(ierr); 1092 ierr = DMDestroy(&bdm);CHKERRQ(ierr); 1093 ierr = DMPlexExtrude(botdm, fac[2], upp[2] - low[2], ordExt, interpolate, dm);CHKERRQ(ierr); 1094 if (low[2] != 0.0) { 1095 Vec v; 1096 PetscScalar *x; 1097 PetscInt cDim, n; 1098 1099 ierr = DMGetCoordinatesLocal(*dm, &v);CHKERRQ(ierr); 1100 ierr = VecGetBlockSize(v, &cDim);CHKERRQ(ierr); 1101 ierr = VecGetLocalSize(v, &n);CHKERRQ(ierr); 1102 ierr = VecGetArray(v, &x);CHKERRQ(ierr); 1103 x += cDim; 1104 for (i=0; i<n; i+=cDim) x[i] += low[2]; 1105 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1106 ierr = DMSetCoordinatesLocal(*dm, v);CHKERRQ(ierr); 1107 } 1108 ierr = DMDestroy(&botdm);CHKERRQ(ierr); 1109 PetscFunctionReturn(0); 1110 } 1111 1112 /*@ 1113 DMPlexExtrude - Creates a (d+1)-D mesh by extruding a d-D mesh in the normal direction using prismatic cells. 1114 1115 Collective on idm 1116 1117 Input Parameters: 1118 + idm - The mesh to be extruted 1119 . layers - The number of layers 1120 . height - The height of the extruded layer 1121 . ordExt - If PETSC_TRUE, orders the extruded cells in the height first. Otherwise, orders the cell on the layers first 1122 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 1123 1124 Output Parameter: 1125 . dm - The DM object 1126 1127 Notes: The object created is an hybrid mesh, the vertex ordering in the cone of the cell is that of the prismatic cells 1128 1129 Level: advanced 1130 1131 .keywords: DM, create 1132 .seealso: DMPlexCreateWedgeCylinderMesh(), DMPlexCreateWedgeBoxMesh(), DMPlexSetHybridBounds(), DMSetType(), DMCreate() 1133 @*/ 1134 PetscErrorCode DMPlexExtrude(DM idm, PetscInt layers, PetscReal height, PetscBool ordExt, PetscBool interpolate, DM* dm) 1135 { 1136 PetscScalar *coordsB; 1137 const PetscScalar *coordsA; 1138 PetscReal *normals = NULL; 1139 Vec coordinatesA, coordinatesB; 1140 PetscSection coordSectionA, coordSectionB; 1141 PetscInt dim, cDim, cDimB, c, l, v, coordSize, *newCone; 1142 PetscInt cStart, cEnd, vStart, vEnd, cellV, numCells, numVertices; 1143 PetscErrorCode ierr; 1144 1145 PetscFunctionBegin; 1146 PetscValidHeaderSpecific(idm, DM_CLASSID, 1); 1147 PetscValidLogicalCollectiveInt(idm, layers, 2); 1148 PetscValidLogicalCollectiveReal(idm, height, 3); 1149 PetscValidLogicalCollectiveBool(idm, interpolate, 4); 1150 ierr = DMGetDimension(idm, &dim);CHKERRQ(ierr); 1151 if (dim < 1 || dim > 3) SETERRQ1(PetscObjectComm((PetscObject)idm), PETSC_ERR_SUP, "Support for dimension %D not coded", dim); 1152 1153 ierr = DMPlexGetHeightStratum(idm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1154 ierr = DMPlexGetDepthStratum(idm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1155 numCells = (cEnd - cStart)*layers; 1156 numVertices = (vEnd - vStart)*(layers+1); 1157 ierr = DMCreate(PetscObjectComm((PetscObject)idm), dm);CHKERRQ(ierr); 1158 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1159 ierr = DMSetDimension(*dm, dim+1);CHKERRQ(ierr); 1160 ierr = DMPlexSetChart(*dm, 0, numCells+numVertices);CHKERRQ(ierr); 1161 for (c = cStart, cellV = 0; c < cEnd; ++c) { 1162 PetscInt *closure = NULL; 1163 PetscInt closureSize, numCorners = 0; 1164 1165 ierr = DMPlexGetTransitiveClosure(idm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1166 for (v = 0; v < closureSize*2; v += 2) if ((closure[v] >= vStart) && (closure[v] < vEnd)) numCorners++; 1167 ierr = DMPlexRestoreTransitiveClosure(idm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1168 for (l = 0; l < layers; ++l) { 1169 ierr = DMPlexSetConeSize(*dm, ordExt ? layers*(c - cStart) + l : l*(cEnd - cStart) + c - cStart, 2*numCorners);CHKERRQ(ierr); 1170 } 1171 cellV = PetscMax(numCorners,cellV); 1172 } 1173 ierr = DMPlexSetHybridBounds(*dm, 0, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 1174 ierr = DMSetUp(*dm);CHKERRQ(ierr); 1175 1176 ierr = DMGetCoordinateDim(idm, &cDim);CHKERRQ(ierr); 1177 if (dim != cDim) { 1178 ierr = PetscCalloc1(cDim*(vEnd - vStart), &normals);CHKERRQ(ierr); 1179 } 1180 ierr = PetscMalloc1(3*cellV,&newCone);CHKERRQ(ierr); 1181 for (c = cStart; c < cEnd; ++c) { 1182 PetscInt *closure = NULL; 1183 PetscInt closureSize, numCorners = 0, l; 1184 PetscReal normal[3] = {0, 0, 0}; 1185 1186 if (normals) { 1187 ierr = DMPlexComputeCellGeometryFVM(idm, c, NULL, NULL, normal);CHKERRQ(ierr); 1188 } 1189 ierr = DMPlexGetTransitiveClosure(idm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1190 for (v = 0; v < closureSize*2; v += 2) { 1191 if ((closure[v] >= vStart) && (closure[v] < vEnd)) { 1192 PetscInt d; 1193 1194 newCone[numCorners++] = closure[v] - vStart; 1195 if (normals) { for (d = 0; d < cDim; ++d) normals[cDim*(closure[v]-vStart)+d] += normal[d]; } 1196 } 1197 } 1198 switch (numCorners) { 1199 case 4: /* do nothing */ 1200 case 2: /* do nothing */ 1201 break; 1202 case 3: /* from counter-clockwise to wedge ordering */ 1203 l = newCone[1]; 1204 newCone[1] = newCone[2]; 1205 newCone[2] = l; 1206 break; 1207 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported number of corners: %D", numCorners); 1208 } 1209 ierr = DMPlexRestoreTransitiveClosure(idm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1210 for (l = 0; l < layers; ++l) { 1211 PetscInt i; 1212 1213 for (i = 0; i < numCorners; ++i) { 1214 newCone[ numCorners + i] = ordExt ? (layers+1)*newCone[i] + l + numCells : l*(vEnd - vStart) + newCone[i] + numCells; 1215 newCone[2*numCorners + i] = ordExt ? (layers+1)*newCone[i] + l + 1 + numCells : (l+1)*(vEnd - vStart) + newCone[i] + numCells; 1216 } 1217 ierr = DMPlexSetCone(*dm, ordExt ? layers*(c - cStart) + l : l*(cEnd - cStart) + c - cStart, newCone + numCorners);CHKERRQ(ierr); 1218 } 1219 } 1220 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1221 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1222 ierr = PetscFree(newCone);CHKERRQ(ierr); 1223 1224 cDimB = cDim == dim ? cDim+1 : cDim; 1225 ierr = DMGetCoordinateSection(*dm, &coordSectionB);CHKERRQ(ierr); 1226 ierr = PetscSectionSetNumFields(coordSectionB, 1);CHKERRQ(ierr); 1227 ierr = PetscSectionSetFieldComponents(coordSectionB, 0, cDimB);CHKERRQ(ierr); 1228 ierr = PetscSectionSetChart(coordSectionB, numCells, numCells+numVertices);CHKERRQ(ierr); 1229 for (v = numCells; v < numCells+numVertices; ++v) { 1230 ierr = PetscSectionSetDof(coordSectionB, v, cDimB);CHKERRQ(ierr); 1231 ierr = PetscSectionSetFieldDof(coordSectionB, v, 0, cDimB);CHKERRQ(ierr); 1232 } 1233 ierr = PetscSectionSetUp(coordSectionB);CHKERRQ(ierr); 1234 ierr = PetscSectionGetStorageSize(coordSectionB, &coordSize);CHKERRQ(ierr); 1235 ierr = VecCreate(PETSC_COMM_SELF, &coordinatesB);CHKERRQ(ierr); 1236 ierr = PetscObjectSetName((PetscObject) coordinatesB, "coordinates");CHKERRQ(ierr); 1237 ierr = VecSetSizes(coordinatesB, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1238 ierr = VecSetBlockSize(coordinatesB, cDimB);CHKERRQ(ierr); 1239 ierr = VecSetType(coordinatesB,VECSTANDARD);CHKERRQ(ierr); 1240 1241 ierr = DMGetCoordinateSection(idm, &coordSectionA);CHKERRQ(ierr); 1242 ierr = DMGetCoordinatesLocal(idm, &coordinatesA);CHKERRQ(ierr); 1243 ierr = VecGetArray(coordinatesB, &coordsB);CHKERRQ(ierr); 1244 ierr = VecGetArrayRead(coordinatesA, &coordsA);CHKERRQ(ierr); 1245 for (v = vStart; v < vEnd; ++v) { 1246 const PetscScalar *cptr; 1247 PetscReal ones2[2] = { 0., 1.}, ones3[3] = { 0., 0., 1.}; 1248 PetscReal *normal, norm, h = height/layers; 1249 PetscInt offA, d, cDimA = cDim; 1250 1251 normal = normals ? normals + cDimB*(v - vStart) : (cDim > 1 ? ones3 : ones2); 1252 if (normals) { 1253 for (d = 0, norm = 0.0; d < cDimB; ++d) norm += normal[d]*normal[d]; 1254 for (d = 0; d < cDimB; ++d) normal[d] *= 1./PetscSqrtReal(norm); 1255 } 1256 1257 ierr = PetscSectionGetOffset(coordSectionA, v, &offA);CHKERRQ(ierr); 1258 cptr = coordsA + offA; 1259 for (l = 0; l < layers+1; ++l) { 1260 PetscInt offB, d, newV; 1261 1262 newV = ordExt ? (layers+1)*(v -vStart) + l + numCells : (vEnd -vStart)*l + (v -vStart) + numCells; 1263 ierr = PetscSectionGetOffset(coordSectionB, newV, &offB);CHKERRQ(ierr); 1264 for (d = 0; d < cDimA; ++d) { coordsB[offB+d] = cptr[d]; } 1265 for (d = 0; d < cDimB; ++d) { coordsB[offB+d] += l ? normal[d]*h : 0.0; } 1266 cptr = coordsB + offB; 1267 cDimA = cDimB; 1268 } 1269 } 1270 ierr = VecRestoreArrayRead(coordinatesA, &coordsA);CHKERRQ(ierr); 1271 ierr = VecRestoreArray(coordinatesB, &coordsB);CHKERRQ(ierr); 1272 ierr = DMSetCoordinatesLocal(*dm, coordinatesB);CHKERRQ(ierr); 1273 ierr = VecDestroy(&coordinatesB);CHKERRQ(ierr); 1274 ierr = PetscFree(normals);CHKERRQ(ierr); 1275 if (interpolate) { 1276 DM idm; 1277 1278 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1279 ierr = DMPlexCopyCoordinates(*dm, idm);CHKERRQ(ierr); 1280 ierr = DMDestroy(dm);CHKERRQ(ierr); 1281 *dm = idm; 1282 } 1283 PetscFunctionReturn(0); 1284 } 1285 1286 /*@C 1287 DMPlexSetOptionsPrefix - Sets the prefix used for searching for all DM options in the database. 1288 1289 Logically Collective on DM 1290 1291 Input Parameters: 1292 + dm - the DM context 1293 - prefix - the prefix to prepend to all option names 1294 1295 Notes: 1296 A hyphen (-) must NOT be given at the beginning of the prefix name. 1297 The first character of all runtime options is AUTOMATICALLY the hyphen. 1298 1299 Level: advanced 1300 1301 .seealso: SNESSetFromOptions() 1302 @*/ 1303 PetscErrorCode DMPlexSetOptionsPrefix(DM dm, const char prefix[]) 1304 { 1305 DM_Plex *mesh = (DM_Plex *) dm->data; 1306 PetscErrorCode ierr; 1307 1308 PetscFunctionBegin; 1309 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1310 ierr = PetscObjectSetOptionsPrefix((PetscObject) dm, prefix);CHKERRQ(ierr); 1311 ierr = PetscObjectSetOptionsPrefix((PetscObject) mesh->partitioner, prefix);CHKERRQ(ierr); 1312 PetscFunctionReturn(0); 1313 } 1314 1315 /*@ 1316 DMPlexCreateHexCylinderMesh - Creates a mesh on the tensor product of the unit interval with the circle (cylinder) using hexahedra. 1317 1318 Collective on MPI_Comm 1319 1320 Input Parameters: 1321 + comm - The communicator for the DM object 1322 . numRefine - The number of regular refinements to the basic 5 cell structure 1323 - periodicZ - The boundary type for the Z direction 1324 1325 Output Parameter: 1326 . dm - The DM object 1327 1328 Note: Here is the output numbering looking from the bottom of the cylinder: 1329 $ 17-----14 1330 $ | | 1331 $ | 2 | 1332 $ | | 1333 $ 17-----8-----7-----14 1334 $ | | | | 1335 $ | 3 | 0 | 1 | 1336 $ | | | | 1337 $ 19-----5-----6-----13 1338 $ | | 1339 $ | 4 | 1340 $ | | 1341 $ 19-----13 1342 $ 1343 $ and up through the top 1344 $ 1345 $ 18-----16 1346 $ | | 1347 $ | 2 | 1348 $ | | 1349 $ 18----10----11-----16 1350 $ | | | | 1351 $ | 3 | 0 | 1 | 1352 $ | | | | 1353 $ 20-----9----12-----15 1354 $ | | 1355 $ | 4 | 1356 $ | | 1357 $ 20-----15 1358 1359 Level: beginner 1360 1361 .keywords: DM, create 1362 .seealso: DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 1363 @*/ 1364 PetscErrorCode DMPlexCreateHexCylinderMesh(MPI_Comm comm, PetscInt numRefine, DMBoundaryType periodicZ, DM *dm) 1365 { 1366 const PetscInt dim = 3; 1367 PetscInt numCells, numVertices, r; 1368 PetscMPIInt rank; 1369 PetscErrorCode ierr; 1370 1371 PetscFunctionBegin; 1372 PetscValidPointer(dm, 4); 1373 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 1374 if (numRefine < 0) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Number of refinements %D cannot be negative", numRefine); 1375 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1376 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1377 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1378 /* Create topology */ 1379 { 1380 PetscInt cone[8], c; 1381 1382 numCells = !rank ? 5 : 0; 1383 numVertices = !rank ? 16 : 0; 1384 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1385 numCells *= 3; 1386 numVertices = !rank ? 24 : 0; 1387 } 1388 ierr = DMPlexSetChart(*dm, 0, numCells+numVertices);CHKERRQ(ierr); 1389 for (c = 0; c < numCells; c++) {ierr = DMPlexSetConeSize(*dm, c, 8);CHKERRQ(ierr);} 1390 ierr = DMSetUp(*dm);CHKERRQ(ierr); 1391 if (!rank) { 1392 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1393 cone[0] = 15; cone[1] = 18; cone[2] = 17; cone[3] = 16; 1394 cone[4] = 31; cone[5] = 32; cone[6] = 33; cone[7] = 34; 1395 ierr = DMPlexSetCone(*dm, 0, cone);CHKERRQ(ierr); 1396 cone[0] = 16; cone[1] = 17; cone[2] = 24; cone[3] = 23; 1397 cone[4] = 32; cone[5] = 36; cone[6] = 37; cone[7] = 33; /* 22 25 26 21 */ 1398 ierr = DMPlexSetCone(*dm, 1, cone);CHKERRQ(ierr); 1399 cone[0] = 18; cone[1] = 27; cone[2] = 24; cone[3] = 17; 1400 cone[4] = 34; cone[5] = 33; cone[6] = 37; cone[7] = 38; 1401 ierr = DMPlexSetCone(*dm, 2, cone);CHKERRQ(ierr); 1402 cone[0] = 29; cone[1] = 27; cone[2] = 18; cone[3] = 15; 1403 cone[4] = 35; cone[5] = 31; cone[6] = 34; cone[7] = 38; 1404 ierr = DMPlexSetCone(*dm, 3, cone);CHKERRQ(ierr); 1405 cone[0] = 29; cone[1] = 15; cone[2] = 16; cone[3] = 23; 1406 cone[4] = 35; cone[5] = 36; cone[6] = 32; cone[7] = 31; 1407 ierr = DMPlexSetCone(*dm, 4, cone);CHKERRQ(ierr); 1408 1409 cone[0] = 31; cone[1] = 34; cone[2] = 33; cone[3] = 32; 1410 cone[4] = 19; cone[5] = 22; cone[6] = 21; cone[7] = 20; 1411 ierr = DMPlexSetCone(*dm, 5, cone);CHKERRQ(ierr); 1412 cone[0] = 32; cone[1] = 33; cone[2] = 37; cone[3] = 36; 1413 cone[4] = 22; cone[5] = 25; cone[6] = 26; cone[7] = 21; 1414 ierr = DMPlexSetCone(*dm, 6, cone);CHKERRQ(ierr); 1415 cone[0] = 34; cone[1] = 38; cone[2] = 37; cone[3] = 33; 1416 cone[4] = 20; cone[5] = 21; cone[6] = 26; cone[7] = 28; 1417 ierr = DMPlexSetCone(*dm, 7, cone);CHKERRQ(ierr); 1418 cone[0] = 35; cone[1] = 38; cone[2] = 34; cone[3] = 31; 1419 cone[4] = 30; cone[5] = 19; cone[6] = 20; cone[7] = 28; 1420 ierr = DMPlexSetCone(*dm, 8, cone);CHKERRQ(ierr); 1421 cone[0] = 35; cone[1] = 31; cone[2] = 32; cone[3] = 36; 1422 cone[4] = 30; cone[5] = 25; cone[6] = 22; cone[7] = 19; 1423 ierr = DMPlexSetCone(*dm, 9, cone);CHKERRQ(ierr); 1424 1425 cone[0] = 19; cone[1] = 20; cone[2] = 21; cone[3] = 22; 1426 cone[4] = 15; cone[5] = 16; cone[6] = 17; cone[7] = 18; 1427 ierr = DMPlexSetCone(*dm, 10, cone);CHKERRQ(ierr); 1428 cone[0] = 22; cone[1] = 21; cone[2] = 26; cone[3] = 25; 1429 cone[4] = 16; cone[5] = 23; cone[6] = 24; cone[7] = 17; 1430 ierr = DMPlexSetCone(*dm, 11, cone);CHKERRQ(ierr); 1431 cone[0] = 20; cone[1] = 28; cone[2] = 26; cone[3] = 21; 1432 cone[4] = 18; cone[5] = 17; cone[6] = 24; cone[7] = 27; 1433 ierr = DMPlexSetCone(*dm, 12, cone);CHKERRQ(ierr); 1434 cone[0] = 30; cone[1] = 28; cone[2] = 20; cone[3] = 19; 1435 cone[4] = 29; cone[5] = 15; cone[6] = 18; cone[7] = 27; 1436 ierr = DMPlexSetCone(*dm, 13, cone);CHKERRQ(ierr); 1437 cone[0] = 30; cone[1] = 19; cone[2] = 22; cone[3] = 25; 1438 cone[4] = 29; cone[5] = 23; cone[6] = 16; cone[7] = 15; 1439 ierr = DMPlexSetCone(*dm, 14, cone);CHKERRQ(ierr); 1440 } else { 1441 cone[0] = 5; cone[1] = 8; cone[2] = 7; cone[3] = 6; 1442 cone[4] = 9; cone[5] = 12; cone[6] = 11; cone[7] = 10; 1443 ierr = DMPlexSetCone(*dm, 0, cone);CHKERRQ(ierr); 1444 cone[0] = 6; cone[1] = 7; cone[2] = 14; cone[3] = 13; 1445 cone[4] = 12; cone[5] = 15; cone[6] = 16; cone[7] = 11; 1446 ierr = DMPlexSetCone(*dm, 1, cone);CHKERRQ(ierr); 1447 cone[0] = 8; cone[1] = 17; cone[2] = 14; cone[3] = 7; 1448 cone[4] = 10; cone[5] = 11; cone[6] = 16; cone[7] = 18; 1449 ierr = DMPlexSetCone(*dm, 2, cone);CHKERRQ(ierr); 1450 cone[0] = 19; cone[1] = 17; cone[2] = 8; cone[3] = 5; 1451 cone[4] = 20; cone[5] = 9; cone[6] = 10; cone[7] = 18; 1452 ierr = DMPlexSetCone(*dm, 3, cone);CHKERRQ(ierr); 1453 cone[0] = 19; cone[1] = 5; cone[2] = 6; cone[3] = 13; 1454 cone[4] = 20; cone[5] = 15; cone[6] = 12; cone[7] = 9; 1455 ierr = DMPlexSetCone(*dm, 4, cone);CHKERRQ(ierr); 1456 } 1457 } 1458 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1459 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1460 } 1461 /* Interpolate */ 1462 { 1463 DM idm; 1464 1465 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1466 ierr = DMDestroy(dm);CHKERRQ(ierr); 1467 *dm = idm; 1468 } 1469 /* Create cube geometry */ 1470 { 1471 Vec coordinates; 1472 PetscSection coordSection; 1473 PetscScalar *coords; 1474 PetscInt coordSize, v; 1475 const PetscReal dis = 1.0/PetscSqrtReal(2.0); 1476 const PetscReal ds2 = dis/2.0; 1477 1478 /* Build coordinates */ 1479 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 1480 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1481 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 1482 ierr = PetscSectionSetChart(coordSection, numCells, numCells+numVertices);CHKERRQ(ierr); 1483 for (v = numCells; v < numCells+numVertices; ++v) { 1484 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 1485 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 1486 } 1487 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1488 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1489 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 1490 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1491 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1492 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 1493 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1494 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1495 if (!rank) { 1496 coords[0*dim+0] = -ds2; coords[0*dim+1] = -ds2; coords[0*dim+2] = 0.0; 1497 coords[1*dim+0] = ds2; coords[1*dim+1] = -ds2; coords[1*dim+2] = 0.0; 1498 coords[2*dim+0] = ds2; coords[2*dim+1] = ds2; coords[2*dim+2] = 0.0; 1499 coords[3*dim+0] = -ds2; coords[3*dim+1] = ds2; coords[3*dim+2] = 0.0; 1500 coords[4*dim+0] = -ds2; coords[4*dim+1] = -ds2; coords[4*dim+2] = 1.0; 1501 coords[5*dim+0] = -ds2; coords[5*dim+1] = ds2; coords[5*dim+2] = 1.0; 1502 coords[6*dim+0] = ds2; coords[6*dim+1] = ds2; coords[6*dim+2] = 1.0; 1503 coords[7*dim+0] = ds2; coords[7*dim+1] = -ds2; coords[7*dim+2] = 1.0; 1504 coords[ 8*dim+0] = dis; coords[ 8*dim+1] = -dis; coords[ 8*dim+2] = 0.0; 1505 coords[ 9*dim+0] = dis; coords[ 9*dim+1] = dis; coords[ 9*dim+2] = 0.0; 1506 coords[10*dim+0] = dis; coords[10*dim+1] = -dis; coords[10*dim+2] = 1.0; 1507 coords[11*dim+0] = dis; coords[11*dim+1] = dis; coords[11*dim+2] = 1.0; 1508 coords[12*dim+0] = -dis; coords[12*dim+1] = dis; coords[12*dim+2] = 0.0; 1509 coords[13*dim+0] = -dis; coords[13*dim+1] = dis; coords[13*dim+2] = 1.0; 1510 coords[14*dim+0] = -dis; coords[14*dim+1] = -dis; coords[14*dim+2] = 0.0; 1511 coords[15*dim+0] = -dis; coords[15*dim+1] = -dis; coords[15*dim+2] = 1.0; 1512 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1513 /* 15 31 19 */ coords[16*dim+0] = -ds2; coords[16*dim+1] = -ds2; coords[16*dim+2] = 0.5; 1514 /* 16 32 22 */ coords[17*dim+0] = ds2; coords[17*dim+1] = -ds2; coords[17*dim+2] = 0.5; 1515 /* 17 33 21 */ coords[18*dim+0] = ds2; coords[18*dim+1] = ds2; coords[18*dim+2] = 0.5; 1516 /* 18 34 20 */ coords[19*dim+0] = -ds2; coords[19*dim+1] = ds2; coords[19*dim+2] = 0.5; 1517 /* 29 35 30 */ coords[20*dim+0] = -dis; coords[20*dim+1] = -dis; coords[20*dim+2] = 0.5; 1518 /* 23 36 25 */ coords[21*dim+0] = dis; coords[21*dim+1] = -dis; coords[21*dim+2] = 0.5; 1519 /* 24 37 26 */ coords[22*dim+0] = dis; coords[22*dim+1] = dis; coords[22*dim+2] = 0.5; 1520 /* 27 38 28 */ coords[23*dim+0] = -dis; coords[23*dim+1] = dis; coords[23*dim+2] = 0.5; 1521 } 1522 } 1523 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1524 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 1525 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1526 } 1527 /* Create periodicity */ 1528 if (periodicZ == DM_BOUNDARY_PERIODIC || periodicZ == DM_BOUNDARY_TWIST) { 1529 PetscReal L[3]; 1530 PetscReal maxCell[3]; 1531 DMBoundaryType bdType[3]; 1532 PetscReal lower[3] = {0.0, 0.0, 0.0}; 1533 PetscReal upper[3] = {1.0, 1.0, 1.5}; 1534 PetscInt i, numZCells = 3; 1535 1536 bdType[0] = DM_BOUNDARY_NONE; 1537 bdType[1] = DM_BOUNDARY_NONE; 1538 bdType[2] = periodicZ; 1539 for (i = 0; i < dim; i++) { 1540 L[i] = upper[i] - lower[i]; 1541 maxCell[i] = 1.1 * (L[i] / numZCells); 1542 } 1543 ierr = DMSetPeriodicity(*dm, PETSC_TRUE, maxCell, L, bdType);CHKERRQ(ierr); 1544 } 1545 /* Refine topology */ 1546 for (r = 0; r < numRefine; ++r) { 1547 DM rdm = NULL; 1548 1549 ierr = DMRefine(*dm, comm, &rdm);CHKERRQ(ierr); 1550 ierr = DMDestroy(dm);CHKERRQ(ierr); 1551 *dm = rdm; 1552 } 1553 /* Remap geometry to cylinder 1554 Interior square: Linear interpolation is correct 1555 The other cells all have vertices on rays from the origin. We want to uniformly expand the spacing 1556 such that the last vertex is on the unit circle. So the closest and farthest vertices are at distance 1557 1558 phi = arctan(y/x) 1559 d_close = sqrt(1/8 + 1/4 sin^2(phi)) 1560 d_far = sqrt(1/2 + sin^2(phi)) 1561 1562 so we remap them using 1563 1564 x_new = x_close + (x - x_close) (1 - d_close) / (d_far - d_close) 1565 y_new = y_close + (y - y_close) (1 - d_close) / (d_far - d_close) 1566 1567 If pi/4 < phi < 3pi/4 or -3pi/4 < phi < -pi/4, then we switch x and y. 1568 */ 1569 { 1570 Vec coordinates; 1571 PetscSection coordSection; 1572 PetscScalar *coords; 1573 PetscInt vStart, vEnd, v; 1574 const PetscReal dis = 1.0/PetscSqrtReal(2.0); 1575 const PetscReal ds2 = 0.5*dis; 1576 1577 ierr = DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1578 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 1579 ierr = DMGetCoordinatesLocal(*dm, &coordinates);CHKERRQ(ierr); 1580 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1581 for (v = vStart; v < vEnd; ++v) { 1582 PetscReal phi, sinp, cosp, dc, df, x, y, xc, yc; 1583 PetscInt off; 1584 1585 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr); 1586 if ((PetscAbsScalar(coords[off+0]) <= ds2) && (PetscAbsScalar(coords[off+1]) <= ds2)) continue; 1587 x = PetscRealPart(coords[off]); 1588 y = PetscRealPart(coords[off+1]); 1589 phi = PetscAtan2Real(y, x); 1590 sinp = PetscSinReal(phi); 1591 cosp = PetscCosReal(phi); 1592 if ((PetscAbsReal(phi) > PETSC_PI/4.0) && (PetscAbsReal(phi) < 3.0*PETSC_PI/4.0)) { 1593 dc = PetscAbsReal(ds2/sinp); 1594 df = PetscAbsReal(dis/sinp); 1595 xc = ds2*x/PetscAbsReal(y); 1596 yc = ds2*PetscSignReal(y); 1597 } else { 1598 dc = PetscAbsReal(ds2/cosp); 1599 df = PetscAbsReal(dis/cosp); 1600 xc = ds2*PetscSignReal(x); 1601 yc = ds2*y/PetscAbsReal(x); 1602 } 1603 coords[off+0] = xc + (coords[off+0] - xc)*(1.0 - dc)/(df - dc); 1604 coords[off+1] = yc + (coords[off+1] - yc)*(1.0 - dc)/(df - dc); 1605 } 1606 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1607 if (periodicZ == DM_BOUNDARY_PERIODIC || periodicZ == DM_BOUNDARY_TWIST) { 1608 ierr = DMLocalizeCoordinates(*dm);CHKERRQ(ierr); 1609 } 1610 } 1611 PetscFunctionReturn(0); 1612 } 1613 1614 /*@ 1615 DMPlexCreateWedgeCylinderMesh - Creates a mesh on the tensor product of the unit interval with the circle (cylinder) using wedges. 1616 1617 Collective on MPI_Comm 1618 1619 Input Parameters: 1620 + comm - The communicator for the DM object 1621 . n - The number of wedges around the origin 1622 - interpolate - Create edges and faces 1623 1624 Output Parameter: 1625 . dm - The DM object 1626 1627 Level: beginner 1628 1629 .keywords: DM, create 1630 .seealso: DMPlexCreateHexCylinderMesh(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 1631 @*/ 1632 PetscErrorCode DMPlexCreateWedgeCylinderMesh(MPI_Comm comm, PetscInt n, PetscBool interpolate, DM *dm) 1633 { 1634 const PetscInt dim = 3; 1635 PetscInt numCells, numVertices; 1636 PetscMPIInt rank; 1637 PetscErrorCode ierr; 1638 1639 PetscFunctionBegin; 1640 PetscValidPointer(dm, 3); 1641 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 1642 if (n < 0) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Number of wedges %D cannot be negative", n); 1643 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1644 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1645 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1646 /* Create topology */ 1647 { 1648 PetscInt cone[6], c; 1649 1650 numCells = !rank ? n : 0; 1651 numVertices = !rank ? 2*(n+1) : 0; 1652 ierr = DMPlexSetChart(*dm, 0, numCells+numVertices);CHKERRQ(ierr); 1653 ierr = DMPlexSetHybridBounds(*dm, 0, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 1654 for (c = 0; c < numCells; c++) {ierr = DMPlexSetConeSize(*dm, c, 6);CHKERRQ(ierr);} 1655 ierr = DMSetUp(*dm);CHKERRQ(ierr); 1656 for (c = 0; c < numCells; c++) { 1657 cone[0] = c+n*1; cone[1] = (c+1)%n+n*1; cone[2] = 0+3*n; 1658 cone[3] = c+n*2; cone[4] = (c+1)%n+n*2; cone[5] = 1+3*n; 1659 ierr = DMPlexSetCone(*dm, c, cone);CHKERRQ(ierr); 1660 } 1661 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1662 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1663 } 1664 /* Interpolate */ 1665 if (interpolate) { 1666 DM idm; 1667 1668 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1669 ierr = DMDestroy(dm);CHKERRQ(ierr); 1670 *dm = idm; 1671 } 1672 /* Create cylinder geometry */ 1673 { 1674 Vec coordinates; 1675 PetscSection coordSection; 1676 PetscScalar *coords; 1677 PetscInt coordSize, v, c; 1678 1679 /* Build coordinates */ 1680 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 1681 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1682 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 1683 ierr = PetscSectionSetChart(coordSection, numCells, numCells+numVertices);CHKERRQ(ierr); 1684 for (v = numCells; v < numCells+numVertices; ++v) { 1685 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 1686 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 1687 } 1688 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1689 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1690 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 1691 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1692 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1693 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 1694 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1695 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1696 for (c = 0; c < numCells; c++) { 1697 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; 1698 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; 1699 } 1700 if (!rank) { 1701 coords[(2*n+0)*dim+0] = 0.0; coords[(2*n+0)*dim+1] = 0.0; coords[(2*n+0)*dim+2] = 1.0; 1702 coords[(2*n+1)*dim+0] = 0.0; coords[(2*n+1)*dim+1] = 0.0; coords[(2*n+1)*dim+2] = 0.0; 1703 } 1704 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1705 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 1706 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1707 } 1708 PetscFunctionReturn(0); 1709 } 1710 1711 PETSC_STATIC_INLINE PetscReal DiffNormReal(PetscInt dim, const PetscReal x[], const PetscReal y[]) 1712 { 1713 PetscReal prod = 0.0; 1714 PetscInt i; 1715 for (i = 0; i < dim; ++i) prod += PetscSqr(x[i] - y[i]); 1716 return PetscSqrtReal(prod); 1717 } 1718 PETSC_STATIC_INLINE PetscReal DotReal(PetscInt dim, const PetscReal x[], const PetscReal y[]) 1719 { 1720 PetscReal prod = 0.0; 1721 PetscInt i; 1722 for (i = 0; i < dim; ++i) prod += x[i]*y[i]; 1723 return prod; 1724 } 1725 1726 /*@ 1727 DMPlexCreateSphereMesh - Creates a mesh on the d-dimensional sphere, S^d. 1728 1729 Collective on MPI_Comm 1730 1731 Input Parameters: 1732 . comm - The communicator for the DM object 1733 . dim - The dimension 1734 - simplex - Use simplices, or tensor product cells 1735 1736 Output Parameter: 1737 . dm - The DM object 1738 1739 Level: beginner 1740 1741 .keywords: DM, create 1742 .seealso: DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 1743 @*/ 1744 PetscErrorCode DMPlexCreateSphereMesh(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *dm) 1745 { 1746 const PetscInt embedDim = dim+1; 1747 PetscSection coordSection; 1748 Vec coordinates; 1749 PetscScalar *coords; 1750 PetscReal *coordsIn; 1751 PetscInt numCells, numEdges, numVerts, firstVertex, v, firstEdge, coordSize, d, c, e; 1752 PetscMPIInt rank; 1753 PetscErrorCode ierr; 1754 1755 PetscFunctionBegin; 1756 PetscValidPointer(dm, 4); 1757 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1758 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1759 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1760 ierr = DMSetCoordinateDim(*dm, dim+1);CHKERRQ(ierr); 1761 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) *dm), &rank);CHKERRQ(ierr); 1762 switch (dim) { 1763 case 2: 1764 if (simplex) { 1765 DM idm; 1766 const PetscReal edgeLen = 2.0/(1.0 + PETSC_PHI); 1767 const PetscReal vertex[3] = {0.0, 1.0/(1.0 + PETSC_PHI), PETSC_PHI/(1.0 + PETSC_PHI)}; 1768 const PetscInt degree = 5; 1769 PetscInt s[3] = {1, 1, 1}; 1770 PetscInt cone[3]; 1771 PetscInt *graph, p, i, j, k; 1772 1773 numCells = !rank ? 20 : 0; 1774 numVerts = !rank ? 12 : 0; 1775 firstVertex = numCells; 1776 /* Use icosahedron, which for a unit sphere has coordinates which are all cyclic permutations of 1777 1778 (0, \pm 1/\phi+1, \pm \phi/\phi+1) 1779 1780 where \phi^2 - \phi - 1 = 0, meaning \phi is the golden ratio \frac{1 + \sqrt{5}}{2}. The edge 1781 length is then given by 2/\phi = 2 * 2.73606 = 5.47214. 1782 */ 1783 /* Construct vertices */ 1784 ierr = PetscCalloc1(numVerts * embedDim, &coordsIn);CHKERRQ(ierr); 1785 for (p = 0, i = 0; p < embedDim; ++p) { 1786 for (s[1] = -1; s[1] < 2; s[1] += 2) { 1787 for (s[2] = -1; s[2] < 2; s[2] += 2) { 1788 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[(d+p)%embedDim]*vertex[(d+p)%embedDim]; 1789 ++i; 1790 } 1791 } 1792 } 1793 /* Construct graph */ 1794 ierr = PetscCalloc1(numVerts * numVerts, &graph);CHKERRQ(ierr); 1795 for (i = 0; i < numVerts; ++i) { 1796 for (j = 0, k = 0; j < numVerts; ++j) { 1797 if (PetscAbsReal(DiffNormReal(embedDim, &coordsIn[i*embedDim], &coordsIn[j*embedDim]) - edgeLen) < PETSC_SMALL) {graph[i*numVerts+j] = 1; ++k;} 1798 } 1799 if (k != degree) SETERRQ3(comm, PETSC_ERR_PLIB, "Invalid icosahedron, vertex %D degree %D != %D", i, k, degree); 1800 } 1801 /* Build Topology */ 1802 ierr = DMPlexSetChart(*dm, 0, numCells+numVerts);CHKERRQ(ierr); 1803 for (c = 0; c < numCells; c++) { 1804 ierr = DMPlexSetConeSize(*dm, c, embedDim);CHKERRQ(ierr); 1805 } 1806 ierr = DMSetUp(*dm);CHKERRQ(ierr); /* Allocate space for cones */ 1807 /* Cells */ 1808 for (i = 0, c = 0; i < numVerts; ++i) { 1809 for (j = 0; j < i; ++j) { 1810 for (k = 0; k < j; ++k) { 1811 if (graph[i*numVerts+j] && graph[j*numVerts+k] && graph[k*numVerts+i]) { 1812 cone[0] = firstVertex+i; cone[1] = firstVertex+j; cone[2] = firstVertex+k; 1813 /* Check orientation */ 1814 { 1815 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}}}; 1816 PetscReal normal[3]; 1817 PetscInt e, f; 1818 1819 for (d = 0; d < embedDim; ++d) { 1820 normal[d] = 0.0; 1821 for (e = 0; e < embedDim; ++e) { 1822 for (f = 0; f < embedDim; ++f) { 1823 normal[d] += epsilon[d][e][f]*(coordsIn[j*embedDim+e] - coordsIn[i*embedDim+e])*(coordsIn[k*embedDim+f] - coordsIn[i*embedDim+f]); 1824 } 1825 } 1826 } 1827 if (DotReal(embedDim, normal, &coordsIn[i*embedDim]) < 0) {PetscInt tmp = cone[1]; cone[1] = cone[2]; cone[2] = tmp;} 1828 } 1829 ierr = DMPlexSetCone(*dm, c++, cone);CHKERRQ(ierr); 1830 } 1831 } 1832 } 1833 } 1834 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1835 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1836 ierr = PetscFree(graph);CHKERRQ(ierr); 1837 /* Interpolate mesh */ 1838 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1839 ierr = DMDestroy(dm);CHKERRQ(ierr); 1840 *dm = idm; 1841 } else { 1842 /* 1843 12-21--13 1844 | | 1845 25 4 24 1846 | | 1847 12-25--9-16--8-24--13 1848 | | | | 1849 23 5 17 0 15 3 22 1850 | | | | 1851 10-20--6-14--7-19--11 1852 | | 1853 20 1 19 1854 | | 1855 10-18--11 1856 | | 1857 23 2 22 1858 | | 1859 12-21--13 1860 */ 1861 const PetscReal dist = 1.0/PetscSqrtReal(3.0); 1862 PetscInt cone[4], ornt[4]; 1863 1864 numCells = !rank ? 6 : 0; 1865 numEdges = !rank ? 12 : 0; 1866 numVerts = !rank ? 8 : 0; 1867 firstVertex = numCells; 1868 firstEdge = numCells + numVerts; 1869 /* Build Topology */ 1870 ierr = DMPlexSetChart(*dm, 0, numCells+numEdges+numVerts);CHKERRQ(ierr); 1871 for (c = 0; c < numCells; c++) { 1872 ierr = DMPlexSetConeSize(*dm, c, 4);CHKERRQ(ierr); 1873 } 1874 for (e = firstEdge; e < firstEdge+numEdges; ++e) { 1875 ierr = DMPlexSetConeSize(*dm, e, 2);CHKERRQ(ierr); 1876 } 1877 ierr = DMSetUp(*dm);CHKERRQ(ierr); /* Allocate space for cones */ 1878 /* Cell 0 */ 1879 cone[0] = 14; cone[1] = 15; cone[2] = 16; cone[3] = 17; 1880 ierr = DMPlexSetCone(*dm, 0, cone);CHKERRQ(ierr); 1881 ornt[0] = 0; ornt[1] = 0; ornt[2] = 0; ornt[3] = 0; 1882 ierr = DMPlexSetConeOrientation(*dm, 0, ornt);CHKERRQ(ierr); 1883 /* Cell 1 */ 1884 cone[0] = 18; cone[1] = 19; cone[2] = 14; cone[3] = 20; 1885 ierr = DMPlexSetCone(*dm, 1, cone);CHKERRQ(ierr); 1886 ornt[0] = 0; ornt[1] = 0; ornt[2] = -2; ornt[3] = 0; 1887 ierr = DMPlexSetConeOrientation(*dm, 1, ornt);CHKERRQ(ierr); 1888 /* Cell 2 */ 1889 cone[0] = 21; cone[1] = 22; cone[2] = 18; cone[3] = 23; 1890 ierr = DMPlexSetCone(*dm, 2, cone);CHKERRQ(ierr); 1891 ornt[0] = 0; ornt[1] = 0; ornt[2] = -2; ornt[3] = 0; 1892 ierr = DMPlexSetConeOrientation(*dm, 2, ornt);CHKERRQ(ierr); 1893 /* Cell 3 */ 1894 cone[0] = 19; cone[1] = 22; cone[2] = 24; cone[3] = 15; 1895 ierr = DMPlexSetCone(*dm, 3, cone);CHKERRQ(ierr); 1896 ornt[0] = -2; ornt[1] = -2; ornt[2] = 0; ornt[3] = -2; 1897 ierr = DMPlexSetConeOrientation(*dm, 3, ornt);CHKERRQ(ierr); 1898 /* Cell 4 */ 1899 cone[0] = 16; cone[1] = 24; cone[2] = 21; cone[3] = 25; 1900 ierr = DMPlexSetCone(*dm, 4, cone);CHKERRQ(ierr); 1901 ornt[0] = -2; ornt[1] = -2; ornt[2] = -2; ornt[3] = 0; 1902 ierr = DMPlexSetConeOrientation(*dm, 4, ornt);CHKERRQ(ierr); 1903 /* Cell 5 */ 1904 cone[0] = 20; cone[1] = 17; cone[2] = 25; cone[3] = 23; 1905 ierr = DMPlexSetCone(*dm, 5, cone);CHKERRQ(ierr); 1906 ornt[0] = -2; ornt[1] = -2; ornt[2] = -2; ornt[3] = -2; 1907 ierr = DMPlexSetConeOrientation(*dm, 5, ornt);CHKERRQ(ierr); 1908 /* Edges */ 1909 cone[0] = 6; cone[1] = 7; 1910 ierr = DMPlexSetCone(*dm, 14, cone);CHKERRQ(ierr); 1911 cone[0] = 7; cone[1] = 8; 1912 ierr = DMPlexSetCone(*dm, 15, cone);CHKERRQ(ierr); 1913 cone[0] = 8; cone[1] = 9; 1914 ierr = DMPlexSetCone(*dm, 16, cone);CHKERRQ(ierr); 1915 cone[0] = 9; cone[1] = 6; 1916 ierr = DMPlexSetCone(*dm, 17, cone);CHKERRQ(ierr); 1917 cone[0] = 10; cone[1] = 11; 1918 ierr = DMPlexSetCone(*dm, 18, cone);CHKERRQ(ierr); 1919 cone[0] = 11; cone[1] = 7; 1920 ierr = DMPlexSetCone(*dm, 19, cone);CHKERRQ(ierr); 1921 cone[0] = 6; cone[1] = 10; 1922 ierr = DMPlexSetCone(*dm, 20, cone);CHKERRQ(ierr); 1923 cone[0] = 12; cone[1] = 13; 1924 ierr = DMPlexSetCone(*dm, 21, cone);CHKERRQ(ierr); 1925 cone[0] = 13; cone[1] = 11; 1926 ierr = DMPlexSetCone(*dm, 22, cone);CHKERRQ(ierr); 1927 cone[0] = 10; cone[1] = 12; 1928 ierr = DMPlexSetCone(*dm, 23, cone);CHKERRQ(ierr); 1929 cone[0] = 13; cone[1] = 8; 1930 ierr = DMPlexSetCone(*dm, 24, cone);CHKERRQ(ierr); 1931 cone[0] = 12; cone[1] = 9; 1932 ierr = DMPlexSetCone(*dm, 25, cone);CHKERRQ(ierr); 1933 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1934 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1935 /* Build coordinates */ 1936 ierr = PetscCalloc1(numVerts * embedDim, &coordsIn);CHKERRQ(ierr); 1937 coordsIn[0*embedDim+0] = -dist; coordsIn[0*embedDim+1] = dist; coordsIn[0*embedDim+2] = -dist; 1938 coordsIn[1*embedDim+0] = dist; coordsIn[1*embedDim+1] = dist; coordsIn[1*embedDim+2] = -dist; 1939 coordsIn[2*embedDim+0] = dist; coordsIn[2*embedDim+1] = -dist; coordsIn[2*embedDim+2] = -dist; 1940 coordsIn[3*embedDim+0] = -dist; coordsIn[3*embedDim+1] = -dist; coordsIn[3*embedDim+2] = -dist; 1941 coordsIn[4*embedDim+0] = -dist; coordsIn[4*embedDim+1] = dist; coordsIn[4*embedDim+2] = dist; 1942 coordsIn[5*embedDim+0] = dist; coordsIn[5*embedDim+1] = dist; coordsIn[5*embedDim+2] = dist; 1943 coordsIn[6*embedDim+0] = -dist; coordsIn[6*embedDim+1] = -dist; coordsIn[6*embedDim+2] = dist; 1944 coordsIn[7*embedDim+0] = dist; coordsIn[7*embedDim+1] = -dist; coordsIn[7*embedDim+2] = dist; 1945 } 1946 break; 1947 case 3: 1948 if (simplex) { 1949 DM idm; 1950 const PetscReal edgeLen = 1.0/PETSC_PHI; 1951 const PetscReal vertexA[4] = {0.5, 0.5, 0.5, 0.5}; 1952 const PetscReal vertexB[4] = {1.0, 0.0, 0.0, 0.0}; 1953 const PetscReal vertexC[4] = {0.5, 0.5*PETSC_PHI, 0.5/PETSC_PHI, 0.0}; 1954 const PetscInt degree = 12; 1955 PetscInt s[4] = {1, 1, 1}; 1956 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}, 1957 {2, 0, 1, 3}, {2, 1, 3, 0}, {2, 3, 0, 1}, {3, 0, 2, 1}, {3, 1, 0, 2}, {3, 2, 1, 0}}; 1958 PetscInt cone[4]; 1959 PetscInt *graph, p, i, j, k, l; 1960 1961 numCells = !rank ? 600 : 0; 1962 numVerts = !rank ? 120 : 0; 1963 firstVertex = numCells; 1964 /* Use the 600-cell, which for a unit sphere has coordinates which are 1965 1966 1/2 (\pm 1, \pm 1, \pm 1, \pm 1) 16 1967 (\pm 1, 0, 0, 0) all cyclic permutations 8 1968 1/2 (\pm 1, \pm phi, \pm 1/phi, 0) all even permutations 96 1969 1970 where \phi^2 - \phi - 1 = 0, meaning \phi is the golden ratio \frac{1 + \sqrt{5}}{2}. The edge 1971 length is then given by 1/\phi = 2.73606. 1972 1973 http://buzzard.pugetsound.edu/sage-practice/ch03s03.html 1974 http://mathworld.wolfram.com/600-Cell.html 1975 */ 1976 /* Construct vertices */ 1977 ierr = PetscCalloc1(numVerts * embedDim, &coordsIn);CHKERRQ(ierr); 1978 i = 0; 1979 for (s[0] = -1; s[0] < 2; s[0] += 2) { 1980 for (s[1] = -1; s[1] < 2; s[1] += 2) { 1981 for (s[2] = -1; s[2] < 2; s[2] += 2) { 1982 for (s[3] = -1; s[3] < 2; s[3] += 2) { 1983 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[d]*vertexA[d]; 1984 ++i; 1985 } 1986 } 1987 } 1988 } 1989 for (p = 0; p < embedDim; ++p) { 1990 s[1] = s[2] = s[3] = 1; 1991 for (s[0] = -1; s[0] < 2; s[0] += 2) { 1992 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[(d+p)%embedDim]*vertexB[(d+p)%embedDim]; 1993 ++i; 1994 } 1995 } 1996 for (p = 0; p < 12; ++p) { 1997 s[3] = 1; 1998 for (s[0] = -1; s[0] < 2; s[0] += 2) { 1999 for (s[1] = -1; s[1] < 2; s[1] += 2) { 2000 for (s[2] = -1; s[2] < 2; s[2] += 2) { 2001 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[evenPerm[p][d]]*vertexC[evenPerm[p][d]]; 2002 ++i; 2003 } 2004 } 2005 } 2006 } 2007 if (i != numVerts) SETERRQ2(comm, PETSC_ERR_PLIB, "Invalid 600-cell, vertices %D != %D", i, numVerts); 2008 /* Construct graph */ 2009 ierr = PetscCalloc1(numVerts * numVerts, &graph);CHKERRQ(ierr); 2010 for (i = 0; i < numVerts; ++i) { 2011 for (j = 0, k = 0; j < numVerts; ++j) { 2012 if (PetscAbsReal(DiffNormReal(embedDim, &coordsIn[i*embedDim], &coordsIn[j*embedDim]) - edgeLen) < PETSC_SMALL) {graph[i*numVerts+j] = 1; ++k;} 2013 } 2014 if (k != degree) SETERRQ3(comm, PETSC_ERR_PLIB, "Invalid 600-cell, vertex %D degree %D != %D", i, k, degree); 2015 } 2016 /* Build Topology */ 2017 ierr = DMPlexSetChart(*dm, 0, numCells+numVerts);CHKERRQ(ierr); 2018 for (c = 0; c < numCells; c++) { 2019 ierr = DMPlexSetConeSize(*dm, c, embedDim);CHKERRQ(ierr); 2020 } 2021 ierr = DMSetUp(*dm);CHKERRQ(ierr); /* Allocate space for cones */ 2022 /* Cells */ 2023 for (i = 0, c = 0; i < numVerts; ++i) { 2024 for (j = 0; j < i; ++j) { 2025 for (k = 0; k < j; ++k) { 2026 for (l = 0; l < k; ++l) { 2027 if (graph[i*numVerts+j] && graph[j*numVerts+k] && graph[k*numVerts+i] && 2028 graph[l*numVerts+i] && graph[l*numVerts+j] && graph[l*numVerts+k]) { 2029 cone[0] = firstVertex+i; cone[1] = firstVertex+j; cone[2] = firstVertex+k; cone[3] = firstVertex+l; 2030 /* Check orientation: https://ef.gy/linear-algebra:normal-vectors-in-higher-dimensional-spaces */ 2031 { 2032 const PetscInt epsilon[4][4][4][4] = {{{{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 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, -1}, { 0, 0, 0, 0}, { 0, 1, 0, 0}}, 2035 {{0, 0, 0, 0}, { 0, 0, 1, 0}, { 0, -1, 0, 0}, { 0, 0, 0, 0}}}, 2036 2037 {{{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, -1}, { 0, 0, 1, 0}}, 2038 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 2039 {{0, 0, 0, 1}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, {-1, 0, 0, 0}}, 2040 {{0, 0, -1, 0}, { 0, 0, 0, 0}, { 1, 0, 0, 0}, { 0, 0, 0, 0}}}, 2041 2042 {{{0, 0, 0, 0}, { 0, 0, 0, 1}, { 0, 0, 0, 0}, { 0, -1, 0, 0}}, 2043 {{0, 0, 0, -1}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 1, 0, 0, 0}}, 2044 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 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 2047 {{{0, 0, 0, 0}, { 0, 0, -1, 0}, { 0, 1, 0, 0}, { 0, 0, 0, 0}}, 2048 {{0, 0, 1, 0}, { 0, 0, 0, 0}, {-1, 0, 0, 0}, { 0, 0, 0, 0}}, 2049 {{0, -1, 0, 0}, { 1, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 2050 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}}}; 2051 PetscReal normal[4]; 2052 PetscInt e, f, g; 2053 2054 for (d = 0; d < embedDim; ++d) { 2055 normal[d] = 0.0; 2056 for (e = 0; e < embedDim; ++e) { 2057 for (f = 0; f < embedDim; ++f) { 2058 for (g = 0; g < embedDim; ++g) { 2059 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]); 2060 } 2061 } 2062 } 2063 } 2064 if (DotReal(embedDim, normal, &coordsIn[i*embedDim]) < 0) {PetscInt tmp = cone[1]; cone[1] = cone[2]; cone[2] = tmp;} 2065 } 2066 ierr = DMPlexSetCone(*dm, c++, cone);CHKERRQ(ierr); 2067 } 2068 } 2069 } 2070 } 2071 } 2072 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 2073 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 2074 ierr = PetscFree(graph);CHKERRQ(ierr); 2075 /* Interpolate mesh */ 2076 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 2077 ierr = DMDestroy(dm);CHKERRQ(ierr); 2078 *dm = idm; 2079 break; 2080 } 2081 default: SETERRQ1(comm, PETSC_ERR_SUP, "Unsupported dimension for sphere: %D", dim); 2082 } 2083 /* Create coordinates */ 2084 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 2085 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 2086 ierr = PetscSectionSetFieldComponents(coordSection, 0, embedDim);CHKERRQ(ierr); 2087 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVerts);CHKERRQ(ierr); 2088 for (v = firstVertex; v < firstVertex+numVerts; ++v) { 2089 ierr = PetscSectionSetDof(coordSection, v, embedDim);CHKERRQ(ierr); 2090 ierr = PetscSectionSetFieldDof(coordSection, v, 0, embedDim);CHKERRQ(ierr); 2091 } 2092 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 2093 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 2094 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 2095 ierr = VecSetBlockSize(coordinates, embedDim);CHKERRQ(ierr); 2096 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 2097 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 2098 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 2099 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2100 for (v = 0; v < numVerts; ++v) for (d = 0; d < embedDim; ++d) {coords[v*embedDim+d] = coordsIn[v*embedDim+d];} 2101 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2102 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 2103 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 2104 ierr = PetscFree(coordsIn);CHKERRQ(ierr); 2105 PetscFunctionReturn(0); 2106 } 2107 2108 /* External function declarations here */ 2109 extern PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling); 2110 extern PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, Mat *mat); 2111 extern PetscErrorCode DMCreateMassMatrix_Plex(DM dmCoarse, DM dmFine, Mat *mat); 2112 extern PetscErrorCode DMCreateDefaultSection_Plex(DM dm); 2113 extern PetscErrorCode DMCreateDefaultConstraints_Plex(DM dm); 2114 extern PetscErrorCode DMCreateMatrix_Plex(DM dm, Mat *J); 2115 extern PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm); 2116 extern PetscErrorCode DMCreateCoordinateField_Plex(DM dm, DMField *field); 2117 PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm); 2118 extern PetscErrorCode DMSetUp_Plex(DM dm); 2119 extern PetscErrorCode DMDestroy_Plex(DM dm); 2120 extern PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer); 2121 extern PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer); 2122 extern PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm); 2123 extern PetscErrorCode DMCreateSuperDM_Plex(DM dms[], PetscInt len, IS **is, DM *superdm); 2124 static PetscErrorCode DMInitialize_Plex(DM dm); 2125 2126 /* Replace dm with the contents of dmNew 2127 - Share the DM_Plex structure 2128 - Share the coordinates 2129 - Share the SF 2130 */ 2131 static PetscErrorCode DMPlexReplace_Static(DM dm, DM dmNew) 2132 { 2133 PetscSF sf; 2134 DM coordDM, coarseDM; 2135 Vec coords; 2136 PetscBool isper; 2137 const PetscReal *maxCell, *L; 2138 const DMBoundaryType *bd; 2139 PetscErrorCode ierr; 2140 2141 PetscFunctionBegin; 2142 ierr = DMGetPointSF(dmNew, &sf);CHKERRQ(ierr); 2143 ierr = DMSetPointSF(dm, sf);CHKERRQ(ierr); 2144 ierr = DMGetCoordinateDM(dmNew, &coordDM);CHKERRQ(ierr); 2145 ierr = DMGetCoordinatesLocal(dmNew, &coords);CHKERRQ(ierr); 2146 ierr = DMSetCoordinateDM(dm, coordDM);CHKERRQ(ierr); 2147 ierr = DMSetCoordinatesLocal(dm, coords);CHKERRQ(ierr); 2148 ierr = DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);CHKERRQ(ierr); 2149 ierr = DMSetPeriodicity(dmNew, isper, maxCell, L, bd);CHKERRQ(ierr); 2150 ierr = DMDestroy_Plex(dm);CHKERRQ(ierr); 2151 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 2152 dm->data = dmNew->data; 2153 ((DM_Plex *) dmNew->data)->refct++; 2154 dmNew->labels->refct++; 2155 if (!--(dm->labels->refct)) { 2156 DMLabelLink next = dm->labels->next; 2157 2158 /* destroy the labels */ 2159 while (next) { 2160 DMLabelLink tmp = next->next; 2161 2162 ierr = DMLabelDestroy(&next->label);CHKERRQ(ierr); 2163 ierr = PetscFree(next);CHKERRQ(ierr); 2164 next = tmp; 2165 } 2166 ierr = PetscFree(dm->labels);CHKERRQ(ierr); 2167 } 2168 dm->labels = dmNew->labels; 2169 dm->depthLabel = dmNew->depthLabel; 2170 ierr = DMGetCoarseDM(dmNew,&coarseDM);CHKERRQ(ierr); 2171 ierr = DMSetCoarseDM(dm,coarseDM);CHKERRQ(ierr); 2172 PetscFunctionReturn(0); 2173 } 2174 2175 /* Swap dm with the contents of dmNew 2176 - Swap the DM_Plex structure 2177 - Swap the coordinates 2178 - Swap the point PetscSF 2179 */ 2180 static PetscErrorCode DMPlexSwap_Static(DM dmA, DM dmB) 2181 { 2182 DM coordDMA, coordDMB; 2183 Vec coordsA, coordsB; 2184 PetscSF sfA, sfB; 2185 void *tmp; 2186 DMLabelLinkList listTmp; 2187 DMLabel depthTmp; 2188 PetscErrorCode ierr; 2189 2190 PetscFunctionBegin; 2191 ierr = DMGetPointSF(dmA, &sfA);CHKERRQ(ierr); 2192 ierr = DMGetPointSF(dmB, &sfB);CHKERRQ(ierr); 2193 ierr = PetscObjectReference((PetscObject) sfA);CHKERRQ(ierr); 2194 ierr = DMSetPointSF(dmA, sfB);CHKERRQ(ierr); 2195 ierr = DMSetPointSF(dmB, sfA);CHKERRQ(ierr); 2196 ierr = PetscObjectDereference((PetscObject) sfA);CHKERRQ(ierr); 2197 2198 ierr = DMGetCoordinateDM(dmA, &coordDMA);CHKERRQ(ierr); 2199 ierr = DMGetCoordinateDM(dmB, &coordDMB);CHKERRQ(ierr); 2200 ierr = PetscObjectReference((PetscObject) coordDMA);CHKERRQ(ierr); 2201 ierr = DMSetCoordinateDM(dmA, coordDMB);CHKERRQ(ierr); 2202 ierr = DMSetCoordinateDM(dmB, coordDMA);CHKERRQ(ierr); 2203 ierr = PetscObjectDereference((PetscObject) coordDMA);CHKERRQ(ierr); 2204 2205 ierr = DMGetCoordinatesLocal(dmA, &coordsA);CHKERRQ(ierr); 2206 ierr = DMGetCoordinatesLocal(dmB, &coordsB);CHKERRQ(ierr); 2207 ierr = PetscObjectReference((PetscObject) coordsA);CHKERRQ(ierr); 2208 ierr = DMSetCoordinatesLocal(dmA, coordsB);CHKERRQ(ierr); 2209 ierr = DMSetCoordinatesLocal(dmB, coordsA);CHKERRQ(ierr); 2210 ierr = PetscObjectDereference((PetscObject) coordsA);CHKERRQ(ierr); 2211 2212 tmp = dmA->data; 2213 dmA->data = dmB->data; 2214 dmB->data = tmp; 2215 listTmp = dmA->labels; 2216 dmA->labels = dmB->labels; 2217 dmB->labels = listTmp; 2218 depthTmp = dmA->depthLabel; 2219 dmA->depthLabel = dmB->depthLabel; 2220 dmB->depthLabel = depthTmp; 2221 PetscFunctionReturn(0); 2222 } 2223 2224 PetscErrorCode DMSetFromOptions_NonRefinement_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 2225 { 2226 DM_Plex *mesh = (DM_Plex*) dm->data; 2227 PetscErrorCode ierr; 2228 2229 PetscFunctionBegin; 2230 /* Handle viewing */ 2231 ierr = PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMPlexMatSetClosure", PETSC_FALSE, &mesh->printSetValues, NULL);CHKERRQ(ierr); 2232 ierr = PetscOptionsInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMPlexSNESComputeResidualFEM", 0, &mesh->printFEM, NULL);CHKERRQ(ierr); 2233 ierr = PetscOptionsReal("-dm_plex_print_tol", "Tolerance for FEM output", "DMPlexSNESComputeResidualFEM", mesh->printTol, &mesh->printTol, NULL);CHKERRQ(ierr); 2234 ierr = PetscOptionsInt("-dm_plex_print_l2", "Debug output level all L2 diff computations", "DMComputeL2Diff", 0, &mesh->printL2, NULL);CHKERRQ(ierr); 2235 /* Point Location */ 2236 ierr = PetscOptionsBool("-dm_plex_hash_location", "Use grid hashing for point location", "DMInterpolate", PETSC_FALSE, &mesh->useHashLocation, NULL);CHKERRQ(ierr); 2237 /* Partitioning and distribution */ 2238 ierr = PetscOptionsBool("-dm_plex_partition_balance", "Attempt to evenly divide points on partition boundary between processes", "DMPlexSetPartitionBalance", PETSC_FALSE, &mesh->partitionBalance, NULL);CHKERRQ(ierr); 2239 /* Generation and remeshing */ 2240 ierr = PetscOptionsBool("-dm_plex_remesh_bd", "Allow changes to the boundary on remeshing", "DMAdapt", PETSC_FALSE, &mesh->remeshBd, NULL);CHKERRQ(ierr); 2241 /* Projection behavior */ 2242 ierr = PetscOptionsInt("-dm_plex_max_projection_height", "Maxmimum mesh point height used to project locally", "DMPlexSetMaxProjectionHeight", 0, &mesh->maxProjectionHeight, NULL);CHKERRQ(ierr); 2243 ierr = PetscOptionsBool("-dm_plex_regular_refinement", "Use special nested projection algorithm for regular refinement", "DMPlexSetRegularRefinement", mesh->regularRefinement, &mesh->regularRefinement, NULL);CHKERRQ(ierr); 2244 2245 ierr = PetscPartitionerSetFromOptions(mesh->partitioner);CHKERRQ(ierr); 2246 PetscFunctionReturn(0); 2247 } 2248 2249 static PetscErrorCode DMSetFromOptions_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 2250 { 2251 PetscInt refine = 0, coarsen = 0, r; 2252 PetscBool isHierarchy; 2253 PetscErrorCode ierr; 2254 2255 PetscFunctionBegin; 2256 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2257 ierr = PetscOptionsHead(PetscOptionsObject,"DMPlex Options");CHKERRQ(ierr); 2258 /* Handle DMPlex refinement */ 2259 ierr = PetscOptionsInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL);CHKERRQ(ierr); 2260 ierr = PetscOptionsInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy);CHKERRQ(ierr); 2261 if (refine) {ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE);CHKERRQ(ierr);} 2262 if (refine && isHierarchy) { 2263 DM *dms, coarseDM; 2264 2265 ierr = DMGetCoarseDM(dm, &coarseDM);CHKERRQ(ierr); 2266 ierr = PetscObjectReference((PetscObject)coarseDM);CHKERRQ(ierr); 2267 ierr = PetscMalloc1(refine,&dms);CHKERRQ(ierr); 2268 ierr = DMRefineHierarchy(dm, refine, dms);CHKERRQ(ierr); 2269 /* Total hack since we do not pass in a pointer */ 2270 ierr = DMPlexSwap_Static(dm, dms[refine-1]);CHKERRQ(ierr); 2271 if (refine == 1) { 2272 ierr = DMSetCoarseDM(dm, dms[0]);CHKERRQ(ierr); 2273 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr); 2274 } else { 2275 ierr = DMSetCoarseDM(dm, dms[refine-2]);CHKERRQ(ierr); 2276 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr); 2277 ierr = DMSetCoarseDM(dms[0], dms[refine-1]);CHKERRQ(ierr); 2278 ierr = DMPlexSetRegularRefinement(dms[0], PETSC_TRUE);CHKERRQ(ierr); 2279 } 2280 ierr = DMSetCoarseDM(dms[refine-1], coarseDM);CHKERRQ(ierr); 2281 ierr = PetscObjectDereference((PetscObject)coarseDM);CHKERRQ(ierr); 2282 /* Free DMs */ 2283 for (r = 0; r < refine; ++r) { 2284 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr); 2285 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 2286 } 2287 ierr = PetscFree(dms);CHKERRQ(ierr); 2288 } else { 2289 for (r = 0; r < refine; ++r) { 2290 DM refinedMesh; 2291 2292 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2293 ierr = DMRefine(dm, PetscObjectComm((PetscObject) dm), &refinedMesh);CHKERRQ(ierr); 2294 /* Total hack since we do not pass in a pointer */ 2295 ierr = DMPlexReplace_Static(dm, refinedMesh);CHKERRQ(ierr); 2296 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2297 ierr = DMDestroy(&refinedMesh);CHKERRQ(ierr); 2298 } 2299 } 2300 /* Handle DMPlex coarsening */ 2301 ierr = PetscOptionsInt("-dm_coarsen", "Coarsen the mesh", "DMCreate", coarsen, &coarsen, NULL);CHKERRQ(ierr); 2302 ierr = PetscOptionsInt("-dm_coarsen_hierarchy", "The number of coarsenings", "DMCreate", coarsen, &coarsen, &isHierarchy);CHKERRQ(ierr); 2303 if (coarsen && isHierarchy) { 2304 DM *dms; 2305 2306 ierr = PetscMalloc1(coarsen, &dms);CHKERRQ(ierr); 2307 ierr = DMCoarsenHierarchy(dm, coarsen, dms);CHKERRQ(ierr); 2308 /* Free DMs */ 2309 for (r = 0; r < coarsen; ++r) { 2310 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr); 2311 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 2312 } 2313 ierr = PetscFree(dms);CHKERRQ(ierr); 2314 } else { 2315 for (r = 0; r < coarsen; ++r) { 2316 DM coarseMesh; 2317 2318 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2319 ierr = DMCoarsen(dm, PetscObjectComm((PetscObject) dm), &coarseMesh);CHKERRQ(ierr); 2320 /* Total hack since we do not pass in a pointer */ 2321 ierr = DMPlexReplace_Static(dm, coarseMesh);CHKERRQ(ierr); 2322 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2323 ierr = DMDestroy(&coarseMesh);CHKERRQ(ierr); 2324 } 2325 } 2326 /* Handle */ 2327 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2328 ierr = PetscOptionsTail();CHKERRQ(ierr); 2329 PetscFunctionReturn(0); 2330 } 2331 2332 static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec) 2333 { 2334 PetscErrorCode ierr; 2335 2336 PetscFunctionBegin; 2337 ierr = DMCreateGlobalVector_Section_Private(dm,vec);CHKERRQ(ierr); 2338 /* ierr = VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM);CHKERRQ(ierr); */ 2339 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex);CHKERRQ(ierr); 2340 ierr = VecSetOperation(*vec, VECOP_VIEWNATIVE, (void (*)(void)) VecView_Plex_Native);CHKERRQ(ierr); 2341 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex);CHKERRQ(ierr); 2342 ierr = VecSetOperation(*vec, VECOP_LOADNATIVE, (void (*)(void)) VecLoad_Plex_Native);CHKERRQ(ierr); 2343 PetscFunctionReturn(0); 2344 } 2345 2346 static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec) 2347 { 2348 PetscErrorCode ierr; 2349 2350 PetscFunctionBegin; 2351 ierr = DMCreateLocalVector_Section_Private(dm,vec);CHKERRQ(ierr); 2352 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex_Local);CHKERRQ(ierr); 2353 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex_Local);CHKERRQ(ierr); 2354 PetscFunctionReturn(0); 2355 } 2356 2357 static PetscErrorCode DMGetDimPoints_Plex(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd) 2358 { 2359 PetscInt depth, d; 2360 PetscErrorCode ierr; 2361 2362 PetscFunctionBegin; 2363 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 2364 if (depth == 1) { 2365 ierr = DMGetDimension(dm, &d);CHKERRQ(ierr); 2366 if (dim == 0) {ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr);} 2367 else if (dim == d) {ierr = DMPlexGetDepthStratum(dm, 1, pStart, pEnd);CHKERRQ(ierr);} 2368 else {*pStart = 0; *pEnd = 0;} 2369 } else { 2370 ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr); 2371 } 2372 PetscFunctionReturn(0); 2373 } 2374 2375 static PetscErrorCode DMGetNeighors_Plex(DM dm, PetscInt *nranks, const PetscMPIInt *ranks[]) 2376 { 2377 PetscSF sf; 2378 PetscErrorCode ierr; 2379 2380 PetscFunctionBegin; 2381 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 2382 ierr = PetscSFGetRanks(sf, nranks, ranks, NULL, NULL, NULL);CHKERRQ(ierr); 2383 PetscFunctionReturn(0); 2384 } 2385 2386 static PetscErrorCode DMHasCreateInjection_Plex(DM dm, PetscBool *flg) 2387 { 2388 PetscFunctionBegin; 2389 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2390 PetscValidPointer(flg,2); 2391 *flg = PETSC_TRUE; 2392 PetscFunctionReturn(0); 2393 } 2394 2395 static PetscErrorCode DMInitialize_Plex(DM dm) 2396 { 2397 PetscErrorCode ierr; 2398 2399 PetscFunctionBegin; 2400 dm->ops->view = DMView_Plex; 2401 dm->ops->load = DMLoad_Plex; 2402 dm->ops->setfromoptions = DMSetFromOptions_Plex; 2403 dm->ops->clone = DMClone_Plex; 2404 dm->ops->setup = DMSetUp_Plex; 2405 dm->ops->createdefaultsection = DMCreateDefaultSection_Plex; 2406 dm->ops->createdefaultconstraints = DMCreateDefaultConstraints_Plex; 2407 dm->ops->createglobalvector = DMCreateGlobalVector_Plex; 2408 dm->ops->createlocalvector = DMCreateLocalVector_Plex; 2409 dm->ops->getlocaltoglobalmapping = NULL; 2410 dm->ops->createfieldis = NULL; 2411 dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex; 2412 dm->ops->createcoordinatefield = DMCreateCoordinateField_Plex; 2413 dm->ops->getcoloring = NULL; 2414 dm->ops->creatematrix = DMCreateMatrix_Plex; 2415 dm->ops->createinterpolation = DMCreateInterpolation_Plex; 2416 dm->ops->createmassmatrix = DMCreateMassMatrix_Plex; 2417 dm->ops->getaggregates = NULL; 2418 dm->ops->getinjection = DMCreateInjection_Plex; 2419 dm->ops->hascreateinjection = DMHasCreateInjection_Plex; 2420 dm->ops->refine = DMRefine_Plex; 2421 dm->ops->coarsen = DMCoarsen_Plex; 2422 dm->ops->refinehierarchy = DMRefineHierarchy_Plex; 2423 dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Plex; 2424 dm->ops->adaptlabel = DMAdaptLabel_Plex; 2425 dm->ops->adaptmetric = DMAdaptMetric_Plex; 2426 dm->ops->globaltolocalbegin = NULL; 2427 dm->ops->globaltolocalend = NULL; 2428 dm->ops->localtoglobalbegin = NULL; 2429 dm->ops->localtoglobalend = NULL; 2430 dm->ops->destroy = DMDestroy_Plex; 2431 dm->ops->createsubdm = DMCreateSubDM_Plex; 2432 dm->ops->createsuperdm = DMCreateSuperDM_Plex; 2433 dm->ops->getdimpoints = DMGetDimPoints_Plex; 2434 dm->ops->locatepoints = DMLocatePoints_Plex; 2435 dm->ops->projectfunctionlocal = DMProjectFunctionLocal_Plex; 2436 dm->ops->projectfunctionlabellocal = DMProjectFunctionLabelLocal_Plex; 2437 dm->ops->projectfieldlocal = DMProjectFieldLocal_Plex; 2438 dm->ops->projectfieldlabellocal = DMProjectFieldLabelLocal_Plex; 2439 dm->ops->computel2diff = DMComputeL2Diff_Plex; 2440 dm->ops->computel2gradientdiff = DMComputeL2GradientDiff_Plex; 2441 dm->ops->computel2fielddiff = DMComputeL2FieldDiff_Plex; 2442 dm->ops->getneighbors = DMGetNeighors_Plex; 2443 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMPlexInsertBoundaryValues_C",DMPlexInsertBoundaryValues_Plex);CHKERRQ(ierr); 2444 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMSetUpGLVisViewer_C",DMSetUpGLVisViewer_Plex);CHKERRQ(ierr); 2445 PetscFunctionReturn(0); 2446 } 2447 2448 PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm) 2449 { 2450 DM_Plex *mesh = (DM_Plex *) dm->data; 2451 PetscErrorCode ierr; 2452 2453 PetscFunctionBegin; 2454 mesh->refct++; 2455 (*newdm)->data = mesh; 2456 ierr = PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);CHKERRQ(ierr); 2457 ierr = DMInitialize_Plex(*newdm);CHKERRQ(ierr); 2458 PetscFunctionReturn(0); 2459 } 2460 2461 /*MC 2462 DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram. 2463 In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is 2464 specified by a PetscSection object. Ownership in the global representation is determined by 2465 ownership of the underlying DMPlex points. This is specified by another PetscSection object. 2466 2467 Options Database Keys: 2468 + -dm_view :mesh.tex:ascii_latex - View the mesh in LaTeX/TikZ 2469 . -dm_plex_view_scale <num> - Scale the TikZ 2470 - -dm_plex_print_fem <num> - View FEM assembly information, such as element vectors and matrices 2471 2472 2473 Level: intermediate 2474 2475 .seealso: DMType, DMPlexCreate(), DMCreate(), DMSetType() 2476 M*/ 2477 2478 PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm) 2479 { 2480 DM_Plex *mesh; 2481 PetscInt unit, d; 2482 PetscErrorCode ierr; 2483 2484 PetscFunctionBegin; 2485 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2486 ierr = PetscNewLog(dm,&mesh);CHKERRQ(ierr); 2487 dm->dim = 0; 2488 dm->data = mesh; 2489 2490 mesh->refct = 1; 2491 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);CHKERRQ(ierr); 2492 mesh->maxConeSize = 0; 2493 mesh->cones = NULL; 2494 mesh->coneOrientations = NULL; 2495 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);CHKERRQ(ierr); 2496 mesh->maxSupportSize = 0; 2497 mesh->supports = NULL; 2498 mesh->refinementUniform = PETSC_TRUE; 2499 mesh->refinementLimit = -1.0; 2500 2501 mesh->facesTmp = NULL; 2502 2503 mesh->tetgenOpts = NULL; 2504 mesh->triangleOpts = NULL; 2505 ierr = PetscPartitionerCreate(PetscObjectComm((PetscObject)dm), &mesh->partitioner);CHKERRQ(ierr); 2506 mesh->remeshBd = PETSC_FALSE; 2507 2508 mesh->subpointMap = NULL; 2509 2510 for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0; 2511 2512 mesh->regularRefinement = PETSC_FALSE; 2513 mesh->depthState = -1; 2514 mesh->globalVertexNumbers = NULL; 2515 mesh->globalCellNumbers = NULL; 2516 mesh->anchorSection = NULL; 2517 mesh->anchorIS = NULL; 2518 mesh->createanchors = NULL; 2519 mesh->computeanchormatrix = NULL; 2520 mesh->parentSection = NULL; 2521 mesh->parents = NULL; 2522 mesh->childIDs = NULL; 2523 mesh->childSection = NULL; 2524 mesh->children = NULL; 2525 mesh->referenceTree = NULL; 2526 mesh->getchildsymmetry = NULL; 2527 for (d = 0; d < 8; ++d) mesh->hybridPointMax[d] = PETSC_DETERMINE; 2528 mesh->vtkCellHeight = 0; 2529 mesh->useAnchors = PETSC_FALSE; 2530 2531 mesh->maxProjectionHeight = 0; 2532 2533 mesh->printSetValues = PETSC_FALSE; 2534 mesh->printFEM = 0; 2535 mesh->printTol = 1.0e-10; 2536 2537 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 2538 PetscFunctionReturn(0); 2539 } 2540 2541 /*@ 2542 DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram. 2543 2544 Collective on MPI_Comm 2545 2546 Input Parameter: 2547 . comm - The communicator for the DMPlex object 2548 2549 Output Parameter: 2550 . mesh - The DMPlex object 2551 2552 Level: beginner 2553 2554 .keywords: DMPlex, create 2555 @*/ 2556 PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh) 2557 { 2558 PetscErrorCode ierr; 2559 2560 PetscFunctionBegin; 2561 PetscValidPointer(mesh,2); 2562 ierr = DMCreate(comm, mesh);CHKERRQ(ierr); 2563 ierr = DMSetType(*mesh, DMPLEX);CHKERRQ(ierr); 2564 PetscFunctionReturn(0); 2565 } 2566 2567 /* 2568 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 2569 */ 2570 /* TODO: invertCells and spaceDim arguments could be added also to to DMPlexCreateFromCellListParallel(), DMPlexBuildFromCellList_Internal() and DMPlexCreateFromCellList() */ 2571 PetscErrorCode DMPlexBuildFromCellList_Parallel_Internal(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[], PetscBool invertCells, PetscSF *sfVert) 2572 { 2573 PetscSF sfPoint; 2574 PetscLayout vLayout; 2575 PetscHSetI vhash; 2576 PetscSFNode *remoteVerticesAdj, *vertexLocal, *vertexOwner, *remoteVertex; 2577 const PetscInt *vrange; 2578 PetscInt numVerticesAdj, off = 0, *verticesAdj, numVerticesGhost = 0, *localVertex, *cone, c, p, v, g; 2579 PetscMPIInt rank, size; 2580 PetscErrorCode ierr; 2581 2582 PetscFunctionBegin; 2583 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr); 2584 ierr = MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);CHKERRQ(ierr); 2585 /* Partition vertices */ 2586 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) dm), &vLayout);CHKERRQ(ierr); 2587 ierr = PetscLayoutSetLocalSize(vLayout, numVertices);CHKERRQ(ierr); 2588 ierr = PetscLayoutSetBlockSize(vLayout, 1);CHKERRQ(ierr); 2589 ierr = PetscLayoutSetUp(vLayout);CHKERRQ(ierr); 2590 ierr = PetscLayoutGetRanges(vLayout, &vrange);CHKERRQ(ierr); 2591 /* Count vertices and map them to procs */ 2592 ierr = PetscHSetICreate(&vhash);CHKERRQ(ierr); 2593 for (c = 0; c < numCells; ++c) { 2594 for (p = 0; p < numCorners; ++p) { 2595 ierr = PetscHSetIAdd(vhash, cells[c*numCorners+p]);CHKERRQ(ierr); 2596 } 2597 } 2598 ierr = PetscHSetIGetSize(vhash, &numVerticesAdj);CHKERRQ(ierr); 2599 ierr = PetscMalloc1(numVerticesAdj, &verticesAdj);CHKERRQ(ierr); 2600 ierr = PetscHSetIGetElems(vhash, &off, verticesAdj);CHKERRQ(ierr); 2601 ierr = PetscHSetIDestroy(&vhash);CHKERRQ(ierr); 2602 if (off != numVerticesAdj) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid number of local vertices %D should be %D", off, numVerticesAdj); 2603 ierr = PetscSortInt(numVerticesAdj, verticesAdj);CHKERRQ(ierr); 2604 ierr = PetscMalloc1(numVerticesAdj, &remoteVerticesAdj);CHKERRQ(ierr); 2605 for (v = 0; v < numVerticesAdj; ++v) { 2606 const PetscInt gv = verticesAdj[v]; 2607 PetscInt vrank; 2608 2609 ierr = PetscFindInt(gv, size+1, vrange, &vrank);CHKERRQ(ierr); 2610 vrank = vrank < 0 ? -(vrank+2) : vrank; 2611 remoteVerticesAdj[v].index = gv - vrange[vrank]; 2612 remoteVerticesAdj[v].rank = vrank; 2613 } 2614 /* Create cones */ 2615 ierr = DMPlexSetChart(dm, 0, numCells+numVerticesAdj);CHKERRQ(ierr); 2616 for (c = 0; c < numCells; ++c) {ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr);} 2617 ierr = DMSetUp(dm);CHKERRQ(ierr); 2618 ierr = DMGetWorkArray(dm, numCorners, MPIU_INT, &cone);CHKERRQ(ierr); 2619 for (c = 0; c < numCells; ++c) { 2620 for (p = 0; p < numCorners; ++p) { 2621 const PetscInt gv = cells[c*numCorners+p]; 2622 PetscInt lv; 2623 2624 ierr = PetscFindInt(gv, numVerticesAdj, verticesAdj, &lv);CHKERRQ(ierr); 2625 if (lv < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not find global vertex %D in local connectivity", gv); 2626 cone[p] = lv+numCells; 2627 } 2628 if (invertCells) { ierr = DMPlexInvertCell_Internal(spaceDim, numCorners, cone);CHKERRQ(ierr); } 2629 ierr = DMPlexSetCone(dm, c, cone);CHKERRQ(ierr); 2630 } 2631 ierr = DMRestoreWorkArray(dm, numCorners, MPIU_INT, &cone);CHKERRQ(ierr); 2632 /* Create SF for vertices */ 2633 ierr = PetscSFCreate(PetscObjectComm((PetscObject)dm), sfVert);CHKERRQ(ierr); 2634 ierr = PetscObjectSetName((PetscObject) *sfVert, "Vertex Ownership SF");CHKERRQ(ierr); 2635 ierr = PetscSFSetFromOptions(*sfVert);CHKERRQ(ierr); 2636 ierr = PetscSFSetGraph(*sfVert, numVertices, numVerticesAdj, NULL, PETSC_OWN_POINTER, remoteVerticesAdj, PETSC_OWN_POINTER);CHKERRQ(ierr); 2637 ierr = PetscFree(verticesAdj);CHKERRQ(ierr); 2638 /* Build pointSF */ 2639 ierr = PetscMalloc2(numVerticesAdj, &vertexLocal, numVertices, &vertexOwner);CHKERRQ(ierr); 2640 for (v = 0; v < numVerticesAdj; ++v) {vertexLocal[v].index = v+numCells; vertexLocal[v].rank = rank;} 2641 for (v = 0; v < numVertices; ++v) {vertexOwner[v].index = -1; vertexOwner[v].rank = -1;} 2642 ierr = PetscSFReduceBegin(*sfVert, MPIU_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);CHKERRQ(ierr); 2643 ierr = PetscSFReduceEnd(*sfVert, MPIU_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);CHKERRQ(ierr); 2644 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); 2645 ierr = PetscSFBcastBegin(*sfVert, MPIU_2INT, vertexOwner, vertexLocal);CHKERRQ(ierr); 2646 ierr = PetscSFBcastEnd(*sfVert, MPIU_2INT, vertexOwner, vertexLocal);CHKERRQ(ierr); 2647 for (v = 0; v < numVerticesAdj; ++v) if (vertexLocal[v].rank != rank) ++numVerticesGhost; 2648 ierr = PetscMalloc1(numVerticesGhost, &localVertex);CHKERRQ(ierr); 2649 ierr = PetscMalloc1(numVerticesGhost, &remoteVertex);CHKERRQ(ierr); 2650 for (v = 0, g = 0; v < numVerticesAdj; ++v) { 2651 if (vertexLocal[v].rank != rank) { 2652 localVertex[g] = v+numCells; 2653 remoteVertex[g].index = vertexLocal[v].index; 2654 remoteVertex[g].rank = vertexLocal[v].rank; 2655 ++g; 2656 } 2657 } 2658 ierr = PetscFree2(vertexLocal, vertexOwner);CHKERRQ(ierr); 2659 if (g != numVerticesGhost) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of vertex ghosts %D should be %D", g, numVerticesGhost); 2660 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 2661 ierr = PetscObjectSetName((PetscObject) sfPoint, "point SF");CHKERRQ(ierr); 2662 ierr = PetscSFSetGraph(sfPoint, numCells+numVerticesAdj, numVerticesGhost, localVertex, PETSC_OWN_POINTER, remoteVertex, PETSC_OWN_POINTER);CHKERRQ(ierr); 2663 ierr = PetscLayoutDestroy(&vLayout);CHKERRQ(ierr); 2664 /* Fill in the rest of the topology structure */ 2665 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 2666 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 2667 PetscFunctionReturn(0); 2668 } 2669 2670 /* 2671 This takes as input the coordinates for each owned vertex 2672 */ 2673 PetscErrorCode DMPlexBuildCoordinates_Parallel_Internal(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numV, PetscSF sfVert, const PetscReal vertexCoords[]) 2674 { 2675 PetscSection coordSection; 2676 Vec coordinates; 2677 PetscScalar *coords; 2678 PetscInt numVertices, numVerticesAdj, coordSize, v; 2679 PetscErrorCode ierr; 2680 2681 PetscFunctionBegin; 2682 ierr = DMSetCoordinateDim(dm, spaceDim);CHKERRQ(ierr); 2683 ierr = PetscSFGetGraph(sfVert, &numVertices, &numVerticesAdj, NULL, NULL);CHKERRQ(ierr); 2684 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2685 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 2686 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 2687 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVerticesAdj);CHKERRQ(ierr); 2688 for (v = numCells; v < numCells+numVerticesAdj; ++v) { 2689 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 2690 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 2691 } 2692 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 2693 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 2694 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 2695 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 2696 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 2697 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 2698 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 2699 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2700 { 2701 MPI_Datatype coordtype; 2702 2703 /* Need a temp buffer for coords if we have complex/single */ 2704 ierr = MPI_Type_contiguous(spaceDim, MPIU_SCALAR, &coordtype);CHKERRQ(ierr); 2705 ierr = MPI_Type_commit(&coordtype);CHKERRQ(ierr); 2706 #if defined(PETSC_USE_COMPLEX) 2707 { 2708 PetscScalar *svertexCoords; 2709 PetscInt i; 2710 ierr = PetscMalloc1(numV*spaceDim,&svertexCoords);CHKERRQ(ierr); 2711 for (i=0; i<numV*spaceDim; i++) svertexCoords[i] = vertexCoords[i]; 2712 ierr = PetscSFBcastBegin(sfVert, coordtype, svertexCoords, coords);CHKERRQ(ierr); 2713 ierr = PetscSFBcastEnd(sfVert, coordtype, svertexCoords, coords);CHKERRQ(ierr); 2714 ierr = PetscFree(svertexCoords);CHKERRQ(ierr); 2715 } 2716 #else 2717 ierr = PetscSFBcastBegin(sfVert, coordtype, vertexCoords, coords);CHKERRQ(ierr); 2718 ierr = PetscSFBcastEnd(sfVert, coordtype, vertexCoords, coords);CHKERRQ(ierr); 2719 #endif 2720 ierr = MPI_Type_free(&coordtype);CHKERRQ(ierr); 2721 } 2722 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2723 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 2724 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 2725 PetscFunctionReturn(0); 2726 } 2727 2728 /*@C 2729 DMPlexCreateFromCellListParallel - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM 2730 2731 Input Parameters: 2732 + comm - The communicator 2733 . dim - The topological dimension of the mesh 2734 . numCells - The number of cells owned by this process 2735 . numVertices - The number of vertices owned by this process 2736 . numCorners - The number of vertices for each cell 2737 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 2738 . cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 2739 . spaceDim - The spatial dimension used for coordinates 2740 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 2741 2742 Output Parameter: 2743 + dm - The DM 2744 - vertexSF - Optional, SF describing complete vertex ownership 2745 2746 Note: Two triangles sharing a face 2747 $ 2748 $ 2 2749 $ / | \ 2750 $ / | \ 2751 $ / | \ 2752 $ 0 0 | 1 3 2753 $ \ | / 2754 $ \ | / 2755 $ \ | / 2756 $ 1 2757 would have input 2758 $ numCells = 2, numVertices = 4 2759 $ cells = [0 1 2 1 3 2] 2760 $ 2761 which would result in the DMPlex 2762 $ 2763 $ 4 2764 $ / | \ 2765 $ / | \ 2766 $ / | \ 2767 $ 2 0 | 1 5 2768 $ \ | / 2769 $ \ | / 2770 $ \ | / 2771 $ 3 2772 2773 Level: beginner 2774 2775 .seealso: DMPlexCreateFromCellList(), DMPlexCreateFromDAG(), DMPlexCreate() 2776 @*/ 2777 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) 2778 { 2779 PetscSF sfVert; 2780 PetscErrorCode ierr; 2781 2782 PetscFunctionBegin; 2783 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 2784 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 2785 PetscValidLogicalCollectiveInt(*dm, dim, 2); 2786 PetscValidLogicalCollectiveInt(*dm, spaceDim, 8); 2787 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 2788 ierr = DMPlexBuildFromCellList_Parallel_Internal(*dm, spaceDim, numCells, numVertices, numCorners, cells, PETSC_FALSE, &sfVert);CHKERRQ(ierr); 2789 if (interpolate) { 2790 DM idm; 2791 2792 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 2793 ierr = DMDestroy(dm);CHKERRQ(ierr); 2794 *dm = idm; 2795 } 2796 ierr = DMPlexBuildCoordinates_Parallel_Internal(*dm, spaceDim, numCells, numVertices, sfVert, vertexCoords);CHKERRQ(ierr); 2797 if (vertexSF) *vertexSF = sfVert; 2798 else {ierr = PetscSFDestroy(&sfVert);CHKERRQ(ierr);} 2799 PetscFunctionReturn(0); 2800 } 2801 2802 /* 2803 This takes as input the common mesh generator output, a list of the vertices for each cell 2804 */ 2805 PetscErrorCode DMPlexBuildFromCellList_Internal(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[], PetscBool invertCells) 2806 { 2807 PetscInt *cone, c, p; 2808 PetscErrorCode ierr; 2809 2810 PetscFunctionBegin; 2811 ierr = DMPlexSetChart(dm, 0, numCells+numVertices);CHKERRQ(ierr); 2812 for (c = 0; c < numCells; ++c) { 2813 ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr); 2814 } 2815 ierr = DMSetUp(dm);CHKERRQ(ierr); 2816 ierr = DMGetWorkArray(dm, numCorners, MPIU_INT, &cone);CHKERRQ(ierr); 2817 for (c = 0; c < numCells; ++c) { 2818 for (p = 0; p < numCorners; ++p) { 2819 cone[p] = cells[c*numCorners+p]+numCells; 2820 } 2821 if (invertCells) { ierr = DMPlexInvertCell_Internal(spaceDim, numCorners, cone);CHKERRQ(ierr); } 2822 ierr = DMPlexSetCone(dm, c, cone);CHKERRQ(ierr); 2823 } 2824 ierr = DMRestoreWorkArray(dm, numCorners, MPIU_INT, &cone);CHKERRQ(ierr); 2825 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 2826 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 2827 PetscFunctionReturn(0); 2828 } 2829 2830 /* 2831 This takes as input the coordinates for each vertex 2832 */ 2833 PetscErrorCode DMPlexBuildCoordinates_Internal(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, const double vertexCoords[]) 2834 { 2835 PetscSection coordSection; 2836 Vec coordinates; 2837 DM cdm; 2838 PetscScalar *coords; 2839 PetscInt v, d; 2840 PetscErrorCode ierr; 2841 2842 PetscFunctionBegin; 2843 ierr = DMSetCoordinateDim(dm, spaceDim);CHKERRQ(ierr); 2844 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2845 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 2846 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 2847 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVertices);CHKERRQ(ierr); 2848 for (v = numCells; v < numCells+numVertices; ++v) { 2849 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 2850 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 2851 } 2852 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 2853 2854 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 2855 ierr = DMCreateLocalVector(cdm, &coordinates);CHKERRQ(ierr); 2856 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 2857 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 2858 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2859 for (v = 0; v < numVertices; ++v) { 2860 for (d = 0; d < spaceDim; ++d) { 2861 coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d]; 2862 } 2863 } 2864 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2865 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 2866 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 2867 PetscFunctionReturn(0); 2868 } 2869 2870 /*@C 2871 DMPlexCreateFromCellList - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM 2872 2873 Input Parameters: 2874 + comm - The communicator 2875 . dim - The topological dimension of the mesh 2876 . numCells - The number of cells 2877 . numVertices - The number of vertices 2878 . numCorners - The number of vertices for each cell 2879 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 2880 . cells - An array of numCells*numCorners numbers, the vertices for each cell 2881 . spaceDim - The spatial dimension used for coordinates 2882 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 2883 2884 Output Parameter: 2885 . dm - The DM 2886 2887 Note: Two triangles sharing a face 2888 $ 2889 $ 2 2890 $ / | \ 2891 $ / | \ 2892 $ / | \ 2893 $ 0 0 | 1 3 2894 $ \ | / 2895 $ \ | / 2896 $ \ | / 2897 $ 1 2898 would have input 2899 $ numCells = 2, numVertices = 4 2900 $ cells = [0 1 2 1 3 2] 2901 $ 2902 which would result in the DMPlex 2903 $ 2904 $ 4 2905 $ / | \ 2906 $ / | \ 2907 $ / | \ 2908 $ 2 0 | 1 5 2909 $ \ | / 2910 $ \ | / 2911 $ \ | / 2912 $ 3 2913 2914 Level: beginner 2915 2916 .seealso: DMPlexCreateFromDAG(), DMPlexCreate() 2917 @*/ 2918 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) 2919 { 2920 PetscErrorCode ierr; 2921 2922 PetscFunctionBegin; 2923 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 2924 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 2925 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 2926 ierr = DMPlexBuildFromCellList_Internal(*dm, spaceDim, numCells, numVertices, numCorners, cells, PETSC_FALSE);CHKERRQ(ierr); 2927 if (interpolate) { 2928 DM idm; 2929 2930 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 2931 ierr = DMDestroy(dm);CHKERRQ(ierr); 2932 *dm = idm; 2933 } 2934 ierr = DMPlexBuildCoordinates_Internal(*dm, spaceDim, numCells, numVertices, vertexCoords);CHKERRQ(ierr); 2935 PetscFunctionReturn(0); 2936 } 2937 2938 /*@ 2939 DMPlexCreateFromDAG - This takes as input the adjacency-list representation of the Directed Acyclic Graph (Hasse Diagram) encoding a mesh, and produces a DM 2940 2941 Input Parameters: 2942 + dm - The empty DM object, usually from DMCreate() and DMSetDimension() 2943 . depth - The depth of the DAG 2944 . numPoints - Array of size depth + 1 containing the number of points at each depth 2945 . coneSize - The cone size of each point 2946 . cones - The concatenation of the cone points for each point, the cone list must be oriented correctly for each point 2947 . coneOrientations - The orientation of each cone point 2948 - vertexCoords - An array of numPoints[0]*spacedim numbers representing the coordinates of each vertex, with spacedim the value set via DMSetCoordinateDim() 2949 2950 Output Parameter: 2951 . dm - The DM 2952 2953 Note: Two triangles sharing a face would have input 2954 $ depth = 1, numPoints = [4 2], coneSize = [3 3 0 0 0 0] 2955 $ cones = [2 3 4 3 5 4], coneOrientations = [0 0 0 0 0 0] 2956 $ vertexCoords = [-1.0 0.0 0.0 -1.0 0.0 1.0 1.0 0.0] 2957 $ 2958 which would result in the DMPlex 2959 $ 2960 $ 4 2961 $ / | \ 2962 $ / | \ 2963 $ / | \ 2964 $ 2 0 | 1 5 2965 $ \ | / 2966 $ \ | / 2967 $ \ | / 2968 $ 3 2969 $ 2970 $ Notice that all points are numbered consecutively, unlikely DMPlexCreateFromCellList() 2971 2972 Level: advanced 2973 2974 .seealso: DMPlexCreateFromCellList(), DMPlexCreate() 2975 @*/ 2976 PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[]) 2977 { 2978 Vec coordinates; 2979 PetscSection coordSection; 2980 PetscScalar *coords; 2981 PetscInt coordSize, firstVertex = -1, pStart = 0, pEnd = 0, p, v, dim, dimEmbed, d, off; 2982 PetscErrorCode ierr; 2983 2984 PetscFunctionBegin; 2985 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2986 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 2987 if (dimEmbed < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Embedding dimension %d cannot be less than intrinsic dimension %d",dimEmbed,dim); 2988 for (d = 0; d <= depth; ++d) pEnd += numPoints[d]; 2989 ierr = DMPlexSetChart(dm, pStart, pEnd);CHKERRQ(ierr); 2990 for (p = pStart; p < pEnd; ++p) { 2991 ierr = DMPlexSetConeSize(dm, p, coneSize[p-pStart]);CHKERRQ(ierr); 2992 if (firstVertex < 0 && !coneSize[p - pStart]) { 2993 firstVertex = p - pStart; 2994 } 2995 } 2996 if (firstVertex < 0 && numPoints[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected %d vertices but could not find any", numPoints[0]); 2997 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 2998 for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) { 2999 ierr = DMPlexSetCone(dm, p, &cones[off]);CHKERRQ(ierr); 3000 ierr = DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);CHKERRQ(ierr); 3001 } 3002 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 3003 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 3004 /* Build coordinates */ 3005 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 3006 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 3007 ierr = PetscSectionSetFieldComponents(coordSection, 0, dimEmbed);CHKERRQ(ierr); 3008 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0]);CHKERRQ(ierr); 3009 for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) { 3010 ierr = PetscSectionSetDof(coordSection, v, dimEmbed);CHKERRQ(ierr); 3011 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed);CHKERRQ(ierr); 3012 } 3013 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 3014 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 3015 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 3016 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 3017 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 3018 ierr = VecSetBlockSize(coordinates, dimEmbed);CHKERRQ(ierr); 3019 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 3020 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 3021 for (v = 0; v < numPoints[0]; ++v) { 3022 PetscInt off; 3023 3024 ierr = PetscSectionGetOffset(coordSection, v+firstVertex, &off);CHKERRQ(ierr); 3025 for (d = 0; d < dimEmbed; ++d) { 3026 coords[off+d] = vertexCoords[v*dimEmbed+d]; 3027 } 3028 } 3029 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 3030 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 3031 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 3032 PetscFunctionReturn(0); 3033 } 3034 3035 /*@C 3036 DMPlexCreateCellVertexFromFile - Create a DMPlex mesh from a simple cell-vertex file. 3037 3038 + comm - The MPI communicator 3039 . filename - Name of the .dat file 3040 - interpolate - Create faces and edges in the mesh 3041 3042 Output Parameter: 3043 . dm - The DM object representing the mesh 3044 3045 Note: The format is the simplest possible: 3046 $ Ne 3047 $ v0 v1 ... vk 3048 $ Nv 3049 $ x y z marker 3050 3051 Level: beginner 3052 3053 .seealso: DMPlexCreateFromFile(), DMPlexCreateMedFromFile(), DMPlexCreateGmsh(), DMPlexCreate() 3054 @*/ 3055 PetscErrorCode DMPlexCreateCellVertexFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 3056 { 3057 DMLabel marker; 3058 PetscViewer viewer; 3059 Vec coordinates; 3060 PetscSection coordSection; 3061 PetscScalar *coords; 3062 char line[PETSC_MAX_PATH_LEN]; 3063 PetscInt dim = 3, cdim = 3, coordSize, v, c, d; 3064 PetscMPIInt rank; 3065 int snum, Nv, Nc; 3066 PetscErrorCode ierr; 3067 3068 PetscFunctionBegin; 3069 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 3070 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 3071 ierr = PetscViewerSetType(viewer, PETSCVIEWERASCII);CHKERRQ(ierr); 3072 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 3073 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 3074 if (!rank) { 3075 ierr = PetscViewerRead(viewer, line, 2, NULL, PETSC_STRING);CHKERRQ(ierr); 3076 snum = sscanf(line, "%d %d", &Nc, &Nv); 3077 if (snum != 2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 3078 } else { 3079 Nc = Nv = 0; 3080 } 3081 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 3082 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 3083 ierr = DMPlexSetChart(*dm, 0, Nc+Nv);CHKERRQ(ierr); 3084 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 3085 ierr = DMSetCoordinateDim(*dm, cdim);CHKERRQ(ierr); 3086 /* Read topology */ 3087 if (!rank) { 3088 PetscInt cone[8], corners = 8; 3089 int vbuf[8], v; 3090 3091 for (c = 0; c < Nc; ++c) {ierr = DMPlexSetConeSize(*dm, c, corners);CHKERRQ(ierr);} 3092 ierr = DMSetUp(*dm);CHKERRQ(ierr); 3093 for (c = 0; c < Nc; ++c) { 3094 ierr = PetscViewerRead(viewer, line, corners, NULL, PETSC_STRING);CHKERRQ(ierr); 3095 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]); 3096 if (snum != corners) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 3097 for (v = 0; v < corners; ++v) cone[v] = vbuf[v] + Nc; 3098 /* Hexahedra are inverted */ 3099 { 3100 PetscInt tmp = cone[1]; 3101 cone[1] = cone[3]; 3102 cone[3] = tmp; 3103 } 3104 ierr = DMPlexSetCone(*dm, c, cone);CHKERRQ(ierr); 3105 } 3106 } 3107 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 3108 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 3109 /* Read coordinates */ 3110 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 3111 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 3112 ierr = PetscSectionSetFieldComponents(coordSection, 0, cdim);CHKERRQ(ierr); 3113 ierr = PetscSectionSetChart(coordSection, Nc, Nc + Nv);CHKERRQ(ierr); 3114 for (v = Nc; v < Nc+Nv; ++v) { 3115 ierr = PetscSectionSetDof(coordSection, v, cdim);CHKERRQ(ierr); 3116 ierr = PetscSectionSetFieldDof(coordSection, v, 0, cdim);CHKERRQ(ierr); 3117 } 3118 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 3119 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 3120 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 3121 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 3122 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 3123 ierr = VecSetBlockSize(coordinates, cdim);CHKERRQ(ierr); 3124 ierr = VecSetType(coordinates, VECSTANDARD);CHKERRQ(ierr); 3125 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 3126 if (!rank) { 3127 double x[3]; 3128 int val; 3129 3130 ierr = DMCreateLabel(*dm, "marker");CHKERRQ(ierr); 3131 ierr = DMGetLabel(*dm, "marker", &marker);CHKERRQ(ierr); 3132 for (v = 0; v < Nv; ++v) { 3133 ierr = PetscViewerRead(viewer, line, 4, NULL, PETSC_STRING);CHKERRQ(ierr); 3134 snum = sscanf(line, "%lg %lg %lg %d", &x[0], &x[1], &x[2], &val); 3135 if (snum != 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 3136 for (d = 0; d < cdim; ++d) coords[v*cdim+d] = x[d]; 3137 if (val) {ierr = DMLabelSetValue(marker, v+Nc, val);CHKERRQ(ierr);} 3138 } 3139 } 3140 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 3141 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 3142 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 3143 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3144 if (interpolate) { 3145 DM idm; 3146 DMLabel bdlabel; 3147 3148 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 3149 ierr = DMDestroy(dm);CHKERRQ(ierr); 3150 *dm = idm; 3151 3152 ierr = DMGetLabel(*dm, "marker", &bdlabel);CHKERRQ(ierr); 3153 ierr = DMPlexMarkBoundaryFaces(*dm, PETSC_DETERMINE, bdlabel);CHKERRQ(ierr); 3154 ierr = DMPlexLabelComplete(*dm, bdlabel);CHKERRQ(ierr); 3155 } 3156 PetscFunctionReturn(0); 3157 } 3158 3159 /*@C 3160 DMPlexCreateFromFile - This takes a filename and produces a DM 3161 3162 Input Parameters: 3163 + comm - The communicator 3164 . filename - A file name 3165 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 3166 3167 Output Parameter: 3168 . dm - The DM 3169 3170 Options Database Keys: 3171 . -dm_plex_create_from_hdf5_xdmf - use the PETSC_VIEWER_HDF5_XDMF format for reading HDF5 3172 3173 Level: beginner 3174 3175 .seealso: DMPlexCreateFromDAG(), DMPlexCreateFromCellList(), DMPlexCreate() 3176 @*/ 3177 PetscErrorCode DMPlexCreateFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 3178 { 3179 const char *extGmsh = ".msh"; 3180 const char *extCGNS = ".cgns"; 3181 const char *extExodus = ".exo"; 3182 const char *extGenesis = ".gen"; 3183 const char *extFluent = ".cas"; 3184 const char *extHDF5 = ".h5"; 3185 const char *extMed = ".med"; 3186 const char *extPLY = ".ply"; 3187 const char *extCV = ".dat"; 3188 size_t len; 3189 PetscBool isGmsh, isCGNS, isExodus, isGenesis, isFluent, isHDF5, isMed, isPLY, isCV; 3190 PetscMPIInt rank; 3191 PetscErrorCode ierr; 3192 3193 PetscFunctionBegin; 3194 PetscValidPointer(filename, 2); 3195 PetscValidPointer(dm, 4); 3196 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 3197 ierr = PetscStrlen(filename, &len);CHKERRQ(ierr); 3198 if (!len) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Filename must be a valid path"); 3199 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGmsh, 4, &isGmsh);CHKERRQ(ierr); 3200 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extCGNS, 5, &isCGNS);CHKERRQ(ierr); 3201 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extExodus, 4, &isExodus);CHKERRQ(ierr); 3202 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGenesis, 4, &isGenesis);CHKERRQ(ierr); 3203 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extFluent, 4, &isFluent);CHKERRQ(ierr); 3204 ierr = PetscStrncmp(&filename[PetscMax(0,len-3)], extHDF5, 3, &isHDF5);CHKERRQ(ierr); 3205 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extMed, 4, &isMed);CHKERRQ(ierr); 3206 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extPLY, 4, &isPLY);CHKERRQ(ierr); 3207 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extCV, 4, &isCV);CHKERRQ(ierr); 3208 if (isGmsh) { 3209 ierr = DMPlexCreateGmshFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3210 } else if (isCGNS) { 3211 ierr = DMPlexCreateCGNSFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3212 } else if (isExodus || isGenesis) { 3213 ierr = DMPlexCreateExodusFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3214 } else if (isFluent) { 3215 ierr = DMPlexCreateFluentFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3216 } else if (isHDF5) { 3217 PetscBool load_hdf5_xdmf = PETSC_FALSE; 3218 PetscViewer viewer; 3219 3220 ierr = PetscOptionsGetBool(NULL, NULL, "-dm_plex_create_from_hdf5_xdmf", &load_hdf5_xdmf, NULL);CHKERRQ(ierr); 3221 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 3222 ierr = PetscViewerSetType(viewer, PETSCVIEWERHDF5);CHKERRQ(ierr); 3223 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 3224 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 3225 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 3226 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 3227 if (load_hdf5_xdmf) {ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_XDMF);CHKERRQ(ierr);} 3228 ierr = DMLoad(*dm, viewer);CHKERRQ(ierr); 3229 if (load_hdf5_xdmf) {ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);} 3230 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3231 3232 if (interpolate) { 3233 DM idm; 3234 3235 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 3236 ierr = DMDestroy(dm);CHKERRQ(ierr); 3237 *dm = idm; 3238 } 3239 } else if (isMed) { 3240 ierr = DMPlexCreateMedFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3241 } else if (isPLY) { 3242 ierr = DMPlexCreatePLYFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3243 } else if (isCV) { 3244 ierr = DMPlexCreateCellVertexFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3245 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot load file %s: unrecognized extension", filename); 3246 PetscFunctionReturn(0); 3247 } 3248 3249 /*@ 3250 DMPlexCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell 3251 3252 Collective on comm 3253 3254 Input Parameters: 3255 + comm - The communicator 3256 . dim - The spatial dimension 3257 - simplex - Flag for simplex, otherwise use a tensor-product cell 3258 3259 Output Parameter: 3260 . refdm - The reference cell 3261 3262 Level: intermediate 3263 3264 .keywords: reference cell 3265 .seealso: 3266 @*/ 3267 PetscErrorCode DMPlexCreateReferenceCell(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *refdm) 3268 { 3269 DM rdm; 3270 Vec coords; 3271 PetscErrorCode ierr; 3272 3273 PetscFunctionBeginUser; 3274 ierr = DMCreate(comm, &rdm);CHKERRQ(ierr); 3275 ierr = DMSetType(rdm, DMPLEX);CHKERRQ(ierr); 3276 ierr = DMSetDimension(rdm, dim);CHKERRQ(ierr); 3277 switch (dim) { 3278 case 0: 3279 { 3280 PetscInt numPoints[1] = {1}; 3281 PetscInt coneSize[1] = {0}; 3282 PetscInt cones[1] = {0}; 3283 PetscInt coneOrientations[1] = {0}; 3284 PetscScalar vertexCoords[1] = {0.0}; 3285 3286 ierr = DMPlexCreateFromDAG(rdm, 0, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3287 } 3288 break; 3289 case 1: 3290 { 3291 PetscInt numPoints[2] = {2, 1}; 3292 PetscInt coneSize[3] = {2, 0, 0}; 3293 PetscInt cones[2] = {1, 2}; 3294 PetscInt coneOrientations[2] = {0, 0}; 3295 PetscScalar vertexCoords[2] = {-1.0, 1.0}; 3296 3297 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3298 } 3299 break; 3300 case 2: 3301 if (simplex) { 3302 PetscInt numPoints[2] = {3, 1}; 3303 PetscInt coneSize[4] = {3, 0, 0, 0}; 3304 PetscInt cones[3] = {1, 2, 3}; 3305 PetscInt coneOrientations[3] = {0, 0, 0}; 3306 PetscScalar vertexCoords[6] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0}; 3307 3308 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3309 } else { 3310 PetscInt numPoints[2] = {4, 1}; 3311 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 3312 PetscInt cones[4] = {1, 2, 3, 4}; 3313 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 3314 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0}; 3315 3316 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3317 } 3318 break; 3319 case 3: 3320 if (simplex) { 3321 PetscInt numPoints[2] = {4, 1}; 3322 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 3323 PetscInt cones[4] = {1, 3, 2, 4}; 3324 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 3325 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}; 3326 3327 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3328 } else { 3329 PetscInt numPoints[2] = {8, 1}; 3330 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0}; 3331 PetscInt cones[8] = {1, 4, 3, 2, 5, 6, 7, 8}; 3332 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 3333 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, 3334 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0}; 3335 3336 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3337 } 3338 break; 3339 default: 3340 SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Cannot create reference cell for dimension %d", dim); 3341 } 3342 ierr = DMPlexInterpolate(rdm, refdm);CHKERRQ(ierr); 3343 if (rdm->coordinateDM) { 3344 DM ncdm; 3345 PetscSection cs; 3346 PetscInt pEnd = -1; 3347 3348 ierr = DMGetSection(rdm->coordinateDM, &cs);CHKERRQ(ierr); 3349 if (cs) {ierr = PetscSectionGetChart(cs, NULL, &pEnd);CHKERRQ(ierr);} 3350 if (pEnd >= 0) { 3351 ierr = DMClone(rdm->coordinateDM, &ncdm);CHKERRQ(ierr); 3352 ierr = DMSetSection(ncdm, cs);CHKERRQ(ierr); 3353 ierr = DMSetCoordinateDM(*refdm, ncdm);CHKERRQ(ierr); 3354 ierr = DMDestroy(&ncdm);CHKERRQ(ierr); 3355 } 3356 } 3357 ierr = DMGetCoordinatesLocal(rdm, &coords);CHKERRQ(ierr); 3358 if (coords) { 3359 ierr = DMSetCoordinatesLocal(*refdm, coords);CHKERRQ(ierr); 3360 } else { 3361 ierr = DMGetCoordinates(rdm, &coords);CHKERRQ(ierr); 3362 if (coords) {ierr = DMSetCoordinates(*refdm, coords);CHKERRQ(ierr);} 3363 } 3364 ierr = DMDestroy(&rdm);CHKERRQ(ierr); 3365 PetscFunctionReturn(0); 3366 } 3367