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