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 PetscInt markerTop = 1, faceMarkerTop = 1; 499 PetscInt markerBottom = 1, faceMarkerBottom = 1; 500 PetscInt markerFront = 1, faceMarkerFront = 1; 501 PetscInt markerBack = 1, faceMarkerBack = 1; 502 PetscInt markerRight = 1, faceMarkerRight = 1; 503 PetscInt markerLeft = 1, faceMarkerLeft = 1; 504 PetscInt dim; 505 PetscBool markerSeparate = PETSC_FALSE; 506 PetscMPIInt rank; 507 PetscErrorCode ierr; 508 509 PetscFunctionBegin; 510 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 511 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 512 ierr = DMCreateLabel(dm,"marker");CHKERRQ(ierr); 513 ierr = DMCreateLabel(dm,"Face Sets");CHKERRQ(ierr); 514 switch (dim) { 515 case 2: 516 faceMarkerTop = 3; 517 faceMarkerBottom = 1; 518 faceMarkerRight = 2; 519 faceMarkerLeft = 4; 520 break; 521 case 3: 522 faceMarkerBottom = 1; 523 faceMarkerTop = 2; 524 faceMarkerFront = 3; 525 faceMarkerBack = 4; 526 faceMarkerRight = 5; 527 faceMarkerLeft = 6; 528 break; 529 default: 530 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Dimension %d not supported",dim); 531 break; 532 } 533 ierr = PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);CHKERRQ(ierr); 534 if (markerSeparate) { 535 markerBottom = faceMarkerBottom; 536 markerTop = faceMarkerTop; 537 markerFront = faceMarkerFront; 538 markerBack = faceMarkerBack; 539 markerRight = faceMarkerRight; 540 markerLeft = faceMarkerLeft; 541 } 542 { 543 const PetscInt numXEdges = !rank ? edges[0] : 0; 544 const PetscInt numYEdges = !rank ? edges[1] : 0; 545 const PetscInt numZEdges = !rank ? edges[2] : 0; 546 const PetscInt numXVertices = !rank ? (bdX == DM_BOUNDARY_PERIODIC || bdX == DM_BOUNDARY_TWIST ? edges[0] : edges[0]+1) : 0; 547 const PetscInt numYVertices = !rank ? (bdY == DM_BOUNDARY_PERIODIC || bdY == DM_BOUNDARY_TWIST ? edges[1] : edges[1]+1) : 0; 548 const PetscInt numZVertices = !rank ? (bdZ == DM_BOUNDARY_PERIODIC || bdZ == DM_BOUNDARY_TWIST ? edges[2] : edges[2]+1) : 0; 549 const PetscInt numCells = numXEdges*numYEdges*numZEdges; 550 const PetscInt numXFaces = numYEdges*numZEdges; 551 const PetscInt numYFaces = numXEdges*numZEdges; 552 const PetscInt numZFaces = numXEdges*numYEdges; 553 const PetscInt numTotXFaces = numXVertices*numXFaces; 554 const PetscInt numTotYFaces = numYVertices*numYFaces; 555 const PetscInt numTotZFaces = numZVertices*numZFaces; 556 const PetscInt numFaces = numTotXFaces + numTotYFaces + numTotZFaces; 557 const PetscInt numTotXEdges = numXEdges*numYVertices*numZVertices; 558 const PetscInt numTotYEdges = numYEdges*numXVertices*numZVertices; 559 const PetscInt numTotZEdges = numZEdges*numXVertices*numYVertices; 560 const PetscInt numVertices = numXVertices*numYVertices*numZVertices; 561 const PetscInt numEdges = numTotXEdges + numTotYEdges + numTotZEdges; 562 const PetscInt firstVertex = (dim == 2) ? numFaces : numCells; 563 const PetscInt firstXFace = (dim == 2) ? 0 : numCells + numVertices; 564 const PetscInt firstYFace = firstXFace + numTotXFaces; 565 const PetscInt firstZFace = firstYFace + numTotYFaces; 566 const PetscInt firstXEdge = numCells + numFaces + numVertices; 567 const PetscInt firstYEdge = firstXEdge + numTotXEdges; 568 const PetscInt firstZEdge = firstYEdge + numTotYEdges; 569 Vec coordinates; 570 PetscSection coordSection; 571 PetscScalar *coords; 572 PetscInt coordSize; 573 PetscInt v, vx, vy, vz; 574 PetscInt c, f, fx, fy, fz, e, ex, ey, ez; 575 576 ierr = DMPlexSetChart(dm, 0, numCells+numFaces+numEdges+numVertices);CHKERRQ(ierr); 577 for (c = 0; c < numCells; c++) { 578 ierr = DMPlexSetConeSize(dm, c, 6);CHKERRQ(ierr); 579 } 580 for (f = firstXFace; f < firstXFace+numFaces; ++f) { 581 ierr = DMPlexSetConeSize(dm, f, 4);CHKERRQ(ierr); 582 } 583 for (e = firstXEdge; e < firstXEdge+numEdges; ++e) { 584 ierr = DMPlexSetConeSize(dm, e, 2);CHKERRQ(ierr); 585 } 586 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 587 /* Build cells */ 588 for (fz = 0; fz < numZEdges; ++fz) { 589 for (fy = 0; fy < numYEdges; ++fy) { 590 for (fx = 0; fx < numXEdges; ++fx) { 591 PetscInt cell = (fz*numYEdges + fy)*numXEdges + fx; 592 PetscInt faceB = firstZFace + (fy*numXEdges+fx)*numZVertices + fz; 593 PetscInt faceT = firstZFace + (fy*numXEdges+fx)*numZVertices + ((fz+1)%numZVertices); 594 PetscInt faceF = firstYFace + (fz*numXEdges+fx)*numYVertices + fy; 595 PetscInt faceK = firstYFace + (fz*numXEdges+fx)*numYVertices + ((fy+1)%numYVertices); 596 PetscInt faceL = firstXFace + (fz*numYEdges+fy)*numXVertices + fx; 597 PetscInt faceR = firstXFace + (fz*numYEdges+fy)*numXVertices + ((fx+1)%numXVertices); 598 /* B, T, F, K, R, L */ 599 PetscInt ornt[6] = {-4, 0, 0, -1, 0, -4}; /* ??? */ 600 PetscInt cone[6]; 601 602 /* no boundary twisting in 3D */ 603 cone[0] = faceB; cone[1] = faceT; cone[2] = faceF; cone[3] = faceK; cone[4] = faceR; cone[5] = faceL; 604 ierr = DMPlexSetCone(dm, cell, cone);CHKERRQ(ierr); 605 ierr = DMPlexSetConeOrientation(dm, cell, ornt);CHKERRQ(ierr); 606 } 607 } 608 } 609 /* Build x faces */ 610 for (fz = 0; fz < numZEdges; ++fz) { 611 for (fy = 0; fy < numYEdges; ++fy) { 612 for (fx = 0; fx < numXVertices; ++fx) { 613 PetscInt face = firstXFace + (fz*numYEdges+fy)*numXVertices + fx; 614 PetscInt edgeL = firstZEdge + ( fy* numXVertices+fx)*numZEdges + fz; 615 PetscInt edgeR = firstZEdge + (((fy+1)%numYVertices)*numXVertices+fx)*numZEdges + fz; 616 PetscInt edgeB = firstYEdge + ( fz* numXVertices+fx)*numYEdges + fy; 617 PetscInt edgeT = firstYEdge + (((fz+1)%numZVertices)*numXVertices+fx)*numYEdges + fy; 618 PetscInt ornt[4] = {0, 0, -2, -2}; 619 PetscInt cone[4]; 620 621 if (dim == 3) { 622 /* markers */ 623 if (bdX != DM_BOUNDARY_PERIODIC) { 624 if (fx == numXVertices-1) { 625 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerRight);CHKERRQ(ierr); 626 ierr = DMSetLabelValue(dm, "marker", face, markerRight);CHKERRQ(ierr); 627 } 628 else if (fx == 0) { 629 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerLeft);CHKERRQ(ierr); 630 ierr = DMSetLabelValue(dm, "marker", face, markerLeft);CHKERRQ(ierr); 631 } 632 } 633 } 634 cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL; 635 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr); 636 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr); 637 } 638 } 639 } 640 /* Build y faces */ 641 for (fz = 0; fz < numZEdges; ++fz) { 642 for (fx = 0; fx < numXEdges; ++fx) { 643 for (fy = 0; fy < numYVertices; ++fy) { 644 PetscInt face = firstYFace + (fz*numXEdges+fx)*numYVertices + fy; 645 PetscInt edgeL = firstZEdge + (fy*numXVertices+ fx )*numZEdges + fz; 646 PetscInt edgeR = firstZEdge + (fy*numXVertices+((fx+1)%numXVertices))*numZEdges + fz; 647 PetscInt edgeB = firstXEdge + ( fz *numYVertices+fy)*numXEdges + fx; 648 PetscInt edgeT = firstXEdge + (((fz+1)%numZVertices)*numYVertices+fy)*numXEdges + fx; 649 PetscInt ornt[4] = {0, 0, -2, -2}; 650 PetscInt cone[4]; 651 652 if (dim == 3) { 653 /* markers */ 654 if (bdY != DM_BOUNDARY_PERIODIC) { 655 if (fy == numYVertices-1) { 656 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerBack);CHKERRQ(ierr); 657 ierr = DMSetLabelValue(dm, "marker", face, markerBack);CHKERRQ(ierr); 658 } 659 else if (fy == 0) { 660 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerFront);CHKERRQ(ierr); 661 ierr = DMSetLabelValue(dm, "marker", face, markerFront);CHKERRQ(ierr); 662 } 663 } 664 } 665 cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL; 666 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr); 667 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr); 668 } 669 } 670 } 671 /* Build z faces */ 672 for (fy = 0; fy < numYEdges; ++fy) { 673 for (fx = 0; fx < numXEdges; ++fx) { 674 for (fz = 0; fz < numZVertices; fz++) { 675 PetscInt face = firstZFace + (fy*numXEdges+fx)*numZVertices + fz; 676 PetscInt edgeL = firstYEdge + (fz*numXVertices+ fx )*numYEdges + fy; 677 PetscInt edgeR = firstYEdge + (fz*numXVertices+((fx+1)%numXVertices))*numYEdges + fy; 678 PetscInt edgeB = firstXEdge + (fz*numYVertices+ fy )*numXEdges + fx; 679 PetscInt edgeT = firstXEdge + (fz*numYVertices+((fy+1)%numYVertices))*numXEdges + fx; 680 PetscInt ornt[4] = {0, 0, -2, -2}; 681 PetscInt cone[4]; 682 683 if (dim == 2) { 684 if (bdX == DM_BOUNDARY_TWIST && fx == numXEdges-1) {edgeR += numYEdges-1-2*fy; ornt[1] = -2;} 685 if (bdY == DM_BOUNDARY_TWIST && fy == numYEdges-1) {edgeT += numXEdges-1-2*fx; ornt[2] = 0;} 686 } 687 else { 688 /* markers */ 689 if (bdZ != DM_BOUNDARY_PERIODIC) { 690 if (fz == numZVertices-1) { 691 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerTop);CHKERRQ(ierr); 692 ierr = DMSetLabelValue(dm, "marker", face, markerTop);CHKERRQ(ierr); 693 } 694 else if (fz == 0) { 695 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerBottom);CHKERRQ(ierr); 696 ierr = DMSetLabelValue(dm, "marker", face, markerBottom);CHKERRQ(ierr); 697 } 698 } 699 } 700 cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL; 701 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr); 702 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr); 703 } 704 } 705 } 706 /* Build Z edges*/ 707 for (vy = 0; vy < numYVertices; vy++) { 708 for (vx = 0; vx < numXVertices; vx++) { 709 for (ez = 0; ez < numZEdges; ez++) { 710 const PetscInt edge = firstZEdge + (vy*numXVertices+vx)*numZEdges + ez; 711 const PetscInt vertexB = firstVertex + ( ez *numYVertices+vy)*numXVertices + vx; 712 const PetscInt vertexT = firstVertex + (((ez+1)%numZVertices)*numYVertices+vy)*numXVertices + vx; 713 PetscInt cone[2]; 714 715 if (dim == 3) { 716 if (bdX != DM_BOUNDARY_PERIODIC) { 717 if (vx == numXVertices-1) { 718 ierr = DMSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 719 } 720 else if (vx == 0) { 721 ierr = DMSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 722 } 723 } 724 if (bdY != DM_BOUNDARY_PERIODIC) { 725 if (vy == numYVertices-1) { 726 ierr = DMSetLabelValue(dm, "marker", edge, markerBack);CHKERRQ(ierr); 727 } 728 else if (vy == 0) { 729 ierr = DMSetLabelValue(dm, "marker", edge, markerFront);CHKERRQ(ierr); 730 } 731 } 732 } 733 cone[0] = vertexB; cone[1] = vertexT; 734 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 735 } 736 } 737 } 738 /* Build Y edges*/ 739 for (vz = 0; vz < numZVertices; vz++) { 740 for (vx = 0; vx < numXVertices; vx++) { 741 for (ey = 0; ey < numYEdges; ey++) { 742 const PetscInt nextv = (dim == 2 && bdY == DM_BOUNDARY_TWIST && ey == numYEdges-1) ? (numXVertices-vx-1) : (vz*numYVertices+((ey+1)%numYVertices))*numXVertices + vx; 743 const PetscInt edge = firstYEdge + (vz*numXVertices+vx)*numYEdges + ey; 744 const PetscInt vertexF = firstVertex + (vz*numYVertices+ey)*numXVertices + vx; 745 const PetscInt vertexK = firstVertex + nextv; 746 PetscInt cone[2]; 747 748 cone[0] = vertexF; cone[1] = vertexK; 749 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 750 if (dim == 2) { 751 if ((bdX != DM_BOUNDARY_PERIODIC) && (bdX != DM_BOUNDARY_TWIST)) { 752 if (vx == numXVertices-1) { 753 ierr = DMSetLabelValue(dm, "Face Sets", edge, faceMarkerRight);CHKERRQ(ierr); 754 ierr = DMSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 755 ierr = DMSetLabelValue(dm, "marker", cone[0], markerRight);CHKERRQ(ierr); 756 if (ey == numYEdges-1) { 757 ierr = DMSetLabelValue(dm, "marker", cone[1], markerRight);CHKERRQ(ierr); 758 } 759 } 760 else if (vx == 0) { 761 ierr = DMSetLabelValue(dm, "Face Sets", edge, faceMarkerLeft);CHKERRQ(ierr); 762 ierr = DMSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 763 ierr = DMSetLabelValue(dm, "marker", cone[0], markerLeft);CHKERRQ(ierr); 764 if (ey == numYEdges-1) { 765 ierr = DMSetLabelValue(dm, "marker", cone[1], markerLeft);CHKERRQ(ierr); 766 } 767 } 768 } 769 } 770 else { 771 if (bdX != DM_BOUNDARY_PERIODIC) { 772 if (vx == numXVertices-1) { 773 ierr = DMSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 774 } 775 else if (vx == 0) { 776 ierr = DMSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 777 } 778 } 779 if (bdZ != DM_BOUNDARY_PERIODIC) { 780 if (vz == numZVertices-1) { 781 ierr = DMSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 782 } 783 else if (vz == 0) { 784 ierr = DMSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 785 } 786 } 787 } 788 } 789 } 790 } 791 /* Build X edges*/ 792 for (vz = 0; vz < numZVertices; vz++) { 793 for (vy = 0; vy < numYVertices; vy++) { 794 for (ex = 0; ex < numXEdges; ex++) { 795 const PetscInt nextv = (dim == 2 && bdX == DM_BOUNDARY_TWIST && ex == numXEdges-1) ? (numYVertices-vy-1)*numXVertices : (vz*numYVertices+vy)*numXVertices + (ex+1)%numXVertices; 796 const PetscInt edge = firstXEdge + (vz*numYVertices+vy)*numXEdges + ex; 797 const PetscInt vertexL = firstVertex + (vz*numYVertices+vy)*numXVertices + ex; 798 const PetscInt vertexR = firstVertex + nextv; 799 PetscInt cone[2]; 800 801 cone[0] = vertexL; cone[1] = vertexR; 802 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 803 if (dim == 2) { 804 if ((bdY != DM_BOUNDARY_PERIODIC) && (bdY != DM_BOUNDARY_TWIST)) { 805 if (vy == numYVertices-1) { 806 ierr = DMSetLabelValue(dm, "Face Sets", edge, faceMarkerTop);CHKERRQ(ierr); 807 ierr = DMSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 808 ierr = DMSetLabelValue(dm, "marker", cone[0], markerTop);CHKERRQ(ierr); 809 if (ex == numXEdges-1) { 810 ierr = DMSetLabelValue(dm, "marker", cone[1], markerTop);CHKERRQ(ierr); 811 } 812 } 813 else if (vy == 0) { 814 ierr = DMSetLabelValue(dm, "Face Sets", edge, faceMarkerBottom);CHKERRQ(ierr); 815 ierr = DMSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 816 ierr = DMSetLabelValue(dm, "marker", cone[0], markerBottom);CHKERRQ(ierr); 817 if (ex == numXEdges-1) { 818 ierr = DMSetLabelValue(dm, "marker", cone[1], markerBottom);CHKERRQ(ierr); 819 } 820 } 821 } 822 } 823 else { 824 if (bdY != DM_BOUNDARY_PERIODIC) { 825 if (vy == numYVertices-1) { 826 ierr = DMSetLabelValue(dm, "marker", edge, markerBack);CHKERRQ(ierr); 827 } 828 else if (vy == 0) { 829 ierr = DMSetLabelValue(dm, "marker", edge, markerFront);CHKERRQ(ierr); 830 } 831 } 832 if (bdZ != DM_BOUNDARY_PERIODIC) { 833 if (vz == numZVertices-1) { 834 ierr = DMSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 835 } 836 else if (vz == 0) { 837 ierr = DMSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 838 } 839 } 840 } 841 } 842 } 843 } 844 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 845 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 846 /* Build coordinates */ 847 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 848 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 849 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 850 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVertices);CHKERRQ(ierr); 851 for (v = firstVertex; v < firstVertex+numVertices; ++v) { 852 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 853 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 854 } 855 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 856 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 857 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 858 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 859 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 860 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 861 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 862 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 863 for (vz = 0; vz < numZVertices; ++vz) { 864 for (vy = 0; vy < numYVertices; ++vy) { 865 for (vx = 0; vx < numXVertices; ++vx) { 866 coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+0] = lower[0] + ((upper[0] - lower[0])/numXEdges)*vx; 867 coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+1] = lower[1] + ((upper[1] - lower[1])/numYEdges)*vy; 868 if (dim == 3) { 869 coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+2] = lower[2] + ((upper[2] - lower[2])/numZEdges)*vz; 870 } 871 } 872 } 873 } 874 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 875 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 876 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 877 } 878 PetscFunctionReturn(0); 879 } 880 881 /*@ 882 DMPlexCreateHexBoxMesh - Creates a mesh on the tensor product of unit intervals (box) using hexahedra. 883 884 Collective on MPI_Comm 885 886 Input Parameters: 887 + comm - The communicator for the DM object 888 . dim - The spatial dimension 889 . periodicX - The boundary type for the X direction 890 . periodicY - The boundary type for the Y direction 891 . periodicZ - The boundary type for the Z direction 892 - cells - The number of cells in each direction 893 894 Output Parameter: 895 . dm - The DM object 896 897 Note: Here is the numbering returned for 2 cells in each direction: 898 $ 22--8-23--9--24 899 $ | | | 900 $ 13 2 14 3 15 901 $ | | | 902 $ 19--6-20--7--21 903 $ | | | 904 $ 10 0 11 1 12 905 $ | | | 906 $ 16--4-17--5--18 907 908 Level: beginner 909 910 .keywords: DM, create 911 .seealso: DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 912 @*/ 913 PetscErrorCode DMPlexCreateHexBoxMesh(MPI_Comm comm, PetscInt dim, const PetscInt cells[], DMBoundaryType periodicX, DMBoundaryType periodicY, DMBoundaryType periodicZ, DM *dm) 914 { 915 PetscInt i; 916 PetscErrorCode ierr; 917 918 PetscFunctionBegin; 919 PetscValidPointer(dm, 7); 920 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 921 PetscValidLogicalCollectiveInt(*dm,dim,2); 922 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 923 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 924 switch (dim) { 925 case 2: 926 { 927 PetscReal lower[3] = {0.0, 0.0, 0.0}; 928 PetscReal upper[3] = {1.0, 1.0, 0.0}; 929 930 ierr = DMPlexCreateCubeMesh_Internal(*dm, lower, upper, cells, periodicX, periodicY, DM_BOUNDARY_NONE);CHKERRQ(ierr); 931 if (periodicX == DM_BOUNDARY_PERIODIC || periodicX == DM_BOUNDARY_TWIST || 932 periodicY == DM_BOUNDARY_PERIODIC || periodicY == DM_BOUNDARY_TWIST) { 933 PetscReal L[2]; 934 PetscReal maxCell[2]; 935 DMBoundaryType bdType[2]; 936 937 bdType[0] = periodicX; 938 bdType[1] = periodicY; 939 for (i = 0; i < dim; i++) { 940 L[i] = upper[i] - lower[i]; 941 maxCell[i] = 1.1 * (L[i] / cells[i]); 942 } 943 944 ierr = DMSetPeriodicity(*dm,maxCell,L,bdType);CHKERRQ(ierr); 945 } 946 break; 947 } 948 case 3: 949 { 950 PetscReal lower[3] = {0.0, 0.0, 0.0}; 951 PetscReal upper[3] = {1.0, 1.0, 1.0}; 952 953 ierr = DMPlexCreateCubeMesh_Internal(*dm, lower, upper, cells, periodicX, periodicY, periodicZ);CHKERRQ(ierr); 954 if (periodicX == DM_BOUNDARY_PERIODIC || periodicX == DM_BOUNDARY_TWIST || 955 periodicY == DM_BOUNDARY_PERIODIC || periodicY == DM_BOUNDARY_TWIST || 956 periodicZ == DM_BOUNDARY_PERIODIC || periodicZ == DM_BOUNDARY_TWIST) { 957 PetscReal L[3]; 958 PetscReal maxCell[3]; 959 DMBoundaryType bdType[3]; 960 961 bdType[0] = periodicX; 962 bdType[1] = periodicY; 963 bdType[2] = periodicZ; 964 for (i = 0; i < dim; i++) { 965 L[i] = upper[i] - lower[i]; 966 maxCell[i] = 1.1 * (L[i] / cells[i]); 967 } 968 969 ierr = DMSetPeriodicity(*dm,maxCell,L,bdType);CHKERRQ(ierr); 970 } 971 break; 972 } 973 default: 974 SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim); 975 } 976 PetscFunctionReturn(0); 977 } 978 979 /*@ 980 DMPlexCreateHexCylinderMesh - Creates a mesh on the tensor product of the unit interval with the circle (cylinder) using hexahedra. 981 982 Collective on MPI_Comm 983 984 Input Parameters: 985 + comm - The communicator for the DM object 986 - periodicZ - The boundary type for the Z direction 987 988 Output Parameter: 989 . dm - The DM object 990 991 Note: Here is the output numbering looking from the bottom of the cylinder: 992 $ 17-----14 993 $ | | 994 $ | 2 | 995 $ | | 996 $ 17-----8-----7-----14 997 $ | | | | 998 $ | 3 | 0 | 1 | 999 $ | | | | 1000 $ 19-----5-----6-----13 1001 $ | | 1002 $ | 4 | 1003 $ | | 1004 $ 19-----13 1005 $ 1006 $ and up through the top 1007 $ 1008 $ 18-----16 1009 $ | | 1010 $ | 2 | 1011 $ | | 1012 $ 18----10----11-----16 1013 $ | | | | 1014 $ | 3 | 0 | 1 | 1015 $ | | | | 1016 $ 20-----9----12-----15 1017 $ | | 1018 $ | 4 | 1019 $ | | 1020 $ 20-----15 1021 1022 Level: beginner 1023 1024 .keywords: DM, create 1025 .seealso: DMPlexCreateHexBoxMesh(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 1026 @*/ 1027 #define REAL_CYL 1028 PetscErrorCode DMPlexCreateHexCylinderMesh(MPI_Comm comm, DMBoundaryType periodicZ, DM *dm) 1029 { 1030 const PetscInt dim = 3; 1031 #ifdef REAL_CYL 1032 const PetscInt numCells = 5; 1033 const PetscInt numVertices = 16; 1034 #else 1035 const PetscInt numCells = 2; 1036 const PetscInt numVertices = 12; 1037 #endif 1038 PetscErrorCode ierr; 1039 1040 PetscFunctionBegin; 1041 PetscValidPointer(dm, 4); 1042 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1043 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1044 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1045 /* Create topology */ 1046 { 1047 PetscInt cone[8], c; 1048 1049 ierr = DMPlexSetChart(*dm, 0, numCells+numVertices);CHKERRQ(ierr); 1050 for (c = 0; c < numCells; c++) {ierr = DMPlexSetConeSize(*dm, c, 8);CHKERRQ(ierr);} 1051 ierr = DMSetUp(*dm);CHKERRQ(ierr); 1052 #ifndef REAL_CYL 1053 cone[0] = 2; cone[1] = 3; cone[2] = 4; cone[3] = 5; 1054 cone[4] = 6; cone[5] = 7; cone[6] = 8; cone[7] = 9; 1055 ierr = DMPlexSetCone(*dm, 0, cone);CHKERRQ(ierr); 1056 cone[0] = 3; cone[1] = 10; cone[2] = 11; cone[3] = 4; 1057 cone[4] = 9; cone[5] = 8; cone[6] = 13; cone[7] = 12; 1058 ierr = DMPlexSetCone(*dm, 1, cone);CHKERRQ(ierr); 1059 #else 1060 cone[0] = 5; cone[1] = 6; cone[2] = 7; cone[3] = 8; 1061 cone[4] = 9; cone[5] = 10; cone[6] = 11; cone[7] = 12; 1062 ierr = DMPlexSetCone(*dm, 0, cone);CHKERRQ(ierr); 1063 cone[0] = 6; cone[1] = 13; cone[2] = 14; cone[3] = 7; 1064 cone[4] = 12; cone[5] = 11; cone[6] = 16; cone[7] = 15; 1065 ierr = DMPlexSetCone(*dm, 1, cone);CHKERRQ(ierr); 1066 cone[0] = 8; cone[1] = 7; cone[2] = 14; cone[3] = 17; 1067 cone[4] = 10; cone[5] = 18; cone[6] = 16; cone[7] = 11; 1068 ierr = DMPlexSetCone(*dm, 2, cone);CHKERRQ(ierr); 1069 cone[0] = 19; cone[1] = 5; cone[2] = 8; cone[3] = 17; 1070 cone[4] = 20; cone[5] = 18; cone[6] = 10; cone[7] = 9; 1071 ierr = DMPlexSetCone(*dm, 3, cone);CHKERRQ(ierr); 1072 cone[0] = 19; cone[1] = 13; cone[2] = 6; cone[3] = 5; 1073 cone[4] = 20; cone[5] = 9; cone[6] = 12; cone[7] = 15; 1074 ierr = DMPlexSetCone(*dm, 4, cone);CHKERRQ(ierr); 1075 #endif 1076 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1077 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1078 } 1079 /* Interpolate */ 1080 { 1081 DM idm = NULL; 1082 1083 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1084 ierr = DMDestroy(dm);CHKERRQ(ierr); 1085 *dm = idm; 1086 } 1087 /* Create geometry */ 1088 { 1089 Vec coordinates; 1090 PetscSection coordSection; 1091 PetscScalar *coords; 1092 PetscInt coordSize, v; 1093 const PetscReal dis = 1.0/PetscSqrtReal(2.0); 1094 const PetscReal ds2 = dis/2.0; 1095 1096 /* Build coordinates */ 1097 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 1098 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1099 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 1100 ierr = PetscSectionSetChart(coordSection, numCells, numCells+numVertices);CHKERRQ(ierr); 1101 for (v = numCells; v < numCells+numVertices; ++v) { 1102 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 1103 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 1104 } 1105 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1106 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1107 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 1108 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1109 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1110 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 1111 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1112 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1113 coords[0*dim+0] = -ds2; coords[0*dim+1] = -ds2; coords[0*dim+2] = 0.0; 1114 coords[1*dim+0] = ds2; coords[1*dim+1] = -ds2; coords[1*dim+2] = 0.0; 1115 coords[2*dim+0] = ds2; coords[2*dim+1] = ds2; coords[2*dim+2] = 0.0; 1116 coords[3*dim+0] = -ds2; coords[3*dim+1] = ds2; coords[3*dim+2] = 0.0; 1117 coords[4*dim+0] = -ds2; coords[4*dim+1] = -ds2; coords[4*dim+2] = 1.0; 1118 coords[5*dim+0] = -ds2; coords[5*dim+1] = ds2; coords[5*dim+2] = 1.0; 1119 coords[6*dim+0] = ds2; coords[6*dim+1] = ds2; coords[6*dim+2] = 1.0; 1120 coords[7*dim+0] = ds2; coords[7*dim+1] = -ds2; coords[7*dim+2] = 1.0; 1121 coords[ 8*dim+0] = dis; coords[ 8*dim+1] = -dis; coords[ 8*dim+2] = 0.0; 1122 coords[ 9*dim+0] = dis; coords[ 9*dim+1] = dis; coords[ 9*dim+2] = 0.0; 1123 coords[10*dim+0] = dis; coords[10*dim+1] = -dis; coords[10*dim+2] = 1.0; 1124 coords[11*dim+0] = dis; coords[11*dim+1] = dis; coords[11*dim+2] = 1.0; 1125 #ifdef REAL_CYL 1126 coords[12*dim+0] = -dis; coords[12*dim+1] = dis; coords[12*dim+2] = 0.0; 1127 coords[13*dim+0] = -dis; coords[13*dim+1] = dis; coords[13*dim+2] = 1.0; 1128 coords[14*dim+0] = -dis; coords[14*dim+1] = -dis; coords[14*dim+2] = 0.0; 1129 coords[15*dim+0] = -dis; coords[15*dim+1] = -dis; coords[15*dim+2] = 1.0; 1130 #endif 1131 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1132 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 1133 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1134 } 1135 /* Create periodicity */ 1136 if (periodicZ == DM_BOUNDARY_PERIODIC || periodicZ == DM_BOUNDARY_TWIST) { 1137 PetscReal L[3]; 1138 PetscReal maxCell[3]; 1139 DMBoundaryType bdType[3]; 1140 PetscReal lower[3] = {0.0, 0.0, 0.0}; 1141 PetscReal upper[3] = {1.0, 1.0, 1.0}; 1142 PetscInt i, numZCells = 3; 1143 1144 bdType[0] = DM_BOUNDARY_NONE; 1145 bdType[1] = DM_BOUNDARY_NONE; 1146 bdType[2] = periodicZ; 1147 for (i = 0; i < dim; i++) { 1148 L[i] = upper[i] - lower[i]; 1149 maxCell[i] = 1.1 * (L[i] / numZCells); 1150 } 1151 ierr = DMSetPeriodicity(*dm, maxCell, L, bdType);CHKERRQ(ierr); 1152 } 1153 PetscFunctionReturn(0); 1154 } 1155 1156 /* External function declarations here */ 1157 extern PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling); 1158 extern PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, Mat *mat); 1159 extern PetscErrorCode DMCreateDefaultSection_Plex(DM dm); 1160 extern PetscErrorCode DMCreateDefaultConstraints_Plex(DM dm); 1161 extern PetscErrorCode DMCreateMatrix_Plex(DM dm, Mat *J); 1162 extern PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm); 1163 PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm); 1164 extern PetscErrorCode DMSetUp_Plex(DM dm); 1165 extern PetscErrorCode DMDestroy_Plex(DM dm); 1166 extern PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer); 1167 extern PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer); 1168 extern PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm); 1169 1170 /* Replace dm with the contents of dmNew 1171 - Share the DM_Plex structure 1172 - Share the coordinates 1173 - Share the SF 1174 */ 1175 static PetscErrorCode DMPlexReplace_Static(DM dm, DM dmNew) 1176 { 1177 PetscSF sf; 1178 DM coordDM, coarseDM; 1179 Vec coords; 1180 const PetscReal *maxCell, *L; 1181 const DMBoundaryType *bd; 1182 PetscErrorCode ierr; 1183 1184 PetscFunctionBegin; 1185 ierr = DMGetPointSF(dmNew, &sf);CHKERRQ(ierr); 1186 ierr = DMSetPointSF(dm, sf);CHKERRQ(ierr); 1187 ierr = DMGetCoordinateDM(dmNew, &coordDM);CHKERRQ(ierr); 1188 ierr = DMGetCoordinatesLocal(dmNew, &coords);CHKERRQ(ierr); 1189 ierr = DMSetCoordinateDM(dm, coordDM);CHKERRQ(ierr); 1190 ierr = DMSetCoordinatesLocal(dm, coords);CHKERRQ(ierr); 1191 ierr = DMGetPeriodicity(dm, &maxCell, &L, &bd);CHKERRQ(ierr); 1192 if (L) {ierr = DMSetPeriodicity(dmNew, maxCell, L, bd);CHKERRQ(ierr);} 1193 ierr = DMDestroy_Plex(dm);CHKERRQ(ierr); 1194 dm->data = dmNew->data; 1195 ((DM_Plex *) dmNew->data)->refct++; 1196 dmNew->labels->refct++; 1197 if (!--(dm->labels->refct)) { 1198 DMLabelLink next = dm->labels->next; 1199 1200 /* destroy the labels */ 1201 while (next) { 1202 DMLabelLink tmp = next->next; 1203 1204 ierr = DMLabelDestroy(&next->label);CHKERRQ(ierr); 1205 ierr = PetscFree(next);CHKERRQ(ierr); 1206 next = tmp; 1207 } 1208 ierr = PetscFree(dm->labels);CHKERRQ(ierr); 1209 } 1210 dm->labels = dmNew->labels; 1211 dm->depthLabel = dmNew->depthLabel; 1212 ierr = DMGetCoarseDM(dmNew,&coarseDM);CHKERRQ(ierr); 1213 ierr = DMSetCoarseDM(dm,coarseDM);CHKERRQ(ierr); 1214 PetscFunctionReturn(0); 1215 } 1216 1217 /* Swap dm with the contents of dmNew 1218 - Swap the DM_Plex structure 1219 - Swap the coordinates 1220 - Swap the point PetscSF 1221 */ 1222 static PetscErrorCode DMPlexSwap_Static(DM dmA, DM dmB) 1223 { 1224 DM coordDMA, coordDMB; 1225 Vec coordsA, coordsB; 1226 PetscSF sfA, sfB; 1227 void *tmp; 1228 DMLabelLinkList listTmp; 1229 DMLabel depthTmp; 1230 PetscErrorCode ierr; 1231 1232 PetscFunctionBegin; 1233 ierr = DMGetPointSF(dmA, &sfA);CHKERRQ(ierr); 1234 ierr = DMGetPointSF(dmB, &sfB);CHKERRQ(ierr); 1235 ierr = PetscObjectReference((PetscObject) sfA);CHKERRQ(ierr); 1236 ierr = DMSetPointSF(dmA, sfB);CHKERRQ(ierr); 1237 ierr = DMSetPointSF(dmB, sfA);CHKERRQ(ierr); 1238 ierr = PetscObjectDereference((PetscObject) sfA);CHKERRQ(ierr); 1239 1240 ierr = DMGetCoordinateDM(dmA, &coordDMA);CHKERRQ(ierr); 1241 ierr = DMGetCoordinateDM(dmB, &coordDMB);CHKERRQ(ierr); 1242 ierr = PetscObjectReference((PetscObject) coordDMA);CHKERRQ(ierr); 1243 ierr = DMSetCoordinateDM(dmA, coordDMB);CHKERRQ(ierr); 1244 ierr = DMSetCoordinateDM(dmB, coordDMA);CHKERRQ(ierr); 1245 ierr = PetscObjectDereference((PetscObject) coordDMA);CHKERRQ(ierr); 1246 1247 ierr = DMGetCoordinatesLocal(dmA, &coordsA);CHKERRQ(ierr); 1248 ierr = DMGetCoordinatesLocal(dmB, &coordsB);CHKERRQ(ierr); 1249 ierr = PetscObjectReference((PetscObject) coordsA);CHKERRQ(ierr); 1250 ierr = DMSetCoordinatesLocal(dmA, coordsB);CHKERRQ(ierr); 1251 ierr = DMSetCoordinatesLocal(dmB, coordsA);CHKERRQ(ierr); 1252 ierr = PetscObjectDereference((PetscObject) coordsA);CHKERRQ(ierr); 1253 1254 tmp = dmA->data; 1255 dmA->data = dmB->data; 1256 dmB->data = tmp; 1257 listTmp = dmA->labels; 1258 dmA->labels = dmB->labels; 1259 dmB->labels = listTmp; 1260 depthTmp = dmA->depthLabel; 1261 dmA->depthLabel = dmB->depthLabel; 1262 dmB->depthLabel = depthTmp; 1263 PetscFunctionReturn(0); 1264 } 1265 1266 PetscErrorCode DMSetFromOptions_NonRefinement_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 1267 { 1268 DM_Plex *mesh = (DM_Plex*) dm->data; 1269 PetscErrorCode ierr; 1270 1271 PetscFunctionBegin; 1272 /* Handle viewing */ 1273 ierr = PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMView", PETSC_FALSE, &mesh->printSetValues, NULL);CHKERRQ(ierr); 1274 ierr = PetscOptionsInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMView", 0, &mesh->printFEM, NULL);CHKERRQ(ierr); 1275 ierr = PetscOptionsReal("-dm_plex_print_tol", "Tolerance for FEM output", "DMView", mesh->printTol, &mesh->printTol, NULL);CHKERRQ(ierr); 1276 /* Point Location */ 1277 ierr = PetscOptionsBool("-dm_plex_hash_location", "Use grid hashing for point location", "DMView", PETSC_FALSE, &mesh->useHashLocation, NULL);CHKERRQ(ierr); 1278 /* Generation and remeshing */ 1279 ierr = PetscOptionsBool("-dm_plex_remesh_bd", "Allow changes to the boundary on remeshing", "DMView", PETSC_FALSE, &mesh->remeshBd, NULL);CHKERRQ(ierr); 1280 /* Projection behavior */ 1281 ierr = PetscOptionsInt("-dm_plex_max_projection_height", "Maxmimum mesh point height used to project locally", "DMPlexSetMaxProjectionHeight", 0, &mesh->maxProjectionHeight, NULL);CHKERRQ(ierr); 1282 ierr = PetscOptionsBool("-dm_plex_regular_refinement", "Use special nested projection algorithm for regular refinement", "DMPlexSetRegularRefinement", mesh->regularRefinement, &mesh->regularRefinement, NULL);CHKERRQ(ierr); 1283 PetscFunctionReturn(0); 1284 } 1285 1286 static PetscErrorCode DMSetFromOptions_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 1287 { 1288 PetscInt refine = 0, coarsen = 0, r; 1289 PetscBool isHierarchy; 1290 PetscErrorCode ierr; 1291 1292 PetscFunctionBegin; 1293 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1294 ierr = PetscOptionsHead(PetscOptionsObject,"DMPlex Options");CHKERRQ(ierr); 1295 /* Handle DMPlex refinement */ 1296 ierr = PetscOptionsInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL);CHKERRQ(ierr); 1297 ierr = PetscOptionsInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy);CHKERRQ(ierr); 1298 if (refine) {ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE);CHKERRQ(ierr);} 1299 if (refine && isHierarchy) { 1300 DM *dms, coarseDM; 1301 1302 ierr = DMGetCoarseDM(dm, &coarseDM);CHKERRQ(ierr); 1303 ierr = PetscObjectReference((PetscObject)coarseDM);CHKERRQ(ierr); 1304 ierr = PetscMalloc1(refine,&dms);CHKERRQ(ierr); 1305 ierr = DMRefineHierarchy(dm, refine, dms);CHKERRQ(ierr); 1306 /* Total hack since we do not pass in a pointer */ 1307 ierr = DMPlexSwap_Static(dm, dms[refine-1]);CHKERRQ(ierr); 1308 if (refine == 1) { 1309 ierr = DMSetCoarseDM(dm, dms[0]);CHKERRQ(ierr); 1310 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr); 1311 } else { 1312 ierr = DMSetCoarseDM(dm, dms[refine-2]);CHKERRQ(ierr); 1313 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr); 1314 ierr = DMSetCoarseDM(dms[0], dms[refine-1]);CHKERRQ(ierr); 1315 ierr = DMPlexSetRegularRefinement(dms[0], PETSC_TRUE);CHKERRQ(ierr); 1316 } 1317 ierr = DMSetCoarseDM(dms[refine-1], coarseDM);CHKERRQ(ierr); 1318 ierr = PetscObjectDereference((PetscObject)coarseDM);CHKERRQ(ierr); 1319 /* Free DMs */ 1320 for (r = 0; r < refine; ++r) { 1321 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr); 1322 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 1323 } 1324 ierr = PetscFree(dms);CHKERRQ(ierr); 1325 } else { 1326 for (r = 0; r < refine; ++r) { 1327 DM refinedMesh; 1328 1329 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 1330 ierr = DMRefine(dm, PetscObjectComm((PetscObject) dm), &refinedMesh);CHKERRQ(ierr); 1331 /* Total hack since we do not pass in a pointer */ 1332 ierr = DMPlexReplace_Static(dm, refinedMesh);CHKERRQ(ierr); 1333 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 1334 ierr = DMDestroy(&refinedMesh);CHKERRQ(ierr); 1335 } 1336 } 1337 /* Handle DMPlex coarsening */ 1338 ierr = PetscOptionsInt("-dm_coarsen", "Coarsen the mesh", "DMCreate", coarsen, &coarsen, NULL);CHKERRQ(ierr); 1339 ierr = PetscOptionsInt("-dm_coarsen_hierarchy", "The number of coarsenings", "DMCreate", coarsen, &coarsen, &isHierarchy);CHKERRQ(ierr); 1340 if (coarsen && isHierarchy) { 1341 DM *dms; 1342 1343 ierr = PetscMalloc1(coarsen, &dms);CHKERRQ(ierr); 1344 ierr = DMCoarsenHierarchy(dm, coarsen, dms);CHKERRQ(ierr); 1345 /* Free DMs */ 1346 for (r = 0; r < coarsen; ++r) { 1347 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr); 1348 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 1349 } 1350 ierr = PetscFree(dms);CHKERRQ(ierr); 1351 } else { 1352 for (r = 0; r < coarsen; ++r) { 1353 DM coarseMesh; 1354 1355 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 1356 ierr = DMCoarsen(dm, PetscObjectComm((PetscObject) dm), &coarseMesh);CHKERRQ(ierr); 1357 /* Total hack since we do not pass in a pointer */ 1358 ierr = DMPlexReplace_Static(dm, coarseMesh);CHKERRQ(ierr); 1359 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 1360 ierr = DMDestroy(&coarseMesh);CHKERRQ(ierr); 1361 } 1362 } 1363 /* Handle */ 1364 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 1365 ierr = PetscOptionsTail();CHKERRQ(ierr); 1366 PetscFunctionReturn(0); 1367 } 1368 1369 static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec) 1370 { 1371 PetscErrorCode ierr; 1372 1373 PetscFunctionBegin; 1374 ierr = DMCreateGlobalVector_Section_Private(dm,vec);CHKERRQ(ierr); 1375 /* ierr = VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM);CHKERRQ(ierr); */ 1376 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex);CHKERRQ(ierr); 1377 ierr = VecSetOperation(*vec, VECOP_VIEWNATIVE, (void (*)(void)) VecView_Plex_Native);CHKERRQ(ierr); 1378 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex);CHKERRQ(ierr); 1379 ierr = VecSetOperation(*vec, VECOP_LOADNATIVE, (void (*)(void)) VecLoad_Plex_Native);CHKERRQ(ierr); 1380 PetscFunctionReturn(0); 1381 } 1382 1383 static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec) 1384 { 1385 PetscErrorCode ierr; 1386 1387 PetscFunctionBegin; 1388 ierr = DMCreateLocalVector_Section_Private(dm,vec);CHKERRQ(ierr); 1389 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex_Local);CHKERRQ(ierr); 1390 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex_Local);CHKERRQ(ierr); 1391 PetscFunctionReturn(0); 1392 } 1393 1394 static PetscErrorCode DMGetDimPoints_Plex(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd) 1395 { 1396 PetscInt depth, d; 1397 PetscErrorCode ierr; 1398 1399 PetscFunctionBegin; 1400 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1401 if (depth == 1) { 1402 ierr = DMGetDimension(dm, &d);CHKERRQ(ierr); 1403 if (dim == 0) {ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr);} 1404 else if (dim == d) {ierr = DMPlexGetDepthStratum(dm, 1, pStart, pEnd);CHKERRQ(ierr);} 1405 else {*pStart = 0; *pEnd = 0;} 1406 } else { 1407 ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr); 1408 } 1409 PetscFunctionReturn(0); 1410 } 1411 1412 static PetscErrorCode DMGetNeighors_Plex(DM dm, PetscInt *nranks, const PetscMPIInt *ranks[]) 1413 { 1414 PetscSF sf; 1415 PetscErrorCode ierr; 1416 1417 PetscFunctionBegin; 1418 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 1419 ierr = PetscSFGetRanks(sf, nranks, ranks, NULL, NULL, NULL);CHKERRQ(ierr); 1420 PetscFunctionReturn(0); 1421 } 1422 1423 static PetscErrorCode DMInitialize_Plex(DM dm) 1424 { 1425 PetscErrorCode ierr; 1426 1427 PetscFunctionBegin; 1428 dm->ops->view = DMView_Plex; 1429 dm->ops->load = DMLoad_Plex; 1430 dm->ops->setfromoptions = DMSetFromOptions_Plex; 1431 dm->ops->clone = DMClone_Plex; 1432 dm->ops->setup = DMSetUp_Plex; 1433 dm->ops->createdefaultsection = DMCreateDefaultSection_Plex; 1434 dm->ops->createdefaultconstraints = DMCreateDefaultConstraints_Plex; 1435 dm->ops->createglobalvector = DMCreateGlobalVector_Plex; 1436 dm->ops->createlocalvector = DMCreateLocalVector_Plex; 1437 dm->ops->getlocaltoglobalmapping = NULL; 1438 dm->ops->createfieldis = NULL; 1439 dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex; 1440 dm->ops->getcoloring = NULL; 1441 dm->ops->creatematrix = DMCreateMatrix_Plex; 1442 dm->ops->createinterpolation = DMCreateInterpolation_Plex; 1443 dm->ops->getaggregates = NULL; 1444 dm->ops->getinjection = DMCreateInjection_Plex; 1445 dm->ops->refine = DMRefine_Plex; 1446 dm->ops->coarsen = DMCoarsen_Plex; 1447 dm->ops->refinehierarchy = DMRefineHierarchy_Plex; 1448 dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Plex; 1449 dm->ops->globaltolocalbegin = NULL; 1450 dm->ops->globaltolocalend = NULL; 1451 dm->ops->localtoglobalbegin = NULL; 1452 dm->ops->localtoglobalend = NULL; 1453 dm->ops->destroy = DMDestroy_Plex; 1454 dm->ops->createsubdm = DMCreateSubDM_Plex; 1455 dm->ops->getdimpoints = DMGetDimPoints_Plex; 1456 dm->ops->locatepoints = DMLocatePoints_Plex; 1457 dm->ops->projectfunctionlocal = DMProjectFunctionLocal_Plex; 1458 dm->ops->projectfunctionlabellocal = DMProjectFunctionLabelLocal_Plex; 1459 dm->ops->projectfieldlocal = DMProjectFieldLocal_Plex; 1460 dm->ops->projectfieldlabellocal = DMProjectFieldLabelLocal_Plex; 1461 dm->ops->computel2diff = DMComputeL2Diff_Plex; 1462 dm->ops->computel2gradientdiff = DMComputeL2GradientDiff_Plex; 1463 dm->ops->computel2fielddiff = DMComputeL2FieldDiff_Plex; 1464 dm->ops->getneighbors = DMGetNeighors_Plex; 1465 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMAdaptLabel_C",DMAdaptLabel_Plex);CHKERRQ(ierr); 1466 PetscFunctionReturn(0); 1467 } 1468 1469 PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm) 1470 { 1471 DM_Plex *mesh = (DM_Plex *) dm->data; 1472 PetscErrorCode ierr; 1473 1474 PetscFunctionBegin; 1475 mesh->refct++; 1476 (*newdm)->data = mesh; 1477 ierr = PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);CHKERRQ(ierr); 1478 ierr = DMInitialize_Plex(*newdm);CHKERRQ(ierr); 1479 PetscFunctionReturn(0); 1480 } 1481 1482 /*MC 1483 DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram. 1484 In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is 1485 specified by a PetscSection object. Ownership in the global representation is determined by 1486 ownership of the underlying DMPlex points. This is specified by another PetscSection object. 1487 1488 Level: intermediate 1489 1490 .seealso: DMType, DMPlexCreate(), DMCreate(), DMSetType() 1491 M*/ 1492 1493 PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm) 1494 { 1495 DM_Plex *mesh; 1496 PetscInt unit, d; 1497 PetscErrorCode ierr; 1498 1499 PetscFunctionBegin; 1500 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1501 ierr = PetscNewLog(dm,&mesh);CHKERRQ(ierr); 1502 dm->dim = 0; 1503 dm->data = mesh; 1504 1505 mesh->refct = 1; 1506 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);CHKERRQ(ierr); 1507 mesh->maxConeSize = 0; 1508 mesh->cones = NULL; 1509 mesh->coneOrientations = NULL; 1510 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);CHKERRQ(ierr); 1511 mesh->maxSupportSize = 0; 1512 mesh->supports = NULL; 1513 mesh->refinementUniform = PETSC_TRUE; 1514 mesh->refinementLimit = -1.0; 1515 1516 mesh->facesTmp = NULL; 1517 1518 mesh->tetgenOpts = NULL; 1519 mesh->triangleOpts = NULL; 1520 ierr = PetscPartitionerCreate(PetscObjectComm((PetscObject)dm), &mesh->partitioner);CHKERRQ(ierr); 1521 ierr = PetscPartitionerSetTypeFromOptions_Internal(mesh->partitioner);CHKERRQ(ierr); 1522 mesh->remeshBd = PETSC_FALSE; 1523 1524 mesh->subpointMap = NULL; 1525 1526 for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0; 1527 1528 mesh->regularRefinement = PETSC_FALSE; 1529 mesh->depthState = -1; 1530 mesh->globalVertexNumbers = NULL; 1531 mesh->globalCellNumbers = NULL; 1532 mesh->anchorSection = NULL; 1533 mesh->anchorIS = NULL; 1534 mesh->createanchors = NULL; 1535 mesh->computeanchormatrix = NULL; 1536 mesh->parentSection = NULL; 1537 mesh->parents = NULL; 1538 mesh->childIDs = NULL; 1539 mesh->childSection = NULL; 1540 mesh->children = NULL; 1541 mesh->referenceTree = NULL; 1542 mesh->getchildsymmetry = NULL; 1543 for (d = 0; d < 8; ++d) mesh->hybridPointMax[d] = PETSC_DETERMINE; 1544 mesh->vtkCellHeight = 0; 1545 mesh->useCone = PETSC_FALSE; 1546 mesh->useClosure = PETSC_TRUE; 1547 mesh->useAnchors = PETSC_FALSE; 1548 1549 mesh->maxProjectionHeight = 0; 1550 1551 mesh->printSetValues = PETSC_FALSE; 1552 mesh->printFEM = 0; 1553 mesh->printTol = 1.0e-10; 1554 1555 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 1556 PetscFunctionReturn(0); 1557 } 1558 1559 /*@ 1560 DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram. 1561 1562 Collective on MPI_Comm 1563 1564 Input Parameter: 1565 . comm - The communicator for the DMPlex object 1566 1567 Output Parameter: 1568 . mesh - The DMPlex object 1569 1570 Level: beginner 1571 1572 .keywords: DMPlex, create 1573 @*/ 1574 PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh) 1575 { 1576 PetscErrorCode ierr; 1577 1578 PetscFunctionBegin; 1579 PetscValidPointer(mesh,2); 1580 ierr = DMCreate(comm, mesh);CHKERRQ(ierr); 1581 ierr = DMSetType(*mesh, DMPLEX);CHKERRQ(ierr); 1582 PetscFunctionReturn(0); 1583 } 1584 1585 /* 1586 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 1587 */ 1588 static PetscErrorCode DMPlexBuildFromCellList_Parallel_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[], PetscSF *sfVert) 1589 { 1590 PetscSF sfPoint; 1591 PetscLayout vLayout; 1592 PetscHashI vhash; 1593 PetscSFNode *remoteVerticesAdj, *vertexLocal, *vertexOwner, *remoteVertex; 1594 const PetscInt *vrange; 1595 PetscInt numVerticesAdj, off, *verticesAdj, numVerticesGhost = 0, *localVertex, *cone, c, p, v, g; 1596 PETSC_UNUSED PetscHashIIter ret, iter; 1597 PetscMPIInt rank, size; 1598 PetscErrorCode ierr; 1599 1600 PetscFunctionBegin; 1601 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr); 1602 ierr = MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);CHKERRQ(ierr); 1603 /* Partition vertices */ 1604 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) dm), &vLayout);CHKERRQ(ierr); 1605 ierr = PetscLayoutSetLocalSize(vLayout, numVertices);CHKERRQ(ierr); 1606 ierr = PetscLayoutSetBlockSize(vLayout, 1);CHKERRQ(ierr); 1607 ierr = PetscLayoutSetUp(vLayout);CHKERRQ(ierr); 1608 ierr = PetscLayoutGetRanges(vLayout, &vrange);CHKERRQ(ierr); 1609 /* Count vertices and map them to procs */ 1610 PetscHashICreate(vhash); 1611 for (c = 0; c < numCells; ++c) { 1612 for (p = 0; p < numCorners; ++p) { 1613 PetscHashIPut(vhash, cells[c*numCorners+p], ret, iter); 1614 } 1615 } 1616 PetscHashISize(vhash, numVerticesAdj); 1617 ierr = PetscMalloc1(numVerticesAdj, &verticesAdj);CHKERRQ(ierr); 1618 off = 0; ierr = PetscHashIGetKeys(vhash, &off, verticesAdj);CHKERRQ(ierr); 1619 if (off != numVerticesAdj) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid number of local vertices %D should be %D", off, numVerticesAdj); 1620 ierr = PetscSortInt(numVerticesAdj, verticesAdj);CHKERRQ(ierr); 1621 ierr = PetscMalloc1(numVerticesAdj, &remoteVerticesAdj);CHKERRQ(ierr); 1622 for (v = 0; v < numVerticesAdj; ++v) { 1623 const PetscInt gv = verticesAdj[v]; 1624 PetscInt vrank; 1625 1626 ierr = PetscFindInt(gv, size+1, vrange, &vrank);CHKERRQ(ierr); 1627 vrank = vrank < 0 ? -(vrank+2) : vrank; 1628 remoteVerticesAdj[v].index = gv - vrange[vrank]; 1629 remoteVerticesAdj[v].rank = vrank; 1630 } 1631 /* Create cones */ 1632 ierr = DMPlexSetChart(dm, 0, numCells+numVerticesAdj);CHKERRQ(ierr); 1633 for (c = 0; c < numCells; ++c) {ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr);} 1634 ierr = DMSetUp(dm);CHKERRQ(ierr); 1635 ierr = DMGetWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 1636 for (c = 0; c < numCells; ++c) { 1637 for (p = 0; p < numCorners; ++p) { 1638 const PetscInt gv = cells[c*numCorners+p]; 1639 PetscInt lv; 1640 1641 ierr = PetscFindInt(gv, numVerticesAdj, verticesAdj, &lv);CHKERRQ(ierr); 1642 if (lv < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not find global vertex %D in local connectivity", gv); 1643 cone[p] = lv+numCells; 1644 } 1645 ierr = DMPlexSetCone(dm, c, cone);CHKERRQ(ierr); 1646 } 1647 ierr = DMRestoreWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 1648 /* Create SF for vertices */ 1649 ierr = PetscSFCreate(PetscObjectComm((PetscObject)dm), sfVert);CHKERRQ(ierr); 1650 ierr = PetscObjectSetName((PetscObject) *sfVert, "Vertex Ownership SF");CHKERRQ(ierr); 1651 ierr = PetscSFSetFromOptions(*sfVert);CHKERRQ(ierr); 1652 ierr = PetscSFSetGraph(*sfVert, numVertices, numVerticesAdj, NULL, PETSC_OWN_POINTER, remoteVerticesAdj, PETSC_OWN_POINTER);CHKERRQ(ierr); 1653 ierr = PetscFree(verticesAdj);CHKERRQ(ierr); 1654 /* Build pointSF */ 1655 ierr = PetscMalloc2(numVerticesAdj, &vertexLocal, numVertices, &vertexOwner);CHKERRQ(ierr); 1656 for (v = 0; v < numVerticesAdj; ++v) {vertexLocal[v].index = v+numCells; vertexLocal[v].rank = rank;} 1657 for (v = 0; v < numVertices; ++v) {vertexOwner[v].index = -1; vertexOwner[v].rank = -1;} 1658 ierr = PetscSFReduceBegin(*sfVert, MPI_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);CHKERRQ(ierr); 1659 ierr = PetscSFReduceEnd(*sfVert, MPI_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);CHKERRQ(ierr); 1660 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); 1661 ierr = PetscSFBcastBegin(*sfVert, MPI_2INT, vertexOwner, vertexLocal);CHKERRQ(ierr); 1662 ierr = PetscSFBcastEnd(*sfVert, MPI_2INT, vertexOwner, vertexLocal);CHKERRQ(ierr); 1663 for (v = 0; v < numVerticesAdj; ++v) if (vertexLocal[v].rank != rank) ++numVerticesGhost; 1664 ierr = PetscMalloc1(numVerticesGhost, &localVertex);CHKERRQ(ierr); 1665 ierr = PetscMalloc1(numVerticesGhost, &remoteVertex);CHKERRQ(ierr); 1666 for (v = 0, g = 0; v < numVerticesAdj; ++v) { 1667 if (vertexLocal[v].rank != rank) { 1668 localVertex[g] = v+numCells; 1669 remoteVertex[g].index = vertexLocal[v].index; 1670 remoteVertex[g].rank = vertexLocal[v].rank; 1671 ++g; 1672 } 1673 } 1674 ierr = PetscFree2(vertexLocal, vertexOwner);CHKERRQ(ierr); 1675 if (g != numVerticesGhost) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of vertex ghosts %D should be %D", g, numVerticesGhost); 1676 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 1677 ierr = PetscObjectSetName((PetscObject) sfPoint, "point SF");CHKERRQ(ierr); 1678 ierr = PetscSFSetGraph(sfPoint, numCells+numVerticesAdj, numVerticesGhost, localVertex, PETSC_OWN_POINTER, remoteVertex, PETSC_OWN_POINTER);CHKERRQ(ierr); 1679 ierr = PetscLayoutDestroy(&vLayout);CHKERRQ(ierr); 1680 PetscHashIDestroy(vhash); 1681 /* Fill in the rest of the topology structure */ 1682 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 1683 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 1684 PetscFunctionReturn(0); 1685 } 1686 1687 /* 1688 This takes as input the coordinates for each owned vertex 1689 */ 1690 static PetscErrorCode DMPlexBuildCoordinates_Parallel_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numV, PetscSF sfVert, const PetscReal vertexCoords[]) 1691 { 1692 PetscSection coordSection; 1693 Vec coordinates; 1694 PetscScalar *coords; 1695 PetscInt numVertices, numVerticesAdj, coordSize, v; 1696 PetscErrorCode ierr; 1697 1698 PetscFunctionBegin; 1699 ierr = PetscSFGetGraph(sfVert, &numVertices, &numVerticesAdj, NULL, NULL);CHKERRQ(ierr); 1700 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1701 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1702 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 1703 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVerticesAdj);CHKERRQ(ierr); 1704 for (v = numCells; v < numCells+numVerticesAdj; ++v) { 1705 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 1706 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 1707 } 1708 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1709 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1710 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 1711 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 1712 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1713 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1714 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1715 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1716 { 1717 MPI_Datatype coordtype; 1718 1719 /* Need a temp buffer for coords if we have complex/single */ 1720 ierr = MPI_Type_contiguous(spaceDim, MPIU_SCALAR, &coordtype);CHKERRQ(ierr); 1721 ierr = MPI_Type_commit(&coordtype);CHKERRQ(ierr); 1722 #if defined(PETSC_USE_COMPLEX) 1723 { 1724 PetscScalar *svertexCoords; 1725 PetscInt i; 1726 ierr = PetscMalloc1(numV*spaceDim,&svertexCoords);CHKERRQ(ierr); 1727 for (i=0; i<numV*spaceDim; i++) svertexCoords[i] = vertexCoords[i]; 1728 ierr = PetscSFBcastBegin(sfVert, coordtype, svertexCoords, coords);CHKERRQ(ierr); 1729 ierr = PetscSFBcastEnd(sfVert, coordtype, svertexCoords, coords);CHKERRQ(ierr); 1730 ierr = PetscFree(svertexCoords);CHKERRQ(ierr); 1731 } 1732 #else 1733 ierr = PetscSFBcastBegin(sfVert, coordtype, vertexCoords, coords);CHKERRQ(ierr); 1734 ierr = PetscSFBcastEnd(sfVert, coordtype, vertexCoords, coords);CHKERRQ(ierr); 1735 #endif 1736 ierr = MPI_Type_free(&coordtype);CHKERRQ(ierr); 1737 } 1738 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1739 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 1740 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1741 PetscFunctionReturn(0); 1742 } 1743 1744 /*@C 1745 DMPlexCreateFromCellListParallel - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM 1746 1747 Input Parameters: 1748 + comm - The communicator 1749 . dim - The topological dimension of the mesh 1750 . numCells - The number of cells owned by this process 1751 . numVertices - The number of vertices owned by this process 1752 . numCorners - The number of vertices for each cell 1753 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 1754 . cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 1755 . spaceDim - The spatial dimension used for coordinates 1756 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 1757 1758 Output Parameter: 1759 + dm - The DM 1760 - vertexSF - Optional, SF describing complete vertex ownership 1761 1762 Note: Two triangles sharing a face 1763 $ 1764 $ 2 1765 $ / | \ 1766 $ / | \ 1767 $ / | \ 1768 $ 0 0 | 1 3 1769 $ \ | / 1770 $ \ | / 1771 $ \ | / 1772 $ 1 1773 would have input 1774 $ numCells = 2, numVertices = 4 1775 $ cells = [0 1 2 1 3 2] 1776 $ 1777 which would result in the DMPlex 1778 $ 1779 $ 4 1780 $ / | \ 1781 $ / | \ 1782 $ / | \ 1783 $ 2 0 | 1 5 1784 $ \ | / 1785 $ \ | / 1786 $ \ | / 1787 $ 3 1788 1789 Level: beginner 1790 1791 .seealso: DMPlexCreateFromCellList(), DMPlexCreateFromDAG(), DMPlexCreate() 1792 @*/ 1793 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) 1794 { 1795 PetscSF sfVert; 1796 PetscErrorCode ierr; 1797 1798 PetscFunctionBegin; 1799 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1800 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1801 PetscValidLogicalCollectiveInt(*dm, dim, 2); 1802 PetscValidLogicalCollectiveInt(*dm, spaceDim, 8); 1803 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1804 ierr = DMPlexBuildFromCellList_Parallel_Private(*dm, numCells, numVertices, numCorners, cells, &sfVert);CHKERRQ(ierr); 1805 if (interpolate) { 1806 DM idm = NULL; 1807 1808 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1809 ierr = DMDestroy(dm);CHKERRQ(ierr); 1810 *dm = idm; 1811 } 1812 ierr = DMPlexBuildCoordinates_Parallel_Private(*dm, spaceDim, numCells, numVertices,sfVert, vertexCoords);CHKERRQ(ierr); 1813 if (vertexSF) *vertexSF = sfVert; 1814 else {ierr = PetscSFDestroy(&sfVert);CHKERRQ(ierr);} 1815 PetscFunctionReturn(0); 1816 } 1817 1818 /* 1819 This takes as input the common mesh generator output, a list of the vertices for each cell 1820 */ 1821 static PetscErrorCode DMPlexBuildFromCellList_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[]) 1822 { 1823 PetscInt *cone, c, p; 1824 PetscErrorCode ierr; 1825 1826 PetscFunctionBegin; 1827 ierr = DMPlexSetChart(dm, 0, numCells+numVertices);CHKERRQ(ierr); 1828 for (c = 0; c < numCells; ++c) { 1829 ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr); 1830 } 1831 ierr = DMSetUp(dm);CHKERRQ(ierr); 1832 ierr = DMGetWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 1833 for (c = 0; c < numCells; ++c) { 1834 for (p = 0; p < numCorners; ++p) { 1835 cone[p] = cells[c*numCorners+p]+numCells; 1836 } 1837 ierr = DMPlexSetCone(dm, c, cone);CHKERRQ(ierr); 1838 } 1839 ierr = DMRestoreWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 1840 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 1841 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 1842 PetscFunctionReturn(0); 1843 } 1844 1845 /* 1846 This takes as input the coordinates for each vertex 1847 */ 1848 static PetscErrorCode DMPlexBuildCoordinates_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, const double vertexCoords[]) 1849 { 1850 PetscSection coordSection; 1851 Vec coordinates; 1852 DM cdm; 1853 PetscScalar *coords; 1854 PetscInt v, d; 1855 PetscErrorCode ierr; 1856 1857 PetscFunctionBegin; 1858 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1859 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1860 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 1861 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVertices);CHKERRQ(ierr); 1862 for (v = numCells; v < numCells+numVertices; ++v) { 1863 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 1864 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 1865 } 1866 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1867 1868 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 1869 ierr = DMCreateLocalVector(cdm, &coordinates);CHKERRQ(ierr); 1870 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 1871 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1872 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1873 for (v = 0; v < numVertices; ++v) { 1874 for (d = 0; d < spaceDim; ++d) { 1875 coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d]; 1876 } 1877 } 1878 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1879 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 1880 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1881 PetscFunctionReturn(0); 1882 } 1883 1884 /*@C 1885 DMPlexCreateFromCellList - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM 1886 1887 Input Parameters: 1888 + comm - The communicator 1889 . dim - The topological dimension of the mesh 1890 . numCells - The number of cells 1891 . numVertices - The number of vertices 1892 . numCorners - The number of vertices for each cell 1893 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 1894 . cells - An array of numCells*numCorners numbers, the vertices for each cell 1895 . spaceDim - The spatial dimension used for coordinates 1896 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 1897 1898 Output Parameter: 1899 . dm - The DM 1900 1901 Note: Two triangles sharing a face 1902 $ 1903 $ 2 1904 $ / | \ 1905 $ / | \ 1906 $ / | \ 1907 $ 0 0 | 1 3 1908 $ \ | / 1909 $ \ | / 1910 $ \ | / 1911 $ 1 1912 would have input 1913 $ numCells = 2, numVertices = 4 1914 $ cells = [0 1 2 1 3 2] 1915 $ 1916 which would result in the DMPlex 1917 $ 1918 $ 4 1919 $ / | \ 1920 $ / | \ 1921 $ / | \ 1922 $ 2 0 | 1 5 1923 $ \ | / 1924 $ \ | / 1925 $ \ | / 1926 $ 3 1927 1928 Level: beginner 1929 1930 .seealso: DMPlexCreateFromDAG(), DMPlexCreate() 1931 @*/ 1932 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) 1933 { 1934 PetscErrorCode ierr; 1935 1936 PetscFunctionBegin; 1937 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1938 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1939 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1940 ierr = DMPlexBuildFromCellList_Private(*dm, numCells, numVertices, numCorners, cells);CHKERRQ(ierr); 1941 if (interpolate) { 1942 DM idm = NULL; 1943 1944 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1945 ierr = DMDestroy(dm);CHKERRQ(ierr); 1946 *dm = idm; 1947 } 1948 ierr = DMPlexBuildCoordinates_Private(*dm, spaceDim, numCells, numVertices, vertexCoords);CHKERRQ(ierr); 1949 PetscFunctionReturn(0); 1950 } 1951 1952 /*@ 1953 DMPlexCreateFromDAG - This takes as input the adjacency-list representation of the Directed Acyclic Graph (Hasse Diagram) encoding a mesh, and produces a DM 1954 1955 Input Parameters: 1956 + dm - The empty DM object, usually from DMCreate() and DMSetDimension() 1957 . depth - The depth of the DAG 1958 . numPoints - The number of points at each depth 1959 . coneSize - The cone size of each point 1960 . cones - The concatenation of the cone points for each point, the cone list must be oriented correctly for each point 1961 . coneOrientations - The orientation of each cone point 1962 - vertexCoords - An array of numVertices*dim numbers, the coordinates of each vertex 1963 1964 Output Parameter: 1965 . dm - The DM 1966 1967 Note: Two triangles sharing a face would have input 1968 $ depth = 1, numPoints = [4 2], coneSize = [3 3 0 0 0 0] 1969 $ cones = [2 3 4 3 5 4], coneOrientations = [0 0 0 0 0 0] 1970 $ vertexCoords = [-1.0 0.0 0.0 -1.0 0.0 1.0 1.0 0.0] 1971 $ 1972 which would result in the DMPlex 1973 $ 1974 $ 4 1975 $ / | \ 1976 $ / | \ 1977 $ / | \ 1978 $ 2 0 | 1 5 1979 $ \ | / 1980 $ \ | / 1981 $ \ | / 1982 $ 3 1983 $ 1984 $ Notice that all points are numbered consecutively, unlikely DMPlexCreateFromCellList() 1985 1986 Level: advanced 1987 1988 .seealso: DMPlexCreateFromCellList(), DMPlexCreate() 1989 @*/ 1990 PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[]) 1991 { 1992 Vec coordinates; 1993 PetscSection coordSection; 1994 PetscScalar *coords; 1995 PetscInt coordSize, firstVertex = -1, pStart = 0, pEnd = 0, p, v, dim, dimEmbed, d, off; 1996 PetscErrorCode ierr; 1997 1998 PetscFunctionBegin; 1999 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2000 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 2001 if (dimEmbed < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Embedding dimension %d cannot be less than intrinsic dimension %d",dimEmbed,dim); 2002 for (d = 0; d <= depth; ++d) pEnd += numPoints[d]; 2003 ierr = DMPlexSetChart(dm, pStart, pEnd);CHKERRQ(ierr); 2004 for (p = pStart; p < pEnd; ++p) { 2005 ierr = DMPlexSetConeSize(dm, p, coneSize[p-pStart]);CHKERRQ(ierr); 2006 if (firstVertex < 0 && !coneSize[p - pStart]) { 2007 firstVertex = p - pStart; 2008 } 2009 } 2010 if (firstVertex < 0 && numPoints[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected %d vertices but could not find any", numPoints[0]); 2011 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 2012 for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) { 2013 ierr = DMPlexSetCone(dm, p, &cones[off]);CHKERRQ(ierr); 2014 ierr = DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);CHKERRQ(ierr); 2015 } 2016 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 2017 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 2018 /* Build coordinates */ 2019 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2020 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 2021 ierr = PetscSectionSetFieldComponents(coordSection, 0, dimEmbed);CHKERRQ(ierr); 2022 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0]);CHKERRQ(ierr); 2023 for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) { 2024 ierr = PetscSectionSetDof(coordSection, v, dimEmbed);CHKERRQ(ierr); 2025 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed);CHKERRQ(ierr); 2026 } 2027 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 2028 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 2029 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 2030 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 2031 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 2032 ierr = VecSetBlockSize(coordinates, dimEmbed);CHKERRQ(ierr); 2033 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 2034 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2035 for (v = 0; v < numPoints[0]; ++v) { 2036 PetscInt off; 2037 2038 ierr = PetscSectionGetOffset(coordSection, v+firstVertex, &off);CHKERRQ(ierr); 2039 for (d = 0; d < dimEmbed; ++d) { 2040 coords[off+d] = vertexCoords[v*dimEmbed+d]; 2041 } 2042 } 2043 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2044 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 2045 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 2046 PetscFunctionReturn(0); 2047 } 2048 2049 /*@C 2050 DMPlexCreateFromFile - This takes a filename and produces a DM 2051 2052 Input Parameters: 2053 + comm - The communicator 2054 . filename - A file name 2055 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 2056 2057 Output Parameter: 2058 . dm - The DM 2059 2060 Level: beginner 2061 2062 .seealso: DMPlexCreateFromDAG(), DMPlexCreateFromCellList(), DMPlexCreate() 2063 @*/ 2064 PetscErrorCode DMPlexCreateFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 2065 { 2066 const char *extGmsh = ".msh"; 2067 const char *extCGNS = ".cgns"; 2068 const char *extExodus = ".exo"; 2069 const char *extFluent = ".cas"; 2070 const char *extHDF5 = ".h5"; 2071 const char *extMed = ".med"; 2072 size_t len; 2073 PetscBool isGmsh, isCGNS, isExodus, isFluent, isHDF5, isMed; 2074 PetscMPIInt rank; 2075 PetscErrorCode ierr; 2076 2077 PetscFunctionBegin; 2078 PetscValidPointer(filename, 2); 2079 PetscValidPointer(dm, 4); 2080 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 2081 ierr = PetscStrlen(filename, &len);CHKERRQ(ierr); 2082 if (!len) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Filename must be a valid path"); 2083 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGmsh, 4, &isGmsh);CHKERRQ(ierr); 2084 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extCGNS, 5, &isCGNS);CHKERRQ(ierr); 2085 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extExodus, 4, &isExodus);CHKERRQ(ierr); 2086 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extFluent, 4, &isFluent);CHKERRQ(ierr); 2087 ierr = PetscStrncmp(&filename[PetscMax(0,len-3)], extHDF5, 3, &isHDF5);CHKERRQ(ierr); 2088 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extMed, 4, &isMed);CHKERRQ(ierr); 2089 if (isGmsh) { 2090 ierr = DMPlexCreateGmshFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2091 } else if (isCGNS) { 2092 ierr = DMPlexCreateCGNSFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2093 } else if (isExodus) { 2094 ierr = DMPlexCreateExodusFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2095 } else if (isFluent) { 2096 ierr = DMPlexCreateFluentFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2097 } else if (isHDF5) { 2098 PetscViewer viewer; 2099 2100 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 2101 ierr = PetscViewerSetType(viewer, PETSCVIEWERHDF5);CHKERRQ(ierr); 2102 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 2103 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 2104 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 2105 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 2106 ierr = DMLoad(*dm, viewer);CHKERRQ(ierr); 2107 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 2108 } else if (isMed) { 2109 ierr = DMPlexCreateMedFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2110 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot load file %s: unrecognized extension", filename); 2111 PetscFunctionReturn(0); 2112 } 2113 2114 /*@ 2115 DMPlexCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell 2116 2117 Collective on comm 2118 2119 Input Parameters: 2120 + comm - The communicator 2121 . dim - The spatial dimension 2122 - simplex - Flag for simplex, otherwise use a tensor-product cell 2123 2124 Output Parameter: 2125 . refdm - The reference cell 2126 2127 Level: intermediate 2128 2129 .keywords: reference cell 2130 .seealso: 2131 @*/ 2132 PetscErrorCode DMPlexCreateReferenceCell(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *refdm) 2133 { 2134 DM rdm; 2135 Vec coords; 2136 PetscErrorCode ierr; 2137 2138 PetscFunctionBeginUser; 2139 ierr = DMCreate(comm, &rdm);CHKERRQ(ierr); 2140 ierr = DMSetType(rdm, DMPLEX);CHKERRQ(ierr); 2141 ierr = DMSetDimension(rdm, dim);CHKERRQ(ierr); 2142 switch (dim) { 2143 case 0: 2144 { 2145 PetscInt numPoints[1] = {1}; 2146 PetscInt coneSize[1] = {0}; 2147 PetscInt cones[1] = {0}; 2148 PetscInt coneOrientations[1] = {0}; 2149 PetscScalar vertexCoords[1] = {0.0}; 2150 2151 ierr = DMPlexCreateFromDAG(rdm, 0, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2152 } 2153 break; 2154 case 1: 2155 { 2156 PetscInt numPoints[2] = {2, 1}; 2157 PetscInt coneSize[3] = {2, 0, 0}; 2158 PetscInt cones[2] = {1, 2}; 2159 PetscInt coneOrientations[2] = {0, 0}; 2160 PetscScalar vertexCoords[2] = {-1.0, 1.0}; 2161 2162 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2163 } 2164 break; 2165 case 2: 2166 if (simplex) { 2167 PetscInt numPoints[2] = {3, 1}; 2168 PetscInt coneSize[4] = {3, 0, 0, 0}; 2169 PetscInt cones[3] = {1, 2, 3}; 2170 PetscInt coneOrientations[3] = {0, 0, 0}; 2171 PetscScalar vertexCoords[6] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0}; 2172 2173 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2174 } else { 2175 PetscInt numPoints[2] = {4, 1}; 2176 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 2177 PetscInt cones[4] = {1, 2, 3, 4}; 2178 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 2179 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0}; 2180 2181 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2182 } 2183 break; 2184 case 3: 2185 if (simplex) { 2186 PetscInt numPoints[2] = {4, 1}; 2187 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 2188 PetscInt cones[4] = {1, 3, 2, 4}; 2189 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 2190 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}; 2191 2192 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2193 } else { 2194 PetscInt numPoints[2] = {8, 1}; 2195 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0}; 2196 PetscInt cones[8] = {1, 4, 3, 2, 5, 6, 7, 8}; 2197 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 2198 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, 2199 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0}; 2200 2201 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2202 } 2203 break; 2204 default: 2205 SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Cannot create reference cell for dimension %d", dim); 2206 } 2207 *refdm = NULL; 2208 ierr = DMPlexInterpolate(rdm, refdm);CHKERRQ(ierr); 2209 if (rdm->coordinateDM) { 2210 DM ncdm; 2211 PetscSection cs; 2212 PetscInt pEnd = -1; 2213 2214 ierr = DMGetDefaultSection(rdm->coordinateDM, &cs);CHKERRQ(ierr); 2215 if (cs) {ierr = PetscSectionGetChart(cs, NULL, &pEnd);CHKERRQ(ierr);} 2216 if (pEnd >= 0) { 2217 ierr = DMClone(rdm->coordinateDM, &ncdm);CHKERRQ(ierr); 2218 ierr = DMSetDefaultSection(ncdm, cs);CHKERRQ(ierr); 2219 ierr = DMSetCoordinateDM(*refdm, ncdm);CHKERRQ(ierr); 2220 ierr = DMDestroy(&ncdm);CHKERRQ(ierr); 2221 } 2222 } 2223 ierr = DMGetCoordinatesLocal(rdm, &coords);CHKERRQ(ierr); 2224 if (coords) { 2225 ierr = DMSetCoordinatesLocal(*refdm, coords);CHKERRQ(ierr); 2226 } else { 2227 ierr = DMGetCoordinates(rdm, &coords);CHKERRQ(ierr); 2228 if (coords) {ierr = DMSetCoordinates(*refdm, coords);CHKERRQ(ierr);} 2229 } 2230 ierr = DMDestroy(&rdm);CHKERRQ(ierr); 2231 PetscFunctionReturn(0); 2232 } 2233