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