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