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