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