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