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