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