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 PetscInt edges[3]; 930 931 edges[0] = cells[0]; 932 edges[1] = cells[1]; 933 edges[2] = 0; 934 935 ierr = DMPlexCreateCubeMesh_Internal(*dm, lower, upper, edges, periodicX, periodicY, DM_BOUNDARY_NONE);CHKERRQ(ierr); 936 if (periodicX == DM_BOUNDARY_PERIODIC || periodicX == DM_BOUNDARY_TWIST || 937 periodicY == DM_BOUNDARY_PERIODIC || periodicY == DM_BOUNDARY_TWIST) { 938 PetscReal L[2]; 939 PetscReal maxCell[2]; 940 DMBoundaryType bdType[2]; 941 942 bdType[0] = periodicX; 943 bdType[1] = periodicY; 944 for (i = 0; i < dim; i++) { 945 L[i] = upper[i] - lower[i]; 946 maxCell[i] = 1.1 * (L[i] / PetscMax(1,cells[i])); 947 } 948 949 ierr = DMSetPeriodicity(*dm,maxCell,L,bdType);CHKERRQ(ierr); 950 } 951 break; 952 } 953 case 3: 954 { 955 PetscReal lower[3] = {0.0, 0.0, 0.0}; 956 PetscReal upper[3] = {1.0, 1.0, 1.0}; 957 958 ierr = DMPlexCreateCubeMesh_Internal(*dm, lower, upper, cells, periodicX, periodicY, periodicZ);CHKERRQ(ierr); 959 if (periodicX == DM_BOUNDARY_PERIODIC || periodicX == DM_BOUNDARY_TWIST || 960 periodicY == DM_BOUNDARY_PERIODIC || periodicY == DM_BOUNDARY_TWIST || 961 periodicZ == DM_BOUNDARY_PERIODIC || periodicZ == DM_BOUNDARY_TWIST) { 962 PetscReal L[3]; 963 PetscReal maxCell[3]; 964 DMBoundaryType bdType[3]; 965 966 bdType[0] = periodicX; 967 bdType[1] = periodicY; 968 bdType[2] = periodicZ; 969 for (i = 0; i < dim; i++) { 970 L[i] = upper[i] - lower[i]; 971 maxCell[i] = 1.1 * (L[i] / cells[i]); 972 } 973 974 ierr = DMSetPeriodicity(*dm,maxCell,L,bdType);CHKERRQ(ierr); 975 } 976 break; 977 } 978 default: 979 SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim); 980 } 981 PetscFunctionReturn(0); 982 } 983 984 /*@ 985 DMPlexCreateHexCylinderMesh - Creates a mesh on the tensor product of the unit interval with the circle (cylinder) using hexahedra. 986 987 Collective on MPI_Comm 988 989 Input Parameters: 990 + comm - The communicator for the DM object 991 . numRefine - The number of regular refinements to the basic 5 cell structure 992 - periodicZ - The boundary type for the Z direction 993 994 Output Parameter: 995 . dm - The DM object 996 997 Note: Here is the output numbering looking from the bottom of the cylinder: 998 $ 17-----14 999 $ | | 1000 $ | 2 | 1001 $ | | 1002 $ 17-----8-----7-----14 1003 $ | | | | 1004 $ | 3 | 0 | 1 | 1005 $ | | | | 1006 $ 19-----5-----6-----13 1007 $ | | 1008 $ | 4 | 1009 $ | | 1010 $ 19-----13 1011 $ 1012 $ and up through the top 1013 $ 1014 $ 18-----16 1015 $ | | 1016 $ | 2 | 1017 $ | | 1018 $ 18----10----11-----16 1019 $ | | | | 1020 $ | 3 | 0 | 1 | 1021 $ | | | | 1022 $ 20-----9----12-----15 1023 $ | | 1024 $ | 4 | 1025 $ | | 1026 $ 20-----15 1027 1028 Level: beginner 1029 1030 .keywords: DM, create 1031 .seealso: DMPlexCreateHexBoxMesh(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 1032 @*/ 1033 PetscErrorCode DMPlexCreateHexCylinderMesh(MPI_Comm comm, PetscInt numRefine, DMBoundaryType periodicZ, DM *dm) 1034 { 1035 const PetscInt dim = 3; 1036 PetscInt numCells, numVertices, r; 1037 PetscErrorCode ierr; 1038 1039 PetscFunctionBegin; 1040 PetscValidPointer(dm, 4); 1041 if (numRefine < 0) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Number of refinements %D cannot be negative", numRefine); 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 numCells = 5; 1050 numVertices = 16; 1051 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1052 numCells *= 2; 1053 } 1054 ierr = DMPlexSetChart(*dm, 0, numCells+numVertices);CHKERRQ(ierr); 1055 for (c = 0; c < numCells; c++) {ierr = DMPlexSetConeSize(*dm, c, 8);CHKERRQ(ierr);} 1056 ierr = DMSetUp(*dm);CHKERRQ(ierr); 1057 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1058 cone[0] = 10; cone[1] = 11; cone[2] = 12; cone[3] = 13; 1059 cone[4] = 14; cone[5] = 15; cone[6] = 16; cone[7] = 17; 1060 ierr = DMPlexSetCone(*dm, 0, cone);CHKERRQ(ierr); 1061 cone[0] = 11; cone[1] = 18; cone[2] = 19; cone[3] = 12; 1062 cone[4] = 17; cone[5] = 16; cone[6] = 21; cone[7] = 20; 1063 ierr = DMPlexSetCone(*dm, 1, cone);CHKERRQ(ierr); 1064 cone[0] = 13; cone[1] = 12; cone[2] = 19; cone[3] = 22; 1065 cone[4] = 15; cone[5] = 23; cone[6] = 21; cone[7] = 16; 1066 ierr = DMPlexSetCone(*dm, 2, cone);CHKERRQ(ierr); 1067 cone[0] = 24; cone[1] = 10; cone[2] = 13; cone[3] = 22; 1068 cone[4] = 25; cone[5] = 23; cone[6] = 15; cone[7] = 14; 1069 ierr = DMPlexSetCone(*dm, 3, cone);CHKERRQ(ierr); 1070 cone[0] = 24; cone[1] = 18; cone[2] = 11; cone[3] = 10; 1071 cone[4] = 25; cone[5] = 14; cone[6] = 17; cone[7] = 20; 1072 ierr = DMPlexSetCone(*dm, 4, cone);CHKERRQ(ierr); 1073 1074 cone[0] = 14; cone[1] = 17; cone[2] = 16; cone[3] = 15; 1075 cone[4] = 10; cone[5] = 13; cone[6] = 12; cone[7] = 11; 1076 ierr = DMPlexSetCone(*dm, 5, cone);CHKERRQ(ierr); 1077 cone[0] = 17; cone[1] = 20; cone[2] = 21; cone[3] = 16; 1078 cone[4] = 11; cone[5] = 12; cone[6] = 19; cone[7] = 18; 1079 ierr = DMPlexSetCone(*dm, 6, cone);CHKERRQ(ierr); 1080 cone[0] = 15; cone[1] = 16; cone[2] = 21; cone[3] = 23; 1081 cone[4] = 13; cone[5] = 22; cone[6] = 19; cone[7] = 12; 1082 ierr = DMPlexSetCone(*dm, 7, cone);CHKERRQ(ierr); 1083 cone[0] = 25; cone[1] = 14; cone[2] = 15; cone[3] = 23; 1084 cone[4] = 24; cone[5] = 22; cone[6] = 13; cone[7] = 10; 1085 ierr = DMPlexSetCone(*dm, 8, cone);CHKERRQ(ierr); 1086 cone[0] = 25; cone[1] = 20; cone[2] = 17; cone[3] = 14; 1087 cone[4] = 24; cone[5] = 10; cone[6] = 11; cone[7] = 18; 1088 ierr = DMPlexSetCone(*dm, 9, cone);CHKERRQ(ierr); 1089 } else { 1090 cone[0] = 5; cone[1] = 6; cone[2] = 7; cone[3] = 8; 1091 cone[4] = 9; cone[5] = 10; cone[6] = 11; cone[7] = 12; 1092 ierr = DMPlexSetCone(*dm, 0, cone);CHKERRQ(ierr); 1093 cone[0] = 6; cone[1] = 13; cone[2] = 14; cone[3] = 7; 1094 cone[4] = 12; cone[5] = 11; cone[6] = 16; cone[7] = 15; 1095 ierr = DMPlexSetCone(*dm, 1, cone);CHKERRQ(ierr); 1096 cone[0] = 8; cone[1] = 7; cone[2] = 14; cone[3] = 17; 1097 cone[4] = 10; cone[5] = 18; cone[6] = 16; cone[7] = 11; 1098 ierr = DMPlexSetCone(*dm, 2, cone);CHKERRQ(ierr); 1099 cone[0] = 19; cone[1] = 5; cone[2] = 8; cone[3] = 17; 1100 cone[4] = 20; cone[5] = 18; cone[6] = 10; cone[7] = 9; 1101 ierr = DMPlexSetCone(*dm, 3, cone);CHKERRQ(ierr); 1102 cone[0] = 19; cone[1] = 13; cone[2] = 6; cone[3] = 5; 1103 cone[4] = 20; cone[5] = 9; cone[6] = 12; cone[7] = 15; 1104 ierr = DMPlexSetCone(*dm, 4, cone);CHKERRQ(ierr); 1105 } 1106 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1107 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1108 } 1109 /* Interpolate */ 1110 { 1111 DM idm = NULL; 1112 1113 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1114 ierr = DMDestroy(dm);CHKERRQ(ierr); 1115 *dm = idm; 1116 } 1117 /* Create cube geometry */ 1118 { 1119 Vec coordinates; 1120 PetscSection coordSection; 1121 PetscScalar *coords; 1122 PetscInt coordSize, v; 1123 const PetscReal dis = 1.0/PetscSqrtReal(2.0); 1124 const PetscReal ds2 = dis/2.0; 1125 1126 /* Build coordinates */ 1127 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 1128 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1129 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 1130 ierr = PetscSectionSetChart(coordSection, numCells, numCells+numVertices);CHKERRQ(ierr); 1131 for (v = numCells; v < numCells+numVertices; ++v) { 1132 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 1133 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 1134 } 1135 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1136 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1137 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 1138 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1139 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1140 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 1141 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1142 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1143 coords[0*dim+0] = -ds2; coords[0*dim+1] = -ds2; coords[0*dim+2] = 0.0; 1144 coords[1*dim+0] = ds2; coords[1*dim+1] = -ds2; coords[1*dim+2] = 0.0; 1145 coords[2*dim+0] = ds2; coords[2*dim+1] = ds2; coords[2*dim+2] = 0.0; 1146 coords[3*dim+0] = -ds2; coords[3*dim+1] = ds2; coords[3*dim+2] = 0.0; 1147 coords[4*dim+0] = -ds2; coords[4*dim+1] = -ds2; coords[4*dim+2] = 1.0; 1148 coords[5*dim+0] = -ds2; coords[5*dim+1] = ds2; coords[5*dim+2] = 1.0; 1149 coords[6*dim+0] = ds2; coords[6*dim+1] = ds2; coords[6*dim+2] = 1.0; 1150 coords[7*dim+0] = ds2; coords[7*dim+1] = -ds2; coords[7*dim+2] = 1.0; 1151 coords[ 8*dim+0] = dis; coords[ 8*dim+1] = -dis; coords[ 8*dim+2] = 0.0; 1152 coords[ 9*dim+0] = dis; coords[ 9*dim+1] = dis; coords[ 9*dim+2] = 0.0; 1153 coords[10*dim+0] = dis; coords[10*dim+1] = -dis; coords[10*dim+2] = 1.0; 1154 coords[11*dim+0] = dis; coords[11*dim+1] = dis; coords[11*dim+2] = 1.0; 1155 coords[12*dim+0] = -dis; coords[12*dim+1] = dis; coords[12*dim+2] = 0.0; 1156 coords[13*dim+0] = -dis; coords[13*dim+1] = dis; coords[13*dim+2] = 1.0; 1157 coords[14*dim+0] = -dis; coords[14*dim+1] = -dis; coords[14*dim+2] = 0.0; 1158 coords[15*dim+0] = -dis; coords[15*dim+1] = -dis; coords[15*dim+2] = 1.0; 1159 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1160 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 1161 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1162 } 1163 /* Create periodicity */ 1164 if (periodicZ == DM_BOUNDARY_PERIODIC || periodicZ == DM_BOUNDARY_TWIST) { 1165 PetscReal L[3]; 1166 PetscReal maxCell[3]; 1167 DMBoundaryType bdType[3]; 1168 PetscReal lower[3] = {0.0, 0.0, 0.0}; 1169 PetscReal upper[3] = {1.0, 1.0, 2.0}; 1170 PetscInt i, numZCells = PetscPowInt(2, numRefine); 1171 1172 bdType[0] = DM_BOUNDARY_NONE; 1173 bdType[1] = DM_BOUNDARY_NONE; 1174 bdType[2] = periodicZ; 1175 for (i = 0; i < dim; i++) { 1176 L[i] = upper[i] - lower[i]; 1177 maxCell[i] = 1.1 * (L[i] / numZCells); 1178 } 1179 ierr = DMSetPeriodicity(*dm, maxCell, L, bdType);CHKERRQ(ierr); 1180 } 1181 /* Refine topology */ 1182 for (r = 0; r < numRefine; ++r) { 1183 DM rdm = NULL; 1184 1185 ierr = DMRefine(*dm, comm, &rdm);CHKERRQ(ierr); 1186 ierr = DMDestroy(dm);CHKERRQ(ierr); 1187 *dm = rdm; 1188 } 1189 /* Remap geometry to cylinder 1190 Interior square: Linear interpolation is correct 1191 The other cells all have vertices on rays from the origin. We want to uniformly expand the spacing 1192 such that the last vertex is on the unit circle. So the closest and farthest vertices are at distance 1193 1194 phi = arctan(y/x) 1195 d_close = sqrt(1/8 + 1/4 sin^2(phi)) 1196 d_far = sqrt(1/2 + sin^2(phi)) 1197 1198 so we remap them using 1199 1200 x_new = x_close + (x - x_close) (1 - d_close) / (d_far - d_close) 1201 y_new = y_close + (y - y_close) (1 - d_close) / (d_far - d_close) 1202 1203 If pi/4 < phi < 3pi/4 or -3pi/4 < phi < -pi/4, then we switch x and y. 1204 */ 1205 { 1206 Vec coordinates; 1207 PetscSection coordSection; 1208 PetscScalar *coords; 1209 PetscInt vStart, vEnd, v; 1210 const PetscReal dis = 1.0/PetscSqrtReal(2.0); 1211 const PetscReal ds2 = 0.5*dis; 1212 1213 ierr = DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1214 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 1215 ierr = DMGetCoordinatesLocal(*dm, &coordinates);CHKERRQ(ierr); 1216 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1217 for (v = vStart; v < vEnd; ++v) { 1218 PetscReal phi, sinp, cosp, dc, df, x, y, xc, yc; 1219 PetscInt off; 1220 1221 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr); 1222 if ((PetscAbsScalar(coords[off+0]) <= ds2) && (PetscAbsScalar(coords[off+1]) <= ds2)) continue; 1223 x = PetscRealPart(coords[off]); 1224 y = PetscRealPart(coords[off+1]); 1225 phi = PetscAtan2Real(y, x); 1226 sinp = PetscSinReal(phi); 1227 cosp = PetscCosReal(phi); 1228 if ((PetscAbsReal(phi) > PETSC_PI/4.0) && (PetscAbsReal(phi) < 3.0*PETSC_PI/4.0)) { 1229 dc = PetscAbsReal(ds2/sinp); 1230 df = PetscAbsReal(dis/sinp); 1231 xc = ds2*x/PetscAbsScalar(y); 1232 yc = ds2*PetscSignReal(y); 1233 } else { 1234 dc = PetscAbsReal(ds2/cosp); 1235 df = PetscAbsReal(dis/cosp); 1236 xc = ds2*PetscSignReal(x); 1237 yc = ds2*y/PetscAbsScalar(x); 1238 } 1239 coords[off+0] = xc + (coords[off+0] - xc)*(1.0 - dc)/(df - dc); 1240 coords[off+1] = yc + (coords[off+1] - yc)*(1.0 - dc)/(df - dc); 1241 } 1242 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1243 if (periodicZ == DM_BOUNDARY_PERIODIC || periodicZ == DM_BOUNDARY_TWIST) { 1244 ierr = DMLocalizeCoordinates(*dm);CHKERRQ(ierr); 1245 } 1246 } 1247 PetscFunctionReturn(0); 1248 } 1249 1250 /* External function declarations here */ 1251 extern PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling); 1252 extern PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, Mat *mat); 1253 extern PetscErrorCode DMCreateDefaultSection_Plex(DM dm); 1254 extern PetscErrorCode DMCreateDefaultConstraints_Plex(DM dm); 1255 extern PetscErrorCode DMCreateMatrix_Plex(DM dm, Mat *J); 1256 extern PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm); 1257 PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm); 1258 extern PetscErrorCode DMSetUp_Plex(DM dm); 1259 extern PetscErrorCode DMDestroy_Plex(DM dm); 1260 extern PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer); 1261 extern PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer); 1262 extern PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm); 1263 1264 /* Replace dm with the contents of dmNew 1265 - Share the DM_Plex structure 1266 - Share the coordinates 1267 - Share the SF 1268 */ 1269 static PetscErrorCode DMPlexReplace_Static(DM dm, DM dmNew) 1270 { 1271 PetscSF sf; 1272 DM coordDM, coarseDM; 1273 Vec coords; 1274 const PetscReal *maxCell, *L; 1275 const DMBoundaryType *bd; 1276 PetscErrorCode ierr; 1277 1278 PetscFunctionBegin; 1279 ierr = DMGetPointSF(dmNew, &sf);CHKERRQ(ierr); 1280 ierr = DMSetPointSF(dm, sf);CHKERRQ(ierr); 1281 ierr = DMGetCoordinateDM(dmNew, &coordDM);CHKERRQ(ierr); 1282 ierr = DMGetCoordinatesLocal(dmNew, &coords);CHKERRQ(ierr); 1283 ierr = DMSetCoordinateDM(dm, coordDM);CHKERRQ(ierr); 1284 ierr = DMSetCoordinatesLocal(dm, coords);CHKERRQ(ierr); 1285 ierr = DMGetPeriodicity(dm, &maxCell, &L, &bd);CHKERRQ(ierr); 1286 if (L) {ierr = DMSetPeriodicity(dmNew, maxCell, L, bd);CHKERRQ(ierr);} 1287 ierr = DMDestroy_Plex(dm);CHKERRQ(ierr); 1288 dm->data = dmNew->data; 1289 ((DM_Plex *) dmNew->data)->refct++; 1290 dmNew->labels->refct++; 1291 if (!--(dm->labels->refct)) { 1292 DMLabelLink next = dm->labels->next; 1293 1294 /* destroy the labels */ 1295 while (next) { 1296 DMLabelLink tmp = next->next; 1297 1298 ierr = DMLabelDestroy(&next->label);CHKERRQ(ierr); 1299 ierr = PetscFree(next);CHKERRQ(ierr); 1300 next = tmp; 1301 } 1302 ierr = PetscFree(dm->labels);CHKERRQ(ierr); 1303 } 1304 dm->labels = dmNew->labels; 1305 dm->depthLabel = dmNew->depthLabel; 1306 ierr = DMGetCoarseDM(dmNew,&coarseDM);CHKERRQ(ierr); 1307 ierr = DMSetCoarseDM(dm,coarseDM);CHKERRQ(ierr); 1308 PetscFunctionReturn(0); 1309 } 1310 1311 /* Swap dm with the contents of dmNew 1312 - Swap the DM_Plex structure 1313 - Swap the coordinates 1314 - Swap the point PetscSF 1315 */ 1316 static PetscErrorCode DMPlexSwap_Static(DM dmA, DM dmB) 1317 { 1318 DM coordDMA, coordDMB; 1319 Vec coordsA, coordsB; 1320 PetscSF sfA, sfB; 1321 void *tmp; 1322 DMLabelLinkList listTmp; 1323 DMLabel depthTmp; 1324 PetscErrorCode ierr; 1325 1326 PetscFunctionBegin; 1327 ierr = DMGetPointSF(dmA, &sfA);CHKERRQ(ierr); 1328 ierr = DMGetPointSF(dmB, &sfB);CHKERRQ(ierr); 1329 ierr = PetscObjectReference((PetscObject) sfA);CHKERRQ(ierr); 1330 ierr = DMSetPointSF(dmA, sfB);CHKERRQ(ierr); 1331 ierr = DMSetPointSF(dmB, sfA);CHKERRQ(ierr); 1332 ierr = PetscObjectDereference((PetscObject) sfA);CHKERRQ(ierr); 1333 1334 ierr = DMGetCoordinateDM(dmA, &coordDMA);CHKERRQ(ierr); 1335 ierr = DMGetCoordinateDM(dmB, &coordDMB);CHKERRQ(ierr); 1336 ierr = PetscObjectReference((PetscObject) coordDMA);CHKERRQ(ierr); 1337 ierr = DMSetCoordinateDM(dmA, coordDMB);CHKERRQ(ierr); 1338 ierr = DMSetCoordinateDM(dmB, coordDMA);CHKERRQ(ierr); 1339 ierr = PetscObjectDereference((PetscObject) coordDMA);CHKERRQ(ierr); 1340 1341 ierr = DMGetCoordinatesLocal(dmA, &coordsA);CHKERRQ(ierr); 1342 ierr = DMGetCoordinatesLocal(dmB, &coordsB);CHKERRQ(ierr); 1343 ierr = PetscObjectReference((PetscObject) coordsA);CHKERRQ(ierr); 1344 ierr = DMSetCoordinatesLocal(dmA, coordsB);CHKERRQ(ierr); 1345 ierr = DMSetCoordinatesLocal(dmB, coordsA);CHKERRQ(ierr); 1346 ierr = PetscObjectDereference((PetscObject) coordsA);CHKERRQ(ierr); 1347 1348 tmp = dmA->data; 1349 dmA->data = dmB->data; 1350 dmB->data = tmp; 1351 listTmp = dmA->labels; 1352 dmA->labels = dmB->labels; 1353 dmB->labels = listTmp; 1354 depthTmp = dmA->depthLabel; 1355 dmA->depthLabel = dmB->depthLabel; 1356 dmB->depthLabel = depthTmp; 1357 PetscFunctionReturn(0); 1358 } 1359 1360 PetscErrorCode DMSetFromOptions_NonRefinement_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 1361 { 1362 DM_Plex *mesh = (DM_Plex*) dm->data; 1363 PetscErrorCode ierr; 1364 1365 PetscFunctionBegin; 1366 /* Handle viewing */ 1367 ierr = PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMView", PETSC_FALSE, &mesh->printSetValues, NULL);CHKERRQ(ierr); 1368 ierr = PetscOptionsInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMView", 0, &mesh->printFEM, NULL);CHKERRQ(ierr); 1369 ierr = PetscOptionsReal("-dm_plex_print_tol", "Tolerance for FEM output", "DMView", mesh->printTol, &mesh->printTol, NULL);CHKERRQ(ierr); 1370 /* Point Location */ 1371 ierr = PetscOptionsBool("-dm_plex_hash_location", "Use grid hashing for point location", "DMView", PETSC_FALSE, &mesh->useHashLocation, NULL);CHKERRQ(ierr); 1372 /* Generation and remeshing */ 1373 ierr = PetscOptionsBool("-dm_plex_remesh_bd", "Allow changes to the boundary on remeshing", "DMView", PETSC_FALSE, &mesh->remeshBd, NULL);CHKERRQ(ierr); 1374 /* Projection behavior */ 1375 ierr = PetscOptionsInt("-dm_plex_max_projection_height", "Maxmimum mesh point height used to project locally", "DMPlexSetMaxProjectionHeight", 0, &mesh->maxProjectionHeight, NULL);CHKERRQ(ierr); 1376 ierr = PetscOptionsBool("-dm_plex_regular_refinement", "Use special nested projection algorithm for regular refinement", "DMPlexSetRegularRefinement", mesh->regularRefinement, &mesh->regularRefinement, NULL);CHKERRQ(ierr); 1377 PetscFunctionReturn(0); 1378 } 1379 1380 static PetscErrorCode DMSetFromOptions_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 1381 { 1382 PetscInt refine = 0, coarsen = 0, r; 1383 PetscBool isHierarchy; 1384 PetscErrorCode ierr; 1385 1386 PetscFunctionBegin; 1387 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1388 ierr = PetscOptionsHead(PetscOptionsObject,"DMPlex Options");CHKERRQ(ierr); 1389 /* Handle DMPlex refinement */ 1390 ierr = PetscOptionsInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL);CHKERRQ(ierr); 1391 ierr = PetscOptionsInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy);CHKERRQ(ierr); 1392 if (refine) {ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE);CHKERRQ(ierr);} 1393 if (refine && isHierarchy) { 1394 DM *dms, coarseDM; 1395 1396 ierr = DMGetCoarseDM(dm, &coarseDM);CHKERRQ(ierr); 1397 ierr = PetscObjectReference((PetscObject)coarseDM);CHKERRQ(ierr); 1398 ierr = PetscMalloc1(refine,&dms);CHKERRQ(ierr); 1399 ierr = DMRefineHierarchy(dm, refine, dms);CHKERRQ(ierr); 1400 /* Total hack since we do not pass in a pointer */ 1401 ierr = DMPlexSwap_Static(dm, dms[refine-1]);CHKERRQ(ierr); 1402 if (refine == 1) { 1403 ierr = DMSetCoarseDM(dm, dms[0]);CHKERRQ(ierr); 1404 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr); 1405 } else { 1406 ierr = DMSetCoarseDM(dm, dms[refine-2]);CHKERRQ(ierr); 1407 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr); 1408 ierr = DMSetCoarseDM(dms[0], dms[refine-1]);CHKERRQ(ierr); 1409 ierr = DMPlexSetRegularRefinement(dms[0], PETSC_TRUE);CHKERRQ(ierr); 1410 } 1411 ierr = DMSetCoarseDM(dms[refine-1], coarseDM);CHKERRQ(ierr); 1412 ierr = PetscObjectDereference((PetscObject)coarseDM);CHKERRQ(ierr); 1413 /* Free DMs */ 1414 for (r = 0; r < refine; ++r) { 1415 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr); 1416 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 1417 } 1418 ierr = PetscFree(dms);CHKERRQ(ierr); 1419 } else { 1420 for (r = 0; r < refine; ++r) { 1421 DM refinedMesh; 1422 1423 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 1424 ierr = DMRefine(dm, PetscObjectComm((PetscObject) dm), &refinedMesh);CHKERRQ(ierr); 1425 /* Total hack since we do not pass in a pointer */ 1426 ierr = DMPlexReplace_Static(dm, refinedMesh);CHKERRQ(ierr); 1427 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 1428 ierr = DMDestroy(&refinedMesh);CHKERRQ(ierr); 1429 } 1430 } 1431 /* Handle DMPlex coarsening */ 1432 ierr = PetscOptionsInt("-dm_coarsen", "Coarsen the mesh", "DMCreate", coarsen, &coarsen, NULL);CHKERRQ(ierr); 1433 ierr = PetscOptionsInt("-dm_coarsen_hierarchy", "The number of coarsenings", "DMCreate", coarsen, &coarsen, &isHierarchy);CHKERRQ(ierr); 1434 if (coarsen && isHierarchy) { 1435 DM *dms; 1436 1437 ierr = PetscMalloc1(coarsen, &dms);CHKERRQ(ierr); 1438 ierr = DMCoarsenHierarchy(dm, coarsen, dms);CHKERRQ(ierr); 1439 /* Free DMs */ 1440 for (r = 0; r < coarsen; ++r) { 1441 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr); 1442 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 1443 } 1444 ierr = PetscFree(dms);CHKERRQ(ierr); 1445 } else { 1446 for (r = 0; r < coarsen; ++r) { 1447 DM coarseMesh; 1448 1449 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 1450 ierr = DMCoarsen(dm, PetscObjectComm((PetscObject) dm), &coarseMesh);CHKERRQ(ierr); 1451 /* Total hack since we do not pass in a pointer */ 1452 ierr = DMPlexReplace_Static(dm, coarseMesh);CHKERRQ(ierr); 1453 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 1454 ierr = DMDestroy(&coarseMesh);CHKERRQ(ierr); 1455 } 1456 } 1457 /* Handle */ 1458 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 1459 ierr = PetscOptionsTail();CHKERRQ(ierr); 1460 PetscFunctionReturn(0); 1461 } 1462 1463 static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec) 1464 { 1465 PetscErrorCode ierr; 1466 1467 PetscFunctionBegin; 1468 ierr = DMCreateGlobalVector_Section_Private(dm,vec);CHKERRQ(ierr); 1469 /* ierr = VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM);CHKERRQ(ierr); */ 1470 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex);CHKERRQ(ierr); 1471 ierr = VecSetOperation(*vec, VECOP_VIEWNATIVE, (void (*)(void)) VecView_Plex_Native);CHKERRQ(ierr); 1472 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex);CHKERRQ(ierr); 1473 ierr = VecSetOperation(*vec, VECOP_LOADNATIVE, (void (*)(void)) VecLoad_Plex_Native);CHKERRQ(ierr); 1474 PetscFunctionReturn(0); 1475 } 1476 1477 static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec) 1478 { 1479 PetscErrorCode ierr; 1480 1481 PetscFunctionBegin; 1482 ierr = DMCreateLocalVector_Section_Private(dm,vec);CHKERRQ(ierr); 1483 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex_Local);CHKERRQ(ierr); 1484 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex_Local);CHKERRQ(ierr); 1485 PetscFunctionReturn(0); 1486 } 1487 1488 static PetscErrorCode DMGetDimPoints_Plex(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd) 1489 { 1490 PetscInt depth, d; 1491 PetscErrorCode ierr; 1492 1493 PetscFunctionBegin; 1494 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1495 if (depth == 1) { 1496 ierr = DMGetDimension(dm, &d);CHKERRQ(ierr); 1497 if (dim == 0) {ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr);} 1498 else if (dim == d) {ierr = DMPlexGetDepthStratum(dm, 1, pStart, pEnd);CHKERRQ(ierr);} 1499 else {*pStart = 0; *pEnd = 0;} 1500 } else { 1501 ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr); 1502 } 1503 PetscFunctionReturn(0); 1504 } 1505 1506 static PetscErrorCode DMGetNeighors_Plex(DM dm, PetscInt *nranks, const PetscMPIInt *ranks[]) 1507 { 1508 PetscSF sf; 1509 PetscErrorCode ierr; 1510 1511 PetscFunctionBegin; 1512 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 1513 ierr = PetscSFGetRanks(sf, nranks, ranks, NULL, NULL, NULL);CHKERRQ(ierr); 1514 PetscFunctionReturn(0); 1515 } 1516 1517 static PetscErrorCode DMInitialize_Plex(DM dm) 1518 { 1519 PetscErrorCode ierr; 1520 1521 PetscFunctionBegin; 1522 dm->ops->view = DMView_Plex; 1523 dm->ops->load = DMLoad_Plex; 1524 dm->ops->setfromoptions = DMSetFromOptions_Plex; 1525 dm->ops->clone = DMClone_Plex; 1526 dm->ops->setup = DMSetUp_Plex; 1527 dm->ops->createdefaultsection = DMCreateDefaultSection_Plex; 1528 dm->ops->createdefaultconstraints = DMCreateDefaultConstraints_Plex; 1529 dm->ops->createglobalvector = DMCreateGlobalVector_Plex; 1530 dm->ops->createlocalvector = DMCreateLocalVector_Plex; 1531 dm->ops->getlocaltoglobalmapping = NULL; 1532 dm->ops->createfieldis = NULL; 1533 dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex; 1534 dm->ops->getcoloring = NULL; 1535 dm->ops->creatematrix = DMCreateMatrix_Plex; 1536 dm->ops->createinterpolation = DMCreateInterpolation_Plex; 1537 dm->ops->getaggregates = NULL; 1538 dm->ops->getinjection = DMCreateInjection_Plex; 1539 dm->ops->refine = DMRefine_Plex; 1540 dm->ops->coarsen = DMCoarsen_Plex; 1541 dm->ops->refinehierarchy = DMRefineHierarchy_Plex; 1542 dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Plex; 1543 dm->ops->globaltolocalbegin = NULL; 1544 dm->ops->globaltolocalend = NULL; 1545 dm->ops->localtoglobalbegin = NULL; 1546 dm->ops->localtoglobalend = NULL; 1547 dm->ops->destroy = DMDestroy_Plex; 1548 dm->ops->createsubdm = DMCreateSubDM_Plex; 1549 dm->ops->getdimpoints = DMGetDimPoints_Plex; 1550 dm->ops->locatepoints = DMLocatePoints_Plex; 1551 dm->ops->projectfunctionlocal = DMProjectFunctionLocal_Plex; 1552 dm->ops->projectfunctionlabellocal = DMProjectFunctionLabelLocal_Plex; 1553 dm->ops->projectfieldlocal = DMProjectFieldLocal_Plex; 1554 dm->ops->projectfieldlabellocal = DMProjectFieldLabelLocal_Plex; 1555 dm->ops->computel2diff = DMComputeL2Diff_Plex; 1556 dm->ops->computel2gradientdiff = DMComputeL2GradientDiff_Plex; 1557 dm->ops->computel2fielddiff = DMComputeL2FieldDiff_Plex; 1558 dm->ops->getneighbors = DMGetNeighors_Plex; 1559 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMAdaptLabel_C",DMAdaptLabel_Plex);CHKERRQ(ierr); 1560 PetscFunctionReturn(0); 1561 } 1562 1563 PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm) 1564 { 1565 DM_Plex *mesh = (DM_Plex *) dm->data; 1566 PetscErrorCode ierr; 1567 1568 PetscFunctionBegin; 1569 mesh->refct++; 1570 (*newdm)->data = mesh; 1571 ierr = PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);CHKERRQ(ierr); 1572 ierr = DMInitialize_Plex(*newdm);CHKERRQ(ierr); 1573 PetscFunctionReturn(0); 1574 } 1575 1576 /*MC 1577 DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram. 1578 In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is 1579 specified by a PetscSection object. Ownership in the global representation is determined by 1580 ownership of the underlying DMPlex points. This is specified by another PetscSection object. 1581 1582 Level: intermediate 1583 1584 .seealso: DMType, DMPlexCreate(), DMCreate(), DMSetType() 1585 M*/ 1586 1587 PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm) 1588 { 1589 DM_Plex *mesh; 1590 PetscInt unit, d; 1591 PetscErrorCode ierr; 1592 1593 PetscFunctionBegin; 1594 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1595 ierr = PetscNewLog(dm,&mesh);CHKERRQ(ierr); 1596 dm->dim = 0; 1597 dm->data = mesh; 1598 1599 mesh->refct = 1; 1600 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);CHKERRQ(ierr); 1601 mesh->maxConeSize = 0; 1602 mesh->cones = NULL; 1603 mesh->coneOrientations = NULL; 1604 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);CHKERRQ(ierr); 1605 mesh->maxSupportSize = 0; 1606 mesh->supports = NULL; 1607 mesh->refinementUniform = PETSC_TRUE; 1608 mesh->refinementLimit = -1.0; 1609 1610 mesh->facesTmp = NULL; 1611 1612 mesh->tetgenOpts = NULL; 1613 mesh->triangleOpts = NULL; 1614 ierr = PetscPartitionerCreate(PetscObjectComm((PetscObject)dm), &mesh->partitioner);CHKERRQ(ierr); 1615 ierr = PetscPartitionerSetTypeFromOptions_Internal(mesh->partitioner);CHKERRQ(ierr); 1616 mesh->remeshBd = PETSC_FALSE; 1617 1618 mesh->subpointMap = NULL; 1619 1620 for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0; 1621 1622 mesh->regularRefinement = PETSC_FALSE; 1623 mesh->depthState = -1; 1624 mesh->globalVertexNumbers = NULL; 1625 mesh->globalCellNumbers = NULL; 1626 mesh->anchorSection = NULL; 1627 mesh->anchorIS = NULL; 1628 mesh->createanchors = NULL; 1629 mesh->computeanchormatrix = NULL; 1630 mesh->parentSection = NULL; 1631 mesh->parents = NULL; 1632 mesh->childIDs = NULL; 1633 mesh->childSection = NULL; 1634 mesh->children = NULL; 1635 mesh->referenceTree = NULL; 1636 mesh->getchildsymmetry = NULL; 1637 for (d = 0; d < 8; ++d) mesh->hybridPointMax[d] = PETSC_DETERMINE; 1638 mesh->vtkCellHeight = 0; 1639 mesh->useCone = PETSC_FALSE; 1640 mesh->useClosure = PETSC_TRUE; 1641 mesh->useAnchors = PETSC_FALSE; 1642 1643 mesh->maxProjectionHeight = 0; 1644 1645 mesh->printSetValues = PETSC_FALSE; 1646 mesh->printFEM = 0; 1647 mesh->printTol = 1.0e-10; 1648 1649 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 1650 PetscFunctionReturn(0); 1651 } 1652 1653 /*@ 1654 DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram. 1655 1656 Collective on MPI_Comm 1657 1658 Input Parameter: 1659 . comm - The communicator for the DMPlex object 1660 1661 Output Parameter: 1662 . mesh - The DMPlex object 1663 1664 Level: beginner 1665 1666 .keywords: DMPlex, create 1667 @*/ 1668 PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh) 1669 { 1670 PetscErrorCode ierr; 1671 1672 PetscFunctionBegin; 1673 PetscValidPointer(mesh,2); 1674 ierr = DMCreate(comm, mesh);CHKERRQ(ierr); 1675 ierr = DMSetType(*mesh, DMPLEX);CHKERRQ(ierr); 1676 PetscFunctionReturn(0); 1677 } 1678 1679 /* 1680 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 1681 */ 1682 static PetscErrorCode DMPlexBuildFromCellList_Parallel_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[], PetscSF *sfVert) 1683 { 1684 PetscSF sfPoint; 1685 PetscLayout vLayout; 1686 PetscHashI vhash; 1687 PetscSFNode *remoteVerticesAdj, *vertexLocal, *vertexOwner, *remoteVertex; 1688 const PetscInt *vrange; 1689 PetscInt numVerticesAdj, off, *verticesAdj, numVerticesGhost = 0, *localVertex, *cone, c, p, v, g; 1690 PETSC_UNUSED PetscHashIIter ret, iter; 1691 PetscMPIInt rank, size; 1692 PetscErrorCode ierr; 1693 1694 PetscFunctionBegin; 1695 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr); 1696 ierr = MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);CHKERRQ(ierr); 1697 /* Partition vertices */ 1698 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) dm), &vLayout);CHKERRQ(ierr); 1699 ierr = PetscLayoutSetLocalSize(vLayout, numVertices);CHKERRQ(ierr); 1700 ierr = PetscLayoutSetBlockSize(vLayout, 1);CHKERRQ(ierr); 1701 ierr = PetscLayoutSetUp(vLayout);CHKERRQ(ierr); 1702 ierr = PetscLayoutGetRanges(vLayout, &vrange);CHKERRQ(ierr); 1703 /* Count vertices and map them to procs */ 1704 PetscHashICreate(vhash); 1705 for (c = 0; c < numCells; ++c) { 1706 for (p = 0; p < numCorners; ++p) { 1707 PetscHashIPut(vhash, cells[c*numCorners+p], ret, iter); 1708 } 1709 } 1710 PetscHashISize(vhash, numVerticesAdj); 1711 ierr = PetscMalloc1(numVerticesAdj, &verticesAdj);CHKERRQ(ierr); 1712 off = 0; ierr = PetscHashIGetKeys(vhash, &off, verticesAdj);CHKERRQ(ierr); 1713 if (off != numVerticesAdj) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid number of local vertices %D should be %D", off, numVerticesAdj); 1714 ierr = PetscSortInt(numVerticesAdj, verticesAdj);CHKERRQ(ierr); 1715 ierr = PetscMalloc1(numVerticesAdj, &remoteVerticesAdj);CHKERRQ(ierr); 1716 for (v = 0; v < numVerticesAdj; ++v) { 1717 const PetscInt gv = verticesAdj[v]; 1718 PetscInt vrank; 1719 1720 ierr = PetscFindInt(gv, size+1, vrange, &vrank);CHKERRQ(ierr); 1721 vrank = vrank < 0 ? -(vrank+2) : vrank; 1722 remoteVerticesAdj[v].index = gv - vrange[vrank]; 1723 remoteVerticesAdj[v].rank = vrank; 1724 } 1725 /* Create cones */ 1726 ierr = DMPlexSetChart(dm, 0, numCells+numVerticesAdj);CHKERRQ(ierr); 1727 for (c = 0; c < numCells; ++c) {ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr);} 1728 ierr = DMSetUp(dm);CHKERRQ(ierr); 1729 ierr = DMGetWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 1730 for (c = 0; c < numCells; ++c) { 1731 for (p = 0; p < numCorners; ++p) { 1732 const PetscInt gv = cells[c*numCorners+p]; 1733 PetscInt lv; 1734 1735 ierr = PetscFindInt(gv, numVerticesAdj, verticesAdj, &lv);CHKERRQ(ierr); 1736 if (lv < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not find global vertex %D in local connectivity", gv); 1737 cone[p] = lv+numCells; 1738 } 1739 ierr = DMPlexSetCone(dm, c, cone);CHKERRQ(ierr); 1740 } 1741 ierr = DMRestoreWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 1742 /* Create SF for vertices */ 1743 ierr = PetscSFCreate(PetscObjectComm((PetscObject)dm), sfVert);CHKERRQ(ierr); 1744 ierr = PetscObjectSetName((PetscObject) *sfVert, "Vertex Ownership SF");CHKERRQ(ierr); 1745 ierr = PetscSFSetFromOptions(*sfVert);CHKERRQ(ierr); 1746 ierr = PetscSFSetGraph(*sfVert, numVertices, numVerticesAdj, NULL, PETSC_OWN_POINTER, remoteVerticesAdj, PETSC_OWN_POINTER);CHKERRQ(ierr); 1747 ierr = PetscFree(verticesAdj);CHKERRQ(ierr); 1748 /* Build pointSF */ 1749 ierr = PetscMalloc2(numVerticesAdj, &vertexLocal, numVertices, &vertexOwner);CHKERRQ(ierr); 1750 for (v = 0; v < numVerticesAdj; ++v) {vertexLocal[v].index = v+numCells; vertexLocal[v].rank = rank;} 1751 for (v = 0; v < numVertices; ++v) {vertexOwner[v].index = -1; vertexOwner[v].rank = -1;} 1752 ierr = PetscSFReduceBegin(*sfVert, MPIU_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);CHKERRQ(ierr); 1753 ierr = PetscSFReduceEnd(*sfVert, MPIU_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);CHKERRQ(ierr); 1754 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); 1755 ierr = PetscSFBcastBegin(*sfVert, MPIU_2INT, vertexOwner, vertexLocal);CHKERRQ(ierr); 1756 ierr = PetscSFBcastEnd(*sfVert, MPIU_2INT, vertexOwner, vertexLocal);CHKERRQ(ierr); 1757 for (v = 0; v < numVerticesAdj; ++v) if (vertexLocal[v].rank != rank) ++numVerticesGhost; 1758 ierr = PetscMalloc1(numVerticesGhost, &localVertex);CHKERRQ(ierr); 1759 ierr = PetscMalloc1(numVerticesGhost, &remoteVertex);CHKERRQ(ierr); 1760 for (v = 0, g = 0; v < numVerticesAdj; ++v) { 1761 if (vertexLocal[v].rank != rank) { 1762 localVertex[g] = v+numCells; 1763 remoteVertex[g].index = vertexLocal[v].index; 1764 remoteVertex[g].rank = vertexLocal[v].rank; 1765 ++g; 1766 } 1767 } 1768 ierr = PetscFree2(vertexLocal, vertexOwner);CHKERRQ(ierr); 1769 if (g != numVerticesGhost) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of vertex ghosts %D should be %D", g, numVerticesGhost); 1770 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 1771 ierr = PetscObjectSetName((PetscObject) sfPoint, "point SF");CHKERRQ(ierr); 1772 ierr = PetscSFSetGraph(sfPoint, numCells+numVerticesAdj, numVerticesGhost, localVertex, PETSC_OWN_POINTER, remoteVertex, PETSC_OWN_POINTER);CHKERRQ(ierr); 1773 ierr = PetscLayoutDestroy(&vLayout);CHKERRQ(ierr); 1774 PetscHashIDestroy(vhash); 1775 /* Fill in the rest of the topology structure */ 1776 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 1777 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 1778 PetscFunctionReturn(0); 1779 } 1780 1781 /* 1782 This takes as input the coordinates for each owned vertex 1783 */ 1784 static PetscErrorCode DMPlexBuildCoordinates_Parallel_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numV, PetscSF sfVert, const PetscReal vertexCoords[]) 1785 { 1786 PetscSection coordSection; 1787 Vec coordinates; 1788 PetscScalar *coords; 1789 PetscInt numVertices, numVerticesAdj, coordSize, v; 1790 PetscErrorCode ierr; 1791 1792 PetscFunctionBegin; 1793 ierr = PetscSFGetGraph(sfVert, &numVertices, &numVerticesAdj, NULL, NULL);CHKERRQ(ierr); 1794 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1795 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1796 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 1797 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVerticesAdj);CHKERRQ(ierr); 1798 for (v = numCells; v < numCells+numVerticesAdj; ++v) { 1799 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 1800 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 1801 } 1802 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1803 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1804 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 1805 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 1806 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1807 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1808 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1809 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1810 { 1811 MPI_Datatype coordtype; 1812 1813 /* Need a temp buffer for coords if we have complex/single */ 1814 ierr = MPI_Type_contiguous(spaceDim, MPIU_SCALAR, &coordtype);CHKERRQ(ierr); 1815 ierr = MPI_Type_commit(&coordtype);CHKERRQ(ierr); 1816 #if defined(PETSC_USE_COMPLEX) 1817 { 1818 PetscScalar *svertexCoords; 1819 PetscInt i; 1820 ierr = PetscMalloc1(numV*spaceDim,&svertexCoords);CHKERRQ(ierr); 1821 for (i=0; i<numV*spaceDim; i++) svertexCoords[i] = vertexCoords[i]; 1822 ierr = PetscSFBcastBegin(sfVert, coordtype, svertexCoords, coords);CHKERRQ(ierr); 1823 ierr = PetscSFBcastEnd(sfVert, coordtype, svertexCoords, coords);CHKERRQ(ierr); 1824 ierr = PetscFree(svertexCoords);CHKERRQ(ierr); 1825 } 1826 #else 1827 ierr = PetscSFBcastBegin(sfVert, coordtype, vertexCoords, coords);CHKERRQ(ierr); 1828 ierr = PetscSFBcastEnd(sfVert, coordtype, vertexCoords, coords);CHKERRQ(ierr); 1829 #endif 1830 ierr = MPI_Type_free(&coordtype);CHKERRQ(ierr); 1831 } 1832 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1833 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 1834 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1835 PetscFunctionReturn(0); 1836 } 1837 1838 /*@C 1839 DMPlexCreateFromCellListParallel - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM 1840 1841 Input Parameters: 1842 + comm - The communicator 1843 . dim - The topological dimension of the mesh 1844 . numCells - The number of cells owned by this process 1845 . numVertices - The number of vertices owned by this process 1846 . numCorners - The number of vertices for each cell 1847 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 1848 . cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 1849 . spaceDim - The spatial dimension used for coordinates 1850 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 1851 1852 Output Parameter: 1853 + dm - The DM 1854 - vertexSF - Optional, SF describing complete vertex ownership 1855 1856 Note: Two triangles sharing a face 1857 $ 1858 $ 2 1859 $ / | \ 1860 $ / | \ 1861 $ / | \ 1862 $ 0 0 | 1 3 1863 $ \ | / 1864 $ \ | / 1865 $ \ | / 1866 $ 1 1867 would have input 1868 $ numCells = 2, numVertices = 4 1869 $ cells = [0 1 2 1 3 2] 1870 $ 1871 which would result in the DMPlex 1872 $ 1873 $ 4 1874 $ / | \ 1875 $ / | \ 1876 $ / | \ 1877 $ 2 0 | 1 5 1878 $ \ | / 1879 $ \ | / 1880 $ \ | / 1881 $ 3 1882 1883 Level: beginner 1884 1885 .seealso: DMPlexCreateFromCellList(), DMPlexCreateFromDAG(), DMPlexCreate() 1886 @*/ 1887 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) 1888 { 1889 PetscSF sfVert; 1890 PetscErrorCode ierr; 1891 1892 PetscFunctionBegin; 1893 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1894 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1895 PetscValidLogicalCollectiveInt(*dm, dim, 2); 1896 PetscValidLogicalCollectiveInt(*dm, spaceDim, 8); 1897 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1898 ierr = DMPlexBuildFromCellList_Parallel_Private(*dm, numCells, numVertices, numCorners, cells, &sfVert);CHKERRQ(ierr); 1899 if (interpolate) { 1900 DM idm = NULL; 1901 1902 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1903 ierr = DMDestroy(dm);CHKERRQ(ierr); 1904 *dm = idm; 1905 } 1906 ierr = DMPlexBuildCoordinates_Parallel_Private(*dm, spaceDim, numCells, numVertices,sfVert, vertexCoords);CHKERRQ(ierr); 1907 if (vertexSF) *vertexSF = sfVert; 1908 else {ierr = PetscSFDestroy(&sfVert);CHKERRQ(ierr);} 1909 PetscFunctionReturn(0); 1910 } 1911 1912 /* 1913 This takes as input the common mesh generator output, a list of the vertices for each cell 1914 */ 1915 static PetscErrorCode DMPlexBuildFromCellList_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[]) 1916 { 1917 PetscInt *cone, c, p; 1918 PetscErrorCode ierr; 1919 1920 PetscFunctionBegin; 1921 ierr = DMPlexSetChart(dm, 0, numCells+numVertices);CHKERRQ(ierr); 1922 for (c = 0; c < numCells; ++c) { 1923 ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr); 1924 } 1925 ierr = DMSetUp(dm);CHKERRQ(ierr); 1926 ierr = DMGetWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 1927 for (c = 0; c < numCells; ++c) { 1928 for (p = 0; p < numCorners; ++p) { 1929 cone[p] = cells[c*numCorners+p]+numCells; 1930 } 1931 ierr = DMPlexSetCone(dm, c, cone);CHKERRQ(ierr); 1932 } 1933 ierr = DMRestoreWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 1934 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 1935 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 1936 PetscFunctionReturn(0); 1937 } 1938 1939 /* 1940 This takes as input the coordinates for each vertex 1941 */ 1942 static PetscErrorCode DMPlexBuildCoordinates_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, const double vertexCoords[]) 1943 { 1944 PetscSection coordSection; 1945 Vec coordinates; 1946 DM cdm; 1947 PetscScalar *coords; 1948 PetscInt v, d; 1949 PetscErrorCode ierr; 1950 1951 PetscFunctionBegin; 1952 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1953 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1954 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 1955 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVertices);CHKERRQ(ierr); 1956 for (v = numCells; v < numCells+numVertices; ++v) { 1957 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 1958 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 1959 } 1960 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1961 1962 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 1963 ierr = DMCreateLocalVector(cdm, &coordinates);CHKERRQ(ierr); 1964 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 1965 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1966 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1967 for (v = 0; v < numVertices; ++v) { 1968 for (d = 0; d < spaceDim; ++d) { 1969 coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d]; 1970 } 1971 } 1972 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1973 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 1974 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1975 PetscFunctionReturn(0); 1976 } 1977 1978 /*@C 1979 DMPlexCreateFromCellList - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM 1980 1981 Input Parameters: 1982 + comm - The communicator 1983 . dim - The topological dimension of the mesh 1984 . numCells - The number of cells 1985 . numVertices - The number of vertices 1986 . numCorners - The number of vertices for each cell 1987 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 1988 . cells - An array of numCells*numCorners numbers, the vertices for each cell 1989 . spaceDim - The spatial dimension used for coordinates 1990 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 1991 1992 Output Parameter: 1993 . dm - The DM 1994 1995 Note: Two triangles sharing a face 1996 $ 1997 $ 2 1998 $ / | \ 1999 $ / | \ 2000 $ / | \ 2001 $ 0 0 | 1 3 2002 $ \ | / 2003 $ \ | / 2004 $ \ | / 2005 $ 1 2006 would have input 2007 $ numCells = 2, numVertices = 4 2008 $ cells = [0 1 2 1 3 2] 2009 $ 2010 which would result in the DMPlex 2011 $ 2012 $ 4 2013 $ / | \ 2014 $ / | \ 2015 $ / | \ 2016 $ 2 0 | 1 5 2017 $ \ | / 2018 $ \ | / 2019 $ \ | / 2020 $ 3 2021 2022 Level: beginner 2023 2024 .seealso: DMPlexCreateFromDAG(), DMPlexCreate() 2025 @*/ 2026 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) 2027 { 2028 PetscErrorCode ierr; 2029 2030 PetscFunctionBegin; 2031 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 2032 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 2033 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 2034 ierr = DMPlexBuildFromCellList_Private(*dm, numCells, numVertices, numCorners, cells);CHKERRQ(ierr); 2035 if (interpolate) { 2036 DM idm = NULL; 2037 2038 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 2039 ierr = DMDestroy(dm);CHKERRQ(ierr); 2040 *dm = idm; 2041 } 2042 ierr = DMPlexBuildCoordinates_Private(*dm, spaceDim, numCells, numVertices, vertexCoords);CHKERRQ(ierr); 2043 PetscFunctionReturn(0); 2044 } 2045 2046 /*@ 2047 DMPlexCreateFromDAG - This takes as input the adjacency-list representation of the Directed Acyclic Graph (Hasse Diagram) encoding a mesh, and produces a DM 2048 2049 Input Parameters: 2050 + dm - The empty DM object, usually from DMCreate() and DMSetDimension() 2051 . depth - The depth of the DAG 2052 . numPoints - The number of points at each depth 2053 . coneSize - The cone size of each point 2054 . cones - The concatenation of the cone points for each point, the cone list must be oriented correctly for each point 2055 . coneOrientations - The orientation of each cone point 2056 - vertexCoords - An array of numVertices*dim numbers, the coordinates of each vertex 2057 2058 Output Parameter: 2059 . dm - The DM 2060 2061 Note: Two triangles sharing a face would have input 2062 $ depth = 1, numPoints = [4 2], coneSize = [3 3 0 0 0 0] 2063 $ cones = [2 3 4 3 5 4], coneOrientations = [0 0 0 0 0 0] 2064 $ vertexCoords = [-1.0 0.0 0.0 -1.0 0.0 1.0 1.0 0.0] 2065 $ 2066 which would result in the DMPlex 2067 $ 2068 $ 4 2069 $ / | \ 2070 $ / | \ 2071 $ / | \ 2072 $ 2 0 | 1 5 2073 $ \ | / 2074 $ \ | / 2075 $ \ | / 2076 $ 3 2077 $ 2078 $ Notice that all points are numbered consecutively, unlikely DMPlexCreateFromCellList() 2079 2080 Level: advanced 2081 2082 .seealso: DMPlexCreateFromCellList(), DMPlexCreate() 2083 @*/ 2084 PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[]) 2085 { 2086 Vec coordinates; 2087 PetscSection coordSection; 2088 PetscScalar *coords; 2089 PetscInt coordSize, firstVertex = -1, pStart = 0, pEnd = 0, p, v, dim, dimEmbed, d, off; 2090 PetscErrorCode ierr; 2091 2092 PetscFunctionBegin; 2093 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2094 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 2095 if (dimEmbed < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Embedding dimension %d cannot be less than intrinsic dimension %d",dimEmbed,dim); 2096 for (d = 0; d <= depth; ++d) pEnd += numPoints[d]; 2097 ierr = DMPlexSetChart(dm, pStart, pEnd);CHKERRQ(ierr); 2098 for (p = pStart; p < pEnd; ++p) { 2099 ierr = DMPlexSetConeSize(dm, p, coneSize[p-pStart]);CHKERRQ(ierr); 2100 if (firstVertex < 0 && !coneSize[p - pStart]) { 2101 firstVertex = p - pStart; 2102 } 2103 } 2104 if (firstVertex < 0 && numPoints[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected %d vertices but could not find any", numPoints[0]); 2105 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 2106 for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) { 2107 ierr = DMPlexSetCone(dm, p, &cones[off]);CHKERRQ(ierr); 2108 ierr = DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);CHKERRQ(ierr); 2109 } 2110 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 2111 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 2112 /* Build coordinates */ 2113 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2114 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 2115 ierr = PetscSectionSetFieldComponents(coordSection, 0, dimEmbed);CHKERRQ(ierr); 2116 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0]);CHKERRQ(ierr); 2117 for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) { 2118 ierr = PetscSectionSetDof(coordSection, v, dimEmbed);CHKERRQ(ierr); 2119 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed);CHKERRQ(ierr); 2120 } 2121 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 2122 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 2123 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 2124 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 2125 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 2126 ierr = VecSetBlockSize(coordinates, dimEmbed);CHKERRQ(ierr); 2127 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 2128 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2129 for (v = 0; v < numPoints[0]; ++v) { 2130 PetscInt off; 2131 2132 ierr = PetscSectionGetOffset(coordSection, v+firstVertex, &off);CHKERRQ(ierr); 2133 for (d = 0; d < dimEmbed; ++d) { 2134 coords[off+d] = vertexCoords[v*dimEmbed+d]; 2135 } 2136 } 2137 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2138 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 2139 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 2140 PetscFunctionReturn(0); 2141 } 2142 2143 /*@C 2144 DMPlexCreateFromFile - This takes a filename and produces a DM 2145 2146 Input Parameters: 2147 + comm - The communicator 2148 . filename - A file name 2149 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 2150 2151 Output Parameter: 2152 . dm - The DM 2153 2154 Level: beginner 2155 2156 .seealso: DMPlexCreateFromDAG(), DMPlexCreateFromCellList(), DMPlexCreate() 2157 @*/ 2158 PetscErrorCode DMPlexCreateFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 2159 { 2160 const char *extGmsh = ".msh"; 2161 const char *extCGNS = ".cgns"; 2162 const char *extExodus = ".exo"; 2163 const char *extFluent = ".cas"; 2164 const char *extHDF5 = ".h5"; 2165 const char *extMed = ".med"; 2166 size_t len; 2167 PetscBool isGmsh, isCGNS, isExodus, isFluent, isHDF5, isMed; 2168 PetscMPIInt rank; 2169 PetscErrorCode ierr; 2170 2171 PetscFunctionBegin; 2172 PetscValidPointer(filename, 2); 2173 PetscValidPointer(dm, 4); 2174 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 2175 ierr = PetscStrlen(filename, &len);CHKERRQ(ierr); 2176 if (!len) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Filename must be a valid path"); 2177 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGmsh, 4, &isGmsh);CHKERRQ(ierr); 2178 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extCGNS, 5, &isCGNS);CHKERRQ(ierr); 2179 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extExodus, 4, &isExodus);CHKERRQ(ierr); 2180 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extFluent, 4, &isFluent);CHKERRQ(ierr); 2181 ierr = PetscStrncmp(&filename[PetscMax(0,len-3)], extHDF5, 3, &isHDF5);CHKERRQ(ierr); 2182 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extMed, 4, &isMed);CHKERRQ(ierr); 2183 if (isGmsh) { 2184 ierr = DMPlexCreateGmshFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2185 } else if (isCGNS) { 2186 ierr = DMPlexCreateCGNSFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2187 } else if (isExodus) { 2188 ierr = DMPlexCreateExodusFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2189 } else if (isFluent) { 2190 ierr = DMPlexCreateFluentFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2191 } else if (isHDF5) { 2192 PetscViewer viewer; 2193 2194 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 2195 ierr = PetscViewerSetType(viewer, PETSCVIEWERHDF5);CHKERRQ(ierr); 2196 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 2197 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 2198 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 2199 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 2200 ierr = DMLoad(*dm, viewer);CHKERRQ(ierr); 2201 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 2202 } else if (isMed) { 2203 ierr = DMPlexCreateMedFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2204 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot load file %s: unrecognized extension", filename); 2205 PetscFunctionReturn(0); 2206 } 2207 2208 /*@ 2209 DMPlexCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell 2210 2211 Collective on comm 2212 2213 Input Parameters: 2214 + comm - The communicator 2215 . dim - The spatial dimension 2216 - simplex - Flag for simplex, otherwise use a tensor-product cell 2217 2218 Output Parameter: 2219 . refdm - The reference cell 2220 2221 Level: intermediate 2222 2223 .keywords: reference cell 2224 .seealso: 2225 @*/ 2226 PetscErrorCode DMPlexCreateReferenceCell(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *refdm) 2227 { 2228 DM rdm; 2229 Vec coords; 2230 PetscErrorCode ierr; 2231 2232 PetscFunctionBeginUser; 2233 ierr = DMCreate(comm, &rdm);CHKERRQ(ierr); 2234 ierr = DMSetType(rdm, DMPLEX);CHKERRQ(ierr); 2235 ierr = DMSetDimension(rdm, dim);CHKERRQ(ierr); 2236 switch (dim) { 2237 case 0: 2238 { 2239 PetscInt numPoints[1] = {1}; 2240 PetscInt coneSize[1] = {0}; 2241 PetscInt cones[1] = {0}; 2242 PetscInt coneOrientations[1] = {0}; 2243 PetscScalar vertexCoords[1] = {0.0}; 2244 2245 ierr = DMPlexCreateFromDAG(rdm, 0, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2246 } 2247 break; 2248 case 1: 2249 { 2250 PetscInt numPoints[2] = {2, 1}; 2251 PetscInt coneSize[3] = {2, 0, 0}; 2252 PetscInt cones[2] = {1, 2}; 2253 PetscInt coneOrientations[2] = {0, 0}; 2254 PetscScalar vertexCoords[2] = {-1.0, 1.0}; 2255 2256 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2257 } 2258 break; 2259 case 2: 2260 if (simplex) { 2261 PetscInt numPoints[2] = {3, 1}; 2262 PetscInt coneSize[4] = {3, 0, 0, 0}; 2263 PetscInt cones[3] = {1, 2, 3}; 2264 PetscInt coneOrientations[3] = {0, 0, 0}; 2265 PetscScalar vertexCoords[6] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0}; 2266 2267 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2268 } else { 2269 PetscInt numPoints[2] = {4, 1}; 2270 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 2271 PetscInt cones[4] = {1, 2, 3, 4}; 2272 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 2273 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0}; 2274 2275 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2276 } 2277 break; 2278 case 3: 2279 if (simplex) { 2280 PetscInt numPoints[2] = {4, 1}; 2281 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 2282 PetscInt cones[4] = {1, 3, 2, 4}; 2283 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 2284 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}; 2285 2286 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2287 } else { 2288 PetscInt numPoints[2] = {8, 1}; 2289 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0}; 2290 PetscInt cones[8] = {1, 4, 3, 2, 5, 6, 7, 8}; 2291 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 2292 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, 2293 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0}; 2294 2295 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2296 } 2297 break; 2298 default: 2299 SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Cannot create reference cell for dimension %d", dim); 2300 } 2301 *refdm = NULL; 2302 ierr = DMPlexInterpolate(rdm, refdm);CHKERRQ(ierr); 2303 if (rdm->coordinateDM) { 2304 DM ncdm; 2305 PetscSection cs; 2306 PetscInt pEnd = -1; 2307 2308 ierr = DMGetDefaultSection(rdm->coordinateDM, &cs);CHKERRQ(ierr); 2309 if (cs) {ierr = PetscSectionGetChart(cs, NULL, &pEnd);CHKERRQ(ierr);} 2310 if (pEnd >= 0) { 2311 ierr = DMClone(rdm->coordinateDM, &ncdm);CHKERRQ(ierr); 2312 ierr = DMSetDefaultSection(ncdm, cs);CHKERRQ(ierr); 2313 ierr = DMSetCoordinateDM(*refdm, ncdm);CHKERRQ(ierr); 2314 ierr = DMDestroy(&ncdm);CHKERRQ(ierr); 2315 } 2316 } 2317 ierr = DMGetCoordinatesLocal(rdm, &coords);CHKERRQ(ierr); 2318 if (coords) { 2319 ierr = DMSetCoordinatesLocal(*refdm, coords);CHKERRQ(ierr); 2320 } else { 2321 ierr = DMGetCoordinates(rdm, &coords);CHKERRQ(ierr); 2322 if (coords) {ierr = DMSetCoordinates(*refdm, coords);CHKERRQ(ierr);} 2323 } 2324 ierr = DMDestroy(&rdm);CHKERRQ(ierr); 2325 PetscFunctionReturn(0); 2326 } 2327