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 n = 3; 1057 ierr = PetscOptionsGetIntArray(NULL, NULL, "-dm_plex_box_faces", fac, &n, &flg);CHKERRQ(ierr); 1058 for (i = 0; i < dim; ++i) fac[i] = faces ? faces[i] : (flg && i < n ? fac[i] : (dim == 1 ? 1 : 4-dim)); 1059 if (lower) for (i = 0; i < dim; ++i) low[i] = lower[i]; 1060 if (upper) for (i = 0; i < dim; ++i) upp[i] = upper[i]; 1061 if (periodicity) for (i = 0; i < dim; ++i) bdt[i] = periodicity[i]; 1062 /* Allow bounds to be specified from the command line */ 1063 n = 3; 1064 ierr = PetscOptionsGetRealArray(NULL, NULL, "-dm_plex_box_lower", low, &n, &flg);CHKERRQ(ierr); 1065 if (flg && (n != dim)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Lower box point had %D values, should have been %D", n, dim); 1066 n = 3; 1067 ierr = PetscOptionsGetRealArray(NULL, NULL, "-dm_plex_box_upper", upp, &n, &flg);CHKERRQ(ierr); 1068 if (flg && (n != dim)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Upper box point had %D values, should have been %D", n, dim); 1069 1070 if (dim == 1) {ierr = DMPlexCreateLineMesh_Internal(comm, fac[0], low[0], upp[0], bdt[0], dm);CHKERRQ(ierr);} 1071 else if (simplex) {ierr = DMPlexCreateBoxMesh_Simplex_Internal(comm, dim, fac, low, upp, bdt, interpolate, dm);CHKERRQ(ierr);} 1072 else {ierr = DMPlexCreateBoxMesh_Tensor_Internal(comm, dim, fac, low, upp, bdt, interpolate, dm);CHKERRQ(ierr);} 1073 PetscFunctionReturn(0); 1074 } 1075 1076 /*@ 1077 DMPlexCreateWedgeBoxMesh - Creates a 3-D mesh tesselating the (x,y) plane and extruding in the third direction using wedge cells. 1078 1079 Collective on MPI_Comm 1080 1081 Input Parameters: 1082 + comm - The communicator for the DM object 1083 . faces - Number of faces per dimension, or NULL for (1, 1, 1) 1084 . lower - The lower left corner, or NULL for (0, 0, 0) 1085 . upper - The upper right corner, or NULL for (1, 1, 1) 1086 . periodicity - The boundary type for the X,Y,Z direction, or NULL for DM_BOUNDARY_NONE 1087 . ordExt - If PETSC_TRUE, orders the extruded cells in the height first. Otherwise, orders the cell on the layers first 1088 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 1089 1090 Output Parameter: 1091 . dm - The DM object 1092 1093 Level: beginner 1094 1095 .keywords: DM, create 1096 .seealso: DMPlexCreateHexCylinderMesh(), DMPlexCreateWedgeCylinderMesh(), DMPlexExtrude(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 1097 @*/ 1098 PetscErrorCode DMPlexCreateWedgeBoxMesh(MPI_Comm comm, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool ordExt, PetscBool interpolate, DM *dm) 1099 { 1100 DM bdm, botdm; 1101 PetscInt i; 1102 PetscInt fac[3] = {0, 0, 0}; 1103 PetscReal low[3] = {0, 0, 0}; 1104 PetscReal upp[3] = {1, 1, 1}; 1105 DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE}; 1106 PetscErrorCode ierr; 1107 1108 PetscFunctionBegin; 1109 for (i = 0; i < 3; ++i) fac[i] = faces ? (faces[i] > 0 ? faces[i] : 1) : 1; 1110 if (lower) for (i = 0; i < 3; ++i) low[i] = lower[i]; 1111 if (upper) for (i = 0; i < 3; ++i) upp[i] = upper[i]; 1112 if (periodicity) for (i = 0; i < 3; ++i) bdt[i] = periodicity[i]; 1113 for (i = 0; i < 3; ++i) if (bdt[i] != DM_BOUNDARY_NONE) SETERRQ(comm, PETSC_ERR_SUP, "Periodicity not yet supported"); 1114 1115 ierr = DMCreate(comm, &bdm);CHKERRQ(ierr); 1116 ierr = DMSetType(bdm, DMPLEX);CHKERRQ(ierr); 1117 ierr = DMSetDimension(bdm, 1);CHKERRQ(ierr); 1118 ierr = DMSetCoordinateDim(bdm, 2);CHKERRQ(ierr); 1119 ierr = DMPlexCreateSquareBoundary(bdm, low, upp, fac);CHKERRQ(ierr); 1120 ierr = DMPlexGenerate(bdm, NULL, PETSC_FALSE, &botdm);CHKERRQ(ierr); 1121 ierr = DMDestroy(&bdm);CHKERRQ(ierr); 1122 ierr = DMPlexExtrude(botdm, fac[2], upp[2] - low[2], ordExt, interpolate, dm);CHKERRQ(ierr); 1123 if (low[2] != 0.0) { 1124 Vec v; 1125 PetscScalar *x; 1126 PetscInt cDim, n; 1127 1128 ierr = DMGetCoordinatesLocal(*dm, &v);CHKERRQ(ierr); 1129 ierr = VecGetBlockSize(v, &cDim);CHKERRQ(ierr); 1130 ierr = VecGetLocalSize(v, &n);CHKERRQ(ierr); 1131 ierr = VecGetArray(v, &x);CHKERRQ(ierr); 1132 x += cDim; 1133 for (i=0; i<n; i+=cDim) x[i] += low[2]; 1134 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1135 ierr = DMSetCoordinatesLocal(*dm, v);CHKERRQ(ierr); 1136 } 1137 ierr = DMDestroy(&botdm);CHKERRQ(ierr); 1138 PetscFunctionReturn(0); 1139 } 1140 1141 /*@ 1142 DMPlexExtrude - Creates a (d+1)-D mesh by extruding a d-D mesh in the normal direction using prismatic cells. 1143 1144 Collective on idm 1145 1146 Input Parameters: 1147 + idm - The mesh to be extruted 1148 . layers - The number of layers 1149 . height - The height of the extruded layer 1150 . ordExt - If PETSC_TRUE, orders the extruded cells in the height first. Otherwise, orders the cell on the layers first 1151 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 1152 1153 Output Parameter: 1154 . dm - The DM object 1155 1156 Notes: The object created is an hybrid mesh, the vertex ordering in the cone of the cell is that of the prismatic cells 1157 1158 Level: advanced 1159 1160 .keywords: DM, create 1161 .seealso: DMPlexCreateWedgeCylinderMesh(), DMPlexCreateWedgeBoxMesh(), DMPlexSetHybridBounds(), DMSetType(), DMCreate() 1162 @*/ 1163 PetscErrorCode DMPlexExtrude(DM idm, PetscInt layers, PetscReal height, PetscBool ordExt, PetscBool interpolate, DM* dm) 1164 { 1165 PetscScalar *coordsB; 1166 const PetscScalar *coordsA; 1167 PetscReal *normals = NULL; 1168 Vec coordinatesA, coordinatesB; 1169 PetscSection coordSectionA, coordSectionB; 1170 PetscInt dim, cDim, cDimB, c, l, v, coordSize, *newCone; 1171 PetscInt cStart, cEnd, vStart, vEnd, cellV, numCells, numVertices; 1172 PetscErrorCode ierr; 1173 1174 PetscFunctionBegin; 1175 PetscValidHeaderSpecific(idm, DM_CLASSID, 1); 1176 PetscValidLogicalCollectiveInt(idm, layers, 2); 1177 PetscValidLogicalCollectiveReal(idm, height, 3); 1178 PetscValidLogicalCollectiveBool(idm, interpolate, 4); 1179 ierr = DMGetDimension(idm, &dim);CHKERRQ(ierr); 1180 if (dim < 1 || dim > 3) SETERRQ1(PetscObjectComm((PetscObject)idm), PETSC_ERR_SUP, "Support for dimension %D not coded", dim); 1181 1182 ierr = DMPlexGetHeightStratum(idm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1183 ierr = DMPlexGetDepthStratum(idm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1184 numCells = (cEnd - cStart)*layers; 1185 numVertices = (vEnd - vStart)*(layers+1); 1186 ierr = DMCreate(PetscObjectComm((PetscObject)idm), dm);CHKERRQ(ierr); 1187 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1188 ierr = DMSetDimension(*dm, dim+1);CHKERRQ(ierr); 1189 ierr = DMPlexSetChart(*dm, 0, numCells+numVertices);CHKERRQ(ierr); 1190 for (c = cStart, cellV = 0; c < cEnd; ++c) { 1191 PetscInt *closure = NULL; 1192 PetscInt closureSize, numCorners = 0; 1193 1194 ierr = DMPlexGetTransitiveClosure(idm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1195 for (v = 0; v < closureSize*2; v += 2) if ((closure[v] >= vStart) && (closure[v] < vEnd)) numCorners++; 1196 ierr = DMPlexRestoreTransitiveClosure(idm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1197 for (l = 0; l < layers; ++l) { 1198 ierr = DMPlexSetConeSize(*dm, ordExt ? layers*(c - cStart) + l : l*(cEnd - cStart) + c - cStart, 2*numCorners);CHKERRQ(ierr); 1199 } 1200 cellV = PetscMax(numCorners,cellV); 1201 } 1202 ierr = DMPlexSetHybridBounds(*dm, 0, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 1203 ierr = DMSetUp(*dm);CHKERRQ(ierr); 1204 1205 ierr = DMGetCoordinateDim(idm, &cDim);CHKERRQ(ierr); 1206 if (dim != cDim) { 1207 ierr = PetscCalloc1(cDim*(vEnd - vStart), &normals);CHKERRQ(ierr); 1208 } 1209 ierr = PetscMalloc1(3*cellV,&newCone);CHKERRQ(ierr); 1210 for (c = cStart; c < cEnd; ++c) { 1211 PetscInt *closure = NULL; 1212 PetscInt closureSize, numCorners = 0, l; 1213 PetscReal normal[3] = {0, 0, 0}; 1214 1215 if (normals) { 1216 ierr = DMPlexComputeCellGeometryFVM(idm, c, NULL, NULL, normal);CHKERRQ(ierr); 1217 } 1218 ierr = DMPlexGetTransitiveClosure(idm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1219 for (v = 0; v < closureSize*2; v += 2) { 1220 if ((closure[v] >= vStart) && (closure[v] < vEnd)) { 1221 PetscInt d; 1222 1223 newCone[numCorners++] = closure[v] - vStart; 1224 if (normals) { for (d = 0; d < cDim; ++d) normals[cDim*(closure[v]-vStart)+d] += normal[d]; } 1225 } 1226 } 1227 switch (numCorners) { 1228 case 4: /* do nothing */ 1229 case 2: /* do nothing */ 1230 break; 1231 case 3: /* from counter-clockwise to wedge ordering */ 1232 l = newCone[1]; 1233 newCone[1] = newCone[2]; 1234 newCone[2] = l; 1235 break; 1236 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported number of corners: %D", numCorners); 1237 } 1238 ierr = DMPlexRestoreTransitiveClosure(idm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1239 for (l = 0; l < layers; ++l) { 1240 PetscInt i; 1241 1242 for (i = 0; i < numCorners; ++i) { 1243 newCone[ numCorners + i] = ordExt ? (layers+1)*newCone[i] + l + numCells : l*(vEnd - vStart) + newCone[i] + numCells; 1244 newCone[2*numCorners + i] = ordExt ? (layers+1)*newCone[i] + l + 1 + numCells : (l+1)*(vEnd - vStart) + newCone[i] + numCells; 1245 } 1246 ierr = DMPlexSetCone(*dm, ordExt ? layers*(c - cStart) + l : l*(cEnd - cStart) + c - cStart, newCone + numCorners);CHKERRQ(ierr); 1247 } 1248 } 1249 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1250 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1251 ierr = PetscFree(newCone);CHKERRQ(ierr); 1252 1253 cDimB = cDim == dim ? cDim+1 : cDim; 1254 ierr = DMGetCoordinateSection(*dm, &coordSectionB);CHKERRQ(ierr); 1255 ierr = PetscSectionSetNumFields(coordSectionB, 1);CHKERRQ(ierr); 1256 ierr = PetscSectionSetFieldComponents(coordSectionB, 0, cDimB);CHKERRQ(ierr); 1257 ierr = PetscSectionSetChart(coordSectionB, numCells, numCells+numVertices);CHKERRQ(ierr); 1258 for (v = numCells; v < numCells+numVertices; ++v) { 1259 ierr = PetscSectionSetDof(coordSectionB, v, cDimB);CHKERRQ(ierr); 1260 ierr = PetscSectionSetFieldDof(coordSectionB, v, 0, cDimB);CHKERRQ(ierr); 1261 } 1262 ierr = PetscSectionSetUp(coordSectionB);CHKERRQ(ierr); 1263 ierr = PetscSectionGetStorageSize(coordSectionB, &coordSize);CHKERRQ(ierr); 1264 ierr = VecCreate(PETSC_COMM_SELF, &coordinatesB);CHKERRQ(ierr); 1265 ierr = PetscObjectSetName((PetscObject) coordinatesB, "coordinates");CHKERRQ(ierr); 1266 ierr = VecSetSizes(coordinatesB, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1267 ierr = VecSetBlockSize(coordinatesB, cDimB);CHKERRQ(ierr); 1268 ierr = VecSetType(coordinatesB,VECSTANDARD);CHKERRQ(ierr); 1269 1270 ierr = DMGetCoordinateSection(idm, &coordSectionA);CHKERRQ(ierr); 1271 ierr = DMGetCoordinatesLocal(idm, &coordinatesA);CHKERRQ(ierr); 1272 ierr = VecGetArray(coordinatesB, &coordsB);CHKERRQ(ierr); 1273 ierr = VecGetArrayRead(coordinatesA, &coordsA);CHKERRQ(ierr); 1274 for (v = vStart; v < vEnd; ++v) { 1275 const PetscScalar *cptr; 1276 PetscReal ones2[2] = { 0., 1.}, ones3[3] = { 0., 0., 1.}; 1277 PetscReal *normal, norm, h = height/layers; 1278 PetscInt offA, d, cDimA = cDim; 1279 1280 normal = normals ? normals + cDimB*(v - vStart) : (cDim > 1 ? ones3 : ones2); 1281 if (normals) { 1282 for (d = 0, norm = 0.0; d < cDimB; ++d) norm += normal[d]*normal[d]; 1283 for (d = 0; d < cDimB; ++d) normal[d] *= 1./PetscSqrtReal(norm); 1284 } 1285 1286 ierr = PetscSectionGetOffset(coordSectionA, v, &offA);CHKERRQ(ierr); 1287 cptr = coordsA + offA; 1288 for (l = 0; l < layers+1; ++l) { 1289 PetscInt offB, d, newV; 1290 1291 newV = ordExt ? (layers+1)*(v -vStart) + l + numCells : (vEnd -vStart)*l + (v -vStart) + numCells; 1292 ierr = PetscSectionGetOffset(coordSectionB, newV, &offB);CHKERRQ(ierr); 1293 for (d = 0; d < cDimA; ++d) { coordsB[offB+d] = cptr[d]; } 1294 for (d = 0; d < cDimB; ++d) { coordsB[offB+d] += l ? normal[d]*h : 0.0; } 1295 cptr = coordsB + offB; 1296 cDimA = cDimB; 1297 } 1298 } 1299 ierr = VecRestoreArrayRead(coordinatesA, &coordsA);CHKERRQ(ierr); 1300 ierr = VecRestoreArray(coordinatesB, &coordsB);CHKERRQ(ierr); 1301 ierr = DMSetCoordinatesLocal(*dm, coordinatesB);CHKERRQ(ierr); 1302 ierr = VecDestroy(&coordinatesB);CHKERRQ(ierr); 1303 ierr = PetscFree(normals);CHKERRQ(ierr); 1304 if (interpolate) { 1305 DM idm; 1306 1307 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1308 ierr = DMPlexCopyCoordinates(*dm, idm);CHKERRQ(ierr); 1309 ierr = DMDestroy(dm);CHKERRQ(ierr); 1310 *dm = idm; 1311 } 1312 PetscFunctionReturn(0); 1313 } 1314 1315 /*@C 1316 DMPlexSetOptionsPrefix - Sets the prefix used for searching for all DM options in the database. 1317 1318 Logically Collective on DM 1319 1320 Input Parameters: 1321 + dm - the DM context 1322 - prefix - the prefix to prepend to all option names 1323 1324 Notes: 1325 A hyphen (-) must NOT be given at the beginning of the prefix name. 1326 The first character of all runtime options is AUTOMATICALLY the hyphen. 1327 1328 Level: advanced 1329 1330 .seealso: SNESSetFromOptions() 1331 @*/ 1332 PetscErrorCode DMPlexSetOptionsPrefix(DM dm, const char prefix[]) 1333 { 1334 DM_Plex *mesh = (DM_Plex *) dm->data; 1335 PetscErrorCode ierr; 1336 1337 PetscFunctionBegin; 1338 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1339 ierr = PetscObjectSetOptionsPrefix((PetscObject) dm, prefix);CHKERRQ(ierr); 1340 ierr = PetscObjectSetOptionsPrefix((PetscObject) mesh->partitioner, prefix);CHKERRQ(ierr); 1341 PetscFunctionReturn(0); 1342 } 1343 1344 /*@ 1345 DMPlexCreateHexCylinderMesh - Creates a mesh on the tensor product of the unit interval with the circle (cylinder) using hexahedra. 1346 1347 Collective on MPI_Comm 1348 1349 Input Parameters: 1350 + comm - The communicator for the DM object 1351 . numRefine - The number of regular refinements to the basic 5 cell structure 1352 - periodicZ - The boundary type for the Z direction 1353 1354 Output Parameter: 1355 . dm - The DM object 1356 1357 Note: Here is the output numbering looking from the bottom of the cylinder: 1358 $ 17-----14 1359 $ | | 1360 $ | 2 | 1361 $ | | 1362 $ 17-----8-----7-----14 1363 $ | | | | 1364 $ | 3 | 0 | 1 | 1365 $ | | | | 1366 $ 19-----5-----6-----13 1367 $ | | 1368 $ | 4 | 1369 $ | | 1370 $ 19-----13 1371 $ 1372 $ and up through the top 1373 $ 1374 $ 18-----16 1375 $ | | 1376 $ | 2 | 1377 $ | | 1378 $ 18----10----11-----16 1379 $ | | | | 1380 $ | 3 | 0 | 1 | 1381 $ | | | | 1382 $ 20-----9----12-----15 1383 $ | | 1384 $ | 4 | 1385 $ | | 1386 $ 20-----15 1387 1388 Level: beginner 1389 1390 .keywords: DM, create 1391 .seealso: DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 1392 @*/ 1393 PetscErrorCode DMPlexCreateHexCylinderMesh(MPI_Comm comm, PetscInt numRefine, DMBoundaryType periodicZ, DM *dm) 1394 { 1395 const PetscInt dim = 3; 1396 PetscInt numCells, numVertices, r; 1397 PetscMPIInt rank; 1398 PetscErrorCode ierr; 1399 1400 PetscFunctionBegin; 1401 PetscValidPointer(dm, 4); 1402 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 1403 if (numRefine < 0) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Number of refinements %D cannot be negative", numRefine); 1404 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1405 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1406 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1407 /* Create topology */ 1408 { 1409 PetscInt cone[8], c; 1410 1411 numCells = !rank ? 5 : 0; 1412 numVertices = !rank ? 16 : 0; 1413 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1414 numCells *= 3; 1415 numVertices = !rank ? 24 : 0; 1416 } 1417 ierr = DMPlexSetChart(*dm, 0, numCells+numVertices);CHKERRQ(ierr); 1418 for (c = 0; c < numCells; c++) {ierr = DMPlexSetConeSize(*dm, c, 8);CHKERRQ(ierr);} 1419 ierr = DMSetUp(*dm);CHKERRQ(ierr); 1420 if (!rank) { 1421 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1422 cone[0] = 15; cone[1] = 18; cone[2] = 17; cone[3] = 16; 1423 cone[4] = 31; cone[5] = 32; cone[6] = 33; cone[7] = 34; 1424 ierr = DMPlexSetCone(*dm, 0, cone);CHKERRQ(ierr); 1425 cone[0] = 16; cone[1] = 17; cone[2] = 24; cone[3] = 23; 1426 cone[4] = 32; cone[5] = 36; cone[6] = 37; cone[7] = 33; /* 22 25 26 21 */ 1427 ierr = DMPlexSetCone(*dm, 1, cone);CHKERRQ(ierr); 1428 cone[0] = 18; cone[1] = 27; cone[2] = 24; cone[3] = 17; 1429 cone[4] = 34; cone[5] = 33; cone[6] = 37; cone[7] = 38; 1430 ierr = DMPlexSetCone(*dm, 2, cone);CHKERRQ(ierr); 1431 cone[0] = 29; cone[1] = 27; cone[2] = 18; cone[3] = 15; 1432 cone[4] = 35; cone[5] = 31; cone[6] = 34; cone[7] = 38; 1433 ierr = DMPlexSetCone(*dm, 3, cone);CHKERRQ(ierr); 1434 cone[0] = 29; cone[1] = 15; cone[2] = 16; cone[3] = 23; 1435 cone[4] = 35; cone[5] = 36; cone[6] = 32; cone[7] = 31; 1436 ierr = DMPlexSetCone(*dm, 4, cone);CHKERRQ(ierr); 1437 1438 cone[0] = 31; cone[1] = 34; cone[2] = 33; cone[3] = 32; 1439 cone[4] = 19; cone[5] = 22; cone[6] = 21; cone[7] = 20; 1440 ierr = DMPlexSetCone(*dm, 5, cone);CHKERRQ(ierr); 1441 cone[0] = 32; cone[1] = 33; cone[2] = 37; cone[3] = 36; 1442 cone[4] = 22; cone[5] = 25; cone[6] = 26; cone[7] = 21; 1443 ierr = DMPlexSetCone(*dm, 6, cone);CHKERRQ(ierr); 1444 cone[0] = 34; cone[1] = 38; cone[2] = 37; cone[3] = 33; 1445 cone[4] = 20; cone[5] = 21; cone[6] = 26; cone[7] = 28; 1446 ierr = DMPlexSetCone(*dm, 7, cone);CHKERRQ(ierr); 1447 cone[0] = 35; cone[1] = 38; cone[2] = 34; cone[3] = 31; 1448 cone[4] = 30; cone[5] = 19; cone[6] = 20; cone[7] = 28; 1449 ierr = DMPlexSetCone(*dm, 8, cone);CHKERRQ(ierr); 1450 cone[0] = 35; cone[1] = 31; cone[2] = 32; cone[3] = 36; 1451 cone[4] = 30; cone[5] = 25; cone[6] = 22; cone[7] = 19; 1452 ierr = DMPlexSetCone(*dm, 9, cone);CHKERRQ(ierr); 1453 1454 cone[0] = 19; cone[1] = 20; cone[2] = 21; cone[3] = 22; 1455 cone[4] = 15; cone[5] = 16; cone[6] = 17; cone[7] = 18; 1456 ierr = DMPlexSetCone(*dm, 10, cone);CHKERRQ(ierr); 1457 cone[0] = 22; cone[1] = 21; cone[2] = 26; cone[3] = 25; 1458 cone[4] = 16; cone[5] = 23; cone[6] = 24; cone[7] = 17; 1459 ierr = DMPlexSetCone(*dm, 11, cone);CHKERRQ(ierr); 1460 cone[0] = 20; cone[1] = 28; cone[2] = 26; cone[3] = 21; 1461 cone[4] = 18; cone[5] = 17; cone[6] = 24; cone[7] = 27; 1462 ierr = DMPlexSetCone(*dm, 12, cone);CHKERRQ(ierr); 1463 cone[0] = 30; cone[1] = 28; cone[2] = 20; cone[3] = 19; 1464 cone[4] = 29; cone[5] = 15; cone[6] = 18; cone[7] = 27; 1465 ierr = DMPlexSetCone(*dm, 13, cone);CHKERRQ(ierr); 1466 cone[0] = 30; cone[1] = 19; cone[2] = 22; cone[3] = 25; 1467 cone[4] = 29; cone[5] = 23; cone[6] = 16; cone[7] = 15; 1468 ierr = DMPlexSetCone(*dm, 14, cone);CHKERRQ(ierr); 1469 } else { 1470 cone[0] = 5; cone[1] = 8; cone[2] = 7; cone[3] = 6; 1471 cone[4] = 9; cone[5] = 12; cone[6] = 11; cone[7] = 10; 1472 ierr = DMPlexSetCone(*dm, 0, cone);CHKERRQ(ierr); 1473 cone[0] = 6; cone[1] = 7; cone[2] = 14; cone[3] = 13; 1474 cone[4] = 12; cone[5] = 15; cone[6] = 16; cone[7] = 11; 1475 ierr = DMPlexSetCone(*dm, 1, cone);CHKERRQ(ierr); 1476 cone[0] = 8; cone[1] = 17; cone[2] = 14; cone[3] = 7; 1477 cone[4] = 10; cone[5] = 11; cone[6] = 16; cone[7] = 18; 1478 ierr = DMPlexSetCone(*dm, 2, cone);CHKERRQ(ierr); 1479 cone[0] = 19; cone[1] = 17; cone[2] = 8; cone[3] = 5; 1480 cone[4] = 20; cone[5] = 9; cone[6] = 10; cone[7] = 18; 1481 ierr = DMPlexSetCone(*dm, 3, cone);CHKERRQ(ierr); 1482 cone[0] = 19; cone[1] = 5; cone[2] = 6; cone[3] = 13; 1483 cone[4] = 20; cone[5] = 15; cone[6] = 12; cone[7] = 9; 1484 ierr = DMPlexSetCone(*dm, 4, cone);CHKERRQ(ierr); 1485 } 1486 } 1487 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1488 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1489 } 1490 /* Interpolate */ 1491 { 1492 DM idm; 1493 1494 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1495 ierr = DMDestroy(dm);CHKERRQ(ierr); 1496 *dm = idm; 1497 } 1498 /* Create cube geometry */ 1499 { 1500 Vec coordinates; 1501 PetscSection coordSection; 1502 PetscScalar *coords; 1503 PetscInt coordSize, v; 1504 const PetscReal dis = 1.0/PetscSqrtReal(2.0); 1505 const PetscReal ds2 = dis/2.0; 1506 1507 /* Build coordinates */ 1508 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 1509 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1510 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 1511 ierr = PetscSectionSetChart(coordSection, numCells, numCells+numVertices);CHKERRQ(ierr); 1512 for (v = numCells; v < numCells+numVertices; ++v) { 1513 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 1514 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 1515 } 1516 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1517 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1518 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 1519 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1520 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1521 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 1522 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1523 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1524 if (!rank) { 1525 coords[0*dim+0] = -ds2; coords[0*dim+1] = -ds2; coords[0*dim+2] = 0.0; 1526 coords[1*dim+0] = ds2; coords[1*dim+1] = -ds2; coords[1*dim+2] = 0.0; 1527 coords[2*dim+0] = ds2; coords[2*dim+1] = ds2; coords[2*dim+2] = 0.0; 1528 coords[3*dim+0] = -ds2; coords[3*dim+1] = ds2; coords[3*dim+2] = 0.0; 1529 coords[4*dim+0] = -ds2; coords[4*dim+1] = -ds2; coords[4*dim+2] = 1.0; 1530 coords[5*dim+0] = -ds2; coords[5*dim+1] = ds2; coords[5*dim+2] = 1.0; 1531 coords[6*dim+0] = ds2; coords[6*dim+1] = ds2; coords[6*dim+2] = 1.0; 1532 coords[7*dim+0] = ds2; coords[7*dim+1] = -ds2; coords[7*dim+2] = 1.0; 1533 coords[ 8*dim+0] = dis; coords[ 8*dim+1] = -dis; coords[ 8*dim+2] = 0.0; 1534 coords[ 9*dim+0] = dis; coords[ 9*dim+1] = dis; coords[ 9*dim+2] = 0.0; 1535 coords[10*dim+0] = dis; coords[10*dim+1] = -dis; coords[10*dim+2] = 1.0; 1536 coords[11*dim+0] = dis; coords[11*dim+1] = dis; coords[11*dim+2] = 1.0; 1537 coords[12*dim+0] = -dis; coords[12*dim+1] = dis; coords[12*dim+2] = 0.0; 1538 coords[13*dim+0] = -dis; coords[13*dim+1] = dis; coords[13*dim+2] = 1.0; 1539 coords[14*dim+0] = -dis; coords[14*dim+1] = -dis; coords[14*dim+2] = 0.0; 1540 coords[15*dim+0] = -dis; coords[15*dim+1] = -dis; coords[15*dim+2] = 1.0; 1541 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1542 /* 15 31 19 */ coords[16*dim+0] = -ds2; coords[16*dim+1] = -ds2; coords[16*dim+2] = 0.5; 1543 /* 16 32 22 */ coords[17*dim+0] = ds2; coords[17*dim+1] = -ds2; coords[17*dim+2] = 0.5; 1544 /* 17 33 21 */ coords[18*dim+0] = ds2; coords[18*dim+1] = ds2; coords[18*dim+2] = 0.5; 1545 /* 18 34 20 */ coords[19*dim+0] = -ds2; coords[19*dim+1] = ds2; coords[19*dim+2] = 0.5; 1546 /* 29 35 30 */ coords[20*dim+0] = -dis; coords[20*dim+1] = -dis; coords[20*dim+2] = 0.5; 1547 /* 23 36 25 */ coords[21*dim+0] = dis; coords[21*dim+1] = -dis; coords[21*dim+2] = 0.5; 1548 /* 24 37 26 */ coords[22*dim+0] = dis; coords[22*dim+1] = dis; coords[22*dim+2] = 0.5; 1549 /* 27 38 28 */ coords[23*dim+0] = -dis; coords[23*dim+1] = dis; coords[23*dim+2] = 0.5; 1550 } 1551 } 1552 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1553 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 1554 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1555 } 1556 /* Create periodicity */ 1557 if (periodicZ == DM_BOUNDARY_PERIODIC || periodicZ == DM_BOUNDARY_TWIST) { 1558 PetscReal L[3]; 1559 PetscReal maxCell[3]; 1560 DMBoundaryType bdType[3]; 1561 PetscReal lower[3] = {0.0, 0.0, 0.0}; 1562 PetscReal upper[3] = {1.0, 1.0, 1.5}; 1563 PetscInt i, numZCells = 3; 1564 1565 bdType[0] = DM_BOUNDARY_NONE; 1566 bdType[1] = DM_BOUNDARY_NONE; 1567 bdType[2] = periodicZ; 1568 for (i = 0; i < dim; i++) { 1569 L[i] = upper[i] - lower[i]; 1570 maxCell[i] = 1.1 * (L[i] / numZCells); 1571 } 1572 ierr = DMSetPeriodicity(*dm, PETSC_TRUE, maxCell, L, bdType);CHKERRQ(ierr); 1573 } 1574 /* Refine topology */ 1575 for (r = 0; r < numRefine; ++r) { 1576 DM rdm = NULL; 1577 1578 ierr = DMRefine(*dm, comm, &rdm);CHKERRQ(ierr); 1579 ierr = DMDestroy(dm);CHKERRQ(ierr); 1580 *dm = rdm; 1581 } 1582 /* Remap geometry to cylinder 1583 Interior square: Linear interpolation is correct 1584 The other cells all have vertices on rays from the origin. We want to uniformly expand the spacing 1585 such that the last vertex is on the unit circle. So the closest and farthest vertices are at distance 1586 1587 phi = arctan(y/x) 1588 d_close = sqrt(1/8 + 1/4 sin^2(phi)) 1589 d_far = sqrt(1/2 + sin^2(phi)) 1590 1591 so we remap them using 1592 1593 x_new = x_close + (x - x_close) (1 - d_close) / (d_far - d_close) 1594 y_new = y_close + (y - y_close) (1 - d_close) / (d_far - d_close) 1595 1596 If pi/4 < phi < 3pi/4 or -3pi/4 < phi < -pi/4, then we switch x and y. 1597 */ 1598 { 1599 Vec coordinates; 1600 PetscSection coordSection; 1601 PetscScalar *coords; 1602 PetscInt vStart, vEnd, v; 1603 const PetscReal dis = 1.0/PetscSqrtReal(2.0); 1604 const PetscReal ds2 = 0.5*dis; 1605 1606 ierr = DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1607 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 1608 ierr = DMGetCoordinatesLocal(*dm, &coordinates);CHKERRQ(ierr); 1609 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1610 for (v = vStart; v < vEnd; ++v) { 1611 PetscReal phi, sinp, cosp, dc, df, x, y, xc, yc; 1612 PetscInt off; 1613 1614 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr); 1615 if ((PetscAbsScalar(coords[off+0]) <= ds2) && (PetscAbsScalar(coords[off+1]) <= ds2)) continue; 1616 x = PetscRealPart(coords[off]); 1617 y = PetscRealPart(coords[off+1]); 1618 phi = PetscAtan2Real(y, x); 1619 sinp = PetscSinReal(phi); 1620 cosp = PetscCosReal(phi); 1621 if ((PetscAbsReal(phi) > PETSC_PI/4.0) && (PetscAbsReal(phi) < 3.0*PETSC_PI/4.0)) { 1622 dc = PetscAbsReal(ds2/sinp); 1623 df = PetscAbsReal(dis/sinp); 1624 xc = ds2*x/PetscAbsReal(y); 1625 yc = ds2*PetscSignReal(y); 1626 } else { 1627 dc = PetscAbsReal(ds2/cosp); 1628 df = PetscAbsReal(dis/cosp); 1629 xc = ds2*PetscSignReal(x); 1630 yc = ds2*y/PetscAbsReal(x); 1631 } 1632 coords[off+0] = xc + (coords[off+0] - xc)*(1.0 - dc)/(df - dc); 1633 coords[off+1] = yc + (coords[off+1] - yc)*(1.0 - dc)/(df - dc); 1634 } 1635 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1636 if (periodicZ == DM_BOUNDARY_PERIODIC || periodicZ == DM_BOUNDARY_TWIST) { 1637 ierr = DMLocalizeCoordinates(*dm);CHKERRQ(ierr); 1638 } 1639 } 1640 PetscFunctionReturn(0); 1641 } 1642 1643 /*@ 1644 DMPlexCreateWedgeCylinderMesh - Creates a mesh on the tensor product of the unit interval with the circle (cylinder) using wedges. 1645 1646 Collective on MPI_Comm 1647 1648 Input Parameters: 1649 + comm - The communicator for the DM object 1650 . n - The number of wedges around the origin 1651 - interpolate - Create edges and faces 1652 1653 Output Parameter: 1654 . dm - The DM object 1655 1656 Level: beginner 1657 1658 .keywords: DM, create 1659 .seealso: DMPlexCreateHexCylinderMesh(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 1660 @*/ 1661 PetscErrorCode DMPlexCreateWedgeCylinderMesh(MPI_Comm comm, PetscInt n, PetscBool interpolate, DM *dm) 1662 { 1663 const PetscInt dim = 3; 1664 PetscInt numCells, numVertices; 1665 PetscMPIInt rank; 1666 PetscErrorCode ierr; 1667 1668 PetscFunctionBegin; 1669 PetscValidPointer(dm, 3); 1670 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 1671 if (n < 0) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Number of wedges %D cannot be negative", n); 1672 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1673 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1674 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1675 /* Create topology */ 1676 { 1677 PetscInt cone[6], c; 1678 1679 numCells = !rank ? n : 0; 1680 numVertices = !rank ? 2*(n+1) : 0; 1681 ierr = DMPlexSetChart(*dm, 0, numCells+numVertices);CHKERRQ(ierr); 1682 ierr = DMPlexSetHybridBounds(*dm, 0, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 1683 for (c = 0; c < numCells; c++) {ierr = DMPlexSetConeSize(*dm, c, 6);CHKERRQ(ierr);} 1684 ierr = DMSetUp(*dm);CHKERRQ(ierr); 1685 for (c = 0; c < numCells; c++) { 1686 cone[0] = c+n*1; cone[1] = (c+1)%n+n*1; cone[2] = 0+3*n; 1687 cone[3] = c+n*2; cone[4] = (c+1)%n+n*2; cone[5] = 1+3*n; 1688 ierr = DMPlexSetCone(*dm, c, cone);CHKERRQ(ierr); 1689 } 1690 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1691 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1692 } 1693 /* Interpolate */ 1694 if (interpolate) { 1695 DM idm; 1696 1697 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1698 ierr = DMDestroy(dm);CHKERRQ(ierr); 1699 *dm = idm; 1700 } 1701 /* Create cylinder geometry */ 1702 { 1703 Vec coordinates; 1704 PetscSection coordSection; 1705 PetscScalar *coords; 1706 PetscInt coordSize, v, c; 1707 1708 /* Build coordinates */ 1709 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 1710 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1711 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 1712 ierr = PetscSectionSetChart(coordSection, numCells, numCells+numVertices);CHKERRQ(ierr); 1713 for (v = numCells; v < numCells+numVertices; ++v) { 1714 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 1715 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 1716 } 1717 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1718 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1719 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 1720 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1721 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1722 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 1723 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1724 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1725 for (c = 0; c < numCells; c++) { 1726 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; 1727 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; 1728 } 1729 if (!rank) { 1730 coords[(2*n+0)*dim+0] = 0.0; coords[(2*n+0)*dim+1] = 0.0; coords[(2*n+0)*dim+2] = 1.0; 1731 coords[(2*n+1)*dim+0] = 0.0; coords[(2*n+1)*dim+1] = 0.0; coords[(2*n+1)*dim+2] = 0.0; 1732 } 1733 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1734 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 1735 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1736 } 1737 PetscFunctionReturn(0); 1738 } 1739 1740 PETSC_STATIC_INLINE PetscReal DiffNormReal(PetscInt dim, const PetscReal x[], const PetscReal y[]) 1741 { 1742 PetscReal prod = 0.0; 1743 PetscInt i; 1744 for (i = 0; i < dim; ++i) prod += PetscSqr(x[i] - y[i]); 1745 return PetscSqrtReal(prod); 1746 } 1747 PETSC_STATIC_INLINE PetscReal DotReal(PetscInt dim, const PetscReal x[], const PetscReal y[]) 1748 { 1749 PetscReal prod = 0.0; 1750 PetscInt i; 1751 for (i = 0; i < dim; ++i) prod += x[i]*y[i]; 1752 return prod; 1753 } 1754 1755 /*@ 1756 DMPlexCreateSphereMesh - Creates a mesh on the d-dimensional sphere, S^d. 1757 1758 Collective on MPI_Comm 1759 1760 Input Parameters: 1761 . comm - The communicator for the DM object 1762 . dim - The dimension 1763 - simplex - Use simplices, or tensor product cells 1764 1765 Output Parameter: 1766 . dm - The DM object 1767 1768 Level: beginner 1769 1770 .keywords: DM, create 1771 .seealso: DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 1772 @*/ 1773 PetscErrorCode DMPlexCreateSphereMesh(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *dm) 1774 { 1775 const PetscInt embedDim = dim+1; 1776 PetscSection coordSection; 1777 Vec coordinates; 1778 PetscScalar *coords; 1779 PetscReal *coordsIn; 1780 PetscInt numCells, numEdges, numVerts, firstVertex, v, firstEdge, coordSize, d, c, e; 1781 PetscMPIInt rank; 1782 PetscErrorCode ierr; 1783 1784 PetscFunctionBegin; 1785 PetscValidPointer(dm, 4); 1786 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1787 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1788 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1789 ierr = DMSetCoordinateDim(*dm, dim+1);CHKERRQ(ierr); 1790 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) *dm), &rank);CHKERRQ(ierr); 1791 switch (dim) { 1792 case 2: 1793 if (simplex) { 1794 DM idm; 1795 const PetscReal edgeLen = 2.0/(1.0 + PETSC_PHI); 1796 const PetscReal vertex[3] = {0.0, 1.0/(1.0 + PETSC_PHI), PETSC_PHI/(1.0 + PETSC_PHI)}; 1797 const PetscInt degree = 5; 1798 PetscInt s[3] = {1, 1, 1}; 1799 PetscInt cone[3]; 1800 PetscInt *graph, p, i, j, k; 1801 1802 numCells = !rank ? 20 : 0; 1803 numVerts = !rank ? 12 : 0; 1804 firstVertex = numCells; 1805 /* Use icosahedron, which for a unit sphere has coordinates which are all cyclic permutations of 1806 1807 (0, \pm 1/\phi+1, \pm \phi/\phi+1) 1808 1809 where \phi^2 - \phi - 1 = 0, meaning \phi is the golden ratio \frac{1 + \sqrt{5}}{2}. The edge 1810 length is then given by 2/\phi = 2 * 2.73606 = 5.47214. 1811 */ 1812 /* Construct vertices */ 1813 ierr = PetscCalloc1(numVerts * embedDim, &coordsIn);CHKERRQ(ierr); 1814 for (p = 0, i = 0; p < embedDim; ++p) { 1815 for (s[1] = -1; s[1] < 2; s[1] += 2) { 1816 for (s[2] = -1; s[2] < 2; s[2] += 2) { 1817 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[(d+p)%embedDim]*vertex[(d+p)%embedDim]; 1818 ++i; 1819 } 1820 } 1821 } 1822 /* Construct graph */ 1823 ierr = PetscCalloc1(numVerts * numVerts, &graph);CHKERRQ(ierr); 1824 for (i = 0; i < numVerts; ++i) { 1825 for (j = 0, k = 0; j < numVerts; ++j) { 1826 if (PetscAbsReal(DiffNormReal(embedDim, &coordsIn[i*embedDim], &coordsIn[j*embedDim]) - edgeLen) < PETSC_SMALL) {graph[i*numVerts+j] = 1; ++k;} 1827 } 1828 if (k != degree) SETERRQ3(comm, PETSC_ERR_PLIB, "Invalid icosahedron, vertex %D degree %D != %D", i, k, degree); 1829 } 1830 /* Build Topology */ 1831 ierr = DMPlexSetChart(*dm, 0, numCells+numVerts);CHKERRQ(ierr); 1832 for (c = 0; c < numCells; c++) { 1833 ierr = DMPlexSetConeSize(*dm, c, embedDim);CHKERRQ(ierr); 1834 } 1835 ierr = DMSetUp(*dm);CHKERRQ(ierr); /* Allocate space for cones */ 1836 /* Cells */ 1837 for (i = 0, c = 0; i < numVerts; ++i) { 1838 for (j = 0; j < i; ++j) { 1839 for (k = 0; k < j; ++k) { 1840 if (graph[i*numVerts+j] && graph[j*numVerts+k] && graph[k*numVerts+i]) { 1841 cone[0] = firstVertex+i; cone[1] = firstVertex+j; cone[2] = firstVertex+k; 1842 /* Check orientation */ 1843 { 1844 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}}}; 1845 PetscReal normal[3]; 1846 PetscInt e, f; 1847 1848 for (d = 0; d < embedDim; ++d) { 1849 normal[d] = 0.0; 1850 for (e = 0; e < embedDim; ++e) { 1851 for (f = 0; f < embedDim; ++f) { 1852 normal[d] += epsilon[d][e][f]*(coordsIn[j*embedDim+e] - coordsIn[i*embedDim+e])*(coordsIn[k*embedDim+f] - coordsIn[i*embedDim+f]); 1853 } 1854 } 1855 } 1856 if (DotReal(embedDim, normal, &coordsIn[i*embedDim]) < 0) {PetscInt tmp = cone[1]; cone[1] = cone[2]; cone[2] = tmp;} 1857 } 1858 ierr = DMPlexSetCone(*dm, c++, cone);CHKERRQ(ierr); 1859 } 1860 } 1861 } 1862 } 1863 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1864 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1865 ierr = PetscFree(graph);CHKERRQ(ierr); 1866 /* Interpolate mesh */ 1867 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1868 ierr = DMDestroy(dm);CHKERRQ(ierr); 1869 *dm = idm; 1870 } else { 1871 /* 1872 12-21--13 1873 | | 1874 25 4 24 1875 | | 1876 12-25--9-16--8-24--13 1877 | | | | 1878 23 5 17 0 15 3 22 1879 | | | | 1880 10-20--6-14--7-19--11 1881 | | 1882 20 1 19 1883 | | 1884 10-18--11 1885 | | 1886 23 2 22 1887 | | 1888 12-21--13 1889 */ 1890 const PetscReal dist = 1.0/PetscSqrtReal(3.0); 1891 PetscInt cone[4], ornt[4]; 1892 1893 numCells = !rank ? 6 : 0; 1894 numEdges = !rank ? 12 : 0; 1895 numVerts = !rank ? 8 : 0; 1896 firstVertex = numCells; 1897 firstEdge = numCells + numVerts; 1898 /* Build Topology */ 1899 ierr = DMPlexSetChart(*dm, 0, numCells+numEdges+numVerts);CHKERRQ(ierr); 1900 for (c = 0; c < numCells; c++) { 1901 ierr = DMPlexSetConeSize(*dm, c, 4);CHKERRQ(ierr); 1902 } 1903 for (e = firstEdge; e < firstEdge+numEdges; ++e) { 1904 ierr = DMPlexSetConeSize(*dm, e, 2);CHKERRQ(ierr); 1905 } 1906 ierr = DMSetUp(*dm);CHKERRQ(ierr); /* Allocate space for cones */ 1907 /* Cell 0 */ 1908 cone[0] = 14; cone[1] = 15; cone[2] = 16; cone[3] = 17; 1909 ierr = DMPlexSetCone(*dm, 0, cone);CHKERRQ(ierr); 1910 ornt[0] = 0; ornt[1] = 0; ornt[2] = 0; ornt[3] = 0; 1911 ierr = DMPlexSetConeOrientation(*dm, 0, ornt);CHKERRQ(ierr); 1912 /* Cell 1 */ 1913 cone[0] = 18; cone[1] = 19; cone[2] = 14; cone[3] = 20; 1914 ierr = DMPlexSetCone(*dm, 1, cone);CHKERRQ(ierr); 1915 ornt[0] = 0; ornt[1] = 0; ornt[2] = -2; ornt[3] = 0; 1916 ierr = DMPlexSetConeOrientation(*dm, 1, ornt);CHKERRQ(ierr); 1917 /* Cell 2 */ 1918 cone[0] = 21; cone[1] = 22; cone[2] = 18; cone[3] = 23; 1919 ierr = DMPlexSetCone(*dm, 2, cone);CHKERRQ(ierr); 1920 ornt[0] = 0; ornt[1] = 0; ornt[2] = -2; ornt[3] = 0; 1921 ierr = DMPlexSetConeOrientation(*dm, 2, ornt);CHKERRQ(ierr); 1922 /* Cell 3 */ 1923 cone[0] = 19; cone[1] = 22; cone[2] = 24; cone[3] = 15; 1924 ierr = DMPlexSetCone(*dm, 3, cone);CHKERRQ(ierr); 1925 ornt[0] = -2; ornt[1] = -2; ornt[2] = 0; ornt[3] = -2; 1926 ierr = DMPlexSetConeOrientation(*dm, 3, ornt);CHKERRQ(ierr); 1927 /* Cell 4 */ 1928 cone[0] = 16; cone[1] = 24; cone[2] = 21; cone[3] = 25; 1929 ierr = DMPlexSetCone(*dm, 4, cone);CHKERRQ(ierr); 1930 ornt[0] = -2; ornt[1] = -2; ornt[2] = -2; ornt[3] = 0; 1931 ierr = DMPlexSetConeOrientation(*dm, 4, ornt);CHKERRQ(ierr); 1932 /* Cell 5 */ 1933 cone[0] = 20; cone[1] = 17; cone[2] = 25; cone[3] = 23; 1934 ierr = DMPlexSetCone(*dm, 5, cone);CHKERRQ(ierr); 1935 ornt[0] = -2; ornt[1] = -2; ornt[2] = -2; ornt[3] = -2; 1936 ierr = DMPlexSetConeOrientation(*dm, 5, ornt);CHKERRQ(ierr); 1937 /* Edges */ 1938 cone[0] = 6; cone[1] = 7; 1939 ierr = DMPlexSetCone(*dm, 14, cone);CHKERRQ(ierr); 1940 cone[0] = 7; cone[1] = 8; 1941 ierr = DMPlexSetCone(*dm, 15, cone);CHKERRQ(ierr); 1942 cone[0] = 8; cone[1] = 9; 1943 ierr = DMPlexSetCone(*dm, 16, cone);CHKERRQ(ierr); 1944 cone[0] = 9; cone[1] = 6; 1945 ierr = DMPlexSetCone(*dm, 17, cone);CHKERRQ(ierr); 1946 cone[0] = 10; cone[1] = 11; 1947 ierr = DMPlexSetCone(*dm, 18, cone);CHKERRQ(ierr); 1948 cone[0] = 11; cone[1] = 7; 1949 ierr = DMPlexSetCone(*dm, 19, cone);CHKERRQ(ierr); 1950 cone[0] = 6; cone[1] = 10; 1951 ierr = DMPlexSetCone(*dm, 20, cone);CHKERRQ(ierr); 1952 cone[0] = 12; cone[1] = 13; 1953 ierr = DMPlexSetCone(*dm, 21, cone);CHKERRQ(ierr); 1954 cone[0] = 13; cone[1] = 11; 1955 ierr = DMPlexSetCone(*dm, 22, cone);CHKERRQ(ierr); 1956 cone[0] = 10; cone[1] = 12; 1957 ierr = DMPlexSetCone(*dm, 23, cone);CHKERRQ(ierr); 1958 cone[0] = 13; cone[1] = 8; 1959 ierr = DMPlexSetCone(*dm, 24, cone);CHKERRQ(ierr); 1960 cone[0] = 12; cone[1] = 9; 1961 ierr = DMPlexSetCone(*dm, 25, cone);CHKERRQ(ierr); 1962 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1963 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1964 /* Build coordinates */ 1965 ierr = PetscCalloc1(numVerts * embedDim, &coordsIn);CHKERRQ(ierr); 1966 coordsIn[0*embedDim+0] = -dist; coordsIn[0*embedDim+1] = dist; coordsIn[0*embedDim+2] = -dist; 1967 coordsIn[1*embedDim+0] = dist; coordsIn[1*embedDim+1] = dist; coordsIn[1*embedDim+2] = -dist; 1968 coordsIn[2*embedDim+0] = dist; coordsIn[2*embedDim+1] = -dist; coordsIn[2*embedDim+2] = -dist; 1969 coordsIn[3*embedDim+0] = -dist; coordsIn[3*embedDim+1] = -dist; coordsIn[3*embedDim+2] = -dist; 1970 coordsIn[4*embedDim+0] = -dist; coordsIn[4*embedDim+1] = dist; coordsIn[4*embedDim+2] = dist; 1971 coordsIn[5*embedDim+0] = dist; coordsIn[5*embedDim+1] = dist; coordsIn[5*embedDim+2] = dist; 1972 coordsIn[6*embedDim+0] = -dist; coordsIn[6*embedDim+1] = -dist; coordsIn[6*embedDim+2] = dist; 1973 coordsIn[7*embedDim+0] = dist; coordsIn[7*embedDim+1] = -dist; coordsIn[7*embedDim+2] = dist; 1974 } 1975 break; 1976 case 3: 1977 if (simplex) { 1978 DM idm; 1979 const PetscReal edgeLen = 1.0/PETSC_PHI; 1980 const PetscReal vertexA[4] = {0.5, 0.5, 0.5, 0.5}; 1981 const PetscReal vertexB[4] = {1.0, 0.0, 0.0, 0.0}; 1982 const PetscReal vertexC[4] = {0.5, 0.5*PETSC_PHI, 0.5/PETSC_PHI, 0.0}; 1983 const PetscInt degree = 12; 1984 PetscInt s[4] = {1, 1, 1}; 1985 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}, 1986 {2, 0, 1, 3}, {2, 1, 3, 0}, {2, 3, 0, 1}, {3, 0, 2, 1}, {3, 1, 0, 2}, {3, 2, 1, 0}}; 1987 PetscInt cone[4]; 1988 PetscInt *graph, p, i, j, k, l; 1989 1990 numCells = !rank ? 600 : 0; 1991 numVerts = !rank ? 120 : 0; 1992 firstVertex = numCells; 1993 /* Use the 600-cell, which for a unit sphere has coordinates which are 1994 1995 1/2 (\pm 1, \pm 1, \pm 1, \pm 1) 16 1996 (\pm 1, 0, 0, 0) all cyclic permutations 8 1997 1/2 (\pm 1, \pm phi, \pm 1/phi, 0) all even permutations 96 1998 1999 where \phi^2 - \phi - 1 = 0, meaning \phi is the golden ratio \frac{1 + \sqrt{5}}{2}. The edge 2000 length is then given by 1/\phi = 2.73606. 2001 2002 http://buzzard.pugetsound.edu/sage-practice/ch03s03.html 2003 http://mathworld.wolfram.com/600-Cell.html 2004 */ 2005 /* Construct vertices */ 2006 ierr = PetscCalloc1(numVerts * embedDim, &coordsIn);CHKERRQ(ierr); 2007 i = 0; 2008 for (s[0] = -1; s[0] < 2; s[0] += 2) { 2009 for (s[1] = -1; s[1] < 2; s[1] += 2) { 2010 for (s[2] = -1; s[2] < 2; s[2] += 2) { 2011 for (s[3] = -1; s[3] < 2; s[3] += 2) { 2012 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[d]*vertexA[d]; 2013 ++i; 2014 } 2015 } 2016 } 2017 } 2018 for (p = 0; p < embedDim; ++p) { 2019 s[1] = s[2] = s[3] = 1; 2020 for (s[0] = -1; s[0] < 2; s[0] += 2) { 2021 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[(d+p)%embedDim]*vertexB[(d+p)%embedDim]; 2022 ++i; 2023 } 2024 } 2025 for (p = 0; p < 12; ++p) { 2026 s[3] = 1; 2027 for (s[0] = -1; s[0] < 2; s[0] += 2) { 2028 for (s[1] = -1; s[1] < 2; s[1] += 2) { 2029 for (s[2] = -1; s[2] < 2; s[2] += 2) { 2030 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[evenPerm[p][d]]*vertexC[evenPerm[p][d]]; 2031 ++i; 2032 } 2033 } 2034 } 2035 } 2036 if (i != numVerts) SETERRQ2(comm, PETSC_ERR_PLIB, "Invalid 600-cell, vertices %D != %D", i, numVerts); 2037 /* Construct graph */ 2038 ierr = PetscCalloc1(numVerts * numVerts, &graph);CHKERRQ(ierr); 2039 for (i = 0; i < numVerts; ++i) { 2040 for (j = 0, k = 0; j < numVerts; ++j) { 2041 if (PetscAbsReal(DiffNormReal(embedDim, &coordsIn[i*embedDim], &coordsIn[j*embedDim]) - edgeLen) < PETSC_SMALL) {graph[i*numVerts+j] = 1; ++k;} 2042 } 2043 if (k != degree) SETERRQ3(comm, PETSC_ERR_PLIB, "Invalid 600-cell, vertex %D degree %D != %D", i, k, degree); 2044 } 2045 /* Build Topology */ 2046 ierr = DMPlexSetChart(*dm, 0, numCells+numVerts);CHKERRQ(ierr); 2047 for (c = 0; c < numCells; c++) { 2048 ierr = DMPlexSetConeSize(*dm, c, embedDim);CHKERRQ(ierr); 2049 } 2050 ierr = DMSetUp(*dm);CHKERRQ(ierr); /* Allocate space for cones */ 2051 /* Cells */ 2052 for (i = 0, c = 0; i < numVerts; ++i) { 2053 for (j = 0; j < i; ++j) { 2054 for (k = 0; k < j; ++k) { 2055 for (l = 0; l < k; ++l) { 2056 if (graph[i*numVerts+j] && graph[j*numVerts+k] && graph[k*numVerts+i] && 2057 graph[l*numVerts+i] && graph[l*numVerts+j] && graph[l*numVerts+k]) { 2058 cone[0] = firstVertex+i; cone[1] = firstVertex+j; cone[2] = firstVertex+k; cone[3] = firstVertex+l; 2059 /* Check orientation: https://ef.gy/linear-algebra:normal-vectors-in-higher-dimensional-spaces */ 2060 { 2061 const PetscInt epsilon[4][4][4][4] = {{{{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 2062 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 1}, { 0, 0, -1, 0}}, 2063 {{0, 0, 0, 0}, { 0, 0, 0, -1}, { 0, 0, 0, 0}, { 0, 1, 0, 0}}, 2064 {{0, 0, 0, 0}, { 0, 0, 1, 0}, { 0, -1, 0, 0}, { 0, 0, 0, 0}}}, 2065 2066 {{{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, -1}, { 0, 0, 1, 0}}, 2067 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 2068 {{0, 0, 0, 1}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, {-1, 0, 0, 0}}, 2069 {{0, 0, -1, 0}, { 0, 0, 0, 0}, { 1, 0, 0, 0}, { 0, 0, 0, 0}}}, 2070 2071 {{{0, 0, 0, 0}, { 0, 0, 0, 1}, { 0, 0, 0, 0}, { 0, -1, 0, 0}}, 2072 {{0, 0, 0, -1}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 1, 0, 0, 0}}, 2073 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 2074 {{0, 1, 0, 0}, {-1, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}}, 2075 2076 {{{0, 0, 0, 0}, { 0, 0, -1, 0}, { 0, 1, 0, 0}, { 0, 0, 0, 0}}, 2077 {{0, 0, 1, 0}, { 0, 0, 0, 0}, {-1, 0, 0, 0}, { 0, 0, 0, 0}}, 2078 {{0, -1, 0, 0}, { 1, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 2079 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}}}; 2080 PetscReal normal[4]; 2081 PetscInt e, f, g; 2082 2083 for (d = 0; d < embedDim; ++d) { 2084 normal[d] = 0.0; 2085 for (e = 0; e < embedDim; ++e) { 2086 for (f = 0; f < embedDim; ++f) { 2087 for (g = 0; g < embedDim; ++g) { 2088 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]); 2089 } 2090 } 2091 } 2092 } 2093 if (DotReal(embedDim, normal, &coordsIn[i*embedDim]) < 0) {PetscInt tmp = cone[1]; cone[1] = cone[2]; cone[2] = tmp;} 2094 } 2095 ierr = DMPlexSetCone(*dm, c++, cone);CHKERRQ(ierr); 2096 } 2097 } 2098 } 2099 } 2100 } 2101 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 2102 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 2103 ierr = PetscFree(graph);CHKERRQ(ierr); 2104 /* Interpolate mesh */ 2105 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 2106 ierr = DMDestroy(dm);CHKERRQ(ierr); 2107 *dm = idm; 2108 break; 2109 } 2110 default: SETERRQ1(comm, PETSC_ERR_SUP, "Unsupported dimension for sphere: %D", dim); 2111 } 2112 /* Create coordinates */ 2113 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 2114 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 2115 ierr = PetscSectionSetFieldComponents(coordSection, 0, embedDim);CHKERRQ(ierr); 2116 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVerts);CHKERRQ(ierr); 2117 for (v = firstVertex; v < firstVertex+numVerts; ++v) { 2118 ierr = PetscSectionSetDof(coordSection, v, embedDim);CHKERRQ(ierr); 2119 ierr = PetscSectionSetFieldDof(coordSection, v, 0, embedDim);CHKERRQ(ierr); 2120 } 2121 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 2122 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 2123 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 2124 ierr = VecSetBlockSize(coordinates, embedDim);CHKERRQ(ierr); 2125 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 2126 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 2127 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 2128 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2129 for (v = 0; v < numVerts; ++v) for (d = 0; d < embedDim; ++d) {coords[v*embedDim+d] = coordsIn[v*embedDim+d];} 2130 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2131 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 2132 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 2133 ierr = PetscFree(coordsIn);CHKERRQ(ierr); 2134 PetscFunctionReturn(0); 2135 } 2136 2137 /* External function declarations here */ 2138 extern PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling); 2139 extern PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, Mat *mat); 2140 extern PetscErrorCode DMCreateMassMatrix_Plex(DM dmCoarse, DM dmFine, Mat *mat); 2141 extern PetscErrorCode DMCreateDefaultSection_Plex(DM dm); 2142 extern PetscErrorCode DMCreateDefaultConstraints_Plex(DM dm); 2143 extern PetscErrorCode DMCreateMatrix_Plex(DM dm, Mat *J); 2144 extern PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm); 2145 extern PetscErrorCode DMCreateCoordinateField_Plex(DM dm, DMField *field); 2146 PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm); 2147 extern PetscErrorCode DMSetUp_Plex(DM dm); 2148 extern PetscErrorCode DMDestroy_Plex(DM dm); 2149 extern PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer); 2150 extern PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer); 2151 extern PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm); 2152 extern PetscErrorCode DMCreateSuperDM_Plex(DM dms[], PetscInt len, IS **is, DM *superdm); 2153 static PetscErrorCode DMInitialize_Plex(DM dm); 2154 2155 /* Replace dm with the contents of dmNew 2156 - Share the DM_Plex structure 2157 - Share the coordinates 2158 - Share the SF 2159 */ 2160 static PetscErrorCode DMPlexReplace_Static(DM dm, DM dmNew) 2161 { 2162 PetscSF sf; 2163 DM coordDM, coarseDM; 2164 Vec coords; 2165 PetscBool isper; 2166 const PetscReal *maxCell, *L; 2167 const DMBoundaryType *bd; 2168 PetscErrorCode ierr; 2169 2170 PetscFunctionBegin; 2171 ierr = DMGetPointSF(dmNew, &sf);CHKERRQ(ierr); 2172 ierr = DMSetPointSF(dm, sf);CHKERRQ(ierr); 2173 ierr = DMGetCoordinateDM(dmNew, &coordDM);CHKERRQ(ierr); 2174 ierr = DMGetCoordinatesLocal(dmNew, &coords);CHKERRQ(ierr); 2175 ierr = DMSetCoordinateDM(dm, coordDM);CHKERRQ(ierr); 2176 ierr = DMSetCoordinatesLocal(dm, coords);CHKERRQ(ierr); 2177 ierr = DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);CHKERRQ(ierr); 2178 ierr = DMSetPeriodicity(dmNew, isper, maxCell, L, bd);CHKERRQ(ierr); 2179 ierr = DMDestroy_Plex(dm);CHKERRQ(ierr); 2180 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 2181 dm->data = dmNew->data; 2182 ((DM_Plex *) dmNew->data)->refct++; 2183 dmNew->labels->refct++; 2184 if (!--(dm->labels->refct)) { 2185 DMLabelLink next = dm->labels->next; 2186 2187 /* destroy the labels */ 2188 while (next) { 2189 DMLabelLink tmp = next->next; 2190 2191 ierr = DMLabelDestroy(&next->label);CHKERRQ(ierr); 2192 ierr = PetscFree(next);CHKERRQ(ierr); 2193 next = tmp; 2194 } 2195 ierr = PetscFree(dm->labels);CHKERRQ(ierr); 2196 } 2197 dm->labels = dmNew->labels; 2198 dm->depthLabel = dmNew->depthLabel; 2199 ierr = DMGetCoarseDM(dmNew,&coarseDM);CHKERRQ(ierr); 2200 ierr = DMSetCoarseDM(dm,coarseDM);CHKERRQ(ierr); 2201 PetscFunctionReturn(0); 2202 } 2203 2204 /* Swap dm with the contents of dmNew 2205 - Swap the DM_Plex structure 2206 - Swap the coordinates 2207 - Swap the point PetscSF 2208 */ 2209 static PetscErrorCode DMPlexSwap_Static(DM dmA, DM dmB) 2210 { 2211 DM coordDMA, coordDMB; 2212 Vec coordsA, coordsB; 2213 PetscSF sfA, sfB; 2214 void *tmp; 2215 DMLabelLinkList listTmp; 2216 DMLabel depthTmp; 2217 PetscInt tmpI; 2218 PetscErrorCode ierr; 2219 2220 PetscFunctionBegin; 2221 ierr = DMGetPointSF(dmA, &sfA);CHKERRQ(ierr); 2222 ierr = DMGetPointSF(dmB, &sfB);CHKERRQ(ierr); 2223 ierr = PetscObjectReference((PetscObject) sfA);CHKERRQ(ierr); 2224 ierr = DMSetPointSF(dmA, sfB);CHKERRQ(ierr); 2225 ierr = DMSetPointSF(dmB, sfA);CHKERRQ(ierr); 2226 ierr = PetscObjectDereference((PetscObject) sfA);CHKERRQ(ierr); 2227 2228 ierr = DMGetCoordinateDM(dmA, &coordDMA);CHKERRQ(ierr); 2229 ierr = DMGetCoordinateDM(dmB, &coordDMB);CHKERRQ(ierr); 2230 ierr = PetscObjectReference((PetscObject) coordDMA);CHKERRQ(ierr); 2231 ierr = DMSetCoordinateDM(dmA, coordDMB);CHKERRQ(ierr); 2232 ierr = DMSetCoordinateDM(dmB, coordDMA);CHKERRQ(ierr); 2233 ierr = PetscObjectDereference((PetscObject) coordDMA);CHKERRQ(ierr); 2234 2235 ierr = DMGetCoordinatesLocal(dmA, &coordsA);CHKERRQ(ierr); 2236 ierr = DMGetCoordinatesLocal(dmB, &coordsB);CHKERRQ(ierr); 2237 ierr = PetscObjectReference((PetscObject) coordsA);CHKERRQ(ierr); 2238 ierr = DMSetCoordinatesLocal(dmA, coordsB);CHKERRQ(ierr); 2239 ierr = DMSetCoordinatesLocal(dmB, coordsA);CHKERRQ(ierr); 2240 ierr = PetscObjectDereference((PetscObject) coordsA);CHKERRQ(ierr); 2241 2242 tmp = dmA->data; 2243 dmA->data = dmB->data; 2244 dmB->data = tmp; 2245 listTmp = dmA->labels; 2246 dmA->labels = dmB->labels; 2247 dmB->labels = listTmp; 2248 depthTmp = dmA->depthLabel; 2249 dmA->depthLabel = dmB->depthLabel; 2250 dmB->depthLabel = depthTmp; 2251 tmpI = dmA->levelup; 2252 dmA->levelup = dmB->levelup; 2253 dmB->levelup = tmpI; 2254 PetscFunctionReturn(0); 2255 } 2256 2257 static PetscErrorCode DMPlexIsSimplex_Static(DM dm, PetscBool *isSimplex) 2258 { 2259 PetscInt dim, cStart, coneSize; 2260 PetscErrorCode ierr; 2261 2262 PetscFunctionBegin; 2263 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2264 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, NULL);CHKERRQ(ierr); 2265 ierr = DMPlexGetConeSize(dm, cStart, &coneSize);CHKERRQ(ierr); 2266 *isSimplex = coneSize == dim+1 ? PETSC_TRUE : PETSC_FALSE; 2267 PetscFunctionReturn(0); 2268 } 2269 2270 PetscErrorCode DMSetFromOptions_NonRefinement_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 2271 { 2272 DM_Plex *mesh = (DM_Plex*) dm->data; 2273 PetscErrorCode ierr; 2274 2275 PetscFunctionBegin; 2276 /* Handle viewing */ 2277 ierr = PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMPlexMatSetClosure", PETSC_FALSE, &mesh->printSetValues, NULL);CHKERRQ(ierr); 2278 ierr = PetscOptionsInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMPlexSNESComputeResidualFEM", 0, &mesh->printFEM, NULL);CHKERRQ(ierr); 2279 ierr = PetscOptionsReal("-dm_plex_print_tol", "Tolerance for FEM output", "DMPlexSNESComputeResidualFEM", mesh->printTol, &mesh->printTol, NULL);CHKERRQ(ierr); 2280 ierr = PetscOptionsInt("-dm_plex_print_l2", "Debug output level all L2 diff computations", "DMComputeL2Diff", 0, &mesh->printL2, NULL);CHKERRQ(ierr); 2281 /* Point Location */ 2282 ierr = PetscOptionsBool("-dm_plex_hash_location", "Use grid hashing for point location", "DMInterpolate", PETSC_FALSE, &mesh->useHashLocation, NULL);CHKERRQ(ierr); 2283 /* Partitioning and distribution */ 2284 ierr = PetscOptionsBool("-dm_plex_partition_balance", "Attempt to evenly divide points on partition boundary between processes", "DMPlexSetPartitionBalance", PETSC_FALSE, &mesh->partitionBalance, NULL);CHKERRQ(ierr); 2285 /* Generation and remeshing */ 2286 ierr = PetscOptionsBool("-dm_plex_remesh_bd", "Allow changes to the boundary on remeshing", "DMAdapt", PETSC_FALSE, &mesh->remeshBd, NULL);CHKERRQ(ierr); 2287 /* Projection behavior */ 2288 ierr = PetscOptionsInt("-dm_plex_max_projection_height", "Maxmimum mesh point height used to project locally", "DMPlexSetMaxProjectionHeight", 0, &mesh->maxProjectionHeight, NULL);CHKERRQ(ierr); 2289 ierr = PetscOptionsBool("-dm_plex_regular_refinement", "Use special nested projection algorithm for regular refinement", "DMPlexSetRegularRefinement", mesh->regularRefinement, &mesh->regularRefinement, NULL);CHKERRQ(ierr); 2290 /* Checking structure */ 2291 { 2292 const char *cellTypes[] = {"simplex", "tensor", "unknown", "DMPlexCellType", "DM_PLEX_CELLTYPE_", NULL}; 2293 PetscEnum ct; 2294 PetscBool flg = PETSC_FALSE, flg2 = PETSC_FALSE; 2295 2296 ierr = PetscOptionsBool("-dm_plex_check_symmetry", "Check that the adjacency information in the mesh is symmetric", "DMPlexCheckSymmetry", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr); 2297 if (flg && flg2) {ierr = DMPlexCheckSymmetry(dm);CHKERRQ(ierr);} 2298 ierr = PetscOptionsEnum("-dm_plex_check_skeleton", "Check that each cell has the correct number of vertices", "DMPlexCheckSkeleton", cellTypes, (PetscEnum) 0, &ct, &flg);CHKERRQ(ierr); 2299 if (flg) { 2300 PetscBool isSimplex = ct ? PETSC_FALSE : PETSC_TRUE; 2301 2302 if (ct == (PetscEnum) DM_PLEX_CELLTYPE_UNKNOWN) {ierr = DMPlexIsSimplex_Static(dm, &isSimplex);CHKERRQ(ierr);} 2303 ierr = DMPlexCheckSkeleton(dm, isSimplex, 0);CHKERRQ(ierr); 2304 } 2305 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); 2306 if (flg) { 2307 PetscBool isSimplex = ct ? PETSC_FALSE : PETSC_TRUE; 2308 2309 if (ct == (PetscEnum) DM_PLEX_CELLTYPE_UNKNOWN) {ierr = DMPlexIsSimplex_Static(dm, &isSimplex);CHKERRQ(ierr);} 2310 ierr = DMPlexCheckFaces(dm, isSimplex, 0);CHKERRQ(ierr); 2311 } 2312 ierr = PetscOptionsBool("-dm_plex_check_geometry", "Check that cells have positive volume", "DMPlexCheckGeometry", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr); 2313 if (flg && flg2) {ierr = DMPlexCheckGeometry(dm);CHKERRQ(ierr);} 2314 } 2315 2316 ierr = PetscPartitionerSetFromOptions(mesh->partitioner);CHKERRQ(ierr); 2317 PetscFunctionReturn(0); 2318 } 2319 2320 static PetscErrorCode DMSetFromOptions_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 2321 { 2322 PetscInt refine = 0, coarsen = 0, r; 2323 PetscBool isHierarchy; 2324 PetscErrorCode ierr; 2325 2326 PetscFunctionBegin; 2327 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2328 ierr = PetscOptionsHead(PetscOptionsObject,"DMPlex Options");CHKERRQ(ierr); 2329 /* Handle DMPlex refinement */ 2330 ierr = PetscOptionsInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL);CHKERRQ(ierr); 2331 ierr = PetscOptionsInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy);CHKERRQ(ierr); 2332 if (refine) {ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE);CHKERRQ(ierr);} 2333 if (refine && isHierarchy) { 2334 DM *dms, coarseDM; 2335 2336 ierr = DMGetCoarseDM(dm, &coarseDM);CHKERRQ(ierr); 2337 ierr = PetscObjectReference((PetscObject)coarseDM);CHKERRQ(ierr); 2338 ierr = PetscMalloc1(refine,&dms);CHKERRQ(ierr); 2339 ierr = DMRefineHierarchy(dm, refine, dms);CHKERRQ(ierr); 2340 /* Total hack since we do not pass in a pointer */ 2341 ierr = DMPlexSwap_Static(dm, dms[refine-1]);CHKERRQ(ierr); 2342 if (refine == 1) { 2343 ierr = DMSetCoarseDM(dm, dms[0]);CHKERRQ(ierr); 2344 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr); 2345 } else { 2346 ierr = DMSetCoarseDM(dm, dms[refine-2]);CHKERRQ(ierr); 2347 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr); 2348 ierr = DMSetCoarseDM(dms[0], dms[refine-1]);CHKERRQ(ierr); 2349 ierr = DMPlexSetRegularRefinement(dms[0], PETSC_TRUE);CHKERRQ(ierr); 2350 } 2351 ierr = DMSetCoarseDM(dms[refine-1], coarseDM);CHKERRQ(ierr); 2352 ierr = PetscObjectDereference((PetscObject)coarseDM);CHKERRQ(ierr); 2353 /* Free DMs */ 2354 for (r = 0; r < refine; ++r) { 2355 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr); 2356 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 2357 } 2358 ierr = PetscFree(dms);CHKERRQ(ierr); 2359 } else { 2360 for (r = 0; r < refine; ++r) { 2361 DM refinedMesh; 2362 2363 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2364 ierr = DMRefine(dm, PetscObjectComm((PetscObject) dm), &refinedMesh);CHKERRQ(ierr); 2365 /* Total hack since we do not pass in a pointer */ 2366 ierr = DMPlexReplace_Static(dm, refinedMesh);CHKERRQ(ierr); 2367 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2368 ierr = DMDestroy(&refinedMesh);CHKERRQ(ierr); 2369 } 2370 } 2371 /* Handle DMPlex coarsening */ 2372 ierr = PetscOptionsInt("-dm_coarsen", "Coarsen the mesh", "DMCreate", coarsen, &coarsen, NULL);CHKERRQ(ierr); 2373 ierr = PetscOptionsInt("-dm_coarsen_hierarchy", "The number of coarsenings", "DMCreate", coarsen, &coarsen, &isHierarchy);CHKERRQ(ierr); 2374 if (coarsen && isHierarchy) { 2375 DM *dms; 2376 2377 ierr = PetscMalloc1(coarsen, &dms);CHKERRQ(ierr); 2378 ierr = DMCoarsenHierarchy(dm, coarsen, dms);CHKERRQ(ierr); 2379 /* Free DMs */ 2380 for (r = 0; r < coarsen; ++r) { 2381 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr); 2382 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 2383 } 2384 ierr = PetscFree(dms);CHKERRQ(ierr); 2385 } else { 2386 for (r = 0; r < coarsen; ++r) { 2387 DM coarseMesh; 2388 2389 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2390 ierr = DMCoarsen(dm, PetscObjectComm((PetscObject) dm), &coarseMesh);CHKERRQ(ierr); 2391 /* Total hack since we do not pass in a pointer */ 2392 ierr = DMPlexReplace_Static(dm, coarseMesh);CHKERRQ(ierr); 2393 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2394 ierr = DMDestroy(&coarseMesh);CHKERRQ(ierr); 2395 } 2396 } 2397 /* Handle */ 2398 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2399 ierr = PetscOptionsTail();CHKERRQ(ierr); 2400 PetscFunctionReturn(0); 2401 } 2402 2403 static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec) 2404 { 2405 PetscErrorCode ierr; 2406 2407 PetscFunctionBegin; 2408 ierr = DMCreateGlobalVector_Section_Private(dm,vec);CHKERRQ(ierr); 2409 /* ierr = VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM);CHKERRQ(ierr); */ 2410 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex);CHKERRQ(ierr); 2411 ierr = VecSetOperation(*vec, VECOP_VIEWNATIVE, (void (*)(void)) VecView_Plex_Native);CHKERRQ(ierr); 2412 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex);CHKERRQ(ierr); 2413 ierr = VecSetOperation(*vec, VECOP_LOADNATIVE, (void (*)(void)) VecLoad_Plex_Native);CHKERRQ(ierr); 2414 PetscFunctionReturn(0); 2415 } 2416 2417 static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec) 2418 { 2419 PetscErrorCode ierr; 2420 2421 PetscFunctionBegin; 2422 ierr = DMCreateLocalVector_Section_Private(dm,vec);CHKERRQ(ierr); 2423 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex_Local);CHKERRQ(ierr); 2424 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex_Local);CHKERRQ(ierr); 2425 PetscFunctionReturn(0); 2426 } 2427 2428 static PetscErrorCode DMGetDimPoints_Plex(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd) 2429 { 2430 PetscInt depth, d; 2431 PetscErrorCode ierr; 2432 2433 PetscFunctionBegin; 2434 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 2435 if (depth == 1) { 2436 ierr = DMGetDimension(dm, &d);CHKERRQ(ierr); 2437 if (dim == 0) {ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr);} 2438 else if (dim == d) {ierr = DMPlexGetDepthStratum(dm, 1, pStart, pEnd);CHKERRQ(ierr);} 2439 else {*pStart = 0; *pEnd = 0;} 2440 } else { 2441 ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr); 2442 } 2443 PetscFunctionReturn(0); 2444 } 2445 2446 static PetscErrorCode DMGetNeighors_Plex(DM dm, PetscInt *nranks, const PetscMPIInt *ranks[]) 2447 { 2448 PetscSF sf; 2449 PetscErrorCode ierr; 2450 2451 PetscFunctionBegin; 2452 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 2453 ierr = PetscSFGetRanks(sf, nranks, ranks, NULL, NULL, NULL);CHKERRQ(ierr); 2454 PetscFunctionReturn(0); 2455 } 2456 2457 static PetscErrorCode DMHasCreateInjection_Plex(DM dm, PetscBool *flg) 2458 { 2459 PetscFunctionBegin; 2460 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2461 PetscValidPointer(flg,2); 2462 *flg = PETSC_TRUE; 2463 PetscFunctionReturn(0); 2464 } 2465 2466 static PetscErrorCode DMInitialize_Plex(DM dm) 2467 { 2468 PetscErrorCode ierr; 2469 2470 PetscFunctionBegin; 2471 dm->ops->view = DMView_Plex; 2472 dm->ops->load = DMLoad_Plex; 2473 dm->ops->setfromoptions = DMSetFromOptions_Plex; 2474 dm->ops->clone = DMClone_Plex; 2475 dm->ops->setup = DMSetUp_Plex; 2476 dm->ops->createdefaultsection = DMCreateDefaultSection_Plex; 2477 dm->ops->createdefaultconstraints = DMCreateDefaultConstraints_Plex; 2478 dm->ops->createglobalvector = DMCreateGlobalVector_Plex; 2479 dm->ops->createlocalvector = DMCreateLocalVector_Plex; 2480 dm->ops->getlocaltoglobalmapping = NULL; 2481 dm->ops->createfieldis = NULL; 2482 dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex; 2483 dm->ops->createcoordinatefield = DMCreateCoordinateField_Plex; 2484 dm->ops->getcoloring = NULL; 2485 dm->ops->creatematrix = DMCreateMatrix_Plex; 2486 dm->ops->createinterpolation = DMCreateInterpolation_Plex; 2487 dm->ops->createmassmatrix = DMCreateMassMatrix_Plex; 2488 dm->ops->getaggregates = NULL; 2489 dm->ops->getinjection = DMCreateInjection_Plex; 2490 dm->ops->hascreateinjection = DMHasCreateInjection_Plex; 2491 dm->ops->refine = DMRefine_Plex; 2492 dm->ops->coarsen = DMCoarsen_Plex; 2493 dm->ops->refinehierarchy = DMRefineHierarchy_Plex; 2494 dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Plex; 2495 dm->ops->adaptlabel = DMAdaptLabel_Plex; 2496 dm->ops->adaptmetric = DMAdaptMetric_Plex; 2497 dm->ops->globaltolocalbegin = NULL; 2498 dm->ops->globaltolocalend = NULL; 2499 dm->ops->localtoglobalbegin = NULL; 2500 dm->ops->localtoglobalend = NULL; 2501 dm->ops->destroy = DMDestroy_Plex; 2502 dm->ops->createsubdm = DMCreateSubDM_Plex; 2503 dm->ops->createsuperdm = DMCreateSuperDM_Plex; 2504 dm->ops->getdimpoints = DMGetDimPoints_Plex; 2505 dm->ops->locatepoints = DMLocatePoints_Plex; 2506 dm->ops->projectfunctionlocal = DMProjectFunctionLocal_Plex; 2507 dm->ops->projectfunctionlabellocal = DMProjectFunctionLabelLocal_Plex; 2508 dm->ops->projectfieldlocal = DMProjectFieldLocal_Plex; 2509 dm->ops->projectfieldlabellocal = DMProjectFieldLabelLocal_Plex; 2510 dm->ops->computel2diff = DMComputeL2Diff_Plex; 2511 dm->ops->computel2gradientdiff = DMComputeL2GradientDiff_Plex; 2512 dm->ops->computel2fielddiff = DMComputeL2FieldDiff_Plex; 2513 dm->ops->getneighbors = DMGetNeighors_Plex; 2514 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMPlexInsertBoundaryValues_C",DMPlexInsertBoundaryValues_Plex);CHKERRQ(ierr); 2515 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMSetUpGLVisViewer_C",DMSetUpGLVisViewer_Plex);CHKERRQ(ierr); 2516 PetscFunctionReturn(0); 2517 } 2518 2519 PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm) 2520 { 2521 DM_Plex *mesh = (DM_Plex *) dm->data; 2522 PetscErrorCode ierr; 2523 2524 PetscFunctionBegin; 2525 mesh->refct++; 2526 (*newdm)->data = mesh; 2527 ierr = PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);CHKERRQ(ierr); 2528 ierr = DMInitialize_Plex(*newdm);CHKERRQ(ierr); 2529 PetscFunctionReturn(0); 2530 } 2531 2532 /*MC 2533 DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram. 2534 In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is 2535 specified by a PetscSection object. Ownership in the global representation is determined by 2536 ownership of the underlying DMPlex points. This is specified by another PetscSection object. 2537 2538 Options Database Keys: 2539 + -dm_plex_hash_location - Use grid hashing for point location 2540 . -dm_plex_partition_balance - Attempt to evenly divide points on partition boundary between processes 2541 . -dm_plex_remesh_bd - Allow changes to the boundary on remeshing 2542 . -dm_plex_max_projection_height - Maxmimum mesh point height used to project locally 2543 . -dm_plex_regular_refinement - Use special nested projection algorithm for regular refinement 2544 . -dm_plex_check_symmetry - Check that the adjacency information in the mesh is symmetric 2545 . -dm_plex_check_skeleton <celltype> - Check that each cell has the correct number of vertices 2546 . -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 2547 . -dm_plex_check_geometry - Check that cells have positive volume 2548 . -dm_view :mesh.tex:ascii_latex - View the mesh in LaTeX/TikZ 2549 . -dm_plex_view_scale <num> - Scale the TikZ 2550 - -dm_plex_print_fem <num> - View FEM assembly information, such as element vectors and matrices 2551 2552 2553 Level: intermediate 2554 2555 .seealso: DMType, DMPlexCreate(), DMCreate(), DMSetType() 2556 M*/ 2557 2558 PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm) 2559 { 2560 DM_Plex *mesh; 2561 PetscInt unit, d; 2562 PetscErrorCode ierr; 2563 2564 PetscFunctionBegin; 2565 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2566 ierr = PetscNewLog(dm,&mesh);CHKERRQ(ierr); 2567 dm->dim = 0; 2568 dm->data = mesh; 2569 2570 mesh->refct = 1; 2571 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);CHKERRQ(ierr); 2572 mesh->maxConeSize = 0; 2573 mesh->cones = NULL; 2574 mesh->coneOrientations = NULL; 2575 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);CHKERRQ(ierr); 2576 mesh->maxSupportSize = 0; 2577 mesh->supports = NULL; 2578 mesh->refinementUniform = PETSC_TRUE; 2579 mesh->refinementLimit = -1.0; 2580 2581 mesh->facesTmp = NULL; 2582 2583 mesh->tetgenOpts = NULL; 2584 mesh->triangleOpts = NULL; 2585 ierr = PetscPartitionerCreate(PetscObjectComm((PetscObject)dm), &mesh->partitioner);CHKERRQ(ierr); 2586 mesh->remeshBd = PETSC_FALSE; 2587 2588 mesh->subpointMap = NULL; 2589 2590 for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0; 2591 2592 mesh->regularRefinement = PETSC_FALSE; 2593 mesh->depthState = -1; 2594 mesh->globalVertexNumbers = NULL; 2595 mesh->globalCellNumbers = NULL; 2596 mesh->anchorSection = NULL; 2597 mesh->anchorIS = NULL; 2598 mesh->createanchors = NULL; 2599 mesh->computeanchormatrix = NULL; 2600 mesh->parentSection = NULL; 2601 mesh->parents = NULL; 2602 mesh->childIDs = NULL; 2603 mesh->childSection = NULL; 2604 mesh->children = NULL; 2605 mesh->referenceTree = NULL; 2606 mesh->getchildsymmetry = NULL; 2607 for (d = 0; d < 8; ++d) mesh->hybridPointMax[d] = PETSC_DETERMINE; 2608 mesh->vtkCellHeight = 0; 2609 mesh->useAnchors = PETSC_FALSE; 2610 2611 mesh->maxProjectionHeight = 0; 2612 2613 mesh->printSetValues = PETSC_FALSE; 2614 mesh->printFEM = 0; 2615 mesh->printTol = 1.0e-10; 2616 2617 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 2618 PetscFunctionReturn(0); 2619 } 2620 2621 /*@ 2622 DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram. 2623 2624 Collective on MPI_Comm 2625 2626 Input Parameter: 2627 . comm - The communicator for the DMPlex object 2628 2629 Output Parameter: 2630 . mesh - The DMPlex object 2631 2632 Level: beginner 2633 2634 .keywords: DMPlex, create 2635 @*/ 2636 PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh) 2637 { 2638 PetscErrorCode ierr; 2639 2640 PetscFunctionBegin; 2641 PetscValidPointer(mesh,2); 2642 ierr = DMCreate(comm, mesh);CHKERRQ(ierr); 2643 ierr = DMSetType(*mesh, DMPLEX);CHKERRQ(ierr); 2644 PetscFunctionReturn(0); 2645 } 2646 2647 /* 2648 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 2649 */ 2650 /* TODO: invertCells and spaceDim arguments could be added also to to DMPlexCreateFromCellListParallel(), DMPlexBuildFromCellList_Internal() and DMPlexCreateFromCellList() */ 2651 PetscErrorCode DMPlexBuildFromCellList_Parallel_Internal(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[], PetscBool invertCells, PetscSF *sfVert) 2652 { 2653 PetscSF sfPoint; 2654 PetscLayout vLayout; 2655 PetscHSetI vhash; 2656 PetscSFNode *remoteVerticesAdj, *vertexLocal, *vertexOwner, *remoteVertex; 2657 const PetscInt *vrange; 2658 PetscInt numVerticesAdj, off = 0, *verticesAdj, numVerticesGhost = 0, *localVertex, *cone, c, p, v, g; 2659 PetscMPIInt rank, size; 2660 PetscErrorCode ierr; 2661 2662 PetscFunctionBegin; 2663 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr); 2664 ierr = MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);CHKERRQ(ierr); 2665 /* Partition vertices */ 2666 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) dm), &vLayout);CHKERRQ(ierr); 2667 ierr = PetscLayoutSetLocalSize(vLayout, numVertices);CHKERRQ(ierr); 2668 ierr = PetscLayoutSetBlockSize(vLayout, 1);CHKERRQ(ierr); 2669 ierr = PetscLayoutSetUp(vLayout);CHKERRQ(ierr); 2670 ierr = PetscLayoutGetRanges(vLayout, &vrange);CHKERRQ(ierr); 2671 /* Count vertices and map them to procs */ 2672 ierr = PetscHSetICreate(&vhash);CHKERRQ(ierr); 2673 for (c = 0; c < numCells; ++c) { 2674 for (p = 0; p < numCorners; ++p) { 2675 ierr = PetscHSetIAdd(vhash, cells[c*numCorners+p]);CHKERRQ(ierr); 2676 } 2677 } 2678 ierr = PetscHSetIGetSize(vhash, &numVerticesAdj);CHKERRQ(ierr); 2679 ierr = PetscMalloc1(numVerticesAdj, &verticesAdj);CHKERRQ(ierr); 2680 ierr = PetscHSetIGetElems(vhash, &off, verticesAdj);CHKERRQ(ierr); 2681 ierr = PetscHSetIDestroy(&vhash);CHKERRQ(ierr); 2682 if (off != numVerticesAdj) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid number of local vertices %D should be %D", off, numVerticesAdj); 2683 ierr = PetscSortInt(numVerticesAdj, verticesAdj);CHKERRQ(ierr); 2684 ierr = PetscMalloc1(numVerticesAdj, &remoteVerticesAdj);CHKERRQ(ierr); 2685 for (v = 0; v < numVerticesAdj; ++v) { 2686 const PetscInt gv = verticesAdj[v]; 2687 PetscInt vrank; 2688 2689 ierr = PetscFindInt(gv, size+1, vrange, &vrank);CHKERRQ(ierr); 2690 vrank = vrank < 0 ? -(vrank+2) : vrank; 2691 remoteVerticesAdj[v].index = gv - vrange[vrank]; 2692 remoteVerticesAdj[v].rank = vrank; 2693 } 2694 /* Create cones */ 2695 ierr = DMPlexSetChart(dm, 0, numCells+numVerticesAdj);CHKERRQ(ierr); 2696 for (c = 0; c < numCells; ++c) {ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr);} 2697 ierr = DMSetUp(dm);CHKERRQ(ierr); 2698 ierr = DMGetWorkArray(dm, numCorners, MPIU_INT, &cone);CHKERRQ(ierr); 2699 for (c = 0; c < numCells; ++c) { 2700 for (p = 0; p < numCorners; ++p) { 2701 const PetscInt gv = cells[c*numCorners+p]; 2702 PetscInt lv; 2703 2704 ierr = PetscFindInt(gv, numVerticesAdj, verticesAdj, &lv);CHKERRQ(ierr); 2705 if (lv < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not find global vertex %D in local connectivity", gv); 2706 cone[p] = lv+numCells; 2707 } 2708 if (invertCells) { ierr = DMPlexInvertCell_Internal(spaceDim, numCorners, cone);CHKERRQ(ierr); } 2709 ierr = DMPlexSetCone(dm, c, cone);CHKERRQ(ierr); 2710 } 2711 ierr = DMRestoreWorkArray(dm, numCorners, MPIU_INT, &cone);CHKERRQ(ierr); 2712 /* Create SF for vertices */ 2713 ierr = PetscSFCreate(PetscObjectComm((PetscObject)dm), sfVert);CHKERRQ(ierr); 2714 ierr = PetscObjectSetName((PetscObject) *sfVert, "Vertex Ownership SF");CHKERRQ(ierr); 2715 ierr = PetscSFSetFromOptions(*sfVert);CHKERRQ(ierr); 2716 ierr = PetscSFSetGraph(*sfVert, numVertices, numVerticesAdj, NULL, PETSC_OWN_POINTER, remoteVerticesAdj, PETSC_OWN_POINTER);CHKERRQ(ierr); 2717 ierr = PetscFree(verticesAdj);CHKERRQ(ierr); 2718 /* Build pointSF */ 2719 ierr = PetscMalloc2(numVerticesAdj, &vertexLocal, numVertices, &vertexOwner);CHKERRQ(ierr); 2720 for (v = 0; v < numVerticesAdj; ++v) {vertexLocal[v].index = v+numCells; vertexLocal[v].rank = rank;} 2721 for (v = 0; v < numVertices; ++v) {vertexOwner[v].index = -1; vertexOwner[v].rank = -1;} 2722 ierr = PetscSFReduceBegin(*sfVert, MPIU_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);CHKERRQ(ierr); 2723 ierr = PetscSFReduceEnd(*sfVert, MPIU_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);CHKERRQ(ierr); 2724 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); 2725 ierr = PetscSFBcastBegin(*sfVert, MPIU_2INT, vertexOwner, vertexLocal);CHKERRQ(ierr); 2726 ierr = PetscSFBcastEnd(*sfVert, MPIU_2INT, vertexOwner, vertexLocal);CHKERRQ(ierr); 2727 for (v = 0; v < numVerticesAdj; ++v) if (vertexLocal[v].rank != rank) ++numVerticesGhost; 2728 ierr = PetscMalloc1(numVerticesGhost, &localVertex);CHKERRQ(ierr); 2729 ierr = PetscMalloc1(numVerticesGhost, &remoteVertex);CHKERRQ(ierr); 2730 for (v = 0, g = 0; v < numVerticesAdj; ++v) { 2731 if (vertexLocal[v].rank != rank) { 2732 localVertex[g] = v+numCells; 2733 remoteVertex[g].index = vertexLocal[v].index; 2734 remoteVertex[g].rank = vertexLocal[v].rank; 2735 ++g; 2736 } 2737 } 2738 ierr = PetscFree2(vertexLocal, vertexOwner);CHKERRQ(ierr); 2739 if (g != numVerticesGhost) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of vertex ghosts %D should be %D", g, numVerticesGhost); 2740 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 2741 ierr = PetscObjectSetName((PetscObject) sfPoint, "point SF");CHKERRQ(ierr); 2742 ierr = PetscSFSetGraph(sfPoint, numCells+numVerticesAdj, numVerticesGhost, localVertex, PETSC_OWN_POINTER, remoteVertex, PETSC_OWN_POINTER);CHKERRQ(ierr); 2743 ierr = PetscLayoutDestroy(&vLayout);CHKERRQ(ierr); 2744 /* Fill in the rest of the topology structure */ 2745 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 2746 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 2747 PetscFunctionReturn(0); 2748 } 2749 2750 /* 2751 This takes as input the coordinates for each owned vertex 2752 */ 2753 PetscErrorCode DMPlexBuildCoordinates_Parallel_Internal(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numV, PetscSF sfVert, const PetscReal vertexCoords[]) 2754 { 2755 PetscSection coordSection; 2756 Vec coordinates; 2757 PetscScalar *coords; 2758 PetscInt numVertices, numVerticesAdj, coordSize, v; 2759 PetscErrorCode ierr; 2760 2761 PetscFunctionBegin; 2762 ierr = DMSetCoordinateDim(dm, spaceDim);CHKERRQ(ierr); 2763 ierr = PetscSFGetGraph(sfVert, &numVertices, &numVerticesAdj, NULL, NULL);CHKERRQ(ierr); 2764 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2765 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 2766 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 2767 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVerticesAdj);CHKERRQ(ierr); 2768 for (v = numCells; v < numCells+numVerticesAdj; ++v) { 2769 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 2770 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 2771 } 2772 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 2773 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 2774 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 2775 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 2776 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 2777 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 2778 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 2779 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2780 { 2781 MPI_Datatype coordtype; 2782 2783 /* Need a temp buffer for coords if we have complex/single */ 2784 ierr = MPI_Type_contiguous(spaceDim, MPIU_SCALAR, &coordtype);CHKERRQ(ierr); 2785 ierr = MPI_Type_commit(&coordtype);CHKERRQ(ierr); 2786 #if defined(PETSC_USE_COMPLEX) 2787 { 2788 PetscScalar *svertexCoords; 2789 PetscInt i; 2790 ierr = PetscMalloc1(numV*spaceDim,&svertexCoords);CHKERRQ(ierr); 2791 for (i=0; i<numV*spaceDim; i++) svertexCoords[i] = vertexCoords[i]; 2792 ierr = PetscSFBcastBegin(sfVert, coordtype, svertexCoords, coords);CHKERRQ(ierr); 2793 ierr = PetscSFBcastEnd(sfVert, coordtype, svertexCoords, coords);CHKERRQ(ierr); 2794 ierr = PetscFree(svertexCoords);CHKERRQ(ierr); 2795 } 2796 #else 2797 ierr = PetscSFBcastBegin(sfVert, coordtype, vertexCoords, coords);CHKERRQ(ierr); 2798 ierr = PetscSFBcastEnd(sfVert, coordtype, vertexCoords, coords);CHKERRQ(ierr); 2799 #endif 2800 ierr = MPI_Type_free(&coordtype);CHKERRQ(ierr); 2801 } 2802 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2803 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 2804 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 2805 PetscFunctionReturn(0); 2806 } 2807 2808 /*@C 2809 DMPlexCreateFromCellListParallel - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM 2810 2811 Input Parameters: 2812 + comm - The communicator 2813 . dim - The topological dimension of the mesh 2814 . numCells - The number of cells owned by this process 2815 . numVertices - The number of vertices owned by this process 2816 . numCorners - The number of vertices for each cell 2817 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 2818 . cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 2819 . spaceDim - The spatial dimension used for coordinates 2820 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 2821 2822 Output Parameter: 2823 + dm - The DM 2824 - vertexSF - Optional, SF describing complete vertex ownership 2825 2826 Note: Two triangles sharing a face 2827 $ 2828 $ 2 2829 $ / | \ 2830 $ / | \ 2831 $ / | \ 2832 $ 0 0 | 1 3 2833 $ \ | / 2834 $ \ | / 2835 $ \ | / 2836 $ 1 2837 would have input 2838 $ numCells = 2, numVertices = 4 2839 $ cells = [0 1 2 1 3 2] 2840 $ 2841 which would result in the DMPlex 2842 $ 2843 $ 4 2844 $ / | \ 2845 $ / | \ 2846 $ / | \ 2847 $ 2 0 | 1 5 2848 $ \ | / 2849 $ \ | / 2850 $ \ | / 2851 $ 3 2852 2853 Level: beginner 2854 2855 .seealso: DMPlexCreateFromCellList(), DMPlexCreateFromDAG(), DMPlexCreate() 2856 @*/ 2857 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) 2858 { 2859 PetscSF sfVert; 2860 PetscErrorCode ierr; 2861 2862 PetscFunctionBegin; 2863 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 2864 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 2865 PetscValidLogicalCollectiveInt(*dm, dim, 2); 2866 PetscValidLogicalCollectiveInt(*dm, spaceDim, 8); 2867 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 2868 ierr = DMPlexBuildFromCellList_Parallel_Internal(*dm, spaceDim, numCells, numVertices, numCorners, cells, PETSC_FALSE, &sfVert);CHKERRQ(ierr); 2869 if (interpolate) { 2870 DM idm; 2871 2872 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 2873 ierr = DMDestroy(dm);CHKERRQ(ierr); 2874 *dm = idm; 2875 } 2876 ierr = DMPlexBuildCoordinates_Parallel_Internal(*dm, spaceDim, numCells, numVertices, sfVert, vertexCoords);CHKERRQ(ierr); 2877 if (vertexSF) *vertexSF = sfVert; 2878 else {ierr = PetscSFDestroy(&sfVert);CHKERRQ(ierr);} 2879 PetscFunctionReturn(0); 2880 } 2881 2882 /* 2883 This takes as input the common mesh generator output, a list of the vertices for each cell 2884 */ 2885 PetscErrorCode DMPlexBuildFromCellList_Internal(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[], PetscBool invertCells) 2886 { 2887 PetscInt *cone, c, p; 2888 PetscErrorCode ierr; 2889 2890 PetscFunctionBegin; 2891 ierr = DMPlexSetChart(dm, 0, numCells+numVertices);CHKERRQ(ierr); 2892 for (c = 0; c < numCells; ++c) { 2893 ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr); 2894 } 2895 ierr = DMSetUp(dm);CHKERRQ(ierr); 2896 ierr = DMGetWorkArray(dm, numCorners, MPIU_INT, &cone);CHKERRQ(ierr); 2897 for (c = 0; c < numCells; ++c) { 2898 for (p = 0; p < numCorners; ++p) { 2899 cone[p] = cells[c*numCorners+p]+numCells; 2900 } 2901 if (invertCells) { ierr = DMPlexInvertCell_Internal(spaceDim, numCorners, cone);CHKERRQ(ierr); } 2902 ierr = DMPlexSetCone(dm, c, cone);CHKERRQ(ierr); 2903 } 2904 ierr = DMRestoreWorkArray(dm, numCorners, MPIU_INT, &cone);CHKERRQ(ierr); 2905 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 2906 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 2907 PetscFunctionReturn(0); 2908 } 2909 2910 /* 2911 This takes as input the coordinates for each vertex 2912 */ 2913 PetscErrorCode DMPlexBuildCoordinates_Internal(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, const double vertexCoords[]) 2914 { 2915 PetscSection coordSection; 2916 Vec coordinates; 2917 DM cdm; 2918 PetscScalar *coords; 2919 PetscInt v, d; 2920 PetscErrorCode ierr; 2921 2922 PetscFunctionBegin; 2923 ierr = DMSetCoordinateDim(dm, spaceDim);CHKERRQ(ierr); 2924 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2925 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 2926 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 2927 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVertices);CHKERRQ(ierr); 2928 for (v = numCells; v < numCells+numVertices; ++v) { 2929 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 2930 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 2931 } 2932 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 2933 2934 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 2935 ierr = DMCreateLocalVector(cdm, &coordinates);CHKERRQ(ierr); 2936 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 2937 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 2938 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2939 for (v = 0; v < numVertices; ++v) { 2940 for (d = 0; d < spaceDim; ++d) { 2941 coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d]; 2942 } 2943 } 2944 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2945 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 2946 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 2947 PetscFunctionReturn(0); 2948 } 2949 2950 /*@C 2951 DMPlexCreateFromCellList - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM 2952 2953 Input Parameters: 2954 + comm - The communicator 2955 . dim - The topological dimension of the mesh 2956 . numCells - The number of cells 2957 . numVertices - The number of vertices 2958 . numCorners - The number of vertices for each cell 2959 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 2960 . cells - An array of numCells*numCorners numbers, the vertices for each cell 2961 . spaceDim - The spatial dimension used for coordinates 2962 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 2963 2964 Output Parameter: 2965 . dm - The DM 2966 2967 Note: Two triangles sharing a face 2968 $ 2969 $ 2 2970 $ / | \ 2971 $ / | \ 2972 $ / | \ 2973 $ 0 0 | 1 3 2974 $ \ | / 2975 $ \ | / 2976 $ \ | / 2977 $ 1 2978 would have input 2979 $ numCells = 2, numVertices = 4 2980 $ cells = [0 1 2 1 3 2] 2981 $ 2982 which would result in the DMPlex 2983 $ 2984 $ 4 2985 $ / | \ 2986 $ / | \ 2987 $ / | \ 2988 $ 2 0 | 1 5 2989 $ \ | / 2990 $ \ | / 2991 $ \ | / 2992 $ 3 2993 2994 Level: beginner 2995 2996 .seealso: DMPlexCreateFromDAG(), DMPlexCreate() 2997 @*/ 2998 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) 2999 { 3000 PetscErrorCode ierr; 3001 3002 PetscFunctionBegin; 3003 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 3004 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 3005 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 3006 ierr = DMPlexBuildFromCellList_Internal(*dm, spaceDim, numCells, numVertices, numCorners, cells, PETSC_FALSE);CHKERRQ(ierr); 3007 if (interpolate) { 3008 DM idm; 3009 3010 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 3011 ierr = DMDestroy(dm);CHKERRQ(ierr); 3012 *dm = idm; 3013 } 3014 ierr = DMPlexBuildCoordinates_Internal(*dm, spaceDim, numCells, numVertices, vertexCoords);CHKERRQ(ierr); 3015 PetscFunctionReturn(0); 3016 } 3017 3018 /*@ 3019 DMPlexCreateFromDAG - This takes as input the adjacency-list representation of the Directed Acyclic Graph (Hasse Diagram) encoding a mesh, and produces a DM 3020 3021 Input Parameters: 3022 + dm - The empty DM object, usually from DMCreate() and DMSetDimension() 3023 . depth - The depth of the DAG 3024 . numPoints - Array of size depth + 1 containing the number of points at each depth 3025 . coneSize - The cone size of each point 3026 . cones - The concatenation of the cone points for each point, the cone list must be oriented correctly for each point 3027 . coneOrientations - The orientation of each cone point 3028 - vertexCoords - An array of numPoints[0]*spacedim numbers representing the coordinates of each vertex, with spacedim the value set via DMSetCoordinateDim() 3029 3030 Output Parameter: 3031 . dm - The DM 3032 3033 Note: Two triangles sharing a face would have input 3034 $ depth = 1, numPoints = [4 2], coneSize = [3 3 0 0 0 0] 3035 $ cones = [2 3 4 3 5 4], coneOrientations = [0 0 0 0 0 0] 3036 $ vertexCoords = [-1.0 0.0 0.0 -1.0 0.0 1.0 1.0 0.0] 3037 $ 3038 which would result in the DMPlex 3039 $ 3040 $ 4 3041 $ / | \ 3042 $ / | \ 3043 $ / | \ 3044 $ 2 0 | 1 5 3045 $ \ | / 3046 $ \ | / 3047 $ \ | / 3048 $ 3 3049 $ 3050 $ Notice that all points are numbered consecutively, unlikely DMPlexCreateFromCellList() 3051 3052 Level: advanced 3053 3054 .seealso: DMPlexCreateFromCellList(), DMPlexCreate() 3055 @*/ 3056 PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[]) 3057 { 3058 Vec coordinates; 3059 PetscSection coordSection; 3060 PetscScalar *coords; 3061 PetscInt coordSize, firstVertex = -1, pStart = 0, pEnd = 0, p, v, dim, dimEmbed, d, off; 3062 PetscErrorCode ierr; 3063 3064 PetscFunctionBegin; 3065 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3066 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 3067 if (dimEmbed < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Embedding dimension %d cannot be less than intrinsic dimension %d",dimEmbed,dim); 3068 for (d = 0; d <= depth; ++d) pEnd += numPoints[d]; 3069 ierr = DMPlexSetChart(dm, pStart, pEnd);CHKERRQ(ierr); 3070 for (p = pStart; p < pEnd; ++p) { 3071 ierr = DMPlexSetConeSize(dm, p, coneSize[p-pStart]);CHKERRQ(ierr); 3072 if (firstVertex < 0 && !coneSize[p - pStart]) { 3073 firstVertex = p - pStart; 3074 } 3075 } 3076 if (firstVertex < 0 && numPoints[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected %d vertices but could not find any", numPoints[0]); 3077 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 3078 for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) { 3079 ierr = DMPlexSetCone(dm, p, &cones[off]);CHKERRQ(ierr); 3080 ierr = DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);CHKERRQ(ierr); 3081 } 3082 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 3083 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 3084 /* Build coordinates */ 3085 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 3086 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 3087 ierr = PetscSectionSetFieldComponents(coordSection, 0, dimEmbed);CHKERRQ(ierr); 3088 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0]);CHKERRQ(ierr); 3089 for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) { 3090 ierr = PetscSectionSetDof(coordSection, v, dimEmbed);CHKERRQ(ierr); 3091 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed);CHKERRQ(ierr); 3092 } 3093 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 3094 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 3095 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 3096 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 3097 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 3098 ierr = VecSetBlockSize(coordinates, dimEmbed);CHKERRQ(ierr); 3099 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 3100 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 3101 for (v = 0; v < numPoints[0]; ++v) { 3102 PetscInt off; 3103 3104 ierr = PetscSectionGetOffset(coordSection, v+firstVertex, &off);CHKERRQ(ierr); 3105 for (d = 0; d < dimEmbed; ++d) { 3106 coords[off+d] = vertexCoords[v*dimEmbed+d]; 3107 } 3108 } 3109 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 3110 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 3111 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 3112 PetscFunctionReturn(0); 3113 } 3114 3115 /*@C 3116 DMPlexCreateCellVertexFromFile - Create a DMPlex mesh from a simple cell-vertex file. 3117 3118 + comm - The MPI communicator 3119 . filename - Name of the .dat file 3120 - interpolate - Create faces and edges in the mesh 3121 3122 Output Parameter: 3123 . dm - The DM object representing the mesh 3124 3125 Note: The format is the simplest possible: 3126 $ Ne 3127 $ v0 v1 ... vk 3128 $ Nv 3129 $ x y z marker 3130 3131 Level: beginner 3132 3133 .seealso: DMPlexCreateFromFile(), DMPlexCreateMedFromFile(), DMPlexCreateGmsh(), DMPlexCreate() 3134 @*/ 3135 PetscErrorCode DMPlexCreateCellVertexFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 3136 { 3137 DMLabel marker; 3138 PetscViewer viewer; 3139 Vec coordinates; 3140 PetscSection coordSection; 3141 PetscScalar *coords; 3142 char line[PETSC_MAX_PATH_LEN]; 3143 PetscInt dim = 3, cdim = 3, coordSize, v, c, d; 3144 PetscMPIInt rank; 3145 int snum, Nv, Nc; 3146 PetscErrorCode ierr; 3147 3148 PetscFunctionBegin; 3149 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 3150 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 3151 ierr = PetscViewerSetType(viewer, PETSCVIEWERASCII);CHKERRQ(ierr); 3152 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 3153 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 3154 if (!rank) { 3155 ierr = PetscViewerRead(viewer, line, 2, NULL, PETSC_STRING);CHKERRQ(ierr); 3156 snum = sscanf(line, "%d %d", &Nc, &Nv); 3157 if (snum != 2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 3158 } else { 3159 Nc = Nv = 0; 3160 } 3161 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 3162 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 3163 ierr = DMPlexSetChart(*dm, 0, Nc+Nv);CHKERRQ(ierr); 3164 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 3165 ierr = DMSetCoordinateDim(*dm, cdim);CHKERRQ(ierr); 3166 /* Read topology */ 3167 if (!rank) { 3168 PetscInt cone[8], corners = 8; 3169 int vbuf[8], v; 3170 3171 for (c = 0; c < Nc; ++c) {ierr = DMPlexSetConeSize(*dm, c, corners);CHKERRQ(ierr);} 3172 ierr = DMSetUp(*dm);CHKERRQ(ierr); 3173 for (c = 0; c < Nc; ++c) { 3174 ierr = PetscViewerRead(viewer, line, corners, NULL, PETSC_STRING);CHKERRQ(ierr); 3175 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]); 3176 if (snum != corners) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 3177 for (v = 0; v < corners; ++v) cone[v] = vbuf[v] + Nc; 3178 /* Hexahedra are inverted */ 3179 { 3180 PetscInt tmp = cone[1]; 3181 cone[1] = cone[3]; 3182 cone[3] = tmp; 3183 } 3184 ierr = DMPlexSetCone(*dm, c, cone);CHKERRQ(ierr); 3185 } 3186 } 3187 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 3188 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 3189 /* Read coordinates */ 3190 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 3191 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 3192 ierr = PetscSectionSetFieldComponents(coordSection, 0, cdim);CHKERRQ(ierr); 3193 ierr = PetscSectionSetChart(coordSection, Nc, Nc + Nv);CHKERRQ(ierr); 3194 for (v = Nc; v < Nc+Nv; ++v) { 3195 ierr = PetscSectionSetDof(coordSection, v, cdim);CHKERRQ(ierr); 3196 ierr = PetscSectionSetFieldDof(coordSection, v, 0, cdim);CHKERRQ(ierr); 3197 } 3198 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 3199 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 3200 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 3201 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 3202 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 3203 ierr = VecSetBlockSize(coordinates, cdim);CHKERRQ(ierr); 3204 ierr = VecSetType(coordinates, VECSTANDARD);CHKERRQ(ierr); 3205 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 3206 if (!rank) { 3207 double x[3]; 3208 int val; 3209 3210 ierr = DMCreateLabel(*dm, "marker");CHKERRQ(ierr); 3211 ierr = DMGetLabel(*dm, "marker", &marker);CHKERRQ(ierr); 3212 for (v = 0; v < Nv; ++v) { 3213 ierr = PetscViewerRead(viewer, line, 4, NULL, PETSC_STRING);CHKERRQ(ierr); 3214 snum = sscanf(line, "%lg %lg %lg %d", &x[0], &x[1], &x[2], &val); 3215 if (snum != 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 3216 for (d = 0; d < cdim; ++d) coords[v*cdim+d] = x[d]; 3217 if (val) {ierr = DMLabelSetValue(marker, v+Nc, val);CHKERRQ(ierr);} 3218 } 3219 } 3220 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 3221 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 3222 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 3223 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3224 if (interpolate) { 3225 DM idm; 3226 DMLabel bdlabel; 3227 3228 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 3229 ierr = DMDestroy(dm);CHKERRQ(ierr); 3230 *dm = idm; 3231 3232 ierr = DMGetLabel(*dm, "marker", &bdlabel);CHKERRQ(ierr); 3233 ierr = DMPlexMarkBoundaryFaces(*dm, PETSC_DETERMINE, bdlabel);CHKERRQ(ierr); 3234 ierr = DMPlexLabelComplete(*dm, bdlabel);CHKERRQ(ierr); 3235 } 3236 PetscFunctionReturn(0); 3237 } 3238 3239 /*@C 3240 DMPlexCreateFromFile - This takes a filename and produces a DM 3241 3242 Input Parameters: 3243 + comm - The communicator 3244 . filename - A file name 3245 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 3246 3247 Output Parameter: 3248 . dm - The DM 3249 3250 Options Database Keys: 3251 . -dm_plex_create_from_hdf5_xdmf - use the PETSC_VIEWER_HDF5_XDMF format for reading HDF5 3252 3253 Level: beginner 3254 3255 .seealso: DMPlexCreateFromDAG(), DMPlexCreateFromCellList(), DMPlexCreate() 3256 @*/ 3257 PetscErrorCode DMPlexCreateFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 3258 { 3259 const char *extGmsh = ".msh"; 3260 const char *extGmsh2 = ".msh2"; 3261 const char *extGmsh4 = ".msh4"; 3262 const char *extCGNS = ".cgns"; 3263 const char *extExodus = ".exo"; 3264 const char *extGenesis = ".gen"; 3265 const char *extFluent = ".cas"; 3266 const char *extHDF5 = ".h5"; 3267 const char *extMed = ".med"; 3268 const char *extPLY = ".ply"; 3269 const char *extCV = ".dat"; 3270 size_t len; 3271 PetscBool isGmsh, isGmsh2, isGmsh4, isCGNS, isExodus, isGenesis, isFluent, isHDF5, isMed, isPLY, isCV; 3272 PetscMPIInt rank; 3273 PetscErrorCode ierr; 3274 3275 PetscFunctionBegin; 3276 PetscValidPointer(filename, 2); 3277 PetscValidPointer(dm, 4); 3278 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 3279 ierr = PetscStrlen(filename, &len);CHKERRQ(ierr); 3280 if (!len) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Filename must be a valid path"); 3281 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGmsh, 4, &isGmsh);CHKERRQ(ierr); 3282 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extGmsh2, 5, &isGmsh2);CHKERRQ(ierr); 3283 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extGmsh4, 5, &isGmsh4);CHKERRQ(ierr); 3284 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extCGNS, 5, &isCGNS);CHKERRQ(ierr); 3285 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extExodus, 4, &isExodus);CHKERRQ(ierr); 3286 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGenesis, 4, &isGenesis);CHKERRQ(ierr); 3287 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extFluent, 4, &isFluent);CHKERRQ(ierr); 3288 ierr = PetscStrncmp(&filename[PetscMax(0,len-3)], extHDF5, 3, &isHDF5);CHKERRQ(ierr); 3289 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extMed, 4, &isMed);CHKERRQ(ierr); 3290 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extPLY, 4, &isPLY);CHKERRQ(ierr); 3291 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extCV, 4, &isCV);CHKERRQ(ierr); 3292 if (isGmsh || isGmsh2 || isGmsh4) { 3293 ierr = DMPlexCreateGmshFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3294 } else if (isCGNS) { 3295 ierr = DMPlexCreateCGNSFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3296 } else if (isExodus || isGenesis) { 3297 ierr = DMPlexCreateExodusFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3298 } else if (isFluent) { 3299 ierr = DMPlexCreateFluentFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3300 } else if (isHDF5) { 3301 PetscBool load_hdf5_xdmf = PETSC_FALSE; 3302 PetscViewer viewer; 3303 3304 /* PETSC_VIEWER_HDF5_XDMF is used if the filename ends with .xdmf.h5, or if -dm_plex_create_from_hdf5_xdmf option is present */ 3305 ierr = PetscStrncmp(&filename[PetscMax(0,len-8)], ".xdmf", 5, &load_hdf5_xdmf);CHKERRQ(ierr); 3306 ierr = PetscOptionsGetBool(NULL, NULL, "-dm_plex_create_from_hdf5_xdmf", &load_hdf5_xdmf, NULL);CHKERRQ(ierr); 3307 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 3308 ierr = PetscViewerSetType(viewer, PETSCVIEWERHDF5);CHKERRQ(ierr); 3309 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 3310 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 3311 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 3312 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 3313 if (load_hdf5_xdmf) {ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_XDMF);CHKERRQ(ierr);} 3314 ierr = DMLoad(*dm, viewer);CHKERRQ(ierr); 3315 if (load_hdf5_xdmf) {ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);} 3316 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3317 3318 if (interpolate) { 3319 DM idm; 3320 3321 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 3322 ierr = DMDestroy(dm);CHKERRQ(ierr); 3323 *dm = idm; 3324 } 3325 } else if (isMed) { 3326 ierr = DMPlexCreateMedFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3327 } else if (isPLY) { 3328 ierr = DMPlexCreatePLYFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3329 } else if (isCV) { 3330 ierr = DMPlexCreateCellVertexFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3331 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot load file %s: unrecognized extension", filename); 3332 PetscFunctionReturn(0); 3333 } 3334 3335 /*@ 3336 DMPlexCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell 3337 3338 Collective on comm 3339 3340 Input Parameters: 3341 + comm - The communicator 3342 . dim - The spatial dimension 3343 - simplex - Flag for simplex, otherwise use a tensor-product cell 3344 3345 Output Parameter: 3346 . refdm - The reference cell 3347 3348 Level: intermediate 3349 3350 .keywords: reference cell 3351 .seealso: 3352 @*/ 3353 PetscErrorCode DMPlexCreateReferenceCell(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *refdm) 3354 { 3355 DM rdm; 3356 Vec coords; 3357 PetscErrorCode ierr; 3358 3359 PetscFunctionBeginUser; 3360 ierr = DMCreate(comm, &rdm);CHKERRQ(ierr); 3361 ierr = DMSetType(rdm, DMPLEX);CHKERRQ(ierr); 3362 ierr = DMSetDimension(rdm, dim);CHKERRQ(ierr); 3363 switch (dim) { 3364 case 0: 3365 { 3366 PetscInt numPoints[1] = {1}; 3367 PetscInt coneSize[1] = {0}; 3368 PetscInt cones[1] = {0}; 3369 PetscInt coneOrientations[1] = {0}; 3370 PetscScalar vertexCoords[1] = {0.0}; 3371 3372 ierr = DMPlexCreateFromDAG(rdm, 0, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3373 } 3374 break; 3375 case 1: 3376 { 3377 PetscInt numPoints[2] = {2, 1}; 3378 PetscInt coneSize[3] = {2, 0, 0}; 3379 PetscInt cones[2] = {1, 2}; 3380 PetscInt coneOrientations[2] = {0, 0}; 3381 PetscScalar vertexCoords[2] = {-1.0, 1.0}; 3382 3383 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3384 } 3385 break; 3386 case 2: 3387 if (simplex) { 3388 PetscInt numPoints[2] = {3, 1}; 3389 PetscInt coneSize[4] = {3, 0, 0, 0}; 3390 PetscInt cones[3] = {1, 2, 3}; 3391 PetscInt coneOrientations[3] = {0, 0, 0}; 3392 PetscScalar vertexCoords[6] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0}; 3393 3394 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3395 } else { 3396 PetscInt numPoints[2] = {4, 1}; 3397 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 3398 PetscInt cones[4] = {1, 2, 3, 4}; 3399 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 3400 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0}; 3401 3402 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3403 } 3404 break; 3405 case 3: 3406 if (simplex) { 3407 PetscInt numPoints[2] = {4, 1}; 3408 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 3409 PetscInt cones[4] = {1, 3, 2, 4}; 3410 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 3411 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}; 3412 3413 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3414 } else { 3415 PetscInt numPoints[2] = {8, 1}; 3416 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0}; 3417 PetscInt cones[8] = {1, 4, 3, 2, 5, 6, 7, 8}; 3418 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 3419 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, 3420 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0}; 3421 3422 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3423 } 3424 break; 3425 default: 3426 SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Cannot create reference cell for dimension %d", dim); 3427 } 3428 ierr = DMPlexInterpolate(rdm, refdm);CHKERRQ(ierr); 3429 if (rdm->coordinateDM) { 3430 DM ncdm; 3431 PetscSection cs; 3432 PetscInt pEnd = -1; 3433 3434 ierr = DMGetSection(rdm->coordinateDM, &cs);CHKERRQ(ierr); 3435 if (cs) {ierr = PetscSectionGetChart(cs, NULL, &pEnd);CHKERRQ(ierr);} 3436 if (pEnd >= 0) { 3437 ierr = DMClone(rdm->coordinateDM, &ncdm);CHKERRQ(ierr); 3438 ierr = DMCopyDisc(rdm->coordinateDM, ncdm);CHKERRQ(ierr); 3439 ierr = DMSetSection(ncdm, cs);CHKERRQ(ierr); 3440 ierr = DMSetCoordinateDM(*refdm, ncdm);CHKERRQ(ierr); 3441 ierr = DMDestroy(&ncdm);CHKERRQ(ierr); 3442 } 3443 } 3444 ierr = DMGetCoordinatesLocal(rdm, &coords);CHKERRQ(ierr); 3445 if (coords) { 3446 ierr = DMSetCoordinatesLocal(*refdm, coords);CHKERRQ(ierr); 3447 } else { 3448 ierr = DMGetCoordinates(rdm, &coords);CHKERRQ(ierr); 3449 if (coords) {ierr = DMSetCoordinates(*refdm, coords);CHKERRQ(ierr);} 3450 } 3451 ierr = DMDestroy(&rdm);CHKERRQ(ierr); 3452 PetscFunctionReturn(0); 3453 } 3454