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