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