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