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 DMCreateMassMatrix_Plex(DM dmCoarse, DM dmFine, Mat *mat); 1853 extern PetscErrorCode DMCreateDefaultSection_Plex(DM dm); 1854 extern PetscErrorCode DMCreateDefaultConstraints_Plex(DM dm); 1855 extern PetscErrorCode DMCreateMatrix_Plex(DM dm, Mat *J); 1856 extern PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm); 1857 PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm); 1858 extern PetscErrorCode DMSetUp_Plex(DM dm); 1859 extern PetscErrorCode DMDestroy_Plex(DM dm); 1860 extern PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer); 1861 extern PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer); 1862 extern PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm); 1863 static PetscErrorCode DMInitialize_Plex(DM dm); 1864 1865 /* Replace dm with the contents of dmNew 1866 - Share the DM_Plex structure 1867 - Share the coordinates 1868 - Share the SF 1869 */ 1870 static PetscErrorCode DMPlexReplace_Static(DM dm, DM dmNew) 1871 { 1872 PetscSF sf; 1873 DM coordDM, coarseDM; 1874 Vec coords; 1875 PetscBool isper; 1876 const PetscReal *maxCell, *L; 1877 const DMBoundaryType *bd; 1878 PetscErrorCode ierr; 1879 1880 PetscFunctionBegin; 1881 ierr = DMGetPointSF(dmNew, &sf);CHKERRQ(ierr); 1882 ierr = DMSetPointSF(dm, sf);CHKERRQ(ierr); 1883 ierr = DMGetCoordinateDM(dmNew, &coordDM);CHKERRQ(ierr); 1884 ierr = DMGetCoordinatesLocal(dmNew, &coords);CHKERRQ(ierr); 1885 ierr = DMSetCoordinateDM(dm, coordDM);CHKERRQ(ierr); 1886 ierr = DMSetCoordinatesLocal(dm, coords);CHKERRQ(ierr); 1887 ierr = DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);CHKERRQ(ierr); 1888 ierr = DMSetPeriodicity(dmNew, isper, maxCell, L, bd);CHKERRQ(ierr); 1889 ierr = DMDestroy_Plex(dm);CHKERRQ(ierr); 1890 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 1891 dm->data = dmNew->data; 1892 ((DM_Plex *) dmNew->data)->refct++; 1893 dmNew->labels->refct++; 1894 if (!--(dm->labels->refct)) { 1895 DMLabelLink next = dm->labels->next; 1896 1897 /* destroy the labels */ 1898 while (next) { 1899 DMLabelLink tmp = next->next; 1900 1901 ierr = DMLabelDestroy(&next->label);CHKERRQ(ierr); 1902 ierr = PetscFree(next);CHKERRQ(ierr); 1903 next = tmp; 1904 } 1905 ierr = PetscFree(dm->labels);CHKERRQ(ierr); 1906 } 1907 dm->labels = dmNew->labels; 1908 dm->depthLabel = dmNew->depthLabel; 1909 ierr = DMGetCoarseDM(dmNew,&coarseDM);CHKERRQ(ierr); 1910 ierr = DMSetCoarseDM(dm,coarseDM);CHKERRQ(ierr); 1911 PetscFunctionReturn(0); 1912 } 1913 1914 /* Swap dm with the contents of dmNew 1915 - Swap the DM_Plex structure 1916 - Swap the coordinates 1917 - Swap the point PetscSF 1918 */ 1919 static PetscErrorCode DMPlexSwap_Static(DM dmA, DM dmB) 1920 { 1921 DM coordDMA, coordDMB; 1922 Vec coordsA, coordsB; 1923 PetscSF sfA, sfB; 1924 void *tmp; 1925 DMLabelLinkList listTmp; 1926 DMLabel depthTmp; 1927 PetscErrorCode ierr; 1928 1929 PetscFunctionBegin; 1930 ierr = DMGetPointSF(dmA, &sfA);CHKERRQ(ierr); 1931 ierr = DMGetPointSF(dmB, &sfB);CHKERRQ(ierr); 1932 ierr = PetscObjectReference((PetscObject) sfA);CHKERRQ(ierr); 1933 ierr = DMSetPointSF(dmA, sfB);CHKERRQ(ierr); 1934 ierr = DMSetPointSF(dmB, sfA);CHKERRQ(ierr); 1935 ierr = PetscObjectDereference((PetscObject) sfA);CHKERRQ(ierr); 1936 1937 ierr = DMGetCoordinateDM(dmA, &coordDMA);CHKERRQ(ierr); 1938 ierr = DMGetCoordinateDM(dmB, &coordDMB);CHKERRQ(ierr); 1939 ierr = PetscObjectReference((PetscObject) coordDMA);CHKERRQ(ierr); 1940 ierr = DMSetCoordinateDM(dmA, coordDMB);CHKERRQ(ierr); 1941 ierr = DMSetCoordinateDM(dmB, coordDMA);CHKERRQ(ierr); 1942 ierr = PetscObjectDereference((PetscObject) coordDMA);CHKERRQ(ierr); 1943 1944 ierr = DMGetCoordinatesLocal(dmA, &coordsA);CHKERRQ(ierr); 1945 ierr = DMGetCoordinatesLocal(dmB, &coordsB);CHKERRQ(ierr); 1946 ierr = PetscObjectReference((PetscObject) coordsA);CHKERRQ(ierr); 1947 ierr = DMSetCoordinatesLocal(dmA, coordsB);CHKERRQ(ierr); 1948 ierr = DMSetCoordinatesLocal(dmB, coordsA);CHKERRQ(ierr); 1949 ierr = PetscObjectDereference((PetscObject) coordsA);CHKERRQ(ierr); 1950 1951 tmp = dmA->data; 1952 dmA->data = dmB->data; 1953 dmB->data = tmp; 1954 listTmp = dmA->labels; 1955 dmA->labels = dmB->labels; 1956 dmB->labels = listTmp; 1957 depthTmp = dmA->depthLabel; 1958 dmA->depthLabel = dmB->depthLabel; 1959 dmB->depthLabel = depthTmp; 1960 PetscFunctionReturn(0); 1961 } 1962 1963 PetscErrorCode DMSetFromOptions_NonRefinement_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 1964 { 1965 DM_Plex *mesh = (DM_Plex*) dm->data; 1966 PetscErrorCode ierr; 1967 1968 PetscFunctionBegin; 1969 /* Handle viewing */ 1970 ierr = PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMView", PETSC_FALSE, &mesh->printSetValues, NULL);CHKERRQ(ierr); 1971 ierr = PetscOptionsInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMView", 0, &mesh->printFEM, NULL);CHKERRQ(ierr); 1972 ierr = PetscOptionsReal("-dm_plex_print_tol", "Tolerance for FEM output", "DMView", mesh->printTol, &mesh->printTol, NULL);CHKERRQ(ierr); 1973 /* Point Location */ 1974 ierr = PetscOptionsBool("-dm_plex_hash_location", "Use grid hashing for point location", "DMView", PETSC_FALSE, &mesh->useHashLocation, NULL);CHKERRQ(ierr); 1975 /* Generation and remeshing */ 1976 ierr = PetscOptionsBool("-dm_plex_remesh_bd", "Allow changes to the boundary on remeshing", "DMView", PETSC_FALSE, &mesh->remeshBd, NULL);CHKERRQ(ierr); 1977 /* Projection behavior */ 1978 ierr = PetscOptionsInt("-dm_plex_max_projection_height", "Maxmimum mesh point height used to project locally", "DMPlexSetMaxProjectionHeight", 0, &mesh->maxProjectionHeight, NULL);CHKERRQ(ierr); 1979 ierr = PetscOptionsBool("-dm_plex_regular_refinement", "Use special nested projection algorithm for regular refinement", "DMPlexSetRegularRefinement", mesh->regularRefinement, &mesh->regularRefinement, NULL);CHKERRQ(ierr); 1980 1981 ierr = PetscPartitionerSetFromOptions(mesh->partitioner);CHKERRQ(ierr); 1982 PetscFunctionReturn(0); 1983 } 1984 1985 static PetscErrorCode DMSetFromOptions_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 1986 { 1987 PetscInt refine = 0, coarsen = 0, r; 1988 PetscBool isHierarchy; 1989 PetscErrorCode ierr; 1990 1991 PetscFunctionBegin; 1992 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1993 ierr = PetscOptionsHead(PetscOptionsObject,"DMPlex Options");CHKERRQ(ierr); 1994 /* Handle DMPlex refinement */ 1995 ierr = PetscOptionsInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL);CHKERRQ(ierr); 1996 ierr = PetscOptionsInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy);CHKERRQ(ierr); 1997 if (refine) {ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE);CHKERRQ(ierr);} 1998 if (refine && isHierarchy) { 1999 DM *dms, coarseDM; 2000 2001 ierr = DMGetCoarseDM(dm, &coarseDM);CHKERRQ(ierr); 2002 ierr = PetscObjectReference((PetscObject)coarseDM);CHKERRQ(ierr); 2003 ierr = PetscMalloc1(refine,&dms);CHKERRQ(ierr); 2004 ierr = DMRefineHierarchy(dm, refine, dms);CHKERRQ(ierr); 2005 /* Total hack since we do not pass in a pointer */ 2006 ierr = DMPlexSwap_Static(dm, dms[refine-1]);CHKERRQ(ierr); 2007 if (refine == 1) { 2008 ierr = DMSetCoarseDM(dm, dms[0]);CHKERRQ(ierr); 2009 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr); 2010 } else { 2011 ierr = DMSetCoarseDM(dm, dms[refine-2]);CHKERRQ(ierr); 2012 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr); 2013 ierr = DMSetCoarseDM(dms[0], dms[refine-1]);CHKERRQ(ierr); 2014 ierr = DMPlexSetRegularRefinement(dms[0], PETSC_TRUE);CHKERRQ(ierr); 2015 } 2016 ierr = DMSetCoarseDM(dms[refine-1], coarseDM);CHKERRQ(ierr); 2017 ierr = PetscObjectDereference((PetscObject)coarseDM);CHKERRQ(ierr); 2018 /* Free DMs */ 2019 for (r = 0; r < refine; ++r) { 2020 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr); 2021 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 2022 } 2023 ierr = PetscFree(dms);CHKERRQ(ierr); 2024 } else { 2025 for (r = 0; r < refine; ++r) { 2026 DM refinedMesh; 2027 2028 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2029 ierr = DMRefine(dm, PetscObjectComm((PetscObject) dm), &refinedMesh);CHKERRQ(ierr); 2030 /* Total hack since we do not pass in a pointer */ 2031 ierr = DMPlexReplace_Static(dm, refinedMesh);CHKERRQ(ierr); 2032 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2033 ierr = DMDestroy(&refinedMesh);CHKERRQ(ierr); 2034 } 2035 } 2036 /* Handle DMPlex coarsening */ 2037 ierr = PetscOptionsInt("-dm_coarsen", "Coarsen the mesh", "DMCreate", coarsen, &coarsen, NULL);CHKERRQ(ierr); 2038 ierr = PetscOptionsInt("-dm_coarsen_hierarchy", "The number of coarsenings", "DMCreate", coarsen, &coarsen, &isHierarchy);CHKERRQ(ierr); 2039 if (coarsen && isHierarchy) { 2040 DM *dms; 2041 2042 ierr = PetscMalloc1(coarsen, &dms);CHKERRQ(ierr); 2043 ierr = DMCoarsenHierarchy(dm, coarsen, dms);CHKERRQ(ierr); 2044 /* Free DMs */ 2045 for (r = 0; r < coarsen; ++r) { 2046 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr); 2047 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 2048 } 2049 ierr = PetscFree(dms);CHKERRQ(ierr); 2050 } else { 2051 for (r = 0; r < coarsen; ++r) { 2052 DM coarseMesh; 2053 2054 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2055 ierr = DMCoarsen(dm, PetscObjectComm((PetscObject) dm), &coarseMesh);CHKERRQ(ierr); 2056 /* Total hack since we do not pass in a pointer */ 2057 ierr = DMPlexReplace_Static(dm, coarseMesh);CHKERRQ(ierr); 2058 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2059 ierr = DMDestroy(&coarseMesh);CHKERRQ(ierr); 2060 } 2061 } 2062 /* Handle */ 2063 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2064 ierr = PetscOptionsTail();CHKERRQ(ierr); 2065 PetscFunctionReturn(0); 2066 } 2067 2068 static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec) 2069 { 2070 PetscErrorCode ierr; 2071 2072 PetscFunctionBegin; 2073 ierr = DMCreateGlobalVector_Section_Private(dm,vec);CHKERRQ(ierr); 2074 /* ierr = VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM);CHKERRQ(ierr); */ 2075 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex);CHKERRQ(ierr); 2076 ierr = VecSetOperation(*vec, VECOP_VIEWNATIVE, (void (*)(void)) VecView_Plex_Native);CHKERRQ(ierr); 2077 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex);CHKERRQ(ierr); 2078 ierr = VecSetOperation(*vec, VECOP_LOADNATIVE, (void (*)(void)) VecLoad_Plex_Native);CHKERRQ(ierr); 2079 PetscFunctionReturn(0); 2080 } 2081 2082 static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec) 2083 { 2084 PetscErrorCode ierr; 2085 2086 PetscFunctionBegin; 2087 ierr = DMCreateLocalVector_Section_Private(dm,vec);CHKERRQ(ierr); 2088 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex_Local);CHKERRQ(ierr); 2089 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex_Local);CHKERRQ(ierr); 2090 PetscFunctionReturn(0); 2091 } 2092 2093 static PetscErrorCode DMGetDimPoints_Plex(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd) 2094 { 2095 PetscInt depth, d; 2096 PetscErrorCode ierr; 2097 2098 PetscFunctionBegin; 2099 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 2100 if (depth == 1) { 2101 ierr = DMGetDimension(dm, &d);CHKERRQ(ierr); 2102 if (dim == 0) {ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr);} 2103 else if (dim == d) {ierr = DMPlexGetDepthStratum(dm, 1, pStart, pEnd);CHKERRQ(ierr);} 2104 else {*pStart = 0; *pEnd = 0;} 2105 } else { 2106 ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr); 2107 } 2108 PetscFunctionReturn(0); 2109 } 2110 2111 static PetscErrorCode DMGetNeighors_Plex(DM dm, PetscInt *nranks, const PetscMPIInt *ranks[]) 2112 { 2113 PetscSF sf; 2114 PetscErrorCode ierr; 2115 2116 PetscFunctionBegin; 2117 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 2118 ierr = PetscSFGetRanks(sf, nranks, ranks, NULL, NULL, NULL);CHKERRQ(ierr); 2119 PetscFunctionReturn(0); 2120 } 2121 2122 static PetscErrorCode DMInitialize_Plex(DM dm) 2123 { 2124 PetscErrorCode ierr; 2125 2126 PetscFunctionBegin; 2127 dm->ops->view = DMView_Plex; 2128 dm->ops->load = DMLoad_Plex; 2129 dm->ops->setfromoptions = DMSetFromOptions_Plex; 2130 dm->ops->clone = DMClone_Plex; 2131 dm->ops->setup = DMSetUp_Plex; 2132 dm->ops->createdefaultsection = DMCreateDefaultSection_Plex; 2133 dm->ops->createdefaultconstraints = DMCreateDefaultConstraints_Plex; 2134 dm->ops->createglobalvector = DMCreateGlobalVector_Plex; 2135 dm->ops->createlocalvector = DMCreateLocalVector_Plex; 2136 dm->ops->getlocaltoglobalmapping = NULL; 2137 dm->ops->createfieldis = NULL; 2138 dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex; 2139 dm->ops->getcoloring = NULL; 2140 dm->ops->creatematrix = DMCreateMatrix_Plex; 2141 dm->ops->createinterpolation = DMCreateInterpolation_Plex; 2142 dm->ops->createmassmatrix = DMCreateMassMatrix_Plex; 2143 dm->ops->getaggregates = NULL; 2144 dm->ops->getinjection = DMCreateInjection_Plex; 2145 dm->ops->refine = DMRefine_Plex; 2146 dm->ops->coarsen = DMCoarsen_Plex; 2147 dm->ops->refinehierarchy = DMRefineHierarchy_Plex; 2148 dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Plex; 2149 dm->ops->adaptlabel = DMAdaptLabel_Plex; 2150 dm->ops->adaptmetric = DMAdaptMetric_Plex; 2151 dm->ops->globaltolocalbegin = NULL; 2152 dm->ops->globaltolocalend = NULL; 2153 dm->ops->localtoglobalbegin = NULL; 2154 dm->ops->localtoglobalend = NULL; 2155 dm->ops->destroy = DMDestroy_Plex; 2156 dm->ops->createsubdm = DMCreateSubDM_Plex; 2157 dm->ops->getdimpoints = DMGetDimPoints_Plex; 2158 dm->ops->locatepoints = DMLocatePoints_Plex; 2159 dm->ops->projectfunctionlocal = DMProjectFunctionLocal_Plex; 2160 dm->ops->projectfunctionlabellocal = DMProjectFunctionLabelLocal_Plex; 2161 dm->ops->projectfieldlocal = DMProjectFieldLocal_Plex; 2162 dm->ops->projectfieldlabellocal = DMProjectFieldLabelLocal_Plex; 2163 dm->ops->computel2diff = DMComputeL2Diff_Plex; 2164 dm->ops->computel2gradientdiff = DMComputeL2GradientDiff_Plex; 2165 dm->ops->computel2fielddiff = DMComputeL2FieldDiff_Plex; 2166 dm->ops->getneighbors = DMGetNeighors_Plex; 2167 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMPlexInsertBoundaryValues_C",DMPlexInsertBoundaryValues_Plex);CHKERRQ(ierr); 2168 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMSetUpGLVisViewer_C",DMSetUpGLVisViewer_Plex);CHKERRQ(ierr); 2169 PetscFunctionReturn(0); 2170 } 2171 2172 PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm) 2173 { 2174 DM_Plex *mesh = (DM_Plex *) dm->data; 2175 PetscErrorCode ierr; 2176 2177 PetscFunctionBegin; 2178 mesh->refct++; 2179 (*newdm)->data = mesh; 2180 ierr = PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);CHKERRQ(ierr); 2181 ierr = DMInitialize_Plex(*newdm);CHKERRQ(ierr); 2182 PetscFunctionReturn(0); 2183 } 2184 2185 /*MC 2186 DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram. 2187 In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is 2188 specified by a PetscSection object. Ownership in the global representation is determined by 2189 ownership of the underlying DMPlex points. This is specified by another PetscSection object. 2190 2191 Level: intermediate 2192 2193 .seealso: DMType, DMPlexCreate(), DMCreate(), DMSetType() 2194 M*/ 2195 2196 PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm) 2197 { 2198 DM_Plex *mesh; 2199 PetscInt unit, d; 2200 PetscErrorCode ierr; 2201 2202 PetscFunctionBegin; 2203 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2204 ierr = PetscNewLog(dm,&mesh);CHKERRQ(ierr); 2205 dm->dim = 0; 2206 dm->data = mesh; 2207 2208 mesh->refct = 1; 2209 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);CHKERRQ(ierr); 2210 mesh->maxConeSize = 0; 2211 mesh->cones = NULL; 2212 mesh->coneOrientations = NULL; 2213 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);CHKERRQ(ierr); 2214 mesh->maxSupportSize = 0; 2215 mesh->supports = NULL; 2216 mesh->refinementUniform = PETSC_TRUE; 2217 mesh->refinementLimit = -1.0; 2218 2219 mesh->facesTmp = NULL; 2220 2221 mesh->tetgenOpts = NULL; 2222 mesh->triangleOpts = NULL; 2223 ierr = PetscPartitionerCreate(PetscObjectComm((PetscObject)dm), &mesh->partitioner);CHKERRQ(ierr); 2224 mesh->remeshBd = PETSC_FALSE; 2225 2226 mesh->subpointMap = NULL; 2227 2228 for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0; 2229 2230 mesh->regularRefinement = PETSC_FALSE; 2231 mesh->depthState = -1; 2232 mesh->globalVertexNumbers = NULL; 2233 mesh->globalCellNumbers = NULL; 2234 mesh->anchorSection = NULL; 2235 mesh->anchorIS = NULL; 2236 mesh->createanchors = NULL; 2237 mesh->computeanchormatrix = NULL; 2238 mesh->parentSection = NULL; 2239 mesh->parents = NULL; 2240 mesh->childIDs = NULL; 2241 mesh->childSection = NULL; 2242 mesh->children = NULL; 2243 mesh->referenceTree = NULL; 2244 mesh->getchildsymmetry = NULL; 2245 for (d = 0; d < 8; ++d) mesh->hybridPointMax[d] = PETSC_DETERMINE; 2246 mesh->vtkCellHeight = 0; 2247 mesh->useCone = PETSC_FALSE; 2248 mesh->useClosure = PETSC_TRUE; 2249 mesh->useAnchors = PETSC_FALSE; 2250 2251 mesh->maxProjectionHeight = 0; 2252 2253 mesh->printSetValues = PETSC_FALSE; 2254 mesh->printFEM = 0; 2255 mesh->printTol = 1.0e-10; 2256 2257 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 2258 PetscFunctionReturn(0); 2259 } 2260 2261 /*@ 2262 DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram. 2263 2264 Collective on MPI_Comm 2265 2266 Input Parameter: 2267 . comm - The communicator for the DMPlex object 2268 2269 Output Parameter: 2270 . mesh - The DMPlex object 2271 2272 Level: beginner 2273 2274 .keywords: DMPlex, create 2275 @*/ 2276 PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh) 2277 { 2278 PetscErrorCode ierr; 2279 2280 PetscFunctionBegin; 2281 PetscValidPointer(mesh,2); 2282 ierr = DMCreate(comm, mesh);CHKERRQ(ierr); 2283 ierr = DMSetType(*mesh, DMPLEX);CHKERRQ(ierr); 2284 PetscFunctionReturn(0); 2285 } 2286 2287 /* 2288 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 2289 */ 2290 static PetscErrorCode DMPlexBuildFromCellList_Parallel_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[], PetscSF *sfVert) 2291 { 2292 PetscSF sfPoint; 2293 PetscLayout vLayout; 2294 PetscHashI vhash; 2295 PetscSFNode *remoteVerticesAdj, *vertexLocal, *vertexOwner, *remoteVertex; 2296 const PetscInt *vrange; 2297 PetscInt numVerticesAdj, off, *verticesAdj, numVerticesGhost = 0, *localVertex, *cone, c, p, v, g; 2298 PETSC_UNUSED PetscHashIIter ret, iter; 2299 PetscMPIInt rank, size; 2300 PetscErrorCode ierr; 2301 2302 PetscFunctionBegin; 2303 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr); 2304 ierr = MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);CHKERRQ(ierr); 2305 /* Partition vertices */ 2306 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) dm), &vLayout);CHKERRQ(ierr); 2307 ierr = PetscLayoutSetLocalSize(vLayout, numVertices);CHKERRQ(ierr); 2308 ierr = PetscLayoutSetBlockSize(vLayout, 1);CHKERRQ(ierr); 2309 ierr = PetscLayoutSetUp(vLayout);CHKERRQ(ierr); 2310 ierr = PetscLayoutGetRanges(vLayout, &vrange);CHKERRQ(ierr); 2311 /* Count vertices and map them to procs */ 2312 PetscHashICreate(vhash); 2313 for (c = 0; c < numCells; ++c) { 2314 for (p = 0; p < numCorners; ++p) { 2315 PetscHashIPut(vhash, cells[c*numCorners+p], ret, iter); 2316 } 2317 } 2318 PetscHashISize(vhash, numVerticesAdj); 2319 ierr = PetscMalloc1(numVerticesAdj, &verticesAdj);CHKERRQ(ierr); 2320 off = 0; ierr = PetscHashIGetKeys(vhash, &off, verticesAdj);CHKERRQ(ierr); 2321 if (off != numVerticesAdj) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid number of local vertices %D should be %D", off, numVerticesAdj); 2322 ierr = PetscSortInt(numVerticesAdj, verticesAdj);CHKERRQ(ierr); 2323 ierr = PetscMalloc1(numVerticesAdj, &remoteVerticesAdj);CHKERRQ(ierr); 2324 for (v = 0; v < numVerticesAdj; ++v) { 2325 const PetscInt gv = verticesAdj[v]; 2326 PetscInt vrank; 2327 2328 ierr = PetscFindInt(gv, size+1, vrange, &vrank);CHKERRQ(ierr); 2329 vrank = vrank < 0 ? -(vrank+2) : vrank; 2330 remoteVerticesAdj[v].index = gv - vrange[vrank]; 2331 remoteVerticesAdj[v].rank = vrank; 2332 } 2333 /* Create cones */ 2334 ierr = DMPlexSetChart(dm, 0, numCells+numVerticesAdj);CHKERRQ(ierr); 2335 for (c = 0; c < numCells; ++c) {ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr);} 2336 ierr = DMSetUp(dm);CHKERRQ(ierr); 2337 ierr = DMGetWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 2338 for (c = 0; c < numCells; ++c) { 2339 for (p = 0; p < numCorners; ++p) { 2340 const PetscInt gv = cells[c*numCorners+p]; 2341 PetscInt lv; 2342 2343 ierr = PetscFindInt(gv, numVerticesAdj, verticesAdj, &lv);CHKERRQ(ierr); 2344 if (lv < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not find global vertex %D in local connectivity", gv); 2345 cone[p] = lv+numCells; 2346 } 2347 ierr = DMPlexSetCone(dm, c, cone);CHKERRQ(ierr); 2348 } 2349 ierr = DMRestoreWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 2350 /* Create SF for vertices */ 2351 ierr = PetscSFCreate(PetscObjectComm((PetscObject)dm), sfVert);CHKERRQ(ierr); 2352 ierr = PetscObjectSetName((PetscObject) *sfVert, "Vertex Ownership SF");CHKERRQ(ierr); 2353 ierr = PetscSFSetFromOptions(*sfVert);CHKERRQ(ierr); 2354 ierr = PetscSFSetGraph(*sfVert, numVertices, numVerticesAdj, NULL, PETSC_OWN_POINTER, remoteVerticesAdj, PETSC_OWN_POINTER);CHKERRQ(ierr); 2355 ierr = PetscFree(verticesAdj);CHKERRQ(ierr); 2356 /* Build pointSF */ 2357 ierr = PetscMalloc2(numVerticesAdj, &vertexLocal, numVertices, &vertexOwner);CHKERRQ(ierr); 2358 for (v = 0; v < numVerticesAdj; ++v) {vertexLocal[v].index = v+numCells; vertexLocal[v].rank = rank;} 2359 for (v = 0; v < numVertices; ++v) {vertexOwner[v].index = -1; vertexOwner[v].rank = -1;} 2360 ierr = PetscSFReduceBegin(*sfVert, MPIU_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);CHKERRQ(ierr); 2361 ierr = PetscSFReduceEnd(*sfVert, MPIU_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);CHKERRQ(ierr); 2362 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); 2363 ierr = PetscSFBcastBegin(*sfVert, MPIU_2INT, vertexOwner, vertexLocal);CHKERRQ(ierr); 2364 ierr = PetscSFBcastEnd(*sfVert, MPIU_2INT, vertexOwner, vertexLocal);CHKERRQ(ierr); 2365 for (v = 0; v < numVerticesAdj; ++v) if (vertexLocal[v].rank != rank) ++numVerticesGhost; 2366 ierr = PetscMalloc1(numVerticesGhost, &localVertex);CHKERRQ(ierr); 2367 ierr = PetscMalloc1(numVerticesGhost, &remoteVertex);CHKERRQ(ierr); 2368 for (v = 0, g = 0; v < numVerticesAdj; ++v) { 2369 if (vertexLocal[v].rank != rank) { 2370 localVertex[g] = v+numCells; 2371 remoteVertex[g].index = vertexLocal[v].index; 2372 remoteVertex[g].rank = vertexLocal[v].rank; 2373 ++g; 2374 } 2375 } 2376 ierr = PetscFree2(vertexLocal, vertexOwner);CHKERRQ(ierr); 2377 if (g != numVerticesGhost) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of vertex ghosts %D should be %D", g, numVerticesGhost); 2378 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 2379 ierr = PetscObjectSetName((PetscObject) sfPoint, "point SF");CHKERRQ(ierr); 2380 ierr = PetscSFSetGraph(sfPoint, numCells+numVerticesAdj, numVerticesGhost, localVertex, PETSC_OWN_POINTER, remoteVertex, PETSC_OWN_POINTER);CHKERRQ(ierr); 2381 ierr = PetscLayoutDestroy(&vLayout);CHKERRQ(ierr); 2382 PetscHashIDestroy(vhash); 2383 /* Fill in the rest of the topology structure */ 2384 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 2385 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 2386 PetscFunctionReturn(0); 2387 } 2388 2389 /* 2390 This takes as input the coordinates for each owned vertex 2391 */ 2392 static PetscErrorCode DMPlexBuildCoordinates_Parallel_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numV, PetscSF sfVert, const PetscReal vertexCoords[]) 2393 { 2394 PetscSection coordSection; 2395 Vec coordinates; 2396 PetscScalar *coords; 2397 PetscInt numVertices, numVerticesAdj, coordSize, v; 2398 PetscErrorCode ierr; 2399 2400 PetscFunctionBegin; 2401 ierr = DMSetCoordinateDim(dm, spaceDim);CHKERRQ(ierr); 2402 ierr = PetscSFGetGraph(sfVert, &numVertices, &numVerticesAdj, NULL, NULL);CHKERRQ(ierr); 2403 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2404 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 2405 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 2406 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVerticesAdj);CHKERRQ(ierr); 2407 for (v = numCells; v < numCells+numVerticesAdj; ++v) { 2408 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 2409 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 2410 } 2411 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 2412 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 2413 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 2414 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 2415 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 2416 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 2417 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 2418 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2419 { 2420 MPI_Datatype coordtype; 2421 2422 /* Need a temp buffer for coords if we have complex/single */ 2423 ierr = MPI_Type_contiguous(spaceDim, MPIU_SCALAR, &coordtype);CHKERRQ(ierr); 2424 ierr = MPI_Type_commit(&coordtype);CHKERRQ(ierr); 2425 #if defined(PETSC_USE_COMPLEX) 2426 { 2427 PetscScalar *svertexCoords; 2428 PetscInt i; 2429 ierr = PetscMalloc1(numV*spaceDim,&svertexCoords);CHKERRQ(ierr); 2430 for (i=0; i<numV*spaceDim; i++) svertexCoords[i] = vertexCoords[i]; 2431 ierr = PetscSFBcastBegin(sfVert, coordtype, svertexCoords, coords);CHKERRQ(ierr); 2432 ierr = PetscSFBcastEnd(sfVert, coordtype, svertexCoords, coords);CHKERRQ(ierr); 2433 ierr = PetscFree(svertexCoords);CHKERRQ(ierr); 2434 } 2435 #else 2436 ierr = PetscSFBcastBegin(sfVert, coordtype, vertexCoords, coords);CHKERRQ(ierr); 2437 ierr = PetscSFBcastEnd(sfVert, coordtype, vertexCoords, coords);CHKERRQ(ierr); 2438 #endif 2439 ierr = MPI_Type_free(&coordtype);CHKERRQ(ierr); 2440 } 2441 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2442 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 2443 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 2444 PetscFunctionReturn(0); 2445 } 2446 2447 /*@C 2448 DMPlexCreateFromCellListParallel - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM 2449 2450 Input Parameters: 2451 + comm - The communicator 2452 . dim - The topological dimension of the mesh 2453 . numCells - The number of cells owned by this process 2454 . numVertices - The number of vertices owned by this process 2455 . numCorners - The number of vertices for each cell 2456 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 2457 . cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 2458 . spaceDim - The spatial dimension used for coordinates 2459 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 2460 2461 Output Parameter: 2462 + dm - The DM 2463 - vertexSF - Optional, SF describing complete vertex ownership 2464 2465 Note: Two triangles sharing a face 2466 $ 2467 $ 2 2468 $ / | \ 2469 $ / | \ 2470 $ / | \ 2471 $ 0 0 | 1 3 2472 $ \ | / 2473 $ \ | / 2474 $ \ | / 2475 $ 1 2476 would have input 2477 $ numCells = 2, numVertices = 4 2478 $ cells = [0 1 2 1 3 2] 2479 $ 2480 which would result in the DMPlex 2481 $ 2482 $ 4 2483 $ / | \ 2484 $ / | \ 2485 $ / | \ 2486 $ 2 0 | 1 5 2487 $ \ | / 2488 $ \ | / 2489 $ \ | / 2490 $ 3 2491 2492 Level: beginner 2493 2494 .seealso: DMPlexCreateFromCellList(), DMPlexCreateFromDAG(), DMPlexCreate() 2495 @*/ 2496 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) 2497 { 2498 PetscSF sfVert; 2499 PetscErrorCode ierr; 2500 2501 PetscFunctionBegin; 2502 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 2503 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 2504 PetscValidLogicalCollectiveInt(*dm, dim, 2); 2505 PetscValidLogicalCollectiveInt(*dm, spaceDim, 8); 2506 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 2507 ierr = DMPlexBuildFromCellList_Parallel_Private(*dm, numCells, numVertices, numCorners, cells, &sfVert);CHKERRQ(ierr); 2508 if (interpolate) { 2509 DM idm = NULL; 2510 2511 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 2512 ierr = DMDestroy(dm);CHKERRQ(ierr); 2513 *dm = idm; 2514 } 2515 ierr = DMPlexBuildCoordinates_Parallel_Private(*dm, spaceDim, numCells, numVertices,sfVert, vertexCoords);CHKERRQ(ierr); 2516 if (vertexSF) *vertexSF = sfVert; 2517 else {ierr = PetscSFDestroy(&sfVert);CHKERRQ(ierr);} 2518 PetscFunctionReturn(0); 2519 } 2520 2521 /* 2522 This takes as input the common mesh generator output, a list of the vertices for each cell 2523 */ 2524 static PetscErrorCode DMPlexBuildFromCellList_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[]) 2525 { 2526 PetscInt *cone, c, p; 2527 PetscErrorCode ierr; 2528 2529 PetscFunctionBegin; 2530 ierr = DMPlexSetChart(dm, 0, numCells+numVertices);CHKERRQ(ierr); 2531 for (c = 0; c < numCells; ++c) { 2532 ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr); 2533 } 2534 ierr = DMSetUp(dm);CHKERRQ(ierr); 2535 ierr = DMGetWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 2536 for (c = 0; c < numCells; ++c) { 2537 for (p = 0; p < numCorners; ++p) { 2538 cone[p] = cells[c*numCorners+p]+numCells; 2539 } 2540 ierr = DMPlexSetCone(dm, c, cone);CHKERRQ(ierr); 2541 } 2542 ierr = DMRestoreWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 2543 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 2544 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 2545 PetscFunctionReturn(0); 2546 } 2547 2548 /* 2549 This takes as input the coordinates for each vertex 2550 */ 2551 static PetscErrorCode DMPlexBuildCoordinates_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, const double vertexCoords[]) 2552 { 2553 PetscSection coordSection; 2554 Vec coordinates; 2555 DM cdm; 2556 PetscScalar *coords; 2557 PetscInt v, d; 2558 PetscErrorCode ierr; 2559 2560 PetscFunctionBegin; 2561 ierr = DMSetCoordinateDim(dm, spaceDim);CHKERRQ(ierr); 2562 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2563 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 2564 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 2565 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVertices);CHKERRQ(ierr); 2566 for (v = numCells; v < numCells+numVertices; ++v) { 2567 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 2568 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 2569 } 2570 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 2571 2572 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 2573 ierr = DMCreateLocalVector(cdm, &coordinates);CHKERRQ(ierr); 2574 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 2575 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 2576 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2577 for (v = 0; v < numVertices; ++v) { 2578 for (d = 0; d < spaceDim; ++d) { 2579 coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d]; 2580 } 2581 } 2582 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2583 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 2584 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 2585 PetscFunctionReturn(0); 2586 } 2587 2588 /*@C 2589 DMPlexCreateFromCellList - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM 2590 2591 Input Parameters: 2592 + comm - The communicator 2593 . dim - The topological dimension of the mesh 2594 . numCells - The number of cells 2595 . numVertices - The number of vertices 2596 . numCorners - The number of vertices for each cell 2597 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 2598 . cells - An array of numCells*numCorners numbers, the vertices for each cell 2599 . spaceDim - The spatial dimension used for coordinates 2600 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 2601 2602 Output Parameter: 2603 . dm - The DM 2604 2605 Note: Two triangles sharing a face 2606 $ 2607 $ 2 2608 $ / | \ 2609 $ / | \ 2610 $ / | \ 2611 $ 0 0 | 1 3 2612 $ \ | / 2613 $ \ | / 2614 $ \ | / 2615 $ 1 2616 would have input 2617 $ numCells = 2, numVertices = 4 2618 $ cells = [0 1 2 1 3 2] 2619 $ 2620 which would result in the DMPlex 2621 $ 2622 $ 4 2623 $ / | \ 2624 $ / | \ 2625 $ / | \ 2626 $ 2 0 | 1 5 2627 $ \ | / 2628 $ \ | / 2629 $ \ | / 2630 $ 3 2631 2632 Level: beginner 2633 2634 .seealso: DMPlexCreateFromDAG(), DMPlexCreate() 2635 @*/ 2636 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) 2637 { 2638 PetscErrorCode ierr; 2639 2640 PetscFunctionBegin; 2641 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 2642 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 2643 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 2644 ierr = DMPlexBuildFromCellList_Private(*dm, numCells, numVertices, numCorners, cells);CHKERRQ(ierr); 2645 if (interpolate) { 2646 DM idm = NULL; 2647 2648 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 2649 ierr = DMDestroy(dm);CHKERRQ(ierr); 2650 *dm = idm; 2651 } 2652 ierr = DMPlexBuildCoordinates_Private(*dm, spaceDim, numCells, numVertices, vertexCoords);CHKERRQ(ierr); 2653 PetscFunctionReturn(0); 2654 } 2655 2656 /*@ 2657 DMPlexCreateFromDAG - This takes as input the adjacency-list representation of the Directed Acyclic Graph (Hasse Diagram) encoding a mesh, and produces a DM 2658 2659 Input Parameters: 2660 + dm - The empty DM object, usually from DMCreate() and DMSetDimension() 2661 . depth - The depth of the DAG 2662 . numPoints - The number of points at each depth 2663 . coneSize - The cone size of each point 2664 . cones - The concatenation of the cone points for each point, the cone list must be oriented correctly for each point 2665 . coneOrientations - The orientation of each cone point 2666 - vertexCoords - An array of numVertices*dim numbers, the coordinates of each vertex 2667 2668 Output Parameter: 2669 . dm - The DM 2670 2671 Note: Two triangles sharing a face would have input 2672 $ depth = 1, numPoints = [4 2], coneSize = [3 3 0 0 0 0] 2673 $ cones = [2 3 4 3 5 4], coneOrientations = [0 0 0 0 0 0] 2674 $ vertexCoords = [-1.0 0.0 0.0 -1.0 0.0 1.0 1.0 0.0] 2675 $ 2676 which would result in the DMPlex 2677 $ 2678 $ 4 2679 $ / | \ 2680 $ / | \ 2681 $ / | \ 2682 $ 2 0 | 1 5 2683 $ \ | / 2684 $ \ | / 2685 $ \ | / 2686 $ 3 2687 $ 2688 $ Notice that all points are numbered consecutively, unlikely DMPlexCreateFromCellList() 2689 2690 Level: advanced 2691 2692 .seealso: DMPlexCreateFromCellList(), DMPlexCreate() 2693 @*/ 2694 PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[]) 2695 { 2696 Vec coordinates; 2697 PetscSection coordSection; 2698 PetscScalar *coords; 2699 PetscInt coordSize, firstVertex = -1, pStart = 0, pEnd = 0, p, v, dim, dimEmbed, d, off; 2700 PetscErrorCode ierr; 2701 2702 PetscFunctionBegin; 2703 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2704 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 2705 if (dimEmbed < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Embedding dimension %d cannot be less than intrinsic dimension %d",dimEmbed,dim); 2706 for (d = 0; d <= depth; ++d) pEnd += numPoints[d]; 2707 ierr = DMPlexSetChart(dm, pStart, pEnd);CHKERRQ(ierr); 2708 for (p = pStart; p < pEnd; ++p) { 2709 ierr = DMPlexSetConeSize(dm, p, coneSize[p-pStart]);CHKERRQ(ierr); 2710 if (firstVertex < 0 && !coneSize[p - pStart]) { 2711 firstVertex = p - pStart; 2712 } 2713 } 2714 if (firstVertex < 0 && numPoints[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected %d vertices but could not find any", numPoints[0]); 2715 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 2716 for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) { 2717 ierr = DMPlexSetCone(dm, p, &cones[off]);CHKERRQ(ierr); 2718 ierr = DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);CHKERRQ(ierr); 2719 } 2720 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 2721 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 2722 /* Build coordinates */ 2723 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2724 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 2725 ierr = PetscSectionSetFieldComponents(coordSection, 0, dimEmbed);CHKERRQ(ierr); 2726 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0]);CHKERRQ(ierr); 2727 for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) { 2728 ierr = PetscSectionSetDof(coordSection, v, dimEmbed);CHKERRQ(ierr); 2729 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed);CHKERRQ(ierr); 2730 } 2731 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 2732 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 2733 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 2734 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 2735 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 2736 ierr = VecSetBlockSize(coordinates, dimEmbed);CHKERRQ(ierr); 2737 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 2738 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2739 for (v = 0; v < numPoints[0]; ++v) { 2740 PetscInt off; 2741 2742 ierr = PetscSectionGetOffset(coordSection, v+firstVertex, &off);CHKERRQ(ierr); 2743 for (d = 0; d < dimEmbed; ++d) { 2744 coords[off+d] = vertexCoords[v*dimEmbed+d]; 2745 } 2746 } 2747 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2748 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 2749 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 2750 PetscFunctionReturn(0); 2751 } 2752 2753 /*@C 2754 DMPlexCreateFromFile - This takes a filename and produces a DM 2755 2756 Input Parameters: 2757 + comm - The communicator 2758 . filename - A file name 2759 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 2760 2761 Output Parameter: 2762 . dm - The DM 2763 2764 Level: beginner 2765 2766 .seealso: DMPlexCreateFromDAG(), DMPlexCreateFromCellList(), DMPlexCreate() 2767 @*/ 2768 PetscErrorCode DMPlexCreateFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 2769 { 2770 const char *extGmsh = ".msh"; 2771 const char *extCGNS = ".cgns"; 2772 const char *extExodus = ".exo"; 2773 const char *extGenesis = ".gen"; 2774 const char *extFluent = ".cas"; 2775 const char *extHDF5 = ".h5"; 2776 const char *extMed = ".med"; 2777 const char *extPLY = ".ply"; 2778 size_t len; 2779 PetscBool isGmsh, isCGNS, isExodus, isGenesis, isFluent, isHDF5, isMed, isPLY; 2780 PetscMPIInt rank; 2781 PetscErrorCode ierr; 2782 2783 PetscFunctionBegin; 2784 PetscValidPointer(filename, 2); 2785 PetscValidPointer(dm, 4); 2786 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 2787 ierr = PetscStrlen(filename, &len);CHKERRQ(ierr); 2788 if (!len) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Filename must be a valid path"); 2789 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGmsh, 4, &isGmsh);CHKERRQ(ierr); 2790 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extCGNS, 5, &isCGNS);CHKERRQ(ierr); 2791 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extExodus, 4, &isExodus);CHKERRQ(ierr); 2792 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGenesis, 4, &isGenesis);CHKERRQ(ierr); 2793 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extFluent, 4, &isFluent);CHKERRQ(ierr); 2794 ierr = PetscStrncmp(&filename[PetscMax(0,len-3)], extHDF5, 3, &isHDF5);CHKERRQ(ierr); 2795 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extMed, 4, &isMed);CHKERRQ(ierr); 2796 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extPLY, 4, &isPLY);CHKERRQ(ierr); 2797 if (isGmsh) { 2798 ierr = DMPlexCreateGmshFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2799 } else if (isCGNS) { 2800 ierr = DMPlexCreateCGNSFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2801 } else if (isExodus || isGenesis) { 2802 ierr = DMPlexCreateExodusFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2803 } else if (isFluent) { 2804 ierr = DMPlexCreateFluentFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2805 } else if (isHDF5) { 2806 PetscViewer viewer; 2807 2808 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 2809 ierr = PetscViewerSetType(viewer, PETSCVIEWERHDF5);CHKERRQ(ierr); 2810 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 2811 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 2812 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 2813 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 2814 ierr = DMLoad(*dm, viewer);CHKERRQ(ierr); 2815 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 2816 } else if (isMed) { 2817 ierr = DMPlexCreateMedFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2818 } else if (isPLY) { 2819 ierr = DMPlexCreatePLYFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2820 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot load file %s: unrecognized extension", filename); 2821 PetscFunctionReturn(0); 2822 } 2823 2824 /*@ 2825 DMPlexCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell 2826 2827 Collective on comm 2828 2829 Input Parameters: 2830 + comm - The communicator 2831 . dim - The spatial dimension 2832 - simplex - Flag for simplex, otherwise use a tensor-product cell 2833 2834 Output Parameter: 2835 . refdm - The reference cell 2836 2837 Level: intermediate 2838 2839 .keywords: reference cell 2840 .seealso: 2841 @*/ 2842 PetscErrorCode DMPlexCreateReferenceCell(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *refdm) 2843 { 2844 DM rdm; 2845 Vec coords; 2846 PetscErrorCode ierr; 2847 2848 PetscFunctionBeginUser; 2849 ierr = DMCreate(comm, &rdm);CHKERRQ(ierr); 2850 ierr = DMSetType(rdm, DMPLEX);CHKERRQ(ierr); 2851 ierr = DMSetDimension(rdm, dim);CHKERRQ(ierr); 2852 switch (dim) { 2853 case 0: 2854 { 2855 PetscInt numPoints[1] = {1}; 2856 PetscInt coneSize[1] = {0}; 2857 PetscInt cones[1] = {0}; 2858 PetscInt coneOrientations[1] = {0}; 2859 PetscScalar vertexCoords[1] = {0.0}; 2860 2861 ierr = DMPlexCreateFromDAG(rdm, 0, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2862 } 2863 break; 2864 case 1: 2865 { 2866 PetscInt numPoints[2] = {2, 1}; 2867 PetscInt coneSize[3] = {2, 0, 0}; 2868 PetscInt cones[2] = {1, 2}; 2869 PetscInt coneOrientations[2] = {0, 0}; 2870 PetscScalar vertexCoords[2] = {-1.0, 1.0}; 2871 2872 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2873 } 2874 break; 2875 case 2: 2876 if (simplex) { 2877 PetscInt numPoints[2] = {3, 1}; 2878 PetscInt coneSize[4] = {3, 0, 0, 0}; 2879 PetscInt cones[3] = {1, 2, 3}; 2880 PetscInt coneOrientations[3] = {0, 0, 0}; 2881 PetscScalar vertexCoords[6] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0}; 2882 2883 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2884 } else { 2885 PetscInt numPoints[2] = {4, 1}; 2886 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 2887 PetscInt cones[4] = {1, 2, 3, 4}; 2888 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 2889 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0}; 2890 2891 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2892 } 2893 break; 2894 case 3: 2895 if (simplex) { 2896 PetscInt numPoints[2] = {4, 1}; 2897 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 2898 PetscInt cones[4] = {1, 3, 2, 4}; 2899 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 2900 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}; 2901 2902 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2903 } else { 2904 PetscInt numPoints[2] = {8, 1}; 2905 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0}; 2906 PetscInt cones[8] = {1, 4, 3, 2, 5, 6, 7, 8}; 2907 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 2908 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, 2909 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0}; 2910 2911 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2912 } 2913 break; 2914 default: 2915 SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Cannot create reference cell for dimension %d", dim); 2916 } 2917 *refdm = NULL; 2918 ierr = DMPlexInterpolate(rdm, refdm);CHKERRQ(ierr); 2919 if (rdm->coordinateDM) { 2920 DM ncdm; 2921 PetscSection cs; 2922 PetscInt pEnd = -1; 2923 2924 ierr = DMGetDefaultSection(rdm->coordinateDM, &cs);CHKERRQ(ierr); 2925 if (cs) {ierr = PetscSectionGetChart(cs, NULL, &pEnd);CHKERRQ(ierr);} 2926 if (pEnd >= 0) { 2927 ierr = DMClone(rdm->coordinateDM, &ncdm);CHKERRQ(ierr); 2928 ierr = DMSetDefaultSection(ncdm, cs);CHKERRQ(ierr); 2929 ierr = DMSetCoordinateDM(*refdm, ncdm);CHKERRQ(ierr); 2930 ierr = DMDestroy(&ncdm);CHKERRQ(ierr); 2931 } 2932 } 2933 ierr = DMGetCoordinatesLocal(rdm, &coords);CHKERRQ(ierr); 2934 if (coords) { 2935 ierr = DMSetCoordinatesLocal(*refdm, coords);CHKERRQ(ierr); 2936 } else { 2937 ierr = DMGetCoordinates(rdm, &coords);CHKERRQ(ierr); 2938 if (coords) {ierr = DMSetCoordinates(*refdm, coords);CHKERRQ(ierr);} 2939 } 2940 ierr = DMDestroy(&rdm);CHKERRQ(ierr); 2941 PetscFunctionReturn(0); 2942 } 2943