1 #define PETSCDM_DLL 2 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 3 #include <petsc/private/hashseti.h> /*I "petscdmplex.h" I*/ 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; 119 120 ierr = DMPlexInterpolate(*newdm, &idm);CHKERRQ(ierr); 121 ierr = DMDestroy(newdm);CHKERRQ(ierr); 122 *newdm = idm; 123 } 124 { 125 DM refinedMesh = NULL; 126 DM distributedMesh = NULL; 127 128 /* Distribute mesh over processes */ 129 ierr = DMPlexDistribute(*newdm, 0, NULL, &distributedMesh);CHKERRQ(ierr); 130 if (distributedMesh) { 131 ierr = DMDestroy(newdm);CHKERRQ(ierr); 132 *newdm = distributedMesh; 133 } 134 if (refinementUniform) { 135 ierr = DMPlexSetRefinementUniform(*newdm, refinementUniform);CHKERRQ(ierr); 136 ierr = DMRefine(*newdm, comm, &refinedMesh);CHKERRQ(ierr); 137 if (refinedMesh) { 138 ierr = DMDestroy(newdm);CHKERRQ(ierr); 139 *newdm = refinedMesh; 140 } 141 } 142 } 143 PetscFunctionReturn(0); 144 } 145 146 /*@ 147 DMPlexCreateSquareBoundary - Creates a 1D mesh the is the boundary of a square lattice. 148 149 Collective on MPI_Comm 150 151 Input Parameters: 152 + comm - The communicator for the DM object 153 . lower - The lower left corner coordinates 154 . upper - The upper right corner coordinates 155 - edges - The number of cells in each direction 156 157 Output Parameter: 158 . dm - The DM object 159 160 Note: Here is the numbering returned for 2 cells in each direction: 161 $ 18--5-17--4--16 162 $ | | | 163 $ 6 10 3 164 $ | | | 165 $ 19-11-20--9--15 166 $ | | | 167 $ 7 8 2 168 $ | | | 169 $ 12--0-13--1--14 170 171 Level: beginner 172 173 .keywords: DM, create 174 .seealso: DMPlexCreateBoxMesh(), DMPlexCreateCubeBoundary(), DMSetType(), DMCreate() 175 @*/ 176 PetscErrorCode DMPlexCreateSquareBoundary(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[]) 177 { 178 const PetscInt numVertices = (edges[0]+1)*(edges[1]+1); 179 const PetscInt numEdges = edges[0]*(edges[1]+1) + (edges[0]+1)*edges[1]; 180 PetscInt markerTop = 1; 181 PetscInt markerBottom = 1; 182 PetscInt markerRight = 1; 183 PetscInt markerLeft = 1; 184 PetscBool markerSeparate = PETSC_FALSE; 185 Vec coordinates; 186 PetscSection coordSection; 187 PetscScalar *coords; 188 PetscInt coordSize; 189 PetscMPIInt rank; 190 PetscInt v, vx, vy; 191 PetscErrorCode ierr; 192 193 PetscFunctionBegin; 194 ierr = PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);CHKERRQ(ierr); 195 if (markerSeparate) { 196 markerTop = 3; 197 markerBottom = 1; 198 markerRight = 2; 199 markerLeft = 4; 200 } 201 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 202 if (!rank) { 203 PetscInt e, ex, ey; 204 205 ierr = DMPlexSetChart(dm, 0, numEdges+numVertices);CHKERRQ(ierr); 206 for (e = 0; e < numEdges; ++e) { 207 ierr = DMPlexSetConeSize(dm, e, 2);CHKERRQ(ierr); 208 } 209 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 210 for (vx = 0; vx <= edges[0]; vx++) { 211 for (ey = 0; ey < edges[1]; ey++) { 212 PetscInt edge = vx*edges[1] + ey + edges[0]*(edges[1]+1); 213 PetscInt vertex = ey*(edges[0]+1) + vx + numEdges; 214 PetscInt cone[2]; 215 216 cone[0] = vertex; cone[1] = vertex+edges[0]+1; 217 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 218 if (vx == edges[0]) { 219 ierr = DMSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 220 ierr = DMSetLabelValue(dm, "marker", cone[0], markerRight);CHKERRQ(ierr); 221 if (ey == edges[1]-1) { 222 ierr = DMSetLabelValue(dm, "marker", cone[1], markerRight);CHKERRQ(ierr); 223 ierr = DMSetLabelValue(dm, "Face Sets", cone[1], markerRight);CHKERRQ(ierr); 224 } 225 } else if (vx == 0) { 226 ierr = DMSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 227 ierr = DMSetLabelValue(dm, "marker", cone[0], markerLeft);CHKERRQ(ierr); 228 if (ey == edges[1]-1) { 229 ierr = DMSetLabelValue(dm, "marker", cone[1], markerLeft);CHKERRQ(ierr); 230 ierr = DMSetLabelValue(dm, "Face Sets", cone[1], markerLeft);CHKERRQ(ierr); 231 } 232 } 233 } 234 } 235 for (vy = 0; vy <= edges[1]; vy++) { 236 for (ex = 0; ex < edges[0]; ex++) { 237 PetscInt edge = vy*edges[0] + ex; 238 PetscInt vertex = vy*(edges[0]+1) + ex + numEdges; 239 PetscInt cone[2]; 240 241 cone[0] = vertex; cone[1] = vertex+1; 242 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 243 if (vy == edges[1]) { 244 ierr = DMSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 245 ierr = DMSetLabelValue(dm, "marker", cone[0], markerTop);CHKERRQ(ierr); 246 if (ex == edges[0]-1) { 247 ierr = DMSetLabelValue(dm, "marker", cone[1], markerTop);CHKERRQ(ierr); 248 ierr = DMSetLabelValue(dm, "Face Sets", cone[1], markerTop);CHKERRQ(ierr); 249 } 250 } else if (vy == 0) { 251 ierr = DMSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 252 ierr = DMSetLabelValue(dm, "marker", cone[0], markerBottom);CHKERRQ(ierr); 253 if (ex == edges[0]-1) { 254 ierr = DMSetLabelValue(dm, "marker", cone[1], markerBottom);CHKERRQ(ierr); 255 ierr = DMSetLabelValue(dm, "Face Sets", cone[1], markerBottom);CHKERRQ(ierr); 256 } 257 } 258 } 259 } 260 } 261 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 262 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 263 /* Build coordinates */ 264 ierr = DMSetCoordinateDim(dm, 2);CHKERRQ(ierr); 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 that 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 339 /* Side 0 (Top) */ 340 for (vy = 0; vy < faces[1]; vy++) { 341 for (vx = 0; vx < faces[0]; vx++) { 342 voffset = numFaces + vertices[0]*vertices[1]*(vertices[2]-1) + vy*vertices[0] + vx; 343 cone[0] = voffset; cone[1] = voffset+1; cone[2] = voffset+vertices[0]+1; cone[3] = voffset+vertices[0]; 344 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 345 ierr = DMSetLabelValue(dm, "marker", iface, 1);CHKERRQ(ierr); 346 ierr = DMSetLabelValue(dm, "marker", voffset+0, 1);CHKERRQ(ierr); 347 ierr = DMSetLabelValue(dm, "marker", voffset+1, 1);CHKERRQ(ierr); 348 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]+0, 1);CHKERRQ(ierr); 349 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]+1, 1);CHKERRQ(ierr); 350 iface++; 351 } 352 } 353 354 /* Side 1 (Bottom) */ 355 for (vy = 0; vy < faces[1]; vy++) { 356 for (vx = 0; vx < faces[0]; vx++) { 357 voffset = numFaces + vy*(faces[0]+1) + vx; 358 cone[0] = voffset+1; cone[1] = voffset; cone[2] = voffset+vertices[0]; cone[3] = voffset+vertices[0]+1; 359 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 360 ierr = DMSetLabelValue(dm, "marker", iface, 1);CHKERRQ(ierr); 361 ierr = DMSetLabelValue(dm, "marker", voffset+0, 1);CHKERRQ(ierr); 362 ierr = DMSetLabelValue(dm, "marker", voffset+1, 1);CHKERRQ(ierr); 363 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]+0, 1);CHKERRQ(ierr); 364 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]+1, 1);CHKERRQ(ierr); 365 iface++; 366 } 367 } 368 369 /* Side 2 (Front) */ 370 for (vz = 0; vz < faces[2]; vz++) { 371 for (vx = 0; vx < faces[0]; vx++) { 372 voffset = numFaces + vz*vertices[0]*vertices[1] + vx; 373 cone[0] = voffset; cone[1] = voffset+1; cone[2] = voffset+vertices[0]*vertices[1]+1; cone[3] = voffset+vertices[0]*vertices[1]; 374 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 375 ierr = DMSetLabelValue(dm, "marker", iface, 1);CHKERRQ(ierr); 376 ierr = DMSetLabelValue(dm, "marker", voffset+0, 1);CHKERRQ(ierr); 377 ierr = DMSetLabelValue(dm, "marker", voffset+1, 1);CHKERRQ(ierr); 378 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+0, 1);CHKERRQ(ierr); 379 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+1, 1);CHKERRQ(ierr); 380 iface++; 381 } 382 } 383 384 /* Side 3 (Back) */ 385 for (vz = 0; vz < faces[2]; vz++) { 386 for (vx = 0; vx < faces[0]; vx++) { 387 voffset = numFaces + vz*vertices[0]*vertices[1] + vertices[0]*(vertices[1]-1) + vx; 388 cone[0] = voffset+vertices[0]*vertices[1]; cone[1] = voffset+vertices[0]*vertices[1]+1; 389 cone[2] = voffset+1; cone[3] = voffset; 390 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 391 ierr = DMSetLabelValue(dm, "marker", iface, 1);CHKERRQ(ierr); 392 ierr = DMSetLabelValue(dm, "marker", voffset+0, 1);CHKERRQ(ierr); 393 ierr = DMSetLabelValue(dm, "marker", voffset+1, 1);CHKERRQ(ierr); 394 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+0, 1);CHKERRQ(ierr); 395 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+1, 1);CHKERRQ(ierr); 396 iface++; 397 } 398 } 399 400 /* Side 4 (Left) */ 401 for (vz = 0; vz < faces[2]; vz++) { 402 for (vy = 0; vy < faces[1]; vy++) { 403 voffset = numFaces + vz*vertices[0]*vertices[1] + vy*vertices[0]; 404 cone[0] = voffset; cone[1] = voffset+vertices[0]*vertices[1]; 405 cone[2] = voffset+vertices[0]*vertices[1]+vertices[0]; cone[3] = voffset+vertices[0]; 406 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 407 ierr = DMSetLabelValue(dm, "marker", iface, 1);CHKERRQ(ierr); 408 ierr = DMSetLabelValue(dm, "marker", voffset+0, 1);CHKERRQ(ierr); 409 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]+0, 1);CHKERRQ(ierr); 410 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[1]+0, 1);CHKERRQ(ierr); 411 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+vertices[0], 1);CHKERRQ(ierr); 412 iface++; 413 } 414 } 415 416 /* Side 5 (Right) */ 417 for (vz = 0; vz < faces[2]; vz++) { 418 for (vy = 0; vy < faces[1]; vy++) { 419 voffset = numFaces + vz*vertices[0]*vertices[1] + vy*vertices[0] + faces[0]; 420 cone[0] = voffset+vertices[0]*vertices[1]; cone[1] = voffset; 421 cone[2] = voffset+vertices[0]; cone[3] = voffset+vertices[0]*vertices[1]+vertices[0]; 422 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 423 ierr = DMSetLabelValue(dm, "marker", iface, 1);CHKERRQ(ierr); 424 ierr = DMSetLabelValue(dm, "marker", voffset+0, 1);CHKERRQ(ierr); 425 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]+0, 1);CHKERRQ(ierr); 426 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+0, 1);CHKERRQ(ierr); 427 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+vertices[0], 1);CHKERRQ(ierr); 428 iface++; 429 } 430 } 431 } 432 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 433 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 434 /* Build coordinates */ 435 ierr = DMSetCoordinateDim(dm, 3);CHKERRQ(ierr); 436 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 437 ierr = PetscSectionSetChart(coordSection, numFaces, numFaces + numVertices);CHKERRQ(ierr); 438 for (v = numFaces; v < numFaces+numVertices; ++v) { 439 ierr = PetscSectionSetDof(coordSection, v, 3);CHKERRQ(ierr); 440 } 441 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 442 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 443 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 444 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 445 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 446 ierr = VecSetBlockSize(coordinates, 3);CHKERRQ(ierr); 447 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 448 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 449 for (vz = 0; vz <= faces[2]; ++vz) { 450 for (vy = 0; vy <= faces[1]; ++vy) { 451 for (vx = 0; vx <= faces[0]; ++vx) { 452 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+0] = lower[0] + ((upper[0] - lower[0])/faces[0])*vx; 453 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+1] = lower[1] + ((upper[1] - lower[1])/faces[1])*vy; 454 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+2] = lower[2] + ((upper[2] - lower[2])/faces[2])*vz; 455 } 456 } 457 } 458 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 459 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 460 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 461 PetscFunctionReturn(0); 462 } 463 464 static PetscErrorCode DMPlexCreateLineMesh_Internal(MPI_Comm comm,PetscInt segments,PetscReal lower,PetscReal upper,DMBoundaryType bd,DM *dm) 465 { 466 PetscInt i,fStart,fEnd,numCells = 0,numVerts = 0; 467 PetscInt numPoints[2],*coneSize,*cones,*coneOrientations; 468 PetscScalar *vertexCoords; 469 PetscReal L,maxCell; 470 PetscBool markerSeparate = PETSC_FALSE; 471 PetscInt markerLeft = 1, faceMarkerLeft = 1; 472 PetscInt markerRight = 1, faceMarkerRight = 2; 473 PetscBool wrap = (bd == DM_BOUNDARY_PERIODIC || bd == DM_BOUNDARY_TWIST) ? PETSC_TRUE : PETSC_FALSE; 474 PetscMPIInt rank; 475 PetscErrorCode ierr; 476 477 PetscFunctionBegin; 478 PetscValidPointer(dm,4); 479 480 ierr = DMCreate(comm,dm);CHKERRQ(ierr); 481 ierr = DMSetType(*dm,DMPLEX);CHKERRQ(ierr); 482 ierr = DMSetDimension(*dm,1);CHKERRQ(ierr); 483 ierr = DMCreateLabel(*dm,"marker");CHKERRQ(ierr); 484 ierr = DMCreateLabel(*dm,"Face Sets");CHKERRQ(ierr); 485 486 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 487 if (!rank) numCells = segments; 488 if (!rank) numVerts = segments + (wrap ? 0 : 1); 489 490 numPoints[0] = numVerts ; numPoints[1] = numCells; 491 ierr = PetscMalloc4(numCells+numVerts,&coneSize,numCells*2,&cones,numCells+numVerts,&coneOrientations,numVerts,&vertexCoords);CHKERRQ(ierr); 492 ierr = PetscMemzero(coneOrientations,(numCells+numVerts)*sizeof(PetscInt));CHKERRQ(ierr); 493 for (i = 0; i < numCells; ++i) { coneSize[i] = 2; } 494 for (i = 0; i < numVerts; ++i) { coneSize[numCells+i] = 0; } 495 for (i = 0; i < numCells; ++i) { cones[2*i] = numCells + i%numVerts; cones[2*i+1] = numCells + (i+1)%numVerts; } 496 for (i = 0; i < numVerts; ++i) { vertexCoords[i] = lower + (upper-lower)*((PetscReal)i/(PetscReal)numCells); } 497 ierr = DMPlexCreateFromDAG(*dm,1,numPoints,coneSize,cones,coneOrientations,vertexCoords);CHKERRQ(ierr); 498 ierr = PetscFree4(coneSize,cones,coneOrientations,vertexCoords);CHKERRQ(ierr); 499 500 ierr = PetscOptionsGetBool(((PetscObject)*dm)->options,((PetscObject)*dm)->prefix,"-dm_plex_separate_marker",&markerSeparate,NULL);CHKERRQ(ierr); 501 if (markerSeparate) { markerLeft = faceMarkerLeft; markerRight = faceMarkerRight;} 502 if (!wrap && !rank) { 503 ierr = DMPlexGetHeightStratum(*dm,1,&fStart,&fEnd);CHKERRQ(ierr); 504 ierr = DMSetLabelValue(*dm,"marker",fStart,markerLeft);CHKERRQ(ierr); 505 ierr = DMSetLabelValue(*dm,"marker",fEnd-1,markerRight);CHKERRQ(ierr); 506 ierr = DMSetLabelValue(*dm,"Face Sets",fStart,faceMarkerLeft);CHKERRQ(ierr); 507 ierr = DMSetLabelValue(*dm,"Face Sets",fEnd-1,faceMarkerRight);CHKERRQ(ierr); 508 } 509 if (wrap) { 510 L = upper - lower; 511 maxCell = (PetscReal)1.1*(L/(PetscReal)PetscMax(1,segments)); 512 ierr = DMSetPeriodicity(*dm,PETSC_TRUE,&maxCell,&L,&bd);CHKERRQ(ierr); 513 } 514 PetscFunctionReturn(0); 515 } 516 517 static PetscErrorCode DMPlexCreateBoxMesh_Simplex_Internal(MPI_Comm comm, PetscInt dim, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool interpolate, DM *dm) 518 { 519 DM boundary; 520 PetscInt i; 521 PetscErrorCode ierr; 522 523 PetscFunctionBegin; 524 PetscValidPointer(dm, 4); 525 for (i = 0; i < dim; ++i) if (periodicity[i] != DM_BOUNDARY_NONE) SETERRQ(comm, PETSC_ERR_SUP, "Periodicity is not supported for simplex meshes"); 526 ierr = DMCreate(comm, &boundary);CHKERRQ(ierr); 527 PetscValidLogicalCollectiveInt(boundary,dim,2); 528 ierr = DMSetType(boundary, DMPLEX);CHKERRQ(ierr); 529 ierr = DMSetDimension(boundary, dim-1);CHKERRQ(ierr); 530 ierr = DMSetCoordinateDim(boundary, dim);CHKERRQ(ierr); 531 switch (dim) { 532 case 2: ierr = DMPlexCreateSquareBoundary(boundary, lower, upper, faces);CHKERRQ(ierr);break; 533 case 3: ierr = DMPlexCreateCubeBoundary(boundary, lower, upper, faces);CHKERRQ(ierr);break; 534 default: SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim); 535 } 536 ierr = DMPlexGenerate(boundary, NULL, interpolate, dm);CHKERRQ(ierr); 537 ierr = DMDestroy(&boundary);CHKERRQ(ierr); 538 PetscFunctionReturn(0); 539 } 540 541 static PetscErrorCode DMPlexCreateCubeMesh_Internal(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[], DMBoundaryType bdX, DMBoundaryType bdY, DMBoundaryType bdZ) 542 { 543 DMLabel cutLabel = NULL; 544 PetscInt markerTop = 1, faceMarkerTop = 1; 545 PetscInt markerBottom = 1, faceMarkerBottom = 1; 546 PetscInt markerFront = 1, faceMarkerFront = 1; 547 PetscInt markerBack = 1, faceMarkerBack = 1; 548 PetscInt markerRight = 1, faceMarkerRight = 1; 549 PetscInt markerLeft = 1, faceMarkerLeft = 1; 550 PetscInt dim; 551 PetscBool markerSeparate = PETSC_FALSE, cutMarker = PETSC_FALSE; 552 PetscMPIInt rank; 553 PetscErrorCode ierr; 554 555 PetscFunctionBegin; 556 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 557 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 558 ierr = DMCreateLabel(dm,"marker");CHKERRQ(ierr); 559 ierr = DMCreateLabel(dm,"Face Sets");CHKERRQ(ierr); 560 ierr = PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_periodic_cut", &cutMarker, NULL);CHKERRQ(ierr); 561 if (bdX == DM_BOUNDARY_PERIODIC || bdX == DM_BOUNDARY_TWIST || 562 bdY == DM_BOUNDARY_PERIODIC || bdY == DM_BOUNDARY_TWIST || 563 bdZ == DM_BOUNDARY_PERIODIC || bdZ == DM_BOUNDARY_TWIST) { 564 565 if (cutMarker) {ierr = DMCreateLabel(dm, "periodic_cut");CHKERRQ(ierr); ierr = DMGetLabel(dm, "periodic_cut", &cutLabel);CHKERRQ(ierr);} 566 } 567 switch (dim) { 568 case 2: 569 faceMarkerTop = 3; 570 faceMarkerBottom = 1; 571 faceMarkerRight = 2; 572 faceMarkerLeft = 4; 573 break; 574 case 3: 575 faceMarkerBottom = 1; 576 faceMarkerTop = 2; 577 faceMarkerFront = 3; 578 faceMarkerBack = 4; 579 faceMarkerRight = 5; 580 faceMarkerLeft = 6; 581 break; 582 default: 583 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Dimension %d not supported",dim); 584 break; 585 } 586 ierr = PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);CHKERRQ(ierr); 587 if (markerSeparate) { 588 markerBottom = faceMarkerBottom; 589 markerTop = faceMarkerTop; 590 markerFront = faceMarkerFront; 591 markerBack = faceMarkerBack; 592 markerRight = faceMarkerRight; 593 markerLeft = faceMarkerLeft; 594 } 595 { 596 const PetscInt numXEdges = !rank ? edges[0] : 0; 597 const PetscInt numYEdges = !rank ? edges[1] : 0; 598 const PetscInt numZEdges = !rank ? edges[2] : 0; 599 const PetscInt numXVertices = !rank ? (bdX == DM_BOUNDARY_PERIODIC || bdX == DM_BOUNDARY_TWIST ? edges[0] : edges[0]+1) : 0; 600 const PetscInt numYVertices = !rank ? (bdY == DM_BOUNDARY_PERIODIC || bdY == DM_BOUNDARY_TWIST ? edges[1] : edges[1]+1) : 0; 601 const PetscInt numZVertices = !rank ? (bdZ == DM_BOUNDARY_PERIODIC || bdZ == DM_BOUNDARY_TWIST ? edges[2] : edges[2]+1) : 0; 602 const PetscInt numCells = numXEdges*numYEdges*numZEdges; 603 const PetscInt numXFaces = numYEdges*numZEdges; 604 const PetscInt numYFaces = numXEdges*numZEdges; 605 const PetscInt numZFaces = numXEdges*numYEdges; 606 const PetscInt numTotXFaces = numXVertices*numXFaces; 607 const PetscInt numTotYFaces = numYVertices*numYFaces; 608 const PetscInt numTotZFaces = numZVertices*numZFaces; 609 const PetscInt numFaces = numTotXFaces + numTotYFaces + numTotZFaces; 610 const PetscInt numTotXEdges = numXEdges*numYVertices*numZVertices; 611 const PetscInt numTotYEdges = numYEdges*numXVertices*numZVertices; 612 const PetscInt numTotZEdges = numZEdges*numXVertices*numYVertices; 613 const PetscInt numVertices = numXVertices*numYVertices*numZVertices; 614 const PetscInt numEdges = numTotXEdges + numTotYEdges + numTotZEdges; 615 const PetscInt firstVertex = (dim == 2) ? numFaces : numCells; 616 const PetscInt firstXFace = (dim == 2) ? 0 : numCells + numVertices; 617 const PetscInt firstYFace = firstXFace + numTotXFaces; 618 const PetscInt firstZFace = firstYFace + numTotYFaces; 619 const PetscInt firstXEdge = numCells + numFaces + numVertices; 620 const PetscInt firstYEdge = firstXEdge + numTotXEdges; 621 const PetscInt firstZEdge = firstYEdge + numTotYEdges; 622 Vec coordinates; 623 PetscSection coordSection; 624 PetscScalar *coords; 625 PetscInt coordSize; 626 PetscInt v, vx, vy, vz; 627 PetscInt c, f, fx, fy, fz, e, ex, ey, ez; 628 629 ierr = DMPlexSetChart(dm, 0, numCells+numFaces+numEdges+numVertices);CHKERRQ(ierr); 630 for (c = 0; c < numCells; c++) { 631 ierr = DMPlexSetConeSize(dm, c, 6);CHKERRQ(ierr); 632 } 633 for (f = firstXFace; f < firstXFace+numFaces; ++f) { 634 ierr = DMPlexSetConeSize(dm, f, 4);CHKERRQ(ierr); 635 } 636 for (e = firstXEdge; e < firstXEdge+numEdges; ++e) { 637 ierr = DMPlexSetConeSize(dm, e, 2);CHKERRQ(ierr); 638 } 639 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 640 /* Build cells */ 641 for (fz = 0; fz < numZEdges; ++fz) { 642 for (fy = 0; fy < numYEdges; ++fy) { 643 for (fx = 0; fx < numXEdges; ++fx) { 644 PetscInt cell = (fz*numYEdges + fy)*numXEdges + fx; 645 PetscInt faceB = firstZFace + (fy*numXEdges+fx)*numZVertices + fz; 646 PetscInt faceT = firstZFace + (fy*numXEdges+fx)*numZVertices + ((fz+1)%numZVertices); 647 PetscInt faceF = firstYFace + (fz*numXEdges+fx)*numYVertices + fy; 648 PetscInt faceK = firstYFace + (fz*numXEdges+fx)*numYVertices + ((fy+1)%numYVertices); 649 PetscInt faceL = firstXFace + (fz*numYEdges+fy)*numXVertices + fx; 650 PetscInt faceR = firstXFace + (fz*numYEdges+fy)*numXVertices + ((fx+1)%numXVertices); 651 /* B, T, F, K, R, L */ 652 PetscInt ornt[6] = {-4, 0, 0, -1, 0, -4}; /* ??? */ 653 PetscInt cone[6]; 654 655 /* no boundary twisting in 3D */ 656 cone[0] = faceB; cone[1] = faceT; cone[2] = faceF; cone[3] = faceK; cone[4] = faceR; cone[5] = faceL; 657 ierr = DMPlexSetCone(dm, cell, cone);CHKERRQ(ierr); 658 ierr = DMPlexSetConeOrientation(dm, cell, ornt);CHKERRQ(ierr); 659 if (bdX != DM_BOUNDARY_NONE && fx == numXEdges-1 && cutLabel) {ierr = DMLabelSetValue(cutLabel, cell, 2);CHKERRQ(ierr);} 660 if (bdY != DM_BOUNDARY_NONE && fy == numYEdges-1 && cutLabel) {ierr = DMLabelSetValue(cutLabel, cell, 2);CHKERRQ(ierr);} 661 if (bdZ != DM_BOUNDARY_NONE && fz == numZEdges-1 && cutLabel) {ierr = DMLabelSetValue(cutLabel, cell, 2);CHKERRQ(ierr);} 662 } 663 } 664 } 665 /* Build x faces */ 666 for (fz = 0; fz < numZEdges; ++fz) { 667 for (fy = 0; fy < numYEdges; ++fy) { 668 for (fx = 0; fx < numXVertices; ++fx) { 669 PetscInt face = firstXFace + (fz*numYEdges+fy)*numXVertices + fx; 670 PetscInt edgeL = firstZEdge + ( fy* numXVertices+fx)*numZEdges + fz; 671 PetscInt edgeR = firstZEdge + (((fy+1)%numYVertices)*numXVertices+fx)*numZEdges + fz; 672 PetscInt edgeB = firstYEdge + ( fz* numXVertices+fx)*numYEdges + fy; 673 PetscInt edgeT = firstYEdge + (((fz+1)%numZVertices)*numXVertices+fx)*numYEdges + fy; 674 PetscInt ornt[4] = {0, 0, -2, -2}; 675 PetscInt cone[4]; 676 677 if (dim == 3) { 678 /* markers */ 679 if (bdX != DM_BOUNDARY_PERIODIC) { 680 if (fx == numXVertices-1) { 681 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerRight);CHKERRQ(ierr); 682 ierr = DMSetLabelValue(dm, "marker", face, markerRight);CHKERRQ(ierr); 683 } 684 else if (fx == 0) { 685 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerLeft);CHKERRQ(ierr); 686 ierr = DMSetLabelValue(dm, "marker", face, markerLeft);CHKERRQ(ierr); 687 } 688 } 689 } 690 cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL; 691 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr); 692 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr); 693 } 694 } 695 } 696 /* Build y faces */ 697 for (fz = 0; fz < numZEdges; ++fz) { 698 for (fx = 0; fx < numXEdges; ++fx) { 699 for (fy = 0; fy < numYVertices; ++fy) { 700 PetscInt face = firstYFace + (fz*numXEdges+fx)*numYVertices + fy; 701 PetscInt edgeL = firstZEdge + (fy*numXVertices+ fx )*numZEdges + fz; 702 PetscInt edgeR = firstZEdge + (fy*numXVertices+((fx+1)%numXVertices))*numZEdges + fz; 703 PetscInt edgeB = firstXEdge + ( fz *numYVertices+fy)*numXEdges + fx; 704 PetscInt edgeT = firstXEdge + (((fz+1)%numZVertices)*numYVertices+fy)*numXEdges + fx; 705 PetscInt ornt[4] = {0, 0, -2, -2}; 706 PetscInt cone[4]; 707 708 if (dim == 3) { 709 /* markers */ 710 if (bdY != DM_BOUNDARY_PERIODIC) { 711 if (fy == numYVertices-1) { 712 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerBack);CHKERRQ(ierr); 713 ierr = DMSetLabelValue(dm, "marker", face, markerBack);CHKERRQ(ierr); 714 } 715 else if (fy == 0) { 716 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerFront);CHKERRQ(ierr); 717 ierr = DMSetLabelValue(dm, "marker", face, markerFront);CHKERRQ(ierr); 718 } 719 } 720 } 721 cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL; 722 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr); 723 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr); 724 } 725 } 726 } 727 /* Build z faces */ 728 for (fy = 0; fy < numYEdges; ++fy) { 729 for (fx = 0; fx < numXEdges; ++fx) { 730 for (fz = 0; fz < numZVertices; fz++) { 731 PetscInt face = firstZFace + (fy*numXEdges+fx)*numZVertices + fz; 732 PetscInt edgeL = firstYEdge + (fz*numXVertices+ fx )*numYEdges + fy; 733 PetscInt edgeR = firstYEdge + (fz*numXVertices+((fx+1)%numXVertices))*numYEdges + fy; 734 PetscInt edgeB = firstXEdge + (fz*numYVertices+ fy )*numXEdges + fx; 735 PetscInt edgeT = firstXEdge + (fz*numYVertices+((fy+1)%numYVertices))*numXEdges + fx; 736 PetscInt ornt[4] = {0, 0, -2, -2}; 737 PetscInt cone[4]; 738 739 if (dim == 2) { 740 if (bdX == DM_BOUNDARY_TWIST && fx == numXEdges-1) {edgeR += numYEdges-1-2*fy; ornt[1] = -2;} 741 if (bdY == DM_BOUNDARY_TWIST && fy == numYEdges-1) {edgeT += numXEdges-1-2*fx; ornt[2] = 0;} 742 if (bdX != DM_BOUNDARY_NONE && fx == numXEdges-1 && cutLabel) {ierr = DMLabelSetValue(cutLabel, face, 2);CHKERRQ(ierr);} 743 if (bdY != DM_BOUNDARY_NONE && fy == numYEdges-1 && cutLabel) {ierr = DMLabelSetValue(cutLabel, face, 2);CHKERRQ(ierr);} 744 } else { 745 /* markers */ 746 if (bdZ != DM_BOUNDARY_PERIODIC) { 747 if (fz == numZVertices-1) { 748 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerTop);CHKERRQ(ierr); 749 ierr = DMSetLabelValue(dm, "marker", face, markerTop);CHKERRQ(ierr); 750 } 751 else if (fz == 0) { 752 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerBottom);CHKERRQ(ierr); 753 ierr = DMSetLabelValue(dm, "marker", face, markerBottom);CHKERRQ(ierr); 754 } 755 } 756 } 757 cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL; 758 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr); 759 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr); 760 } 761 } 762 } 763 /* Build Z edges*/ 764 for (vy = 0; vy < numYVertices; vy++) { 765 for (vx = 0; vx < numXVertices; vx++) { 766 for (ez = 0; ez < numZEdges; ez++) { 767 const PetscInt edge = firstZEdge + (vy*numXVertices+vx)*numZEdges + ez; 768 const PetscInt vertexB = firstVertex + ( ez *numYVertices+vy)*numXVertices + vx; 769 const PetscInt vertexT = firstVertex + (((ez+1)%numZVertices)*numYVertices+vy)*numXVertices + vx; 770 PetscInt cone[2]; 771 772 if (dim == 3) { 773 if (bdX != DM_BOUNDARY_PERIODIC) { 774 if (vx == numXVertices-1) { 775 ierr = DMSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 776 } 777 else if (vx == 0) { 778 ierr = DMSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 779 } 780 } 781 if (bdY != DM_BOUNDARY_PERIODIC) { 782 if (vy == numYVertices-1) { 783 ierr = DMSetLabelValue(dm, "marker", edge, markerBack);CHKERRQ(ierr); 784 } 785 else if (vy == 0) { 786 ierr = DMSetLabelValue(dm, "marker", edge, markerFront);CHKERRQ(ierr); 787 } 788 } 789 } 790 cone[0] = vertexB; cone[1] = vertexT; 791 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 792 } 793 } 794 } 795 /* Build Y edges*/ 796 for (vz = 0; vz < numZVertices; vz++) { 797 for (vx = 0; vx < numXVertices; vx++) { 798 for (ey = 0; ey < numYEdges; ey++) { 799 const PetscInt nextv = (dim == 2 && bdY == DM_BOUNDARY_TWIST && ey == numYEdges-1) ? (numXVertices-vx-1) : (vz*numYVertices+((ey+1)%numYVertices))*numXVertices + vx; 800 const PetscInt edge = firstYEdge + (vz*numXVertices+vx)*numYEdges + ey; 801 const PetscInt vertexF = firstVertex + (vz*numYVertices+ey)*numXVertices + vx; 802 const PetscInt vertexK = firstVertex + nextv; 803 PetscInt cone[2]; 804 805 cone[0] = vertexF; cone[1] = vertexK; 806 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 807 if (dim == 2) { 808 if ((bdX != DM_BOUNDARY_PERIODIC) && (bdX != DM_BOUNDARY_TWIST)) { 809 if (vx == numXVertices-1) { 810 ierr = DMSetLabelValue(dm, "Face Sets", edge, faceMarkerRight);CHKERRQ(ierr); 811 ierr = DMSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 812 ierr = DMSetLabelValue(dm, "marker", cone[0], markerRight);CHKERRQ(ierr); 813 if (ey == numYEdges-1) { 814 ierr = DMSetLabelValue(dm, "marker", cone[1], markerRight);CHKERRQ(ierr); 815 } 816 } else if (vx == 0) { 817 ierr = DMSetLabelValue(dm, "Face Sets", edge, faceMarkerLeft);CHKERRQ(ierr); 818 ierr = DMSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 819 ierr = DMSetLabelValue(dm, "marker", cone[0], markerLeft);CHKERRQ(ierr); 820 if (ey == numYEdges-1) { 821 ierr = DMSetLabelValue(dm, "marker", cone[1], markerLeft);CHKERRQ(ierr); 822 } 823 } 824 } else { 825 if (vx == 0 && cutLabel) { 826 ierr = DMLabelSetValue(cutLabel, edge, 1);CHKERRQ(ierr); 827 ierr = DMLabelSetValue(cutLabel, cone[0], 1);CHKERRQ(ierr); 828 if (ey == numYEdges-1) { 829 ierr = DMLabelSetValue(cutLabel, cone[1], 1);CHKERRQ(ierr); 830 } 831 } 832 } 833 } else { 834 if (bdX != DM_BOUNDARY_PERIODIC) { 835 if (vx == numXVertices-1) { 836 ierr = DMSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 837 } else if (vx == 0) { 838 ierr = DMSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 839 } 840 } 841 if (bdZ != DM_BOUNDARY_PERIODIC) { 842 if (vz == numZVertices-1) { 843 ierr = DMSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 844 } else if (vz == 0) { 845 ierr = DMSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 846 } 847 } 848 } 849 } 850 } 851 } 852 /* Build X edges*/ 853 for (vz = 0; vz < numZVertices; vz++) { 854 for (vy = 0; vy < numYVertices; vy++) { 855 for (ex = 0; ex < numXEdges; ex++) { 856 const PetscInt nextv = (dim == 2 && bdX == DM_BOUNDARY_TWIST && ex == numXEdges-1) ? (numYVertices-vy-1)*numXVertices : (vz*numYVertices+vy)*numXVertices + (ex+1)%numXVertices; 857 const PetscInt edge = firstXEdge + (vz*numYVertices+vy)*numXEdges + ex; 858 const PetscInt vertexL = firstVertex + (vz*numYVertices+vy)*numXVertices + ex; 859 const PetscInt vertexR = firstVertex + nextv; 860 PetscInt cone[2]; 861 862 cone[0] = vertexL; cone[1] = vertexR; 863 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 864 if (dim == 2) { 865 if ((bdY != DM_BOUNDARY_PERIODIC) && (bdY != DM_BOUNDARY_TWIST)) { 866 if (vy == numYVertices-1) { 867 ierr = DMSetLabelValue(dm, "Face Sets", edge, faceMarkerTop);CHKERRQ(ierr); 868 ierr = DMSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 869 ierr = DMSetLabelValue(dm, "marker", cone[0], markerTop);CHKERRQ(ierr); 870 if (ex == numXEdges-1) { 871 ierr = DMSetLabelValue(dm, "marker", cone[1], markerTop);CHKERRQ(ierr); 872 } 873 } else if (vy == 0) { 874 ierr = DMSetLabelValue(dm, "Face Sets", edge, faceMarkerBottom);CHKERRQ(ierr); 875 ierr = DMSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 876 ierr = DMSetLabelValue(dm, "marker", cone[0], markerBottom);CHKERRQ(ierr); 877 if (ex == numXEdges-1) { 878 ierr = DMSetLabelValue(dm, "marker", cone[1], markerBottom);CHKERRQ(ierr); 879 } 880 } 881 } else { 882 if (vy == 0 && cutLabel) { 883 ierr = DMLabelSetValue(cutLabel, edge, 1);CHKERRQ(ierr); 884 ierr = DMLabelSetValue(cutLabel, cone[0], 1);CHKERRQ(ierr); 885 if (ex == numXEdges-1) { 886 ierr = DMLabelSetValue(cutLabel, cone[1], 1);CHKERRQ(ierr); 887 } 888 } 889 } 890 } else { 891 if (bdY != DM_BOUNDARY_PERIODIC) { 892 if (vy == numYVertices-1) { 893 ierr = DMSetLabelValue(dm, "marker", edge, markerBack);CHKERRQ(ierr); 894 } 895 else if (vy == 0) { 896 ierr = DMSetLabelValue(dm, "marker", edge, markerFront);CHKERRQ(ierr); 897 } 898 } 899 if (bdZ != DM_BOUNDARY_PERIODIC) { 900 if (vz == numZVertices-1) { 901 ierr = DMSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 902 } 903 else if (vz == 0) { 904 ierr = DMSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 905 } 906 } 907 } 908 } 909 } 910 } 911 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 912 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 913 /* Build coordinates */ 914 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 915 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 916 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 917 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVertices);CHKERRQ(ierr); 918 for (v = firstVertex; v < firstVertex+numVertices; ++v) { 919 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 920 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 921 } 922 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 923 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 924 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 925 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 926 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 927 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 928 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 929 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 930 for (vz = 0; vz < numZVertices; ++vz) { 931 for (vy = 0; vy < numYVertices; ++vy) { 932 for (vx = 0; vx < numXVertices; ++vx) { 933 coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+0] = lower[0] + ((upper[0] - lower[0])/numXEdges)*vx; 934 coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+1] = lower[1] + ((upper[1] - lower[1])/numYEdges)*vy; 935 if (dim == 3) { 936 coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+2] = lower[2] + ((upper[2] - lower[2])/numZEdges)*vz; 937 } 938 } 939 } 940 } 941 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 942 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 943 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 944 } 945 PetscFunctionReturn(0); 946 } 947 948 static PetscErrorCode DMPlexCreateBoxMesh_Tensor_Internal(MPI_Comm comm, PetscInt dim, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool interpolate, DM *dm) 949 { 950 PetscInt i; 951 PetscErrorCode ierr; 952 953 PetscFunctionBegin; 954 PetscValidPointer(dm, 7); 955 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 956 PetscValidLogicalCollectiveInt(*dm,dim,2); 957 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 958 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 959 switch (dim) { 960 case 2: {ierr = DMPlexCreateCubeMesh_Internal(*dm, lower, upper, faces, periodicity[0], periodicity[1], DM_BOUNDARY_NONE);CHKERRQ(ierr);break;} 961 case 3: {ierr = DMPlexCreateCubeMesh_Internal(*dm, lower, upper, faces, periodicity[0], periodicity[1], periodicity[2]);CHKERRQ(ierr);break;} 962 default: SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim); 963 } 964 if (periodicity[0] == DM_BOUNDARY_PERIODIC || periodicity[0] == DM_BOUNDARY_TWIST || 965 periodicity[1] == DM_BOUNDARY_PERIODIC || periodicity[1] == DM_BOUNDARY_TWIST || 966 (dim > 2 && (periodicity[2] == DM_BOUNDARY_PERIODIC || periodicity[2] == DM_BOUNDARY_TWIST))) { 967 PetscReal L[3]; 968 PetscReal maxCell[3]; 969 970 for (i = 0; i < dim; i++) { 971 L[i] = upper[i] - lower[i]; 972 maxCell[i] = 1.1 * (L[i] / PetscMax(1,faces[i])); 973 } 974 ierr = DMSetPeriodicity(*dm,PETSC_TRUE,maxCell,L,periodicity);CHKERRQ(ierr); 975 } 976 if (!interpolate) { 977 DM udm; 978 979 ierr = DMPlexUninterpolate(*dm, &udm);CHKERRQ(ierr); 980 ierr = DMPlexCopyCoordinates(*dm, udm);CHKERRQ(ierr); 981 ierr = DMDestroy(dm);CHKERRQ(ierr); 982 *dm = udm; 983 } 984 PetscFunctionReturn(0); 985 } 986 987 /*@C 988 DMPlexCreateBoxMesh - Creates a mesh on the tensor product of unit intervals (box) using simplices or tensor cells (hexahedra). 989 990 Collective on MPI_Comm 991 992 Input Parameters: 993 + comm - The communicator for the DM object 994 . dim - The spatial dimension 995 . simplex - PETSC_TRUE for simplices, PETSC_FALSE for tensor cells 996 . faces - Number of faces per dimension, or NULL for (1,) in 1D and (2, 2) in 2D and (1, 1, 1) in 3D 997 . lower - The lower left corner, or NULL for (0, 0, 0) 998 . upper - The upper right corner, or NULL for (1, 1, 1) 999 . periodicity - The boundary type for the X,Y,Z direction, or NULL for DM_BOUNDARY_NONE 1000 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 1001 1002 Output Parameter: 1003 . dm - The DM object 1004 1005 Options Database Keys: 1006 + -dm_plex_box_lower <x,y,z> - Specify lower-left-bottom coordinates for the box 1007 - -dm_plex_box_upper <x,y,z> - Specify upper-right-top coordinates for the box 1008 1009 Note: Here is the numbering returned for 2 faces in each direction for tensor cells: 1010 $ 10---17---11---18----12 1011 $ | | | 1012 $ | | | 1013 $ 20 2 22 3 24 1014 $ | | | 1015 $ | | | 1016 $ 7---15----8---16----9 1017 $ | | | 1018 $ | | | 1019 $ 19 0 21 1 23 1020 $ | | | 1021 $ | | | 1022 $ 4---13----5---14----6 1023 1024 and for simplicial cells 1025 1026 $ 14----8---15----9----16 1027 $ |\ 5 |\ 7 | 1028 $ | \ | \ | 1029 $ 13 2 14 3 15 1030 $ | 4 \ | 6 \ | 1031 $ | \ | \ | 1032 $ 11----6---12----7----13 1033 $ |\ |\ | 1034 $ | \ 1 | \ 3 | 1035 $ 10 0 11 1 12 1036 $ | 0 \ | 2 \ | 1037 $ | \ | \ | 1038 $ 8----4----9----5----10 1039 1040 Level: beginner 1041 1042 .keywords: DM, create 1043 .seealso: DMPlexCreateFromFile(), DMPlexCreateHexCylinderMesh(), DMSetType(), DMCreate() 1044 @*/ 1045 PetscErrorCode DMPlexCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscBool simplex, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool interpolate, DM *dm) 1046 { 1047 PetscInt fac[3] = {0, 0, 0}; 1048 PetscReal low[3] = {0, 0, 0}; 1049 PetscReal upp[3] = {1, 1, 1}; 1050 DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE}; 1051 PetscInt i, n; 1052 PetscBool flg; 1053 PetscErrorCode ierr; 1054 1055 PetscFunctionBegin; 1056 for (i = 0; i < dim; ++i) fac[i] = faces ? faces[i] : (dim == 1 ? 1 : 4-dim); 1057 if (lower) for (i = 0; i < dim; ++i) low[i] = lower[i]; 1058 if (upper) for (i = 0; i < dim; ++i) upp[i] = upper[i]; 1059 if (periodicity) for (i = 0; i < dim; ++i) bdt[i] = periodicity[i]; 1060 /* Allow bounds to be specified from the command line */ 1061 n = 3; 1062 ierr = PetscOptionsGetRealArray(NULL, NULL, "-dm_plex_box_lower", low, &n, &flg);CHKERRQ(ierr); 1063 if (flg && (n != dim)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Lower box point had %D values, should have been %D", n, dim); 1064 n = 3; 1065 ierr = PetscOptionsGetRealArray(NULL, NULL, "-dm_plex_box_upper", upp, &n, &flg);CHKERRQ(ierr); 1066 if (flg && (n != dim)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Upper box point had %D values, should have been %D", n, dim); 1067 1068 if (dim == 1) {ierr = DMPlexCreateLineMesh_Internal(comm, fac[0], low[0], upp[0], bdt[0], dm);CHKERRQ(ierr);} 1069 else if (simplex) {ierr = DMPlexCreateBoxMesh_Simplex_Internal(comm, dim, fac, low, upp, bdt, interpolate, dm);CHKERRQ(ierr);} 1070 else {ierr = DMPlexCreateBoxMesh_Tensor_Internal(comm, dim, fac, low, upp, bdt, interpolate, dm);CHKERRQ(ierr);} 1071 PetscFunctionReturn(0); 1072 } 1073 1074 /*@ 1075 DMPlexCreateWedgeBoxMesh - Creates a 3-D mesh tesselating the (x,y) plane and extruding in the third direction using wedge cells. 1076 1077 Collective on MPI_Comm 1078 1079 Input Parameters: 1080 + comm - The communicator for the DM object 1081 . faces - Number of faces per dimension, or NULL for (1, 1, 1) 1082 . lower - The lower left corner, or NULL for (0, 0, 0) 1083 . upper - The upper right corner, or NULL for (1, 1, 1) 1084 . periodicity - The boundary type for the X,Y,Z direction, or NULL for DM_BOUNDARY_NONE 1085 . ordExt - If PETSC_TRUE, orders the extruded cells in the height first. Otherwise, orders the cell on the layers first 1086 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 1087 1088 Output Parameter: 1089 . dm - The DM object 1090 1091 Level: beginner 1092 1093 .keywords: DM, create 1094 .seealso: DMPlexCreateHexCylinderMesh(), DMPlexCreateWedgeCylinderMesh(), DMPlexExtrude(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 1095 @*/ 1096 PetscErrorCode DMPlexCreateWedgeBoxMesh(MPI_Comm comm, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool ordExt, PetscBool interpolate, DM *dm) 1097 { 1098 DM bdm, botdm; 1099 PetscInt i; 1100 PetscInt fac[3] = {0, 0, 0}; 1101 PetscReal low[3] = {0, 0, 0}; 1102 PetscReal upp[3] = {1, 1, 1}; 1103 DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE}; 1104 PetscErrorCode ierr; 1105 1106 PetscFunctionBegin; 1107 for (i = 0; i < 3; ++i) fac[i] = faces ? (faces[i] > 0 ? faces[i] : 1) : 1; 1108 if (lower) for (i = 0; i < 3; ++i) low[i] = lower[i]; 1109 if (upper) for (i = 0; i < 3; ++i) upp[i] = upper[i]; 1110 if (periodicity) for (i = 0; i < 3; ++i) bdt[i] = periodicity[i]; 1111 for (i = 0; i < 3; ++i) if (bdt[i] != DM_BOUNDARY_NONE) SETERRQ(comm, PETSC_ERR_SUP, "Periodicity not yet supported"); 1112 1113 ierr = DMCreate(comm, &bdm);CHKERRQ(ierr); 1114 ierr = DMSetType(bdm, DMPLEX);CHKERRQ(ierr); 1115 ierr = DMSetDimension(bdm, 1);CHKERRQ(ierr); 1116 ierr = DMSetCoordinateDim(bdm, 2);CHKERRQ(ierr); 1117 ierr = DMPlexCreateSquareBoundary(bdm, low, upp, fac);CHKERRQ(ierr); 1118 ierr = DMPlexGenerate(bdm, NULL, PETSC_FALSE, &botdm);CHKERRQ(ierr); 1119 ierr = DMDestroy(&bdm);CHKERRQ(ierr); 1120 ierr = DMPlexExtrude(botdm, fac[2], upp[2] - low[2], ordExt, interpolate, dm);CHKERRQ(ierr); 1121 if (low[2] != 0.0) { 1122 Vec v; 1123 PetscScalar *x; 1124 PetscInt cDim, n; 1125 1126 ierr = DMGetCoordinatesLocal(*dm, &v);CHKERRQ(ierr); 1127 ierr = VecGetBlockSize(v, &cDim);CHKERRQ(ierr); 1128 ierr = VecGetLocalSize(v, &n);CHKERRQ(ierr); 1129 ierr = VecGetArray(v, &x);CHKERRQ(ierr); 1130 x += cDim; 1131 for (i=0; i<n; i+=cDim) x[i] += low[2]; 1132 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1133 ierr = DMSetCoordinatesLocal(*dm, v);CHKERRQ(ierr); 1134 } 1135 ierr = DMDestroy(&botdm);CHKERRQ(ierr); 1136 PetscFunctionReturn(0); 1137 } 1138 1139 /*@ 1140 DMPlexExtrude - Creates a (d+1)-D mesh by extruding a d-D mesh in the normal direction using prismatic cells. 1141 1142 Collective on idm 1143 1144 Input Parameters: 1145 + idm - The mesh to be extruted 1146 . layers - The number of layers 1147 . height - The height of the extruded layer 1148 . ordExt - If PETSC_TRUE, orders the extruded cells in the height first. Otherwise, orders the cell on the layers first 1149 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 1150 1151 Output Parameter: 1152 . dm - The DM object 1153 1154 Notes: The object created is an hybrid mesh, the vertex ordering in the cone of the cell is that of the prismatic cells 1155 1156 Level: advanced 1157 1158 .keywords: DM, create 1159 .seealso: DMPlexCreateWedgeCylinderMesh(), DMPlexCreateWedgeBoxMesh(), DMPlexSetHybridBounds(), DMSetType(), DMCreate() 1160 @*/ 1161 PetscErrorCode DMPlexExtrude(DM idm, PetscInt layers, PetscReal height, PetscBool ordExt, PetscBool interpolate, DM* dm) 1162 { 1163 PetscScalar *coordsB; 1164 const PetscScalar *coordsA; 1165 PetscReal *normals = NULL; 1166 Vec coordinatesA, coordinatesB; 1167 PetscSection coordSectionA, coordSectionB; 1168 PetscInt dim, cDim, cDimB, c, l, v, coordSize, *newCone; 1169 PetscInt cStart, cEnd, vStart, vEnd, cellV, numCells, numVertices; 1170 PetscErrorCode ierr; 1171 1172 PetscFunctionBegin; 1173 PetscValidHeaderSpecific(idm, DM_CLASSID, 1); 1174 PetscValidLogicalCollectiveInt(idm, layers, 2); 1175 PetscValidLogicalCollectiveReal(idm, height, 3); 1176 PetscValidLogicalCollectiveBool(idm, interpolate, 4); 1177 ierr = DMGetDimension(idm, &dim);CHKERRQ(ierr); 1178 if (dim < 1 || dim > 3) SETERRQ1(PetscObjectComm((PetscObject)idm), PETSC_ERR_SUP, "Support for dimension %D not coded", dim); 1179 1180 ierr = DMPlexGetHeightStratum(idm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1181 ierr = DMPlexGetDepthStratum(idm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1182 numCells = (cEnd - cStart)*layers; 1183 numVertices = (vEnd - vStart)*(layers+1); 1184 ierr = DMCreate(PetscObjectComm((PetscObject)idm), dm);CHKERRQ(ierr); 1185 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1186 ierr = DMSetDimension(*dm, dim+1);CHKERRQ(ierr); 1187 ierr = DMPlexSetChart(*dm, 0, numCells+numVertices);CHKERRQ(ierr); 1188 for (c = cStart, cellV = 0; c < cEnd; ++c) { 1189 PetscInt *closure = NULL; 1190 PetscInt closureSize, numCorners = 0; 1191 1192 ierr = DMPlexGetTransitiveClosure(idm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1193 for (v = 0; v < closureSize*2; v += 2) if ((closure[v] >= vStart) && (closure[v] < vEnd)) numCorners++; 1194 ierr = DMPlexRestoreTransitiveClosure(idm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1195 for (l = 0; l < layers; ++l) { 1196 ierr = DMPlexSetConeSize(*dm, ordExt ? layers*(c - cStart) + l : l*(cEnd - cStart) + c - cStart, 2*numCorners);CHKERRQ(ierr); 1197 } 1198 cellV = PetscMax(numCorners,cellV); 1199 } 1200 ierr = DMPlexSetHybridBounds(*dm, 0, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 1201 ierr = DMSetUp(*dm);CHKERRQ(ierr); 1202 1203 ierr = DMGetCoordinateDim(idm, &cDim);CHKERRQ(ierr); 1204 if (dim != cDim) { 1205 ierr = PetscCalloc1(cDim*(vEnd - vStart), &normals);CHKERRQ(ierr); 1206 } 1207 ierr = PetscMalloc1(3*cellV,&newCone);CHKERRQ(ierr); 1208 for (c = cStart; c < cEnd; ++c) { 1209 PetscInt *closure = NULL; 1210 PetscInt closureSize, numCorners = 0, l; 1211 PetscReal normal[3] = {0, 0, 0}; 1212 1213 if (normals) { 1214 ierr = DMPlexComputeCellGeometryFVM(idm, c, NULL, NULL, normal);CHKERRQ(ierr); 1215 } 1216 ierr = DMPlexGetTransitiveClosure(idm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1217 for (v = 0; v < closureSize*2; v += 2) { 1218 if ((closure[v] >= vStart) && (closure[v] < vEnd)) { 1219 PetscInt d; 1220 1221 newCone[numCorners++] = closure[v] - vStart; 1222 if (normals) { for (d = 0; d < cDim; ++d) normals[cDim*(closure[v]-vStart)+d] += normal[d]; } 1223 } 1224 } 1225 switch (numCorners) { 1226 case 4: /* do nothing */ 1227 case 2: /* do nothing */ 1228 break; 1229 case 3: /* from counter-clockwise to wedge ordering */ 1230 l = newCone[1]; 1231 newCone[1] = newCone[2]; 1232 newCone[2] = l; 1233 break; 1234 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported number of corners: %D", numCorners); 1235 } 1236 ierr = DMPlexRestoreTransitiveClosure(idm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1237 for (l = 0; l < layers; ++l) { 1238 PetscInt i; 1239 1240 for (i = 0; i < numCorners; ++i) { 1241 newCone[ numCorners + i] = ordExt ? (layers+1)*newCone[i] + l + numCells : l*(vEnd - vStart) + newCone[i] + numCells; 1242 newCone[2*numCorners + i] = ordExt ? (layers+1)*newCone[i] + l + 1 + numCells : (l+1)*(vEnd - vStart) + newCone[i] + numCells; 1243 } 1244 ierr = DMPlexSetCone(*dm, ordExt ? layers*(c - cStart) + l : l*(cEnd - cStart) + c - cStart, newCone + numCorners);CHKERRQ(ierr); 1245 } 1246 } 1247 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1248 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1249 ierr = PetscFree(newCone);CHKERRQ(ierr); 1250 1251 cDimB = cDim == dim ? cDim+1 : cDim; 1252 ierr = DMGetCoordinateSection(*dm, &coordSectionB);CHKERRQ(ierr); 1253 ierr = PetscSectionSetNumFields(coordSectionB, 1);CHKERRQ(ierr); 1254 ierr = PetscSectionSetFieldComponents(coordSectionB, 0, cDimB);CHKERRQ(ierr); 1255 ierr = PetscSectionSetChart(coordSectionB, numCells, numCells+numVertices);CHKERRQ(ierr); 1256 for (v = numCells; v < numCells+numVertices; ++v) { 1257 ierr = PetscSectionSetDof(coordSectionB, v, cDimB);CHKERRQ(ierr); 1258 ierr = PetscSectionSetFieldDof(coordSectionB, v, 0, cDimB);CHKERRQ(ierr); 1259 } 1260 ierr = PetscSectionSetUp(coordSectionB);CHKERRQ(ierr); 1261 ierr = PetscSectionGetStorageSize(coordSectionB, &coordSize);CHKERRQ(ierr); 1262 ierr = VecCreate(PETSC_COMM_SELF, &coordinatesB);CHKERRQ(ierr); 1263 ierr = PetscObjectSetName((PetscObject) coordinatesB, "coordinates");CHKERRQ(ierr); 1264 ierr = VecSetSizes(coordinatesB, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1265 ierr = VecSetBlockSize(coordinatesB, cDimB);CHKERRQ(ierr); 1266 ierr = VecSetType(coordinatesB,VECSTANDARD);CHKERRQ(ierr); 1267 1268 ierr = DMGetCoordinateSection(idm, &coordSectionA);CHKERRQ(ierr); 1269 ierr = DMGetCoordinatesLocal(idm, &coordinatesA);CHKERRQ(ierr); 1270 ierr = VecGetArray(coordinatesB, &coordsB);CHKERRQ(ierr); 1271 ierr = VecGetArrayRead(coordinatesA, &coordsA);CHKERRQ(ierr); 1272 for (v = vStart; v < vEnd; ++v) { 1273 const PetscScalar *cptr; 1274 PetscReal ones2[2] = { 0., 1.}, ones3[3] = { 0., 0., 1.}; 1275 PetscReal *normal, norm, h = height/layers; 1276 PetscInt offA, d, cDimA = cDim; 1277 1278 normal = normals ? normals + cDimB*(v - vStart) : (cDim > 1 ? ones3 : ones2); 1279 if (normals) { 1280 for (d = 0, norm = 0.0; d < cDimB; ++d) norm += normal[d]*normal[d]; 1281 for (d = 0; d < cDimB; ++d) normal[d] *= 1./PetscSqrtReal(norm); 1282 } 1283 1284 ierr = PetscSectionGetOffset(coordSectionA, v, &offA);CHKERRQ(ierr); 1285 cptr = coordsA + offA; 1286 for (l = 0; l < layers+1; ++l) { 1287 PetscInt offB, d, newV; 1288 1289 newV = ordExt ? (layers+1)*(v -vStart) + l + numCells : (vEnd -vStart)*l + (v -vStart) + numCells; 1290 ierr = PetscSectionGetOffset(coordSectionB, newV, &offB);CHKERRQ(ierr); 1291 for (d = 0; d < cDimA; ++d) { coordsB[offB+d] = cptr[d]; } 1292 for (d = 0; d < cDimB; ++d) { coordsB[offB+d] += l ? normal[d]*h : 0.0; } 1293 cptr = coordsB + offB; 1294 cDimA = cDimB; 1295 } 1296 } 1297 ierr = VecRestoreArrayRead(coordinatesA, &coordsA);CHKERRQ(ierr); 1298 ierr = VecRestoreArray(coordinatesB, &coordsB);CHKERRQ(ierr); 1299 ierr = DMSetCoordinatesLocal(*dm, coordinatesB);CHKERRQ(ierr); 1300 ierr = VecDestroy(&coordinatesB);CHKERRQ(ierr); 1301 ierr = PetscFree(normals);CHKERRQ(ierr); 1302 if (interpolate) { 1303 DM idm; 1304 1305 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1306 ierr = DMPlexCopyCoordinates(*dm, idm);CHKERRQ(ierr); 1307 ierr = DMDestroy(dm);CHKERRQ(ierr); 1308 *dm = idm; 1309 } 1310 PetscFunctionReturn(0); 1311 } 1312 1313 /*@C 1314 DMPlexSetOptionsPrefix - Sets the prefix used for searching for all DM options in the database. 1315 1316 Logically Collective on DM 1317 1318 Input Parameters: 1319 + dm - the DM context 1320 - prefix - the prefix to prepend to all option names 1321 1322 Notes: 1323 A hyphen (-) must NOT be given at the beginning of the prefix name. 1324 The first character of all runtime options is AUTOMATICALLY the hyphen. 1325 1326 Level: advanced 1327 1328 .seealso: SNESSetFromOptions() 1329 @*/ 1330 PetscErrorCode DMPlexSetOptionsPrefix(DM dm, const char prefix[]) 1331 { 1332 DM_Plex *mesh = (DM_Plex *) dm->data; 1333 PetscErrorCode ierr; 1334 1335 PetscFunctionBegin; 1336 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1337 ierr = PetscObjectSetOptionsPrefix((PetscObject) dm, prefix);CHKERRQ(ierr); 1338 ierr = PetscObjectSetOptionsPrefix((PetscObject) mesh->partitioner, prefix);CHKERRQ(ierr); 1339 PetscFunctionReturn(0); 1340 } 1341 1342 /*@ 1343 DMPlexCreateHexCylinderMesh - Creates a mesh on the tensor product of the unit interval with the circle (cylinder) using hexahedra. 1344 1345 Collective on MPI_Comm 1346 1347 Input Parameters: 1348 + comm - The communicator for the DM object 1349 . numRefine - The number of regular refinements to the basic 5 cell structure 1350 - periodicZ - The boundary type for the Z direction 1351 1352 Output Parameter: 1353 . dm - The DM object 1354 1355 Note: Here is the output numbering looking from the bottom of the cylinder: 1356 $ 17-----14 1357 $ | | 1358 $ | 2 | 1359 $ | | 1360 $ 17-----8-----7-----14 1361 $ | | | | 1362 $ | 3 | 0 | 1 | 1363 $ | | | | 1364 $ 19-----5-----6-----13 1365 $ | | 1366 $ | 4 | 1367 $ | | 1368 $ 19-----13 1369 $ 1370 $ and up through the top 1371 $ 1372 $ 18-----16 1373 $ | | 1374 $ | 2 | 1375 $ | | 1376 $ 18----10----11-----16 1377 $ | | | | 1378 $ | 3 | 0 | 1 | 1379 $ | | | | 1380 $ 20-----9----12-----15 1381 $ | | 1382 $ | 4 | 1383 $ | | 1384 $ 20-----15 1385 1386 Level: beginner 1387 1388 .keywords: DM, create 1389 .seealso: DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 1390 @*/ 1391 PetscErrorCode DMPlexCreateHexCylinderMesh(MPI_Comm comm, PetscInt numRefine, DMBoundaryType periodicZ, DM *dm) 1392 { 1393 const PetscInt dim = 3; 1394 PetscInt numCells, numVertices, r; 1395 PetscMPIInt rank; 1396 PetscErrorCode ierr; 1397 1398 PetscFunctionBegin; 1399 PetscValidPointer(dm, 4); 1400 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 1401 if (numRefine < 0) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Number of refinements %D cannot be negative", numRefine); 1402 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1403 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1404 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1405 /* Create topology */ 1406 { 1407 PetscInt cone[8], c; 1408 1409 numCells = !rank ? 5 : 0; 1410 numVertices = !rank ? 16 : 0; 1411 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1412 numCells *= 3; 1413 numVertices = !rank ? 24 : 0; 1414 } 1415 ierr = DMPlexSetChart(*dm, 0, numCells+numVertices);CHKERRQ(ierr); 1416 for (c = 0; c < numCells; c++) {ierr = DMPlexSetConeSize(*dm, c, 8);CHKERRQ(ierr);} 1417 ierr = DMSetUp(*dm);CHKERRQ(ierr); 1418 if (!rank) { 1419 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1420 cone[0] = 15; cone[1] = 18; cone[2] = 17; cone[3] = 16; 1421 cone[4] = 31; cone[5] = 32; cone[6] = 33; cone[7] = 34; 1422 ierr = DMPlexSetCone(*dm, 0, cone);CHKERRQ(ierr); 1423 cone[0] = 16; cone[1] = 17; cone[2] = 24; cone[3] = 23; 1424 cone[4] = 32; cone[5] = 36; cone[6] = 37; cone[7] = 33; /* 22 25 26 21 */ 1425 ierr = DMPlexSetCone(*dm, 1, cone);CHKERRQ(ierr); 1426 cone[0] = 18; cone[1] = 27; cone[2] = 24; cone[3] = 17; 1427 cone[4] = 34; cone[5] = 33; cone[6] = 37; cone[7] = 38; 1428 ierr = DMPlexSetCone(*dm, 2, cone);CHKERRQ(ierr); 1429 cone[0] = 29; cone[1] = 27; cone[2] = 18; cone[3] = 15; 1430 cone[4] = 35; cone[5] = 31; cone[6] = 34; cone[7] = 38; 1431 ierr = DMPlexSetCone(*dm, 3, cone);CHKERRQ(ierr); 1432 cone[0] = 29; cone[1] = 15; cone[2] = 16; cone[3] = 23; 1433 cone[4] = 35; cone[5] = 36; cone[6] = 32; cone[7] = 31; 1434 ierr = DMPlexSetCone(*dm, 4, cone);CHKERRQ(ierr); 1435 1436 cone[0] = 31; cone[1] = 34; cone[2] = 33; cone[3] = 32; 1437 cone[4] = 19; cone[5] = 22; cone[6] = 21; cone[7] = 20; 1438 ierr = DMPlexSetCone(*dm, 5, cone);CHKERRQ(ierr); 1439 cone[0] = 32; cone[1] = 33; cone[2] = 37; cone[3] = 36; 1440 cone[4] = 22; cone[5] = 25; cone[6] = 26; cone[7] = 21; 1441 ierr = DMPlexSetCone(*dm, 6, cone);CHKERRQ(ierr); 1442 cone[0] = 34; cone[1] = 38; cone[2] = 37; cone[3] = 33; 1443 cone[4] = 20; cone[5] = 21; cone[6] = 26; cone[7] = 28; 1444 ierr = DMPlexSetCone(*dm, 7, cone);CHKERRQ(ierr); 1445 cone[0] = 35; cone[1] = 38; cone[2] = 34; cone[3] = 31; 1446 cone[4] = 30; cone[5] = 19; cone[6] = 20; cone[7] = 28; 1447 ierr = DMPlexSetCone(*dm, 8, cone);CHKERRQ(ierr); 1448 cone[0] = 35; cone[1] = 31; cone[2] = 32; cone[3] = 36; 1449 cone[4] = 30; cone[5] = 25; cone[6] = 22; cone[7] = 19; 1450 ierr = DMPlexSetCone(*dm, 9, cone);CHKERRQ(ierr); 1451 1452 cone[0] = 19; cone[1] = 20; cone[2] = 21; cone[3] = 22; 1453 cone[4] = 15; cone[5] = 16; cone[6] = 17; cone[7] = 18; 1454 ierr = DMPlexSetCone(*dm, 10, cone);CHKERRQ(ierr); 1455 cone[0] = 22; cone[1] = 21; cone[2] = 26; cone[3] = 25; 1456 cone[4] = 16; cone[5] = 23; cone[6] = 24; cone[7] = 17; 1457 ierr = DMPlexSetCone(*dm, 11, cone);CHKERRQ(ierr); 1458 cone[0] = 20; cone[1] = 28; cone[2] = 26; cone[3] = 21; 1459 cone[4] = 18; cone[5] = 17; cone[6] = 24; cone[7] = 27; 1460 ierr = DMPlexSetCone(*dm, 12, cone);CHKERRQ(ierr); 1461 cone[0] = 30; cone[1] = 28; cone[2] = 20; cone[3] = 19; 1462 cone[4] = 29; cone[5] = 15; cone[6] = 18; cone[7] = 27; 1463 ierr = DMPlexSetCone(*dm, 13, cone);CHKERRQ(ierr); 1464 cone[0] = 30; cone[1] = 19; cone[2] = 22; cone[3] = 25; 1465 cone[4] = 29; cone[5] = 23; cone[6] = 16; cone[7] = 15; 1466 ierr = DMPlexSetCone(*dm, 14, cone);CHKERRQ(ierr); 1467 } else { 1468 cone[0] = 5; cone[1] = 8; cone[2] = 7; cone[3] = 6; 1469 cone[4] = 9; cone[5] = 12; cone[6] = 11; cone[7] = 10; 1470 ierr = DMPlexSetCone(*dm, 0, cone);CHKERRQ(ierr); 1471 cone[0] = 6; cone[1] = 7; cone[2] = 14; cone[3] = 13; 1472 cone[4] = 12; cone[5] = 15; cone[6] = 16; cone[7] = 11; 1473 ierr = DMPlexSetCone(*dm, 1, cone);CHKERRQ(ierr); 1474 cone[0] = 8; cone[1] = 17; cone[2] = 14; cone[3] = 7; 1475 cone[4] = 10; cone[5] = 11; cone[6] = 16; cone[7] = 18; 1476 ierr = DMPlexSetCone(*dm, 2, cone);CHKERRQ(ierr); 1477 cone[0] = 19; cone[1] = 17; cone[2] = 8; cone[3] = 5; 1478 cone[4] = 20; cone[5] = 9; cone[6] = 10; cone[7] = 18; 1479 ierr = DMPlexSetCone(*dm, 3, cone);CHKERRQ(ierr); 1480 cone[0] = 19; cone[1] = 5; cone[2] = 6; cone[3] = 13; 1481 cone[4] = 20; cone[5] = 15; cone[6] = 12; cone[7] = 9; 1482 ierr = DMPlexSetCone(*dm, 4, cone);CHKERRQ(ierr); 1483 } 1484 } 1485 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1486 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1487 } 1488 /* Interpolate */ 1489 { 1490 DM idm; 1491 1492 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1493 ierr = DMDestroy(dm);CHKERRQ(ierr); 1494 *dm = idm; 1495 } 1496 /* Create cube geometry */ 1497 { 1498 Vec coordinates; 1499 PetscSection coordSection; 1500 PetscScalar *coords; 1501 PetscInt coordSize, v; 1502 const PetscReal dis = 1.0/PetscSqrtReal(2.0); 1503 const PetscReal ds2 = dis/2.0; 1504 1505 /* Build coordinates */ 1506 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 1507 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1508 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 1509 ierr = PetscSectionSetChart(coordSection, numCells, numCells+numVertices);CHKERRQ(ierr); 1510 for (v = numCells; v < numCells+numVertices; ++v) { 1511 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 1512 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 1513 } 1514 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1515 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1516 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 1517 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1518 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1519 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 1520 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1521 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1522 if (!rank) { 1523 coords[0*dim+0] = -ds2; coords[0*dim+1] = -ds2; coords[0*dim+2] = 0.0; 1524 coords[1*dim+0] = ds2; coords[1*dim+1] = -ds2; coords[1*dim+2] = 0.0; 1525 coords[2*dim+0] = ds2; coords[2*dim+1] = ds2; coords[2*dim+2] = 0.0; 1526 coords[3*dim+0] = -ds2; coords[3*dim+1] = ds2; coords[3*dim+2] = 0.0; 1527 coords[4*dim+0] = -ds2; coords[4*dim+1] = -ds2; coords[4*dim+2] = 1.0; 1528 coords[5*dim+0] = -ds2; coords[5*dim+1] = ds2; coords[5*dim+2] = 1.0; 1529 coords[6*dim+0] = ds2; coords[6*dim+1] = ds2; coords[6*dim+2] = 1.0; 1530 coords[7*dim+0] = ds2; coords[7*dim+1] = -ds2; coords[7*dim+2] = 1.0; 1531 coords[ 8*dim+0] = dis; coords[ 8*dim+1] = -dis; coords[ 8*dim+2] = 0.0; 1532 coords[ 9*dim+0] = dis; coords[ 9*dim+1] = dis; coords[ 9*dim+2] = 0.0; 1533 coords[10*dim+0] = dis; coords[10*dim+1] = -dis; coords[10*dim+2] = 1.0; 1534 coords[11*dim+0] = dis; coords[11*dim+1] = dis; coords[11*dim+2] = 1.0; 1535 coords[12*dim+0] = -dis; coords[12*dim+1] = dis; coords[12*dim+2] = 0.0; 1536 coords[13*dim+0] = -dis; coords[13*dim+1] = dis; coords[13*dim+2] = 1.0; 1537 coords[14*dim+0] = -dis; coords[14*dim+1] = -dis; coords[14*dim+2] = 0.0; 1538 coords[15*dim+0] = -dis; coords[15*dim+1] = -dis; coords[15*dim+2] = 1.0; 1539 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1540 /* 15 31 19 */ coords[16*dim+0] = -ds2; coords[16*dim+1] = -ds2; coords[16*dim+2] = 0.5; 1541 /* 16 32 22 */ coords[17*dim+0] = ds2; coords[17*dim+1] = -ds2; coords[17*dim+2] = 0.5; 1542 /* 17 33 21 */ coords[18*dim+0] = ds2; coords[18*dim+1] = ds2; coords[18*dim+2] = 0.5; 1543 /* 18 34 20 */ coords[19*dim+0] = -ds2; coords[19*dim+1] = ds2; coords[19*dim+2] = 0.5; 1544 /* 29 35 30 */ coords[20*dim+0] = -dis; coords[20*dim+1] = -dis; coords[20*dim+2] = 0.5; 1545 /* 23 36 25 */ coords[21*dim+0] = dis; coords[21*dim+1] = -dis; coords[21*dim+2] = 0.5; 1546 /* 24 37 26 */ coords[22*dim+0] = dis; coords[22*dim+1] = dis; coords[22*dim+2] = 0.5; 1547 /* 27 38 28 */ coords[23*dim+0] = -dis; coords[23*dim+1] = dis; coords[23*dim+2] = 0.5; 1548 } 1549 } 1550 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1551 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 1552 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1553 } 1554 /* Create periodicity */ 1555 if (periodicZ == DM_BOUNDARY_PERIODIC || periodicZ == DM_BOUNDARY_TWIST) { 1556 PetscReal L[3]; 1557 PetscReal maxCell[3]; 1558 DMBoundaryType bdType[3]; 1559 PetscReal lower[3] = {0.0, 0.0, 0.0}; 1560 PetscReal upper[3] = {1.0, 1.0, 1.5}; 1561 PetscInt i, numZCells = 3; 1562 1563 bdType[0] = DM_BOUNDARY_NONE; 1564 bdType[1] = DM_BOUNDARY_NONE; 1565 bdType[2] = periodicZ; 1566 for (i = 0; i < dim; i++) { 1567 L[i] = upper[i] - lower[i]; 1568 maxCell[i] = 1.1 * (L[i] / numZCells); 1569 } 1570 ierr = DMSetPeriodicity(*dm, PETSC_TRUE, maxCell, L, bdType);CHKERRQ(ierr); 1571 } 1572 /* Refine topology */ 1573 for (r = 0; r < numRefine; ++r) { 1574 DM rdm = NULL; 1575 1576 ierr = DMRefine(*dm, comm, &rdm);CHKERRQ(ierr); 1577 ierr = DMDestroy(dm);CHKERRQ(ierr); 1578 *dm = rdm; 1579 } 1580 /* Remap geometry to cylinder 1581 Interior square: Linear interpolation is correct 1582 The other cells all have vertices on rays from the origin. We want to uniformly expand the spacing 1583 such that the last vertex is on the unit circle. So the closest and farthest vertices are at distance 1584 1585 phi = arctan(y/x) 1586 d_close = sqrt(1/8 + 1/4 sin^2(phi)) 1587 d_far = sqrt(1/2 + sin^2(phi)) 1588 1589 so we remap them using 1590 1591 x_new = x_close + (x - x_close) (1 - d_close) / (d_far - d_close) 1592 y_new = y_close + (y - y_close) (1 - d_close) / (d_far - d_close) 1593 1594 If pi/4 < phi < 3pi/4 or -3pi/4 < phi < -pi/4, then we switch x and y. 1595 */ 1596 { 1597 Vec coordinates; 1598 PetscSection coordSection; 1599 PetscScalar *coords; 1600 PetscInt vStart, vEnd, v; 1601 const PetscReal dis = 1.0/PetscSqrtReal(2.0); 1602 const PetscReal ds2 = 0.5*dis; 1603 1604 ierr = DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1605 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 1606 ierr = DMGetCoordinatesLocal(*dm, &coordinates);CHKERRQ(ierr); 1607 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1608 for (v = vStart; v < vEnd; ++v) { 1609 PetscReal phi, sinp, cosp, dc, df, x, y, xc, yc; 1610 PetscInt off; 1611 1612 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr); 1613 if ((PetscAbsScalar(coords[off+0]) <= ds2) && (PetscAbsScalar(coords[off+1]) <= ds2)) continue; 1614 x = PetscRealPart(coords[off]); 1615 y = PetscRealPart(coords[off+1]); 1616 phi = PetscAtan2Real(y, x); 1617 sinp = PetscSinReal(phi); 1618 cosp = PetscCosReal(phi); 1619 if ((PetscAbsReal(phi) > PETSC_PI/4.0) && (PetscAbsReal(phi) < 3.0*PETSC_PI/4.0)) { 1620 dc = PetscAbsReal(ds2/sinp); 1621 df = PetscAbsReal(dis/sinp); 1622 xc = ds2*x/PetscAbsReal(y); 1623 yc = ds2*PetscSignReal(y); 1624 } else { 1625 dc = PetscAbsReal(ds2/cosp); 1626 df = PetscAbsReal(dis/cosp); 1627 xc = ds2*PetscSignReal(x); 1628 yc = ds2*y/PetscAbsReal(x); 1629 } 1630 coords[off+0] = xc + (coords[off+0] - xc)*(1.0 - dc)/(df - dc); 1631 coords[off+1] = yc + (coords[off+1] - yc)*(1.0 - dc)/(df - dc); 1632 } 1633 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1634 if (periodicZ == DM_BOUNDARY_PERIODIC || periodicZ == DM_BOUNDARY_TWIST) { 1635 ierr = DMLocalizeCoordinates(*dm);CHKERRQ(ierr); 1636 } 1637 } 1638 PetscFunctionReturn(0); 1639 } 1640 1641 /*@ 1642 DMPlexCreateWedgeCylinderMesh - Creates a mesh on the tensor product of the unit interval with the circle (cylinder) using wedges. 1643 1644 Collective on MPI_Comm 1645 1646 Input Parameters: 1647 + comm - The communicator for the DM object 1648 . n - The number of wedges around the origin 1649 - interpolate - Create edges and faces 1650 1651 Output Parameter: 1652 . dm - The DM object 1653 1654 Level: beginner 1655 1656 .keywords: DM, create 1657 .seealso: DMPlexCreateHexCylinderMesh(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 1658 @*/ 1659 PetscErrorCode DMPlexCreateWedgeCylinderMesh(MPI_Comm comm, PetscInt n, PetscBool interpolate, DM *dm) 1660 { 1661 const PetscInt dim = 3; 1662 PetscInt numCells, numVertices; 1663 PetscMPIInt rank; 1664 PetscErrorCode ierr; 1665 1666 PetscFunctionBegin; 1667 PetscValidPointer(dm, 3); 1668 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 1669 if (n < 0) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Number of wedges %D cannot be negative", n); 1670 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1671 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1672 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1673 /* Create topology */ 1674 { 1675 PetscInt cone[6], c; 1676 1677 numCells = !rank ? n : 0; 1678 numVertices = !rank ? 2*(n+1) : 0; 1679 ierr = DMPlexSetChart(*dm, 0, numCells+numVertices);CHKERRQ(ierr); 1680 ierr = DMPlexSetHybridBounds(*dm, 0, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 1681 for (c = 0; c < numCells; c++) {ierr = DMPlexSetConeSize(*dm, c, 6);CHKERRQ(ierr);} 1682 ierr = DMSetUp(*dm);CHKERRQ(ierr); 1683 for (c = 0; c < numCells; c++) { 1684 cone[0] = c+n*1; cone[1] = (c+1)%n+n*1; cone[2] = 0+3*n; 1685 cone[3] = c+n*2; cone[4] = (c+1)%n+n*2; cone[5] = 1+3*n; 1686 ierr = DMPlexSetCone(*dm, c, cone);CHKERRQ(ierr); 1687 } 1688 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1689 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1690 } 1691 /* Interpolate */ 1692 if (interpolate) { 1693 DM idm; 1694 1695 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1696 ierr = DMDestroy(dm);CHKERRQ(ierr); 1697 *dm = idm; 1698 } 1699 /* Create cylinder geometry */ 1700 { 1701 Vec coordinates; 1702 PetscSection coordSection; 1703 PetscScalar *coords; 1704 PetscInt coordSize, v, c; 1705 1706 /* Build coordinates */ 1707 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 1708 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1709 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 1710 ierr = PetscSectionSetChart(coordSection, numCells, numCells+numVertices);CHKERRQ(ierr); 1711 for (v = numCells; v < numCells+numVertices; ++v) { 1712 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 1713 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 1714 } 1715 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1716 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1717 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 1718 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1719 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1720 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 1721 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1722 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1723 for (c = 0; c < numCells; c++) { 1724 coords[(c+0*n)*dim+0] = PetscCosReal(2.0*c*PETSC_PI/n); coords[(c+0*n)*dim+1] = PetscSinReal(2.0*c*PETSC_PI/n); coords[(c+0*n)*dim+2] = 1.0; 1725 coords[(c+1*n)*dim+0] = PetscCosReal(2.0*c*PETSC_PI/n); coords[(c+1*n)*dim+1] = PetscSinReal(2.0*c*PETSC_PI/n); coords[(c+1*n)*dim+2] = 0.0; 1726 } 1727 if (!rank) { 1728 coords[(2*n+0)*dim+0] = 0.0; coords[(2*n+0)*dim+1] = 0.0; coords[(2*n+0)*dim+2] = 1.0; 1729 coords[(2*n+1)*dim+0] = 0.0; coords[(2*n+1)*dim+1] = 0.0; coords[(2*n+1)*dim+2] = 0.0; 1730 } 1731 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1732 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 1733 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1734 } 1735 PetscFunctionReturn(0); 1736 } 1737 1738 PETSC_STATIC_INLINE PetscReal DiffNormReal(PetscInt dim, const PetscReal x[], const PetscReal y[]) 1739 { 1740 PetscReal prod = 0.0; 1741 PetscInt i; 1742 for (i = 0; i < dim; ++i) prod += PetscSqr(x[i] - y[i]); 1743 return PetscSqrtReal(prod); 1744 } 1745 PETSC_STATIC_INLINE PetscReal DotReal(PetscInt dim, const PetscReal x[], const PetscReal y[]) 1746 { 1747 PetscReal prod = 0.0; 1748 PetscInt i; 1749 for (i = 0; i < dim; ++i) prod += x[i]*y[i]; 1750 return prod; 1751 } 1752 1753 /*@ 1754 DMPlexCreateSphereMesh - Creates a mesh on the d-dimensional sphere, S^d. 1755 1756 Collective on MPI_Comm 1757 1758 Input Parameters: 1759 + comm - The communicator for the DM object 1760 . dim - The dimension 1761 - simplex - Use simplices, or tensor product cells 1762 1763 Output Parameter: 1764 . dm - The DM object 1765 1766 Level: beginner 1767 1768 .keywords: DM, create 1769 .seealso: DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 1770 @*/ 1771 PetscErrorCode DMPlexCreateSphereMesh(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *dm) 1772 { 1773 const PetscInt embedDim = dim+1; 1774 PetscSection coordSection; 1775 Vec coordinates; 1776 PetscScalar *coords; 1777 PetscReal *coordsIn; 1778 PetscInt numCells, numEdges, numVerts, firstVertex, v, firstEdge, coordSize, d, c, e; 1779 PetscMPIInt rank; 1780 PetscErrorCode ierr; 1781 1782 PetscFunctionBegin; 1783 PetscValidPointer(dm, 4); 1784 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1785 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1786 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1787 ierr = DMSetCoordinateDim(*dm, dim+1);CHKERRQ(ierr); 1788 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) *dm), &rank);CHKERRQ(ierr); 1789 switch (dim) { 1790 case 2: 1791 if (simplex) { 1792 DM idm; 1793 const PetscReal edgeLen = 2.0/(1.0 + PETSC_PHI); 1794 const PetscReal vertex[3] = {0.0, 1.0/(1.0 + PETSC_PHI), PETSC_PHI/(1.0 + PETSC_PHI)}; 1795 const PetscInt degree = 5; 1796 PetscInt s[3] = {1, 1, 1}; 1797 PetscInt cone[3]; 1798 PetscInt *graph, p, i, j, k; 1799 1800 numCells = !rank ? 20 : 0; 1801 numVerts = !rank ? 12 : 0; 1802 firstVertex = numCells; 1803 /* Use icosahedron, which for a unit sphere has coordinates which are all cyclic permutations of 1804 1805 (0, \pm 1/\phi+1, \pm \phi/\phi+1) 1806 1807 where \phi^2 - \phi - 1 = 0, meaning \phi is the golden ratio \frac{1 + \sqrt{5}}{2}. The edge 1808 length is then given by 2/\phi = 2 * 2.73606 = 5.47214. 1809 */ 1810 /* Construct vertices */ 1811 ierr = PetscCalloc1(numVerts * embedDim, &coordsIn);CHKERRQ(ierr); 1812 for (p = 0, i = 0; p < embedDim; ++p) { 1813 for (s[1] = -1; s[1] < 2; s[1] += 2) { 1814 for (s[2] = -1; s[2] < 2; s[2] += 2) { 1815 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[(d+p)%embedDim]*vertex[(d+p)%embedDim]; 1816 ++i; 1817 } 1818 } 1819 } 1820 /* Construct graph */ 1821 ierr = PetscCalloc1(numVerts * numVerts, &graph);CHKERRQ(ierr); 1822 for (i = 0; i < numVerts; ++i) { 1823 for (j = 0, k = 0; j < numVerts; ++j) { 1824 if (PetscAbsReal(DiffNormReal(embedDim, &coordsIn[i*embedDim], &coordsIn[j*embedDim]) - edgeLen) < PETSC_SMALL) {graph[i*numVerts+j] = 1; ++k;} 1825 } 1826 if (k != degree) SETERRQ3(comm, PETSC_ERR_PLIB, "Invalid icosahedron, vertex %D degree %D != %D", i, k, degree); 1827 } 1828 /* Build Topology */ 1829 ierr = DMPlexSetChart(*dm, 0, numCells+numVerts);CHKERRQ(ierr); 1830 for (c = 0; c < numCells; c++) { 1831 ierr = DMPlexSetConeSize(*dm, c, embedDim);CHKERRQ(ierr); 1832 } 1833 ierr = DMSetUp(*dm);CHKERRQ(ierr); /* Allocate space for cones */ 1834 /* Cells */ 1835 for (i = 0, c = 0; i < numVerts; ++i) { 1836 for (j = 0; j < i; ++j) { 1837 for (k = 0; k < j; ++k) { 1838 if (graph[i*numVerts+j] && graph[j*numVerts+k] && graph[k*numVerts+i]) { 1839 cone[0] = firstVertex+i; cone[1] = firstVertex+j; cone[2] = firstVertex+k; 1840 /* Check orientation */ 1841 { 1842 const PetscInt epsilon[3][3][3] = {{{0, 0, 0}, {0, 0, 1}, {0, -1, 0}}, {{0, 0, -1}, {0, 0, 0}, {1, 0, 0}}, {{0, 1, 0}, {-1, 0, 0}, {0, 0, 0}}}; 1843 PetscReal normal[3]; 1844 PetscInt e, f; 1845 1846 for (d = 0; d < embedDim; ++d) { 1847 normal[d] = 0.0; 1848 for (e = 0; e < embedDim; ++e) { 1849 for (f = 0; f < embedDim; ++f) { 1850 normal[d] += epsilon[d][e][f]*(coordsIn[j*embedDim+e] - coordsIn[i*embedDim+e])*(coordsIn[k*embedDim+f] - coordsIn[i*embedDim+f]); 1851 } 1852 } 1853 } 1854 if (DotReal(embedDim, normal, &coordsIn[i*embedDim]) < 0) {PetscInt tmp = cone[1]; cone[1] = cone[2]; cone[2] = tmp;} 1855 } 1856 ierr = DMPlexSetCone(*dm, c++, cone);CHKERRQ(ierr); 1857 } 1858 } 1859 } 1860 } 1861 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1862 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1863 ierr = PetscFree(graph);CHKERRQ(ierr); 1864 /* Interpolate mesh */ 1865 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1866 ierr = DMDestroy(dm);CHKERRQ(ierr); 1867 *dm = idm; 1868 } else { 1869 /* 1870 12-21--13 1871 | | 1872 25 4 24 1873 | | 1874 12-25--9-16--8-24--13 1875 | | | | 1876 23 5 17 0 15 3 22 1877 | | | | 1878 10-20--6-14--7-19--11 1879 | | 1880 20 1 19 1881 | | 1882 10-18--11 1883 | | 1884 23 2 22 1885 | | 1886 12-21--13 1887 */ 1888 const PetscReal dist = 1.0/PetscSqrtReal(3.0); 1889 PetscInt cone[4], ornt[4]; 1890 1891 numCells = !rank ? 6 : 0; 1892 numEdges = !rank ? 12 : 0; 1893 numVerts = !rank ? 8 : 0; 1894 firstVertex = numCells; 1895 firstEdge = numCells + numVerts; 1896 /* Build Topology */ 1897 ierr = DMPlexSetChart(*dm, 0, numCells+numEdges+numVerts);CHKERRQ(ierr); 1898 for (c = 0; c < numCells; c++) { 1899 ierr = DMPlexSetConeSize(*dm, c, 4);CHKERRQ(ierr); 1900 } 1901 for (e = firstEdge; e < firstEdge+numEdges; ++e) { 1902 ierr = DMPlexSetConeSize(*dm, e, 2);CHKERRQ(ierr); 1903 } 1904 ierr = DMSetUp(*dm);CHKERRQ(ierr); /* Allocate space for cones */ 1905 /* Cell 0 */ 1906 cone[0] = 14; cone[1] = 15; cone[2] = 16; cone[3] = 17; 1907 ierr = DMPlexSetCone(*dm, 0, cone);CHKERRQ(ierr); 1908 ornt[0] = 0; ornt[1] = 0; ornt[2] = 0; ornt[3] = 0; 1909 ierr = DMPlexSetConeOrientation(*dm, 0, ornt);CHKERRQ(ierr); 1910 /* Cell 1 */ 1911 cone[0] = 18; cone[1] = 19; cone[2] = 14; cone[3] = 20; 1912 ierr = DMPlexSetCone(*dm, 1, cone);CHKERRQ(ierr); 1913 ornt[0] = 0; ornt[1] = 0; ornt[2] = -2; ornt[3] = 0; 1914 ierr = DMPlexSetConeOrientation(*dm, 1, ornt);CHKERRQ(ierr); 1915 /* Cell 2 */ 1916 cone[0] = 21; cone[1] = 22; cone[2] = 18; cone[3] = 23; 1917 ierr = DMPlexSetCone(*dm, 2, cone);CHKERRQ(ierr); 1918 ornt[0] = 0; ornt[1] = 0; ornt[2] = -2; ornt[3] = 0; 1919 ierr = DMPlexSetConeOrientation(*dm, 2, ornt);CHKERRQ(ierr); 1920 /* Cell 3 */ 1921 cone[0] = 19; cone[1] = 22; cone[2] = 24; cone[3] = 15; 1922 ierr = DMPlexSetCone(*dm, 3, cone);CHKERRQ(ierr); 1923 ornt[0] = -2; ornt[1] = -2; ornt[2] = 0; ornt[3] = -2; 1924 ierr = DMPlexSetConeOrientation(*dm, 3, ornt);CHKERRQ(ierr); 1925 /* Cell 4 */ 1926 cone[0] = 16; cone[1] = 24; cone[2] = 21; cone[3] = 25; 1927 ierr = DMPlexSetCone(*dm, 4, cone);CHKERRQ(ierr); 1928 ornt[0] = -2; ornt[1] = -2; ornt[2] = -2; ornt[3] = 0; 1929 ierr = DMPlexSetConeOrientation(*dm, 4, ornt);CHKERRQ(ierr); 1930 /* Cell 5 */ 1931 cone[0] = 20; cone[1] = 17; cone[2] = 25; cone[3] = 23; 1932 ierr = DMPlexSetCone(*dm, 5, cone);CHKERRQ(ierr); 1933 ornt[0] = -2; ornt[1] = -2; ornt[2] = -2; ornt[3] = -2; 1934 ierr = DMPlexSetConeOrientation(*dm, 5, ornt);CHKERRQ(ierr); 1935 /* Edges */ 1936 cone[0] = 6; cone[1] = 7; 1937 ierr = DMPlexSetCone(*dm, 14, cone);CHKERRQ(ierr); 1938 cone[0] = 7; cone[1] = 8; 1939 ierr = DMPlexSetCone(*dm, 15, cone);CHKERRQ(ierr); 1940 cone[0] = 8; cone[1] = 9; 1941 ierr = DMPlexSetCone(*dm, 16, cone);CHKERRQ(ierr); 1942 cone[0] = 9; cone[1] = 6; 1943 ierr = DMPlexSetCone(*dm, 17, cone);CHKERRQ(ierr); 1944 cone[0] = 10; cone[1] = 11; 1945 ierr = DMPlexSetCone(*dm, 18, cone);CHKERRQ(ierr); 1946 cone[0] = 11; cone[1] = 7; 1947 ierr = DMPlexSetCone(*dm, 19, cone);CHKERRQ(ierr); 1948 cone[0] = 6; cone[1] = 10; 1949 ierr = DMPlexSetCone(*dm, 20, cone);CHKERRQ(ierr); 1950 cone[0] = 12; cone[1] = 13; 1951 ierr = DMPlexSetCone(*dm, 21, cone);CHKERRQ(ierr); 1952 cone[0] = 13; cone[1] = 11; 1953 ierr = DMPlexSetCone(*dm, 22, cone);CHKERRQ(ierr); 1954 cone[0] = 10; cone[1] = 12; 1955 ierr = DMPlexSetCone(*dm, 23, cone);CHKERRQ(ierr); 1956 cone[0] = 13; cone[1] = 8; 1957 ierr = DMPlexSetCone(*dm, 24, cone);CHKERRQ(ierr); 1958 cone[0] = 12; cone[1] = 9; 1959 ierr = DMPlexSetCone(*dm, 25, cone);CHKERRQ(ierr); 1960 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1961 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1962 /* Build coordinates */ 1963 ierr = PetscCalloc1(numVerts * embedDim, &coordsIn);CHKERRQ(ierr); 1964 coordsIn[0*embedDim+0] = -dist; coordsIn[0*embedDim+1] = dist; coordsIn[0*embedDim+2] = -dist; 1965 coordsIn[1*embedDim+0] = dist; coordsIn[1*embedDim+1] = dist; coordsIn[1*embedDim+2] = -dist; 1966 coordsIn[2*embedDim+0] = dist; coordsIn[2*embedDim+1] = -dist; coordsIn[2*embedDim+2] = -dist; 1967 coordsIn[3*embedDim+0] = -dist; coordsIn[3*embedDim+1] = -dist; coordsIn[3*embedDim+2] = -dist; 1968 coordsIn[4*embedDim+0] = -dist; coordsIn[4*embedDim+1] = dist; coordsIn[4*embedDim+2] = dist; 1969 coordsIn[5*embedDim+0] = dist; coordsIn[5*embedDim+1] = dist; coordsIn[5*embedDim+2] = dist; 1970 coordsIn[6*embedDim+0] = -dist; coordsIn[6*embedDim+1] = -dist; coordsIn[6*embedDim+2] = dist; 1971 coordsIn[7*embedDim+0] = dist; coordsIn[7*embedDim+1] = -dist; coordsIn[7*embedDim+2] = dist; 1972 } 1973 break; 1974 case 3: 1975 if (simplex) { 1976 DM idm; 1977 const PetscReal edgeLen = 1.0/PETSC_PHI; 1978 const PetscReal vertexA[4] = {0.5, 0.5, 0.5, 0.5}; 1979 const PetscReal vertexB[4] = {1.0, 0.0, 0.0, 0.0}; 1980 const PetscReal vertexC[4] = {0.5, 0.5*PETSC_PHI, 0.5/PETSC_PHI, 0.0}; 1981 const PetscInt degree = 12; 1982 PetscInt s[4] = {1, 1, 1}; 1983 PetscInt evenPerm[12][4] = {{0, 1, 2, 3}, {0, 2, 3, 1}, {0, 3, 1, 2}, {1, 0, 3, 2}, {1, 2, 0, 3}, {1, 3, 2, 0}, 1984 {2, 0, 1, 3}, {2, 1, 3, 0}, {2, 3, 0, 1}, {3, 0, 2, 1}, {3, 1, 0, 2}, {3, 2, 1, 0}}; 1985 PetscInt cone[4]; 1986 PetscInt *graph, p, i, j, k, l; 1987 1988 numCells = !rank ? 600 : 0; 1989 numVerts = !rank ? 120 : 0; 1990 firstVertex = numCells; 1991 /* Use the 600-cell, which for a unit sphere has coordinates which are 1992 1993 1/2 (\pm 1, \pm 1, \pm 1, \pm 1) 16 1994 (\pm 1, 0, 0, 0) all cyclic permutations 8 1995 1/2 (\pm 1, \pm phi, \pm 1/phi, 0) all even permutations 96 1996 1997 where \phi^2 - \phi - 1 = 0, meaning \phi is the golden ratio \frac{1 + \sqrt{5}}{2}. The edge 1998 length is then given by 1/\phi = 2.73606. 1999 2000 http://buzzard.pugetsound.edu/sage-practice/ch03s03.html 2001 http://mathworld.wolfram.com/600-Cell.html 2002 */ 2003 /* Construct vertices */ 2004 ierr = PetscCalloc1(numVerts * embedDim, &coordsIn);CHKERRQ(ierr); 2005 i = 0; 2006 for (s[0] = -1; s[0] < 2; s[0] += 2) { 2007 for (s[1] = -1; s[1] < 2; s[1] += 2) { 2008 for (s[2] = -1; s[2] < 2; s[2] += 2) { 2009 for (s[3] = -1; s[3] < 2; s[3] += 2) { 2010 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[d]*vertexA[d]; 2011 ++i; 2012 } 2013 } 2014 } 2015 } 2016 for (p = 0; p < embedDim; ++p) { 2017 s[1] = s[2] = s[3] = 1; 2018 for (s[0] = -1; s[0] < 2; s[0] += 2) { 2019 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[(d+p)%embedDim]*vertexB[(d+p)%embedDim]; 2020 ++i; 2021 } 2022 } 2023 for (p = 0; p < 12; ++p) { 2024 s[3] = 1; 2025 for (s[0] = -1; s[0] < 2; s[0] += 2) { 2026 for (s[1] = -1; s[1] < 2; s[1] += 2) { 2027 for (s[2] = -1; s[2] < 2; s[2] += 2) { 2028 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[evenPerm[p][d]]*vertexC[evenPerm[p][d]]; 2029 ++i; 2030 } 2031 } 2032 } 2033 } 2034 if (i != numVerts) SETERRQ2(comm, PETSC_ERR_PLIB, "Invalid 600-cell, vertices %D != %D", i, numVerts); 2035 /* Construct graph */ 2036 ierr = PetscCalloc1(numVerts * numVerts, &graph);CHKERRQ(ierr); 2037 for (i = 0; i < numVerts; ++i) { 2038 for (j = 0, k = 0; j < numVerts; ++j) { 2039 if (PetscAbsReal(DiffNormReal(embedDim, &coordsIn[i*embedDim], &coordsIn[j*embedDim]) - edgeLen) < PETSC_SMALL) {graph[i*numVerts+j] = 1; ++k;} 2040 } 2041 if (k != degree) SETERRQ3(comm, PETSC_ERR_PLIB, "Invalid 600-cell, vertex %D degree %D != %D", i, k, degree); 2042 } 2043 /* Build Topology */ 2044 ierr = DMPlexSetChart(*dm, 0, numCells+numVerts);CHKERRQ(ierr); 2045 for (c = 0; c < numCells; c++) { 2046 ierr = DMPlexSetConeSize(*dm, c, embedDim);CHKERRQ(ierr); 2047 } 2048 ierr = DMSetUp(*dm);CHKERRQ(ierr); /* Allocate space for cones */ 2049 /* Cells */ 2050 for (i = 0, c = 0; i < numVerts; ++i) { 2051 for (j = 0; j < i; ++j) { 2052 for (k = 0; k < j; ++k) { 2053 for (l = 0; l < k; ++l) { 2054 if (graph[i*numVerts+j] && graph[j*numVerts+k] && graph[k*numVerts+i] && 2055 graph[l*numVerts+i] && graph[l*numVerts+j] && graph[l*numVerts+k]) { 2056 cone[0] = firstVertex+i; cone[1] = firstVertex+j; cone[2] = firstVertex+k; cone[3] = firstVertex+l; 2057 /* Check orientation: https://ef.gy/linear-algebra:normal-vectors-in-higher-dimensional-spaces */ 2058 { 2059 const PetscInt epsilon[4][4][4][4] = {{{{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 2060 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 1}, { 0, 0, -1, 0}}, 2061 {{0, 0, 0, 0}, { 0, 0, 0, -1}, { 0, 0, 0, 0}, { 0, 1, 0, 0}}, 2062 {{0, 0, 0, 0}, { 0, 0, 1, 0}, { 0, -1, 0, 0}, { 0, 0, 0, 0}}}, 2063 2064 {{{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, -1}, { 0, 0, 1, 0}}, 2065 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 2066 {{0, 0, 0, 1}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, {-1, 0, 0, 0}}, 2067 {{0, 0, -1, 0}, { 0, 0, 0, 0}, { 1, 0, 0, 0}, { 0, 0, 0, 0}}}, 2068 2069 {{{0, 0, 0, 0}, { 0, 0, 0, 1}, { 0, 0, 0, 0}, { 0, -1, 0, 0}}, 2070 {{0, 0, 0, -1}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 1, 0, 0, 0}}, 2071 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 2072 {{0, 1, 0, 0}, {-1, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}}, 2073 2074 {{{0, 0, 0, 0}, { 0, 0, -1, 0}, { 0, 1, 0, 0}, { 0, 0, 0, 0}}, 2075 {{0, 0, 1, 0}, { 0, 0, 0, 0}, {-1, 0, 0, 0}, { 0, 0, 0, 0}}, 2076 {{0, -1, 0, 0}, { 1, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 2077 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}}}; 2078 PetscReal normal[4]; 2079 PetscInt e, f, g; 2080 2081 for (d = 0; d < embedDim; ++d) { 2082 normal[d] = 0.0; 2083 for (e = 0; e < embedDim; ++e) { 2084 for (f = 0; f < embedDim; ++f) { 2085 for (g = 0; g < embedDim; ++g) { 2086 normal[d] += epsilon[d][e][f][g]*(coordsIn[j*embedDim+e] - coordsIn[i*embedDim+e])*(coordsIn[k*embedDim+f] - coordsIn[i*embedDim+f])*(coordsIn[l*embedDim+f] - coordsIn[i*embedDim+f]); 2087 } 2088 } 2089 } 2090 } 2091 if (DotReal(embedDim, normal, &coordsIn[i*embedDim]) < 0) {PetscInt tmp = cone[1]; cone[1] = cone[2]; cone[2] = tmp;} 2092 } 2093 ierr = DMPlexSetCone(*dm, c++, cone);CHKERRQ(ierr); 2094 } 2095 } 2096 } 2097 } 2098 } 2099 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 2100 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 2101 ierr = PetscFree(graph);CHKERRQ(ierr); 2102 /* Interpolate mesh */ 2103 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 2104 ierr = DMDestroy(dm);CHKERRQ(ierr); 2105 *dm = idm; 2106 break; 2107 } 2108 default: SETERRQ1(comm, PETSC_ERR_SUP, "Unsupported dimension for sphere: %D", dim); 2109 } 2110 /* Create coordinates */ 2111 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 2112 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 2113 ierr = PetscSectionSetFieldComponents(coordSection, 0, embedDim);CHKERRQ(ierr); 2114 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVerts);CHKERRQ(ierr); 2115 for (v = firstVertex; v < firstVertex+numVerts; ++v) { 2116 ierr = PetscSectionSetDof(coordSection, v, embedDim);CHKERRQ(ierr); 2117 ierr = PetscSectionSetFieldDof(coordSection, v, 0, embedDim);CHKERRQ(ierr); 2118 } 2119 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 2120 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 2121 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 2122 ierr = VecSetBlockSize(coordinates, embedDim);CHKERRQ(ierr); 2123 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 2124 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 2125 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 2126 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2127 for (v = 0; v < numVerts; ++v) for (d = 0; d < embedDim; ++d) {coords[v*embedDim+d] = coordsIn[v*embedDim+d];} 2128 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2129 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 2130 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 2131 ierr = PetscFree(coordsIn);CHKERRQ(ierr); 2132 PetscFunctionReturn(0); 2133 } 2134 2135 /* External function declarations here */ 2136 extern PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling); 2137 extern PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, Mat *mat); 2138 extern PetscErrorCode DMCreateMassMatrix_Plex(DM dmCoarse, DM dmFine, Mat *mat); 2139 extern PetscErrorCode DMCreateDefaultSection_Plex(DM dm); 2140 extern PetscErrorCode DMCreateDefaultConstraints_Plex(DM dm); 2141 extern PetscErrorCode DMCreateMatrix_Plex(DM dm, Mat *J); 2142 extern PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm); 2143 extern PetscErrorCode DMCreateCoordinateField_Plex(DM dm, DMField *field); 2144 PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm); 2145 extern PetscErrorCode DMSetUp_Plex(DM dm); 2146 extern PetscErrorCode DMDestroy_Plex(DM dm); 2147 extern PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer); 2148 extern PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer); 2149 extern PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm); 2150 extern PetscErrorCode DMCreateSuperDM_Plex(DM dms[], PetscInt len, IS **is, DM *superdm); 2151 static PetscErrorCode DMInitialize_Plex(DM dm); 2152 2153 /* Replace dm with the contents of dmNew 2154 - Share the DM_Plex structure 2155 - Share the coordinates 2156 - Share the SF 2157 */ 2158 static PetscErrorCode DMPlexReplace_Static(DM dm, DM dmNew) 2159 { 2160 PetscSF sf; 2161 DM coordDM, coarseDM; 2162 Vec coords; 2163 PetscBool isper; 2164 const PetscReal *maxCell, *L; 2165 const DMBoundaryType *bd; 2166 PetscErrorCode ierr; 2167 2168 PetscFunctionBegin; 2169 ierr = DMGetPointSF(dmNew, &sf);CHKERRQ(ierr); 2170 ierr = DMSetPointSF(dm, sf);CHKERRQ(ierr); 2171 ierr = DMGetCoordinateDM(dmNew, &coordDM);CHKERRQ(ierr); 2172 ierr = DMGetCoordinatesLocal(dmNew, &coords);CHKERRQ(ierr); 2173 ierr = DMSetCoordinateDM(dm, coordDM);CHKERRQ(ierr); 2174 ierr = DMSetCoordinatesLocal(dm, coords);CHKERRQ(ierr); 2175 ierr = DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);CHKERRQ(ierr); 2176 ierr = DMSetPeriodicity(dmNew, isper, maxCell, L, bd);CHKERRQ(ierr); 2177 ierr = DMDestroy_Plex(dm);CHKERRQ(ierr); 2178 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 2179 dm->data = dmNew->data; 2180 ((DM_Plex *) dmNew->data)->refct++; 2181 dmNew->labels->refct++; 2182 if (!--(dm->labels->refct)) { 2183 DMLabelLink next = dm->labels->next; 2184 2185 /* destroy the labels */ 2186 while (next) { 2187 DMLabelLink tmp = next->next; 2188 2189 ierr = DMLabelDestroy(&next->label);CHKERRQ(ierr); 2190 ierr = PetscFree(next);CHKERRQ(ierr); 2191 next = tmp; 2192 } 2193 ierr = PetscFree(dm->labels);CHKERRQ(ierr); 2194 } 2195 dm->labels = dmNew->labels; 2196 dm->depthLabel = dmNew->depthLabel; 2197 ierr = DMGetCoarseDM(dmNew,&coarseDM);CHKERRQ(ierr); 2198 ierr = DMSetCoarseDM(dm,coarseDM);CHKERRQ(ierr); 2199 PetscFunctionReturn(0); 2200 } 2201 2202 /* Swap dm with the contents of dmNew 2203 - Swap the DM_Plex structure 2204 - Swap the coordinates 2205 - Swap the point PetscSF 2206 */ 2207 static PetscErrorCode DMPlexSwap_Static(DM dmA, DM dmB) 2208 { 2209 DM coordDMA, coordDMB; 2210 Vec coordsA, coordsB; 2211 PetscSF sfA, sfB; 2212 void *tmp; 2213 DMLabelLinkList listTmp; 2214 DMLabel depthTmp; 2215 PetscInt tmpI; 2216 PetscErrorCode ierr; 2217 2218 PetscFunctionBegin; 2219 ierr = DMGetPointSF(dmA, &sfA);CHKERRQ(ierr); 2220 ierr = DMGetPointSF(dmB, &sfB);CHKERRQ(ierr); 2221 ierr = PetscObjectReference((PetscObject) sfA);CHKERRQ(ierr); 2222 ierr = DMSetPointSF(dmA, sfB);CHKERRQ(ierr); 2223 ierr = DMSetPointSF(dmB, sfA);CHKERRQ(ierr); 2224 ierr = PetscObjectDereference((PetscObject) sfA);CHKERRQ(ierr); 2225 2226 ierr = DMGetCoordinateDM(dmA, &coordDMA);CHKERRQ(ierr); 2227 ierr = DMGetCoordinateDM(dmB, &coordDMB);CHKERRQ(ierr); 2228 ierr = PetscObjectReference((PetscObject) coordDMA);CHKERRQ(ierr); 2229 ierr = DMSetCoordinateDM(dmA, coordDMB);CHKERRQ(ierr); 2230 ierr = DMSetCoordinateDM(dmB, coordDMA);CHKERRQ(ierr); 2231 ierr = PetscObjectDereference((PetscObject) coordDMA);CHKERRQ(ierr); 2232 2233 ierr = DMGetCoordinatesLocal(dmA, &coordsA);CHKERRQ(ierr); 2234 ierr = DMGetCoordinatesLocal(dmB, &coordsB);CHKERRQ(ierr); 2235 ierr = PetscObjectReference((PetscObject) coordsA);CHKERRQ(ierr); 2236 ierr = DMSetCoordinatesLocal(dmA, coordsB);CHKERRQ(ierr); 2237 ierr = DMSetCoordinatesLocal(dmB, coordsA);CHKERRQ(ierr); 2238 ierr = PetscObjectDereference((PetscObject) coordsA);CHKERRQ(ierr); 2239 2240 tmp = dmA->data; 2241 dmA->data = dmB->data; 2242 dmB->data = tmp; 2243 listTmp = dmA->labels; 2244 dmA->labels = dmB->labels; 2245 dmB->labels = listTmp; 2246 depthTmp = dmA->depthLabel; 2247 dmA->depthLabel = dmB->depthLabel; 2248 dmB->depthLabel = depthTmp; 2249 tmpI = dmA->levelup; 2250 dmA->levelup = dmB->levelup; 2251 dmB->levelup = tmpI; 2252 PetscFunctionReturn(0); 2253 } 2254 2255 static PetscErrorCode DMPlexIsSimplex_Static(DM dm, PetscBool *isSimplex) 2256 { 2257 PetscInt dim, cStart, coneSize; 2258 PetscErrorCode ierr; 2259 2260 PetscFunctionBegin; 2261 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2262 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, NULL);CHKERRQ(ierr); 2263 ierr = DMPlexGetConeSize(dm, cStart, &coneSize);CHKERRQ(ierr); 2264 *isSimplex = coneSize == dim+1 ? PETSC_TRUE : PETSC_FALSE; 2265 PetscFunctionReturn(0); 2266 } 2267 2268 PetscErrorCode DMSetFromOptions_NonRefinement_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 2269 { 2270 DM_Plex *mesh = (DM_Plex*) dm->data; 2271 PetscErrorCode ierr; 2272 2273 PetscFunctionBegin; 2274 /* Handle viewing */ 2275 ierr = PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMPlexMatSetClosure", PETSC_FALSE, &mesh->printSetValues, NULL);CHKERRQ(ierr); 2276 ierr = PetscOptionsInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMPlexSNESComputeResidualFEM", 0, &mesh->printFEM, NULL);CHKERRQ(ierr); 2277 ierr = PetscOptionsReal("-dm_plex_print_tol", "Tolerance for FEM output", "DMPlexSNESComputeResidualFEM", mesh->printTol, &mesh->printTol, NULL);CHKERRQ(ierr); 2278 ierr = PetscOptionsInt("-dm_plex_print_l2", "Debug output level all L2 diff computations", "DMComputeL2Diff", 0, &mesh->printL2, NULL);CHKERRQ(ierr); 2279 /* Point Location */ 2280 ierr = PetscOptionsBool("-dm_plex_hash_location", "Use grid hashing for point location", "DMInterpolate", PETSC_FALSE, &mesh->useHashLocation, NULL);CHKERRQ(ierr); 2281 /* Partitioning and distribution */ 2282 ierr = PetscOptionsBool("-dm_plex_partition_balance", "Attempt to evenly divide points on partition boundary between processes", "DMPlexSetPartitionBalance", PETSC_FALSE, &mesh->partitionBalance, NULL);CHKERRQ(ierr); 2283 /* Generation and remeshing */ 2284 ierr = PetscOptionsBool("-dm_plex_remesh_bd", "Allow changes to the boundary on remeshing", "DMAdapt", PETSC_FALSE, &mesh->remeshBd, NULL);CHKERRQ(ierr); 2285 /* Projection behavior */ 2286 ierr = PetscOptionsInt("-dm_plex_max_projection_height", "Maxmimum mesh point height used to project locally", "DMPlexSetMaxProjectionHeight", 0, &mesh->maxProjectionHeight, NULL);CHKERRQ(ierr); 2287 ierr = PetscOptionsBool("-dm_plex_regular_refinement", "Use special nested projection algorithm for regular refinement", "DMPlexSetRegularRefinement", mesh->regularRefinement, &mesh->regularRefinement, NULL);CHKERRQ(ierr); 2288 /* Checking structure */ 2289 { 2290 const char *cellTypes[] = {"simplex", "tensor", "unknown", "DMPlexCellType", "DM_PLEX_CELLTYPE_", NULL}; 2291 PetscEnum ct; 2292 PetscBool flg = PETSC_FALSE, flg2 = PETSC_FALSE; 2293 2294 ierr = PetscOptionsBool("-dm_plex_check_symmetry", "Check that the adjacency information in the mesh is symmetric", "DMPlexCheckSymmetry", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr); 2295 if (flg && flg2) {ierr = DMPlexCheckSymmetry(dm);CHKERRQ(ierr);} 2296 ierr = PetscOptionsEnum("-dm_plex_check_skeleton", "Check that each cell has the correct number of vertices", "DMPlexCheckSkeleton", cellTypes, (PetscEnum) 0, &ct, &flg);CHKERRQ(ierr); 2297 if (flg) { 2298 PetscBool isSimplex = ct ? PETSC_FALSE : PETSC_TRUE; 2299 2300 if (ct == (PetscEnum) DM_PLEX_CELLTYPE_UNKNOWN) {ierr = DMPlexIsSimplex_Static(dm, &isSimplex);CHKERRQ(ierr);} 2301 ierr = DMPlexCheckSkeleton(dm, isSimplex, 0);CHKERRQ(ierr); 2302 } 2303 ierr = PetscOptionsEnum("-dm_plex_check_faces", "Check that the faces of each cell give a vertex order this is consistent with what we expect from the cell type", "DMPlexCheckFaces", cellTypes, (PetscEnum) 0, &ct, &flg);CHKERRQ(ierr); 2304 if (flg) { 2305 PetscBool isSimplex = ct ? PETSC_FALSE : PETSC_TRUE; 2306 2307 if (ct == (PetscEnum) DM_PLEX_CELLTYPE_UNKNOWN) {ierr = DMPlexIsSimplex_Static(dm, &isSimplex);CHKERRQ(ierr);} 2308 ierr = DMPlexCheckFaces(dm, isSimplex, 0);CHKERRQ(ierr); 2309 } 2310 ierr = PetscOptionsBool("-dm_plex_check_geometry", "Check that cells have positive volume", "DMPlexCheckGeometry", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr); 2311 if (flg && flg2) {ierr = DMPlexCheckGeometry(dm);CHKERRQ(ierr);} 2312 } 2313 2314 ierr = PetscPartitionerSetFromOptions(mesh->partitioner);CHKERRQ(ierr); 2315 PetscFunctionReturn(0); 2316 } 2317 2318 static PetscErrorCode DMSetFromOptions_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 2319 { 2320 PetscInt refine = 0, coarsen = 0, r; 2321 PetscBool isHierarchy; 2322 PetscErrorCode ierr; 2323 2324 PetscFunctionBegin; 2325 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2326 ierr = PetscOptionsHead(PetscOptionsObject,"DMPlex Options");CHKERRQ(ierr); 2327 /* Handle DMPlex refinement */ 2328 ierr = PetscOptionsInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL);CHKERRQ(ierr); 2329 ierr = PetscOptionsInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy);CHKERRQ(ierr); 2330 if (refine) {ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE);CHKERRQ(ierr);} 2331 if (refine && isHierarchy) { 2332 DM *dms, coarseDM; 2333 2334 ierr = DMGetCoarseDM(dm, &coarseDM);CHKERRQ(ierr); 2335 ierr = PetscObjectReference((PetscObject)coarseDM);CHKERRQ(ierr); 2336 ierr = PetscMalloc1(refine,&dms);CHKERRQ(ierr); 2337 ierr = DMRefineHierarchy(dm, refine, dms);CHKERRQ(ierr); 2338 /* Total hack since we do not pass in a pointer */ 2339 ierr = DMPlexSwap_Static(dm, dms[refine-1]);CHKERRQ(ierr); 2340 if (refine == 1) { 2341 ierr = DMSetCoarseDM(dm, dms[0]);CHKERRQ(ierr); 2342 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr); 2343 } else { 2344 ierr = DMSetCoarseDM(dm, dms[refine-2]);CHKERRQ(ierr); 2345 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr); 2346 ierr = DMSetCoarseDM(dms[0], dms[refine-1]);CHKERRQ(ierr); 2347 ierr = DMPlexSetRegularRefinement(dms[0], PETSC_TRUE);CHKERRQ(ierr); 2348 } 2349 ierr = DMSetCoarseDM(dms[refine-1], coarseDM);CHKERRQ(ierr); 2350 ierr = PetscObjectDereference((PetscObject)coarseDM);CHKERRQ(ierr); 2351 /* Free DMs */ 2352 for (r = 0; r < refine; ++r) { 2353 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr); 2354 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 2355 } 2356 ierr = PetscFree(dms);CHKERRQ(ierr); 2357 } else { 2358 for (r = 0; r < refine; ++r) { 2359 DM refinedMesh; 2360 2361 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2362 ierr = DMRefine(dm, PetscObjectComm((PetscObject) dm), &refinedMesh);CHKERRQ(ierr); 2363 /* Total hack since we do not pass in a pointer */ 2364 ierr = DMPlexReplace_Static(dm, refinedMesh);CHKERRQ(ierr); 2365 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2366 ierr = DMDestroy(&refinedMesh);CHKERRQ(ierr); 2367 } 2368 } 2369 /* Handle DMPlex coarsening */ 2370 ierr = PetscOptionsInt("-dm_coarsen", "Coarsen the mesh", "DMCreate", coarsen, &coarsen, NULL);CHKERRQ(ierr); 2371 ierr = PetscOptionsInt("-dm_coarsen_hierarchy", "The number of coarsenings", "DMCreate", coarsen, &coarsen, &isHierarchy);CHKERRQ(ierr); 2372 if (coarsen && isHierarchy) { 2373 DM *dms; 2374 2375 ierr = PetscMalloc1(coarsen, &dms);CHKERRQ(ierr); 2376 ierr = DMCoarsenHierarchy(dm, coarsen, dms);CHKERRQ(ierr); 2377 /* Free DMs */ 2378 for (r = 0; r < coarsen; ++r) { 2379 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr); 2380 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 2381 } 2382 ierr = PetscFree(dms);CHKERRQ(ierr); 2383 } else { 2384 for (r = 0; r < coarsen; ++r) { 2385 DM coarseMesh; 2386 2387 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2388 ierr = DMCoarsen(dm, PetscObjectComm((PetscObject) dm), &coarseMesh);CHKERRQ(ierr); 2389 /* Total hack since we do not pass in a pointer */ 2390 ierr = DMPlexReplace_Static(dm, coarseMesh);CHKERRQ(ierr); 2391 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2392 ierr = DMDestroy(&coarseMesh);CHKERRQ(ierr); 2393 } 2394 } 2395 /* Handle */ 2396 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2397 ierr = PetscOptionsTail();CHKERRQ(ierr); 2398 PetscFunctionReturn(0); 2399 } 2400 2401 static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec) 2402 { 2403 PetscErrorCode ierr; 2404 2405 PetscFunctionBegin; 2406 ierr = DMCreateGlobalVector_Section_Private(dm,vec);CHKERRQ(ierr); 2407 /* ierr = VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM);CHKERRQ(ierr); */ 2408 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex);CHKERRQ(ierr); 2409 ierr = VecSetOperation(*vec, VECOP_VIEWNATIVE, (void (*)(void)) VecView_Plex_Native);CHKERRQ(ierr); 2410 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex);CHKERRQ(ierr); 2411 ierr = VecSetOperation(*vec, VECOP_LOADNATIVE, (void (*)(void)) VecLoad_Plex_Native);CHKERRQ(ierr); 2412 PetscFunctionReturn(0); 2413 } 2414 2415 static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec) 2416 { 2417 PetscErrorCode ierr; 2418 2419 PetscFunctionBegin; 2420 ierr = DMCreateLocalVector_Section_Private(dm,vec);CHKERRQ(ierr); 2421 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex_Local);CHKERRQ(ierr); 2422 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex_Local);CHKERRQ(ierr); 2423 PetscFunctionReturn(0); 2424 } 2425 2426 static PetscErrorCode DMGetDimPoints_Plex(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd) 2427 { 2428 PetscInt depth, d; 2429 PetscErrorCode ierr; 2430 2431 PetscFunctionBegin; 2432 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 2433 if (depth == 1) { 2434 ierr = DMGetDimension(dm, &d);CHKERRQ(ierr); 2435 if (dim == 0) {ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr);} 2436 else if (dim == d) {ierr = DMPlexGetDepthStratum(dm, 1, pStart, pEnd);CHKERRQ(ierr);} 2437 else {*pStart = 0; *pEnd = 0;} 2438 } else { 2439 ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr); 2440 } 2441 PetscFunctionReturn(0); 2442 } 2443 2444 static PetscErrorCode DMGetNeighors_Plex(DM dm, PetscInt *nranks, const PetscMPIInt *ranks[]) 2445 { 2446 PetscSF sf; 2447 PetscErrorCode ierr; 2448 2449 PetscFunctionBegin; 2450 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 2451 ierr = PetscSFGetRanks(sf, nranks, ranks, NULL, NULL, NULL);CHKERRQ(ierr); 2452 PetscFunctionReturn(0); 2453 } 2454 2455 static PetscErrorCode DMHasCreateInjection_Plex(DM dm, PetscBool *flg) 2456 { 2457 PetscFunctionBegin; 2458 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2459 PetscValidPointer(flg,2); 2460 *flg = PETSC_TRUE; 2461 PetscFunctionReturn(0); 2462 } 2463 2464 static PetscErrorCode DMInitialize_Plex(DM dm) 2465 { 2466 PetscErrorCode ierr; 2467 2468 PetscFunctionBegin; 2469 dm->ops->view = DMView_Plex; 2470 dm->ops->load = DMLoad_Plex; 2471 dm->ops->setfromoptions = DMSetFromOptions_Plex; 2472 dm->ops->clone = DMClone_Plex; 2473 dm->ops->setup = DMSetUp_Plex; 2474 dm->ops->createdefaultsection = DMCreateDefaultSection_Plex; 2475 dm->ops->createdefaultconstraints = DMCreateDefaultConstraints_Plex; 2476 dm->ops->createglobalvector = DMCreateGlobalVector_Plex; 2477 dm->ops->createlocalvector = DMCreateLocalVector_Plex; 2478 dm->ops->getlocaltoglobalmapping = NULL; 2479 dm->ops->createfieldis = NULL; 2480 dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex; 2481 dm->ops->createcoordinatefield = DMCreateCoordinateField_Plex; 2482 dm->ops->getcoloring = NULL; 2483 dm->ops->creatematrix = DMCreateMatrix_Plex; 2484 dm->ops->createinterpolation = DMCreateInterpolation_Plex; 2485 dm->ops->createmassmatrix = DMCreateMassMatrix_Plex; 2486 dm->ops->getaggregates = NULL; 2487 dm->ops->getinjection = DMCreateInjection_Plex; 2488 dm->ops->hascreateinjection = DMHasCreateInjection_Plex; 2489 dm->ops->refine = DMRefine_Plex; 2490 dm->ops->coarsen = DMCoarsen_Plex; 2491 dm->ops->refinehierarchy = DMRefineHierarchy_Plex; 2492 dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Plex; 2493 dm->ops->adaptlabel = DMAdaptLabel_Plex; 2494 dm->ops->adaptmetric = DMAdaptMetric_Plex; 2495 dm->ops->globaltolocalbegin = NULL; 2496 dm->ops->globaltolocalend = NULL; 2497 dm->ops->localtoglobalbegin = NULL; 2498 dm->ops->localtoglobalend = NULL; 2499 dm->ops->destroy = DMDestroy_Plex; 2500 dm->ops->createsubdm = DMCreateSubDM_Plex; 2501 dm->ops->createsuperdm = DMCreateSuperDM_Plex; 2502 dm->ops->getdimpoints = DMGetDimPoints_Plex; 2503 dm->ops->locatepoints = DMLocatePoints_Plex; 2504 dm->ops->projectfunctionlocal = DMProjectFunctionLocal_Plex; 2505 dm->ops->projectfunctionlabellocal = DMProjectFunctionLabelLocal_Plex; 2506 dm->ops->projectfieldlocal = DMProjectFieldLocal_Plex; 2507 dm->ops->projectfieldlabellocal = DMProjectFieldLabelLocal_Plex; 2508 dm->ops->computel2diff = DMComputeL2Diff_Plex; 2509 dm->ops->computel2gradientdiff = DMComputeL2GradientDiff_Plex; 2510 dm->ops->computel2fielddiff = DMComputeL2FieldDiff_Plex; 2511 dm->ops->getneighbors = DMGetNeighors_Plex; 2512 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMPlexInsertBoundaryValues_C",DMPlexInsertBoundaryValues_Plex);CHKERRQ(ierr); 2513 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMSetUpGLVisViewer_C",DMSetUpGLVisViewer_Plex);CHKERRQ(ierr); 2514 PetscFunctionReturn(0); 2515 } 2516 2517 PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm) 2518 { 2519 DM_Plex *mesh = (DM_Plex *) dm->data; 2520 PetscErrorCode ierr; 2521 2522 PetscFunctionBegin; 2523 mesh->refct++; 2524 (*newdm)->data = mesh; 2525 ierr = PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);CHKERRQ(ierr); 2526 ierr = DMInitialize_Plex(*newdm);CHKERRQ(ierr); 2527 PetscFunctionReturn(0); 2528 } 2529 2530 /*MC 2531 DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram. 2532 In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is 2533 specified by a PetscSection object. Ownership in the global representation is determined by 2534 ownership of the underlying DMPlex points. This is specified by another PetscSection object. 2535 2536 Options Database Keys: 2537 + -dm_plex_hash_location - Use grid hashing for point location 2538 . -dm_plex_partition_balance - Attempt to evenly divide points on partition boundary between processes 2539 . -dm_plex_remesh_bd - Allow changes to the boundary on remeshing 2540 . -dm_plex_max_projection_height - Maxmimum mesh point height used to project locally 2541 . -dm_plex_regular_refinement - Use special nested projection algorithm for regular refinement 2542 . -dm_plex_check_symmetry - Check that the adjacency information in the mesh is symmetric 2543 . -dm_plex_check_skeleton <celltype> - Check that each cell has the correct number of vertices 2544 . -dm_plex_check_faces <celltype> - Check that the faces of each cell give a vertex order this is consistent with what we expect from the cell type 2545 . -dm_plex_check_geometry - Check that cells have positive volume 2546 . -dm_view :mesh.tex:ascii_latex - View the mesh in LaTeX/TikZ 2547 . -dm_plex_view_scale <num> - Scale the TikZ 2548 - -dm_plex_print_fem <num> - View FEM assembly information, such as element vectors and matrices 2549 2550 2551 Level: intermediate 2552 2553 .seealso: DMType, DMPlexCreate(), DMCreate(), DMSetType() 2554 M*/ 2555 2556 PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm) 2557 { 2558 DM_Plex *mesh; 2559 PetscInt unit, d; 2560 PetscErrorCode ierr; 2561 2562 PetscFunctionBegin; 2563 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2564 ierr = PetscNewLog(dm,&mesh);CHKERRQ(ierr); 2565 dm->dim = 0; 2566 dm->data = mesh; 2567 2568 mesh->refct = 1; 2569 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);CHKERRQ(ierr); 2570 mesh->maxConeSize = 0; 2571 mesh->cones = NULL; 2572 mesh->coneOrientations = NULL; 2573 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);CHKERRQ(ierr); 2574 mesh->maxSupportSize = 0; 2575 mesh->supports = NULL; 2576 mesh->refinementUniform = PETSC_TRUE; 2577 mesh->refinementLimit = -1.0; 2578 2579 mesh->facesTmp = NULL; 2580 2581 mesh->tetgenOpts = NULL; 2582 mesh->triangleOpts = NULL; 2583 ierr = PetscPartitionerCreate(PetscObjectComm((PetscObject)dm), &mesh->partitioner);CHKERRQ(ierr); 2584 mesh->remeshBd = PETSC_FALSE; 2585 2586 mesh->subpointMap = NULL; 2587 2588 for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0; 2589 2590 mesh->regularRefinement = PETSC_FALSE; 2591 mesh->depthState = -1; 2592 mesh->globalVertexNumbers = NULL; 2593 mesh->globalCellNumbers = NULL; 2594 mesh->anchorSection = NULL; 2595 mesh->anchorIS = NULL; 2596 mesh->createanchors = NULL; 2597 mesh->computeanchormatrix = NULL; 2598 mesh->parentSection = NULL; 2599 mesh->parents = NULL; 2600 mesh->childIDs = NULL; 2601 mesh->childSection = NULL; 2602 mesh->children = NULL; 2603 mesh->referenceTree = NULL; 2604 mesh->getchildsymmetry = NULL; 2605 for (d = 0; d < 8; ++d) mesh->hybridPointMax[d] = PETSC_DETERMINE; 2606 mesh->vtkCellHeight = 0; 2607 mesh->useAnchors = PETSC_FALSE; 2608 2609 mesh->maxProjectionHeight = 0; 2610 2611 mesh->printSetValues = PETSC_FALSE; 2612 mesh->printFEM = 0; 2613 mesh->printTol = 1.0e-10; 2614 2615 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 2616 PetscFunctionReturn(0); 2617 } 2618 2619 /*@ 2620 DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram. 2621 2622 Collective on MPI_Comm 2623 2624 Input Parameter: 2625 . comm - The communicator for the DMPlex object 2626 2627 Output Parameter: 2628 . mesh - The DMPlex object 2629 2630 Level: beginner 2631 2632 .keywords: DMPlex, create 2633 @*/ 2634 PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh) 2635 { 2636 PetscErrorCode ierr; 2637 2638 PetscFunctionBegin; 2639 PetscValidPointer(mesh,2); 2640 ierr = DMCreate(comm, mesh);CHKERRQ(ierr); 2641 ierr = DMSetType(*mesh, DMPLEX);CHKERRQ(ierr); 2642 PetscFunctionReturn(0); 2643 } 2644 2645 /* 2646 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 2647 */ 2648 /* TODO: invertCells and spaceDim arguments could be added also to to DMPlexCreateFromCellListParallel(), DMPlexBuildFromCellList_Internal() and DMPlexCreateFromCellList() */ 2649 PetscErrorCode DMPlexBuildFromCellList_Parallel_Internal(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[], PetscBool invertCells, PetscSF *sfVert) 2650 { 2651 PetscSF sfPoint; 2652 PetscLayout vLayout; 2653 PetscHSetI vhash; 2654 PetscSFNode *remoteVerticesAdj, *vertexLocal, *vertexOwner, *remoteVertex; 2655 const PetscInt *vrange; 2656 PetscInt numVerticesAdj, off = 0, *verticesAdj, numVerticesGhost = 0, *localVertex, *cone, c, p, v, g; 2657 PetscMPIInt rank, size; 2658 PetscErrorCode ierr; 2659 2660 PetscFunctionBegin; 2661 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr); 2662 ierr = MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);CHKERRQ(ierr); 2663 /* Partition vertices */ 2664 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) dm), &vLayout);CHKERRQ(ierr); 2665 ierr = PetscLayoutSetLocalSize(vLayout, numVertices);CHKERRQ(ierr); 2666 ierr = PetscLayoutSetBlockSize(vLayout, 1);CHKERRQ(ierr); 2667 ierr = PetscLayoutSetUp(vLayout);CHKERRQ(ierr); 2668 ierr = PetscLayoutGetRanges(vLayout, &vrange);CHKERRQ(ierr); 2669 /* Count vertices and map them to procs */ 2670 ierr = PetscHSetICreate(&vhash);CHKERRQ(ierr); 2671 for (c = 0; c < numCells; ++c) { 2672 for (p = 0; p < numCorners; ++p) { 2673 ierr = PetscHSetIAdd(vhash, cells[c*numCorners+p]);CHKERRQ(ierr); 2674 } 2675 } 2676 ierr = PetscHSetIGetSize(vhash, &numVerticesAdj);CHKERRQ(ierr); 2677 ierr = PetscMalloc1(numVerticesAdj, &verticesAdj);CHKERRQ(ierr); 2678 ierr = PetscHSetIGetElems(vhash, &off, verticesAdj);CHKERRQ(ierr); 2679 ierr = PetscHSetIDestroy(&vhash);CHKERRQ(ierr); 2680 if (off != numVerticesAdj) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid number of local vertices %D should be %D", off, numVerticesAdj); 2681 ierr = PetscSortInt(numVerticesAdj, verticesAdj);CHKERRQ(ierr); 2682 ierr = PetscMalloc1(numVerticesAdj, &remoteVerticesAdj);CHKERRQ(ierr); 2683 for (v = 0; v < numVerticesAdj; ++v) { 2684 const PetscInt gv = verticesAdj[v]; 2685 PetscInt vrank; 2686 2687 ierr = PetscFindInt(gv, size+1, vrange, &vrank);CHKERRQ(ierr); 2688 vrank = vrank < 0 ? -(vrank+2) : vrank; 2689 remoteVerticesAdj[v].index = gv - vrange[vrank]; 2690 remoteVerticesAdj[v].rank = vrank; 2691 } 2692 /* Create cones */ 2693 ierr = DMPlexSetChart(dm, 0, numCells+numVerticesAdj);CHKERRQ(ierr); 2694 for (c = 0; c < numCells; ++c) {ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr);} 2695 ierr = DMSetUp(dm);CHKERRQ(ierr); 2696 ierr = DMGetWorkArray(dm, numCorners, MPIU_INT, &cone);CHKERRQ(ierr); 2697 for (c = 0; c < numCells; ++c) { 2698 for (p = 0; p < numCorners; ++p) { 2699 const PetscInt gv = cells[c*numCorners+p]; 2700 PetscInt lv; 2701 2702 ierr = PetscFindInt(gv, numVerticesAdj, verticesAdj, &lv);CHKERRQ(ierr); 2703 if (lv < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not find global vertex %D in local connectivity", gv); 2704 cone[p] = lv+numCells; 2705 } 2706 if (invertCells) { ierr = DMPlexInvertCell_Internal(spaceDim, numCorners, cone);CHKERRQ(ierr); } 2707 ierr = DMPlexSetCone(dm, c, cone);CHKERRQ(ierr); 2708 } 2709 ierr = DMRestoreWorkArray(dm, numCorners, MPIU_INT, &cone);CHKERRQ(ierr); 2710 /* Create SF for vertices */ 2711 ierr = PetscSFCreate(PetscObjectComm((PetscObject)dm), sfVert);CHKERRQ(ierr); 2712 ierr = PetscObjectSetName((PetscObject) *sfVert, "Vertex Ownership SF");CHKERRQ(ierr); 2713 ierr = PetscSFSetFromOptions(*sfVert);CHKERRQ(ierr); 2714 ierr = PetscSFSetGraph(*sfVert, numVertices, numVerticesAdj, NULL, PETSC_OWN_POINTER, remoteVerticesAdj, PETSC_OWN_POINTER);CHKERRQ(ierr); 2715 ierr = PetscFree(verticesAdj);CHKERRQ(ierr); 2716 /* Build pointSF */ 2717 ierr = PetscMalloc2(numVerticesAdj, &vertexLocal, numVertices, &vertexOwner);CHKERRQ(ierr); 2718 for (v = 0; v < numVerticesAdj; ++v) {vertexLocal[v].index = v+numCells; vertexLocal[v].rank = rank;} 2719 for (v = 0; v < numVertices; ++v) {vertexOwner[v].index = -1; vertexOwner[v].rank = -1;} 2720 ierr = PetscSFReduceBegin(*sfVert, MPIU_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);CHKERRQ(ierr); 2721 ierr = PetscSFReduceEnd(*sfVert, MPIU_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);CHKERRQ(ierr); 2722 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); 2723 ierr = PetscSFBcastBegin(*sfVert, MPIU_2INT, vertexOwner, vertexLocal);CHKERRQ(ierr); 2724 ierr = PetscSFBcastEnd(*sfVert, MPIU_2INT, vertexOwner, vertexLocal);CHKERRQ(ierr); 2725 for (v = 0; v < numVerticesAdj; ++v) if (vertexLocal[v].rank != rank) ++numVerticesGhost; 2726 ierr = PetscMalloc1(numVerticesGhost, &localVertex);CHKERRQ(ierr); 2727 ierr = PetscMalloc1(numVerticesGhost, &remoteVertex);CHKERRQ(ierr); 2728 for (v = 0, g = 0; v < numVerticesAdj; ++v) { 2729 if (vertexLocal[v].rank != rank) { 2730 localVertex[g] = v+numCells; 2731 remoteVertex[g].index = vertexLocal[v].index; 2732 remoteVertex[g].rank = vertexLocal[v].rank; 2733 ++g; 2734 } 2735 } 2736 ierr = PetscFree2(vertexLocal, vertexOwner);CHKERRQ(ierr); 2737 if (g != numVerticesGhost) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of vertex ghosts %D should be %D", g, numVerticesGhost); 2738 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 2739 ierr = PetscObjectSetName((PetscObject) sfPoint, "point SF");CHKERRQ(ierr); 2740 ierr = PetscSFSetGraph(sfPoint, numCells+numVerticesAdj, numVerticesGhost, localVertex, PETSC_OWN_POINTER, remoteVertex, PETSC_OWN_POINTER);CHKERRQ(ierr); 2741 ierr = PetscLayoutDestroy(&vLayout);CHKERRQ(ierr); 2742 /* Fill in the rest of the topology structure */ 2743 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 2744 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 2745 PetscFunctionReturn(0); 2746 } 2747 2748 /* 2749 This takes as input the coordinates for each owned vertex 2750 */ 2751 PetscErrorCode DMPlexBuildCoordinates_Parallel_Internal(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numV, PetscSF sfVert, const PetscReal vertexCoords[]) 2752 { 2753 PetscSection coordSection; 2754 Vec coordinates; 2755 PetscScalar *coords; 2756 PetscInt numVertices, numVerticesAdj, coordSize, v; 2757 PetscErrorCode ierr; 2758 2759 PetscFunctionBegin; 2760 ierr = DMSetCoordinateDim(dm, spaceDim);CHKERRQ(ierr); 2761 ierr = PetscSFGetGraph(sfVert, &numVertices, &numVerticesAdj, NULL, NULL);CHKERRQ(ierr); 2762 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2763 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 2764 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 2765 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVerticesAdj);CHKERRQ(ierr); 2766 for (v = numCells; v < numCells+numVerticesAdj; ++v) { 2767 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 2768 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 2769 } 2770 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 2771 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 2772 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 2773 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 2774 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 2775 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 2776 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 2777 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2778 { 2779 MPI_Datatype coordtype; 2780 2781 /* Need a temp buffer for coords if we have complex/single */ 2782 ierr = MPI_Type_contiguous(spaceDim, MPIU_SCALAR, &coordtype);CHKERRQ(ierr); 2783 ierr = MPI_Type_commit(&coordtype);CHKERRQ(ierr); 2784 #if defined(PETSC_USE_COMPLEX) 2785 { 2786 PetscScalar *svertexCoords; 2787 PetscInt i; 2788 ierr = PetscMalloc1(numV*spaceDim,&svertexCoords);CHKERRQ(ierr); 2789 for (i=0; i<numV*spaceDim; i++) svertexCoords[i] = vertexCoords[i]; 2790 ierr = PetscSFBcastBegin(sfVert, coordtype, svertexCoords, coords);CHKERRQ(ierr); 2791 ierr = PetscSFBcastEnd(sfVert, coordtype, svertexCoords, coords);CHKERRQ(ierr); 2792 ierr = PetscFree(svertexCoords);CHKERRQ(ierr); 2793 } 2794 #else 2795 ierr = PetscSFBcastBegin(sfVert, coordtype, vertexCoords, coords);CHKERRQ(ierr); 2796 ierr = PetscSFBcastEnd(sfVert, coordtype, vertexCoords, coords);CHKERRQ(ierr); 2797 #endif 2798 ierr = MPI_Type_free(&coordtype);CHKERRQ(ierr); 2799 } 2800 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2801 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 2802 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 2803 PetscFunctionReturn(0); 2804 } 2805 2806 /*@C 2807 DMPlexCreateFromCellListParallel - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM 2808 2809 Input Parameters: 2810 + comm - The communicator 2811 . dim - The topological dimension of the mesh 2812 . numCells - The number of cells owned by this process 2813 . numVertices - The number of vertices owned by this process 2814 . numCorners - The number of vertices for each cell 2815 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 2816 . cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 2817 . spaceDim - The spatial dimension used for coordinates 2818 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 2819 2820 Output Parameter: 2821 + dm - The DM 2822 - vertexSF - Optional, SF describing complete vertex ownership 2823 2824 Note: Two triangles sharing a face 2825 $ 2826 $ 2 2827 $ / | \ 2828 $ / | \ 2829 $ / | \ 2830 $ 0 0 | 1 3 2831 $ \ | / 2832 $ \ | / 2833 $ \ | / 2834 $ 1 2835 would have input 2836 $ numCells = 2, numVertices = 4 2837 $ cells = [0 1 2 1 3 2] 2838 $ 2839 which would result in the DMPlex 2840 $ 2841 $ 4 2842 $ / | \ 2843 $ / | \ 2844 $ / | \ 2845 $ 2 0 | 1 5 2846 $ \ | / 2847 $ \ | / 2848 $ \ | / 2849 $ 3 2850 2851 Level: beginner 2852 2853 .seealso: DMPlexCreateFromCellList(), DMPlexCreateFromDAG(), DMPlexCreate() 2854 @*/ 2855 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) 2856 { 2857 PetscSF sfVert; 2858 PetscErrorCode ierr; 2859 2860 PetscFunctionBegin; 2861 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 2862 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 2863 PetscValidLogicalCollectiveInt(*dm, dim, 2); 2864 PetscValidLogicalCollectiveInt(*dm, spaceDim, 8); 2865 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 2866 ierr = DMPlexBuildFromCellList_Parallel_Internal(*dm, spaceDim, numCells, numVertices, numCorners, cells, PETSC_FALSE, &sfVert);CHKERRQ(ierr); 2867 if (interpolate) { 2868 DM idm; 2869 2870 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 2871 ierr = DMDestroy(dm);CHKERRQ(ierr); 2872 *dm = idm; 2873 } 2874 ierr = DMPlexBuildCoordinates_Parallel_Internal(*dm, spaceDim, numCells, numVertices, sfVert, vertexCoords);CHKERRQ(ierr); 2875 if (vertexSF) *vertexSF = sfVert; 2876 else {ierr = PetscSFDestroy(&sfVert);CHKERRQ(ierr);} 2877 PetscFunctionReturn(0); 2878 } 2879 2880 /* 2881 This takes as input the common mesh generator output, a list of the vertices for each cell 2882 */ 2883 PetscErrorCode DMPlexBuildFromCellList_Internal(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[], PetscBool invertCells) 2884 { 2885 PetscInt *cone, c, p; 2886 PetscErrorCode ierr; 2887 2888 PetscFunctionBegin; 2889 ierr = DMPlexSetChart(dm, 0, numCells+numVertices);CHKERRQ(ierr); 2890 for (c = 0; c < numCells; ++c) { 2891 ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr); 2892 } 2893 ierr = DMSetUp(dm);CHKERRQ(ierr); 2894 ierr = DMGetWorkArray(dm, numCorners, MPIU_INT, &cone);CHKERRQ(ierr); 2895 for (c = 0; c < numCells; ++c) { 2896 for (p = 0; p < numCorners; ++p) { 2897 cone[p] = cells[c*numCorners+p]+numCells; 2898 } 2899 if (invertCells) { ierr = DMPlexInvertCell_Internal(spaceDim, numCorners, cone);CHKERRQ(ierr); } 2900 ierr = DMPlexSetCone(dm, c, cone);CHKERRQ(ierr); 2901 } 2902 ierr = DMRestoreWorkArray(dm, numCorners, MPIU_INT, &cone);CHKERRQ(ierr); 2903 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 2904 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 2905 PetscFunctionReturn(0); 2906 } 2907 2908 /* 2909 This takes as input the coordinates for each vertex 2910 */ 2911 PetscErrorCode DMPlexBuildCoordinates_Internal(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, const double vertexCoords[]) 2912 { 2913 PetscSection coordSection; 2914 Vec coordinates; 2915 DM cdm; 2916 PetscScalar *coords; 2917 PetscInt v, d; 2918 PetscErrorCode ierr; 2919 2920 PetscFunctionBegin; 2921 ierr = DMSetCoordinateDim(dm, spaceDim);CHKERRQ(ierr); 2922 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2923 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 2924 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 2925 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVertices);CHKERRQ(ierr); 2926 for (v = numCells; v < numCells+numVertices; ++v) { 2927 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 2928 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 2929 } 2930 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 2931 2932 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 2933 ierr = DMCreateLocalVector(cdm, &coordinates);CHKERRQ(ierr); 2934 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 2935 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 2936 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2937 for (v = 0; v < numVertices; ++v) { 2938 for (d = 0; d < spaceDim; ++d) { 2939 coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d]; 2940 } 2941 } 2942 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2943 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 2944 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 2945 PetscFunctionReturn(0); 2946 } 2947 2948 /*@C 2949 DMPlexCreateFromCellList - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM 2950 2951 Input Parameters: 2952 + comm - The communicator 2953 . dim - The topological dimension of the mesh 2954 . numCells - The number of cells 2955 . numVertices - The number of vertices 2956 . numCorners - The number of vertices for each cell 2957 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 2958 . cells - An array of numCells*numCorners numbers, the vertices for each cell 2959 . spaceDim - The spatial dimension used for coordinates 2960 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 2961 2962 Output Parameter: 2963 . dm - The DM 2964 2965 Note: Two triangles sharing a face 2966 $ 2967 $ 2 2968 $ / | \ 2969 $ / | \ 2970 $ / | \ 2971 $ 0 0 | 1 3 2972 $ \ | / 2973 $ \ | / 2974 $ \ | / 2975 $ 1 2976 would have input 2977 $ numCells = 2, numVertices = 4 2978 $ cells = [0 1 2 1 3 2] 2979 $ 2980 which would result in the DMPlex 2981 $ 2982 $ 4 2983 $ / | \ 2984 $ / | \ 2985 $ / | \ 2986 $ 2 0 | 1 5 2987 $ \ | / 2988 $ \ | / 2989 $ \ | / 2990 $ 3 2991 2992 Level: beginner 2993 2994 .seealso: DMPlexCreateFromDAG(), DMPlexCreate() 2995 @*/ 2996 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) 2997 { 2998 PetscErrorCode ierr; 2999 3000 PetscFunctionBegin; 3001 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 3002 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 3003 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 3004 ierr = DMPlexBuildFromCellList_Internal(*dm, spaceDim, numCells, numVertices, numCorners, cells, PETSC_FALSE);CHKERRQ(ierr); 3005 if (interpolate) { 3006 DM idm; 3007 3008 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 3009 ierr = DMDestroy(dm);CHKERRQ(ierr); 3010 *dm = idm; 3011 } 3012 ierr = DMPlexBuildCoordinates_Internal(*dm, spaceDim, numCells, numVertices, vertexCoords);CHKERRQ(ierr); 3013 PetscFunctionReturn(0); 3014 } 3015 3016 /*@ 3017 DMPlexCreateFromDAG - This takes as input the adjacency-list representation of the Directed Acyclic Graph (Hasse Diagram) encoding a mesh, and produces a DM 3018 3019 Input Parameters: 3020 + dm - The empty DM object, usually from DMCreate() and DMSetDimension() 3021 . depth - The depth of the DAG 3022 . numPoints - Array of size depth + 1 containing the number of points at each depth 3023 . coneSize - The cone size of each point 3024 . cones - The concatenation of the cone points for each point, the cone list must be oriented correctly for each point 3025 . coneOrientations - The orientation of each cone point 3026 - vertexCoords - An array of numPoints[0]*spacedim numbers representing the coordinates of each vertex, with spacedim the value set via DMSetCoordinateDim() 3027 3028 Output Parameter: 3029 . dm - The DM 3030 3031 Note: Two triangles sharing a face would have input 3032 $ depth = 1, numPoints = [4 2], coneSize = [3 3 0 0 0 0] 3033 $ cones = [2 3 4 3 5 4], coneOrientations = [0 0 0 0 0 0] 3034 $ vertexCoords = [-1.0 0.0 0.0 -1.0 0.0 1.0 1.0 0.0] 3035 $ 3036 which would result in the DMPlex 3037 $ 3038 $ 4 3039 $ / | \ 3040 $ / | \ 3041 $ / | \ 3042 $ 2 0 | 1 5 3043 $ \ | / 3044 $ \ | / 3045 $ \ | / 3046 $ 3 3047 $ 3048 $ Notice that all points are numbered consecutively, unlikely DMPlexCreateFromCellList() 3049 3050 Level: advanced 3051 3052 .seealso: DMPlexCreateFromCellList(), DMPlexCreate() 3053 @*/ 3054 PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[]) 3055 { 3056 Vec coordinates; 3057 PetscSection coordSection; 3058 PetscScalar *coords; 3059 PetscInt coordSize, firstVertex = -1, pStart = 0, pEnd = 0, p, v, dim, dimEmbed, d, off; 3060 PetscErrorCode ierr; 3061 3062 PetscFunctionBegin; 3063 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3064 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 3065 if (dimEmbed < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Embedding dimension %d cannot be less than intrinsic dimension %d",dimEmbed,dim); 3066 for (d = 0; d <= depth; ++d) pEnd += numPoints[d]; 3067 ierr = DMPlexSetChart(dm, pStart, pEnd);CHKERRQ(ierr); 3068 for (p = pStart; p < pEnd; ++p) { 3069 ierr = DMPlexSetConeSize(dm, p, coneSize[p-pStart]);CHKERRQ(ierr); 3070 if (firstVertex < 0 && !coneSize[p - pStart]) { 3071 firstVertex = p - pStart; 3072 } 3073 } 3074 if (firstVertex < 0 && numPoints[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected %d vertices but could not find any", numPoints[0]); 3075 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 3076 for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) { 3077 ierr = DMPlexSetCone(dm, p, &cones[off]);CHKERRQ(ierr); 3078 ierr = DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);CHKERRQ(ierr); 3079 } 3080 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 3081 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 3082 /* Build coordinates */ 3083 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 3084 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 3085 ierr = PetscSectionSetFieldComponents(coordSection, 0, dimEmbed);CHKERRQ(ierr); 3086 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0]);CHKERRQ(ierr); 3087 for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) { 3088 ierr = PetscSectionSetDof(coordSection, v, dimEmbed);CHKERRQ(ierr); 3089 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed);CHKERRQ(ierr); 3090 } 3091 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 3092 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 3093 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 3094 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 3095 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 3096 ierr = VecSetBlockSize(coordinates, dimEmbed);CHKERRQ(ierr); 3097 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 3098 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 3099 for (v = 0; v < numPoints[0]; ++v) { 3100 PetscInt off; 3101 3102 ierr = PetscSectionGetOffset(coordSection, v+firstVertex, &off);CHKERRQ(ierr); 3103 for (d = 0; d < dimEmbed; ++d) { 3104 coords[off+d] = vertexCoords[v*dimEmbed+d]; 3105 } 3106 } 3107 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 3108 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 3109 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 3110 PetscFunctionReturn(0); 3111 } 3112 3113 /*@C 3114 DMPlexCreateCellVertexFromFile - Create a DMPlex mesh from a simple cell-vertex file. 3115 3116 + comm - The MPI communicator 3117 . filename - Name of the .dat file 3118 - interpolate - Create faces and edges in the mesh 3119 3120 Output Parameter: 3121 . dm - The DM object representing the mesh 3122 3123 Note: The format is the simplest possible: 3124 $ Ne 3125 $ v0 v1 ... vk 3126 $ Nv 3127 $ x y z marker 3128 3129 Level: beginner 3130 3131 .seealso: DMPlexCreateFromFile(), DMPlexCreateMedFromFile(), DMPlexCreateGmsh(), DMPlexCreate() 3132 @*/ 3133 PetscErrorCode DMPlexCreateCellVertexFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 3134 { 3135 DMLabel marker; 3136 PetscViewer viewer; 3137 Vec coordinates; 3138 PetscSection coordSection; 3139 PetscScalar *coords; 3140 char line[PETSC_MAX_PATH_LEN]; 3141 PetscInt dim = 3, cdim = 3, coordSize, v, c, d; 3142 PetscMPIInt rank; 3143 int snum, Nv, Nc; 3144 PetscErrorCode ierr; 3145 3146 PetscFunctionBegin; 3147 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 3148 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 3149 ierr = PetscViewerSetType(viewer, PETSCVIEWERASCII);CHKERRQ(ierr); 3150 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 3151 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 3152 if (!rank) { 3153 ierr = PetscViewerRead(viewer, line, 2, NULL, PETSC_STRING);CHKERRQ(ierr); 3154 snum = sscanf(line, "%d %d", &Nc, &Nv); 3155 if (snum != 2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 3156 } else { 3157 Nc = Nv = 0; 3158 } 3159 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 3160 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 3161 ierr = DMPlexSetChart(*dm, 0, Nc+Nv);CHKERRQ(ierr); 3162 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 3163 ierr = DMSetCoordinateDim(*dm, cdim);CHKERRQ(ierr); 3164 /* Read topology */ 3165 if (!rank) { 3166 PetscInt cone[8], corners = 8; 3167 int vbuf[8], v; 3168 3169 for (c = 0; c < Nc; ++c) {ierr = DMPlexSetConeSize(*dm, c, corners);CHKERRQ(ierr);} 3170 ierr = DMSetUp(*dm);CHKERRQ(ierr); 3171 for (c = 0; c < Nc; ++c) { 3172 ierr = PetscViewerRead(viewer, line, corners, NULL, PETSC_STRING);CHKERRQ(ierr); 3173 snum = sscanf(line, "%d %d %d %d %d %d %d %d", &vbuf[0], &vbuf[1], &vbuf[2], &vbuf[3], &vbuf[4], &vbuf[5], &vbuf[6], &vbuf[7]); 3174 if (snum != corners) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 3175 for (v = 0; v < corners; ++v) cone[v] = vbuf[v] + Nc; 3176 /* Hexahedra are inverted */ 3177 { 3178 PetscInt tmp = cone[1]; 3179 cone[1] = cone[3]; 3180 cone[3] = tmp; 3181 } 3182 ierr = DMPlexSetCone(*dm, c, cone);CHKERRQ(ierr); 3183 } 3184 } 3185 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 3186 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 3187 /* Read coordinates */ 3188 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 3189 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 3190 ierr = PetscSectionSetFieldComponents(coordSection, 0, cdim);CHKERRQ(ierr); 3191 ierr = PetscSectionSetChart(coordSection, Nc, Nc + Nv);CHKERRQ(ierr); 3192 for (v = Nc; v < Nc+Nv; ++v) { 3193 ierr = PetscSectionSetDof(coordSection, v, cdim);CHKERRQ(ierr); 3194 ierr = PetscSectionSetFieldDof(coordSection, v, 0, cdim);CHKERRQ(ierr); 3195 } 3196 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 3197 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 3198 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 3199 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 3200 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 3201 ierr = VecSetBlockSize(coordinates, cdim);CHKERRQ(ierr); 3202 ierr = VecSetType(coordinates, VECSTANDARD);CHKERRQ(ierr); 3203 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 3204 if (!rank) { 3205 double x[3]; 3206 int val; 3207 3208 ierr = DMCreateLabel(*dm, "marker");CHKERRQ(ierr); 3209 ierr = DMGetLabel(*dm, "marker", &marker);CHKERRQ(ierr); 3210 for (v = 0; v < Nv; ++v) { 3211 ierr = PetscViewerRead(viewer, line, 4, NULL, PETSC_STRING);CHKERRQ(ierr); 3212 snum = sscanf(line, "%lg %lg %lg %d", &x[0], &x[1], &x[2], &val); 3213 if (snum != 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 3214 for (d = 0; d < cdim; ++d) coords[v*cdim+d] = x[d]; 3215 if (val) {ierr = DMLabelSetValue(marker, v+Nc, val);CHKERRQ(ierr);} 3216 } 3217 } 3218 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 3219 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 3220 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 3221 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3222 if (interpolate) { 3223 DM idm; 3224 DMLabel bdlabel; 3225 3226 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 3227 ierr = DMDestroy(dm);CHKERRQ(ierr); 3228 *dm = idm; 3229 3230 ierr = DMGetLabel(*dm, "marker", &bdlabel);CHKERRQ(ierr); 3231 ierr = DMPlexMarkBoundaryFaces(*dm, PETSC_DETERMINE, bdlabel);CHKERRQ(ierr); 3232 ierr = DMPlexLabelComplete(*dm, bdlabel);CHKERRQ(ierr); 3233 } 3234 PetscFunctionReturn(0); 3235 } 3236 3237 /*@C 3238 DMPlexCreateFromFile - This takes a filename and produces a DM 3239 3240 Input Parameters: 3241 + comm - The communicator 3242 . filename - A file name 3243 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 3244 3245 Output Parameter: 3246 . dm - The DM 3247 3248 Options Database Keys: 3249 . -dm_plex_create_from_hdf5_xdmf - use the PETSC_VIEWER_HDF5_XDMF format for reading HDF5 3250 3251 Level: beginner 3252 3253 .seealso: DMPlexCreateFromDAG(), DMPlexCreateFromCellList(), DMPlexCreate() 3254 @*/ 3255 PetscErrorCode DMPlexCreateFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 3256 { 3257 const char *extGmsh = ".msh"; 3258 const char *extGmsh2 = ".msh2"; 3259 const char *extGmsh4 = ".msh4"; 3260 const char *extCGNS = ".cgns"; 3261 const char *extExodus = ".exo"; 3262 const char *extGenesis = ".gen"; 3263 const char *extFluent = ".cas"; 3264 const char *extHDF5 = ".h5"; 3265 const char *extMed = ".med"; 3266 const char *extPLY = ".ply"; 3267 const char *extCV = ".dat"; 3268 size_t len; 3269 PetscBool isGmsh, isGmsh2, isGmsh4, isCGNS, isExodus, isGenesis, isFluent, isHDF5, isMed, isPLY, isCV; 3270 PetscMPIInt rank; 3271 PetscErrorCode ierr; 3272 3273 PetscFunctionBegin; 3274 PetscValidPointer(filename, 2); 3275 PetscValidPointer(dm, 4); 3276 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 3277 ierr = PetscStrlen(filename, &len);CHKERRQ(ierr); 3278 if (!len) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Filename must be a valid path"); 3279 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGmsh, 4, &isGmsh);CHKERRQ(ierr); 3280 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extGmsh2, 5, &isGmsh2);CHKERRQ(ierr); 3281 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extGmsh4, 5, &isGmsh4);CHKERRQ(ierr); 3282 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extCGNS, 5, &isCGNS);CHKERRQ(ierr); 3283 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extExodus, 4, &isExodus);CHKERRQ(ierr); 3284 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGenesis, 4, &isGenesis);CHKERRQ(ierr); 3285 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extFluent, 4, &isFluent);CHKERRQ(ierr); 3286 ierr = PetscStrncmp(&filename[PetscMax(0,len-3)], extHDF5, 3, &isHDF5);CHKERRQ(ierr); 3287 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extMed, 4, &isMed);CHKERRQ(ierr); 3288 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extPLY, 4, &isPLY);CHKERRQ(ierr); 3289 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extCV, 4, &isCV);CHKERRQ(ierr); 3290 if (isGmsh || isGmsh2 || isGmsh4) { 3291 ierr = DMPlexCreateGmshFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3292 } else if (isCGNS) { 3293 ierr = DMPlexCreateCGNSFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3294 } else if (isExodus || isGenesis) { 3295 ierr = DMPlexCreateExodusFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3296 } else if (isFluent) { 3297 ierr = DMPlexCreateFluentFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3298 } else if (isHDF5) { 3299 PetscBool load_hdf5_xdmf = PETSC_FALSE; 3300 PetscViewer viewer; 3301 3302 /* PETSC_VIEWER_HDF5_XDMF is used if the filename ends with .xdmf.h5, or if -dm_plex_create_from_hdf5_xdmf option is present */ 3303 ierr = PetscStrncmp(&filename[PetscMax(0,len-8)], ".xdmf", 5, &load_hdf5_xdmf);CHKERRQ(ierr); 3304 ierr = PetscOptionsGetBool(NULL, NULL, "-dm_plex_create_from_hdf5_xdmf", &load_hdf5_xdmf, NULL);CHKERRQ(ierr); 3305 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 3306 ierr = PetscViewerSetType(viewer, PETSCVIEWERHDF5);CHKERRQ(ierr); 3307 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 3308 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 3309 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 3310 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 3311 if (load_hdf5_xdmf) {ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_XDMF);CHKERRQ(ierr);} 3312 ierr = DMLoad(*dm, viewer);CHKERRQ(ierr); 3313 if (load_hdf5_xdmf) {ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);} 3314 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3315 3316 if (interpolate) { 3317 DM idm; 3318 3319 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 3320 ierr = DMDestroy(dm);CHKERRQ(ierr); 3321 *dm = idm; 3322 } 3323 } else if (isMed) { 3324 ierr = DMPlexCreateMedFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3325 } else if (isPLY) { 3326 ierr = DMPlexCreatePLYFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3327 } else if (isCV) { 3328 ierr = DMPlexCreateCellVertexFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3329 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot load file %s: unrecognized extension", filename); 3330 PetscFunctionReturn(0); 3331 } 3332 3333 /*@ 3334 DMPlexCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell 3335 3336 Collective on comm 3337 3338 Input Parameters: 3339 + comm - The communicator 3340 . dim - The spatial dimension 3341 - simplex - Flag for simplex, otherwise use a tensor-product cell 3342 3343 Output Parameter: 3344 . refdm - The reference cell 3345 3346 Level: intermediate 3347 3348 .keywords: reference cell 3349 .seealso: 3350 @*/ 3351 PetscErrorCode DMPlexCreateReferenceCell(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *refdm) 3352 { 3353 DM rdm; 3354 Vec coords; 3355 PetscErrorCode ierr; 3356 3357 PetscFunctionBeginUser; 3358 ierr = DMCreate(comm, &rdm);CHKERRQ(ierr); 3359 ierr = DMSetType(rdm, DMPLEX);CHKERRQ(ierr); 3360 ierr = DMSetDimension(rdm, dim);CHKERRQ(ierr); 3361 switch (dim) { 3362 case 0: 3363 { 3364 PetscInt numPoints[1] = {1}; 3365 PetscInt coneSize[1] = {0}; 3366 PetscInt cones[1] = {0}; 3367 PetscInt coneOrientations[1] = {0}; 3368 PetscScalar vertexCoords[1] = {0.0}; 3369 3370 ierr = DMPlexCreateFromDAG(rdm, 0, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3371 } 3372 break; 3373 case 1: 3374 { 3375 PetscInt numPoints[2] = {2, 1}; 3376 PetscInt coneSize[3] = {2, 0, 0}; 3377 PetscInt cones[2] = {1, 2}; 3378 PetscInt coneOrientations[2] = {0, 0}; 3379 PetscScalar vertexCoords[2] = {-1.0, 1.0}; 3380 3381 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3382 } 3383 break; 3384 case 2: 3385 if (simplex) { 3386 PetscInt numPoints[2] = {3, 1}; 3387 PetscInt coneSize[4] = {3, 0, 0, 0}; 3388 PetscInt cones[3] = {1, 2, 3}; 3389 PetscInt coneOrientations[3] = {0, 0, 0}; 3390 PetscScalar vertexCoords[6] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0}; 3391 3392 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3393 } else { 3394 PetscInt numPoints[2] = {4, 1}; 3395 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 3396 PetscInt cones[4] = {1, 2, 3, 4}; 3397 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 3398 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0}; 3399 3400 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3401 } 3402 break; 3403 case 3: 3404 if (simplex) { 3405 PetscInt numPoints[2] = {4, 1}; 3406 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 3407 PetscInt cones[4] = {1, 3, 2, 4}; 3408 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 3409 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}; 3410 3411 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3412 } else { 3413 PetscInt numPoints[2] = {8, 1}; 3414 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0}; 3415 PetscInt cones[8] = {1, 4, 3, 2, 5, 6, 7, 8}; 3416 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 3417 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, 3418 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0}; 3419 3420 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3421 } 3422 break; 3423 default: 3424 SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Cannot create reference cell for dimension %d", dim); 3425 } 3426 ierr = DMPlexInterpolate(rdm, refdm);CHKERRQ(ierr); 3427 if (rdm->coordinateDM) { 3428 DM ncdm; 3429 PetscSection cs; 3430 PetscInt pEnd = -1; 3431 3432 ierr = DMGetSection(rdm->coordinateDM, &cs);CHKERRQ(ierr); 3433 if (cs) {ierr = PetscSectionGetChart(cs, NULL, &pEnd);CHKERRQ(ierr);} 3434 if (pEnd >= 0) { 3435 ierr = DMClone(rdm->coordinateDM, &ncdm);CHKERRQ(ierr); 3436 ierr = DMCopyDisc(rdm->coordinateDM, ncdm);CHKERRQ(ierr); 3437 ierr = DMSetSection(ncdm, cs);CHKERRQ(ierr); 3438 ierr = DMSetCoordinateDM(*refdm, ncdm);CHKERRQ(ierr); 3439 ierr = DMDestroy(&ncdm);CHKERRQ(ierr); 3440 } 3441 } 3442 ierr = DMGetCoordinatesLocal(rdm, &coords);CHKERRQ(ierr); 3443 if (coords) { 3444 ierr = DMSetCoordinatesLocal(*refdm, coords);CHKERRQ(ierr); 3445 } else { 3446 ierr = DMGetCoordinates(rdm, &coords);CHKERRQ(ierr); 3447 if (coords) {ierr = DMSetCoordinates(*refdm, coords);CHKERRQ(ierr);} 3448 } 3449 ierr = DMDestroy(&rdm);CHKERRQ(ierr); 3450 PetscFunctionReturn(0); 3451 } 3452