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