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