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