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