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, sfVert; 2877 PetscLayout vLayout; 2878 PetscSFNode *vertexLocal, *vertexOwner, *remoteVertex; 2879 PetscInt numVerticesAdj, *verticesAdj, numVerticesGhost = 0, *localVertex, *cones, c, p, v, g, dim; 2880 PetscMPIInt rank, size; 2881 PetscErrorCode ierr; 2882 2883 PetscFunctionBegin; 2884 PetscValidLogicalCollectiveInt(dm,NVertices,4); 2885 ierr = PetscLogEventBegin(DMPLEX_BuildFromCellList,dm,0,0,0);CHKERRQ(ierr); 2886 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRMPI(ierr); 2887 ierr = MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);CHKERRMPI(ierr); 2888 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2889 /* Get/check global number of vertices */ 2890 { 2891 PetscInt NVerticesInCells, i; 2892 const PetscInt len = numCells * numCorners; 2893 2894 /* NVerticesInCells = max(cells) + 1 */ 2895 NVerticesInCells = PETSC_MIN_INT; 2896 for (i=0; i<len; i++) if (cells[i] > NVerticesInCells) NVerticesInCells = cells[i]; 2897 ++NVerticesInCells; 2898 ierr = MPI_Allreduce(MPI_IN_PLACE, &NVerticesInCells, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject) dm));CHKERRMPI(ierr); 2899 2900 if (numVertices == PETSC_DECIDE && NVertices == PETSC_DECIDE) NVertices = NVerticesInCells; 2901 else if (NVertices != PETSC_DECIDE && NVertices < NVerticesInCells) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Specified global number of vertices %D must be greater than or equal to the number of vertices in cells %D",NVertices,NVerticesInCells); 2902 } 2903 /* Partition vertices */ 2904 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) dm), &vLayout);CHKERRQ(ierr); 2905 ierr = PetscLayoutSetLocalSize(vLayout, numVertices);CHKERRQ(ierr); 2906 ierr = PetscLayoutSetSize(vLayout, NVertices);CHKERRQ(ierr); 2907 ierr = PetscLayoutSetBlockSize(vLayout, 1);CHKERRQ(ierr); 2908 ierr = PetscLayoutSetUp(vLayout);CHKERRQ(ierr); 2909 /* Count locally unique vertices */ 2910 { 2911 PetscHSetI vhash; 2912 PetscInt off = 0; 2913 2914 ierr = PetscHSetICreate(&vhash);CHKERRQ(ierr); 2915 for (c = 0; c < numCells; ++c) { 2916 for (p = 0; p < numCorners; ++p) { 2917 ierr = PetscHSetIAdd(vhash, cells[c*numCorners+p]);CHKERRQ(ierr); 2918 } 2919 } 2920 ierr = PetscHSetIGetSize(vhash, &numVerticesAdj);CHKERRQ(ierr); 2921 ierr = PetscMalloc1(numVerticesAdj, &verticesAdj);CHKERRQ(ierr); 2922 ierr = PetscHSetIGetElems(vhash, &off, verticesAdj);CHKERRQ(ierr); 2923 ierr = PetscHSetIDestroy(&vhash);CHKERRQ(ierr); 2924 if (off != numVerticesAdj) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid number of local vertices %D should be %D", off, numVerticesAdj); 2925 } 2926 ierr = PetscSortInt(numVerticesAdj, verticesAdj);CHKERRQ(ierr); 2927 /* Create cones */ 2928 ierr = DMPlexSetChart(dm, 0, numCells+numVerticesAdj);CHKERRQ(ierr); 2929 for (c = 0; c < numCells; ++c) {ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr);} 2930 ierr = DMSetUp(dm);CHKERRQ(ierr); 2931 ierr = DMPlexGetCones(dm,&cones);CHKERRQ(ierr); 2932 for (c = 0; c < numCells; ++c) { 2933 for (p = 0; p < numCorners; ++p) { 2934 const PetscInt gv = cells[c*numCorners+p]; 2935 PetscInt lv; 2936 2937 /* Positions within verticesAdj form 0-based local vertex numbering; 2938 we need to shift it by numCells to get correct DAG points (cells go first) */ 2939 ierr = PetscFindInt(gv, numVerticesAdj, verticesAdj, &lv);CHKERRQ(ierr); 2940 if (lv < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not find global vertex %D in local connectivity", gv); 2941 cones[c*numCorners+p] = lv+numCells; 2942 } 2943 } 2944 /* Create SF for vertices */ 2945 ierr = PetscSFCreate(PetscObjectComm((PetscObject)dm), &sfVert);CHKERRQ(ierr); 2946 ierr = PetscObjectSetName((PetscObject) sfVert, "Vertex Ownership SF");CHKERRQ(ierr); 2947 ierr = PetscSFSetFromOptions(sfVert);CHKERRQ(ierr); 2948 ierr = PetscSFSetGraphLayout(sfVert, vLayout, numVerticesAdj, NULL, PETSC_OWN_POINTER, verticesAdj);CHKERRQ(ierr); 2949 ierr = PetscFree(verticesAdj);CHKERRQ(ierr); 2950 /* Build pointSF */ 2951 ierr = PetscMalloc2(numVerticesAdj, &vertexLocal, numVertices, &vertexOwner);CHKERRQ(ierr); 2952 for (v = 0; v < numVerticesAdj; ++v) {vertexLocal[v].index = v+numCells; vertexLocal[v].rank = rank;} 2953 for (v = 0; v < numVertices; ++v) {vertexOwner[v].index = -1; vertexOwner[v].rank = -1;} 2954 ierr = PetscSFReduceBegin(sfVert, MPIU_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);CHKERRQ(ierr); 2955 ierr = PetscSFReduceEnd(sfVert, MPIU_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);CHKERRQ(ierr); 2956 for (v = 0; v < numVertices; ++v) if (vertexOwner[v].rank < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Global vertex %D was unclaimed", v + vLayout->rstart); 2957 ierr = PetscSFBcastBegin(sfVert, MPIU_2INT, vertexOwner, vertexLocal,MPI_REPLACE);CHKERRQ(ierr); 2958 ierr = PetscSFBcastEnd(sfVert, MPIU_2INT, vertexOwner, vertexLocal,MPI_REPLACE);CHKERRQ(ierr); 2959 for (v = 0; v < numVerticesAdj; ++v) if (vertexLocal[v].rank != rank) ++numVerticesGhost; 2960 ierr = PetscMalloc1(numVerticesGhost, &localVertex);CHKERRQ(ierr); 2961 ierr = PetscMalloc1(numVerticesGhost, &remoteVertex);CHKERRQ(ierr); 2962 for (v = 0, g = 0; v < numVerticesAdj; ++v) { 2963 if (vertexLocal[v].rank != rank) { 2964 localVertex[g] = v+numCells; 2965 remoteVertex[g].index = vertexLocal[v].index; 2966 remoteVertex[g].rank = vertexLocal[v].rank; 2967 ++g; 2968 } 2969 } 2970 ierr = PetscFree2(vertexLocal, vertexOwner);CHKERRQ(ierr); 2971 if (g != numVerticesGhost) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of vertex ghosts %D should be %D", g, numVerticesGhost); 2972 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 2973 ierr = PetscObjectSetName((PetscObject) sfPoint, "point SF");CHKERRQ(ierr); 2974 ierr = PetscSFSetGraph(sfPoint, numCells+numVerticesAdj, numVerticesGhost, localVertex, PETSC_OWN_POINTER, remoteVertex, PETSC_OWN_POINTER);CHKERRQ(ierr); 2975 ierr = PetscLayoutDestroy(&vLayout);CHKERRQ(ierr); 2976 /* Fill in the rest of the topology structure */ 2977 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 2978 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 2979 ierr = PetscLogEventEnd(DMPLEX_BuildFromCellList,dm,0,0,0);CHKERRQ(ierr); 2980 if (vertexSF) *vertexSF = sfVert; 2981 else {ierr = PetscSFDestroy(&sfVert);CHKERRQ(ierr);} 2982 PetscFunctionReturn(0); 2983 } 2984 2985 /*@C 2986 DMPlexBuildCoordinatesFromCellListParallel - Build DM coordinates from a list of coordinates for each owned vertex (common mesh generator output) 2987 2988 Input Parameters: 2989 + dm - The DM 2990 . spaceDim - The spatial dimension used for coordinates 2991 . sfVert - SF describing complete vertex ownership 2992 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 2993 2994 Level: advanced 2995 2996 Notes: 2997 Not currently supported in Fortran. 2998 2999 .seealso: DMPlexBuildCoordinatesFromCellList(), DMPlexCreateFromCellListParallelPetsc(), DMPlexBuildFromCellListParallel() 3000 @*/ 3001 PetscErrorCode DMPlexBuildCoordinatesFromCellListParallel(DM dm, PetscInt spaceDim, PetscSF sfVert, const PetscReal vertexCoords[]) 3002 { 3003 PetscSection coordSection; 3004 Vec coordinates; 3005 PetscScalar *coords; 3006 PetscInt numVertices, numVerticesAdj, coordSize, v, vStart, vEnd; 3007 PetscErrorCode ierr; 3008 3009 PetscFunctionBegin; 3010 ierr = PetscLogEventBegin(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0);CHKERRQ(ierr); 3011 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 3012 if (vStart < 0 || vEnd < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM is not set up properly. DMPlexBuildFromCellList() should be called first."); 3013 ierr = DMSetCoordinateDim(dm, spaceDim);CHKERRQ(ierr); 3014 ierr = PetscSFGetGraph(sfVert, &numVertices, &numVerticesAdj, NULL, NULL);CHKERRQ(ierr); 3015 if (vEnd - vStart != numVerticesAdj) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Supplied sfVert has wrong number of leaves = %D != %D = vEnd - vStart",numVerticesAdj,vEnd - vStart); 3016 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 3017 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 3018 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 3019 ierr = PetscSectionSetChart(coordSection, vStart, vEnd);CHKERRQ(ierr); 3020 for (v = vStart; v < vEnd; ++v) { 3021 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 3022 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 3023 } 3024 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 3025 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 3026 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 3027 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 3028 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 3029 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 3030 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 3031 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 3032 { 3033 MPI_Datatype coordtype; 3034 3035 /* Need a temp buffer for coords if we have complex/single */ 3036 ierr = MPI_Type_contiguous(spaceDim, MPIU_SCALAR, &coordtype);CHKERRMPI(ierr); 3037 ierr = MPI_Type_commit(&coordtype);CHKERRMPI(ierr); 3038 #if defined(PETSC_USE_COMPLEX) 3039 { 3040 PetscScalar *svertexCoords; 3041 PetscInt i; 3042 ierr = PetscMalloc1(numVertices*spaceDim,&svertexCoords);CHKERRQ(ierr); 3043 for (i=0; i<numVertices*spaceDim; i++) svertexCoords[i] = vertexCoords[i]; 3044 ierr = PetscSFBcastBegin(sfVert, coordtype, svertexCoords, coords,MPI_REPLACE);CHKERRQ(ierr); 3045 ierr = PetscSFBcastEnd(sfVert, coordtype, svertexCoords, coords,MPI_REPLACE);CHKERRQ(ierr); 3046 ierr = PetscFree(svertexCoords);CHKERRQ(ierr); 3047 } 3048 #else 3049 ierr = PetscSFBcastBegin(sfVert, coordtype, vertexCoords, coords,MPI_REPLACE);CHKERRQ(ierr); 3050 ierr = PetscSFBcastEnd(sfVert, coordtype, vertexCoords, coords,MPI_REPLACE);CHKERRQ(ierr); 3051 #endif 3052 ierr = MPI_Type_free(&coordtype);CHKERRMPI(ierr); 3053 } 3054 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 3055 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 3056 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 3057 ierr = PetscLogEventEnd(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0);CHKERRQ(ierr); 3058 PetscFunctionReturn(0); 3059 } 3060 3061 /*@ 3062 DMPlexCreateFromCellListParallelPetsc - Create distributed DMPLEX from a list of vertices for each cell (common mesh generator output) 3063 3064 Input Parameters: 3065 + comm - The communicator 3066 . dim - The topological dimension of the mesh 3067 . numCells - The number of cells owned by this process 3068 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE 3069 . NVertices - The global number of vertices, or PETSC_DECIDE 3070 . numCorners - The number of vertices for each cell 3071 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 3072 . cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 3073 . spaceDim - The spatial dimension used for coordinates 3074 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 3075 3076 Output Parameter: 3077 + dm - The DM 3078 - vertexSF - (Optional) SF describing complete vertex ownership 3079 3080 Notes: 3081 This function is just a convenient sequence of DMCreate(), DMSetType(), DMSetDimension(), 3082 DMPlexBuildFromCellListParallel(), DMPlexInterpolate(), DMPlexBuildCoordinatesFromCellListParallel() 3083 3084 See DMPlexBuildFromCellListParallel() for an example and details about the topology-related parameters. 3085 See DMPlexBuildCoordinatesFromCellListParallel() for details about the geometry-related parameters. 3086 3087 Level: intermediate 3088 3089 .seealso: DMPlexCreateFromCellListPetsc(), DMPlexBuildFromCellListParallel(), DMPlexBuildCoordinatesFromCellListParallel(), DMPlexCreateFromDAG(), DMPlexCreate() 3090 @*/ 3091 PetscErrorCode DMPlexCreateFromCellListParallelPetsc(MPI_Comm comm, PetscInt dim, PetscInt numCells, PetscInt numVertices, PetscInt NVertices, PetscInt numCorners, PetscBool interpolate, const PetscInt cells[], PetscInt spaceDim, const PetscReal vertexCoords[], PetscSF *vertexSF, DM *dm) 3092 { 3093 PetscSF sfVert; 3094 PetscErrorCode ierr; 3095 3096 PetscFunctionBegin; 3097 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 3098 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 3099 PetscValidLogicalCollectiveInt(*dm, dim, 2); 3100 PetscValidLogicalCollectiveInt(*dm, spaceDim, 8); 3101 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 3102 ierr = DMPlexBuildFromCellListParallel(*dm, numCells, numVertices, NVertices, numCorners, cells, &sfVert);CHKERRQ(ierr); 3103 if (interpolate) { 3104 DM idm; 3105 3106 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 3107 ierr = DMDestroy(dm);CHKERRQ(ierr); 3108 *dm = idm; 3109 } 3110 ierr = DMPlexBuildCoordinatesFromCellListParallel(*dm, spaceDim, sfVert, vertexCoords);CHKERRQ(ierr); 3111 if (vertexSF) *vertexSF = sfVert; 3112 else {ierr = PetscSFDestroy(&sfVert);CHKERRQ(ierr);} 3113 PetscFunctionReturn(0); 3114 } 3115 3116 3117 /*@ 3118 DMPlexCreateFromCellListParallel - Deprecated, use DMPlexCreateFromCellListParallelPetsc() 3119 3120 Level: deprecated 3121 3122 .seealso: DMPlexCreateFromCellListParallelPetsc() 3123 @*/ 3124 PetscErrorCode DMPlexCreateFromCellListParallel(MPI_Comm comm, PetscInt dim, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, PetscBool interpolate, const int cells[], PetscInt spaceDim, const PetscReal vertexCoords[], PetscSF *vertexSF, DM *dm) 3125 { 3126 PetscErrorCode ierr; 3127 PetscInt i; 3128 PetscInt *pintCells; 3129 3130 PetscFunctionBegin; 3131 if (sizeof(int) > sizeof(PetscInt)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Size of int %zd greater than size of PetscInt %zd. Reconfigure PETSc --with-64-bit-indices=1", sizeof(int), sizeof(PetscInt)); 3132 if (sizeof(int) == sizeof(PetscInt)) { 3133 pintCells = (PetscInt *) cells; 3134 } else { 3135 ierr = PetscMalloc1(numCells*numCorners, &pintCells);CHKERRQ(ierr); 3136 for (i = 0; i < numCells*numCorners; i++) { 3137 pintCells[i] = (PetscInt) cells[i]; 3138 } 3139 } 3140 ierr = DMPlexCreateFromCellListParallelPetsc(comm, dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, pintCells, spaceDim, vertexCoords, vertexSF, dm);CHKERRQ(ierr); 3141 if (sizeof(int) != sizeof(PetscInt)) { 3142 ierr = PetscFree(pintCells);CHKERRQ(ierr); 3143 } 3144 PetscFunctionReturn(0); 3145 } 3146 3147 /*@C 3148 DMPlexBuildFromCellList - Build DMPLEX topology from a list of vertices for each cell (common mesh generator output) 3149 3150 Input Parameters: 3151 + dm - The DM 3152 . numCells - The number of cells owned by this process 3153 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE 3154 . numCorners - The number of vertices for each cell 3155 - cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 3156 3157 Level: advanced 3158 3159 Notes: 3160 Two triangles sharing a face 3161 $ 3162 $ 2 3163 $ / | \ 3164 $ / | \ 3165 $ / | \ 3166 $ 0 0 | 1 3 3167 $ \ | / 3168 $ \ | / 3169 $ \ | / 3170 $ 1 3171 would have input 3172 $ numCells = 2, numVertices = 4 3173 $ cells = [0 1 2 1 3 2] 3174 $ 3175 which would result in the DMPlex 3176 $ 3177 $ 4 3178 $ / | \ 3179 $ / | \ 3180 $ / | \ 3181 $ 2 0 | 1 5 3182 $ \ | / 3183 $ \ | / 3184 $ \ | / 3185 $ 3 3186 3187 If numVertices is PETSC_DECIDE, it is computed by PETSc as the maximum vertex index in cells + 1. 3188 3189 Not currently supported in Fortran. 3190 3191 .seealso: DMPlexBuildFromCellListParallel(), DMPlexBuildCoordinatesFromCellList(), DMPlexCreateFromCellListPetsc() 3192 @*/ 3193 PetscErrorCode DMPlexBuildFromCellList(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const PetscInt cells[]) 3194 { 3195 PetscInt *cones, c, p, dim; 3196 PetscErrorCode ierr; 3197 3198 PetscFunctionBegin; 3199 ierr = PetscLogEventBegin(DMPLEX_BuildFromCellList,dm,0,0,0);CHKERRQ(ierr); 3200 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3201 /* Get/check global number of vertices */ 3202 { 3203 PetscInt NVerticesInCells, i; 3204 const PetscInt len = numCells * numCorners; 3205 3206 /* NVerticesInCells = max(cells) + 1 */ 3207 NVerticesInCells = PETSC_MIN_INT; 3208 for (i=0; i<len; i++) if (cells[i] > NVerticesInCells) NVerticesInCells = cells[i]; 3209 ++NVerticesInCells; 3210 3211 if (numVertices == PETSC_DECIDE) numVertices = NVerticesInCells; 3212 else if (numVertices < NVerticesInCells) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Specified number of vertices %D must be greater than or equal to the number of vertices in cells %D",numVertices,NVerticesInCells); 3213 } 3214 ierr = DMPlexSetChart(dm, 0, numCells+numVertices);CHKERRQ(ierr); 3215 for (c = 0; c < numCells; ++c) { 3216 ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr); 3217 } 3218 ierr = DMSetUp(dm);CHKERRQ(ierr); 3219 ierr = DMPlexGetCones(dm,&cones);CHKERRQ(ierr); 3220 for (c = 0; c < numCells; ++c) { 3221 for (p = 0; p < numCorners; ++p) { 3222 cones[c*numCorners+p] = cells[c*numCorners+p]+numCells; 3223 } 3224 } 3225 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 3226 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 3227 ierr = PetscLogEventEnd(DMPLEX_BuildFromCellList,dm,0,0,0);CHKERRQ(ierr); 3228 PetscFunctionReturn(0); 3229 } 3230 3231 /*@C 3232 DMPlexBuildCoordinatesFromCellList - Build DM coordinates from a list of coordinates for each owned vertex (common mesh generator output) 3233 3234 Input Parameters: 3235 + dm - The DM 3236 . spaceDim - The spatial dimension used for coordinates 3237 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 3238 3239 Level: advanced 3240 3241 Notes: 3242 Not currently supported in Fortran. 3243 3244 .seealso: DMPlexBuildCoordinatesFromCellListParallel(), DMPlexCreateFromCellListPetsc(), DMPlexBuildFromCellList() 3245 @*/ 3246 PetscErrorCode DMPlexBuildCoordinatesFromCellList(DM dm, PetscInt spaceDim, const PetscReal vertexCoords[]) 3247 { 3248 PetscSection coordSection; 3249 Vec coordinates; 3250 DM cdm; 3251 PetscScalar *coords; 3252 PetscInt v, vStart, vEnd, d; 3253 PetscErrorCode ierr; 3254 3255 PetscFunctionBegin; 3256 ierr = PetscLogEventBegin(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0);CHKERRQ(ierr); 3257 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 3258 if (vStart < 0 || vEnd < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM is not set up properly. DMPlexBuildFromCellList() should be called first."); 3259 ierr = DMSetCoordinateDim(dm, spaceDim);CHKERRQ(ierr); 3260 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 3261 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 3262 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 3263 ierr = PetscSectionSetChart(coordSection, vStart, vEnd);CHKERRQ(ierr); 3264 for (v = vStart; v < vEnd; ++v) { 3265 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 3266 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 3267 } 3268 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 3269 3270 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 3271 ierr = DMCreateLocalVector(cdm, &coordinates);CHKERRQ(ierr); 3272 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 3273 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 3274 ierr = VecGetArrayWrite(coordinates, &coords);CHKERRQ(ierr); 3275 for (v = 0; v < vEnd-vStart; ++v) { 3276 for (d = 0; d < spaceDim; ++d) { 3277 coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d]; 3278 } 3279 } 3280 ierr = VecRestoreArrayWrite(coordinates, &coords);CHKERRQ(ierr); 3281 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 3282 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 3283 ierr = PetscLogEventEnd(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0);CHKERRQ(ierr); 3284 PetscFunctionReturn(0); 3285 } 3286 3287 /*@ 3288 DMPlexCreateFromCellListPetsc - Create DMPLEX from a list of vertices for each cell (common mesh generator output) 3289 3290 Input Parameters: 3291 + comm - The communicator 3292 . dim - The topological dimension of the mesh 3293 . numCells - The number of cells 3294 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE 3295 . numCorners - The number of vertices for each cell 3296 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 3297 . cells - An array of numCells*numCorners numbers, the vertices for each cell 3298 . spaceDim - The spatial dimension used for coordinates 3299 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 3300 3301 Output Parameter: 3302 . dm - The DM 3303 3304 Notes: 3305 This function is just a convenient sequence of DMCreate(), DMSetType(), DMSetDimension(), DMPlexBuildFromCellList(), 3306 DMPlexInterpolate(), DMPlexBuildCoordinatesFromCellList() 3307 3308 See DMPlexBuildFromCellList() for an example and details about the topology-related parameters. 3309 See DMPlexBuildCoordinatesFromCellList() for details about the geometry-related parameters. 3310 3311 Level: intermediate 3312 3313 .seealso: DMPlexCreateFromCellListParallelPetsc(), DMPlexBuildFromCellList(), DMPlexBuildCoordinatesFromCellList(), DMPlexCreateFromDAG(), DMPlexCreate() 3314 @*/ 3315 PetscErrorCode DMPlexCreateFromCellListPetsc(MPI_Comm comm, PetscInt dim, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, PetscBool interpolate, const PetscInt cells[], PetscInt spaceDim, const PetscReal vertexCoords[], DM *dm) 3316 { 3317 PetscErrorCode ierr; 3318 3319 PetscFunctionBegin; 3320 if (!dim) SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "This is not appropriate for 0-dimensional meshes. Consider either creating the DM using DMPlexCreateFromDAG(), by hand, or using DMSwarm."); 3321 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 3322 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 3323 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 3324 ierr = DMPlexBuildFromCellList(*dm, numCells, numVertices, numCorners, cells);CHKERRQ(ierr); 3325 if (interpolate) { 3326 DM idm; 3327 3328 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 3329 ierr = DMDestroy(dm);CHKERRQ(ierr); 3330 *dm = idm; 3331 } 3332 ierr = DMPlexBuildCoordinatesFromCellList(*dm, spaceDim, vertexCoords);CHKERRQ(ierr); 3333 PetscFunctionReturn(0); 3334 } 3335 3336 /*@ 3337 DMPlexCreateFromCellList - Deprecated, use DMPlexCreateFromCellListPetsc() 3338 3339 Level: deprecated 3340 3341 .seealso: DMPlexCreateFromCellListPetsc() 3342 @*/ 3343 PetscErrorCode DMPlexCreateFromCellList(MPI_Comm comm, PetscInt dim, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, PetscBool interpolate, const int cells[], PetscInt spaceDim, const double vertexCoords[], DM *dm) 3344 { 3345 PetscErrorCode ierr; 3346 PetscInt i; 3347 PetscInt *pintCells; 3348 PetscReal *prealVC; 3349 3350 PetscFunctionBegin; 3351 if (sizeof(int) > sizeof(PetscInt)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Size of int %zd greater than size of PetscInt %zd. Reconfigure PETSc --with-64-bit-indices=1", sizeof(int), sizeof(PetscInt)); 3352 if (sizeof(int) == sizeof(PetscInt)) { 3353 pintCells = (PetscInt *) cells; 3354 } else { 3355 ierr = PetscMalloc1(numCells*numCorners, &pintCells);CHKERRQ(ierr); 3356 for (i = 0; i < numCells*numCorners; i++) { 3357 pintCells[i] = (PetscInt) cells[i]; 3358 } 3359 } 3360 if (sizeof(double) > sizeof(PetscReal)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Size of double %zd greater than size of PetscReal %zd. Reconfigure PETSc --with-precision=<higher precision>.", sizeof(double), sizeof(PetscReal)); 3361 if (sizeof(double) == sizeof(PetscReal)) { 3362 prealVC = (PetscReal *) vertexCoords; 3363 } else { 3364 ierr = PetscMalloc1(numVertices*spaceDim, &prealVC);CHKERRQ(ierr); 3365 for (i = 0; i < numVertices*spaceDim; i++) { 3366 prealVC[i] = (PetscReal) vertexCoords[i]; 3367 } 3368 } 3369 ierr = DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, pintCells, spaceDim, prealVC, dm);CHKERRQ(ierr); 3370 if (sizeof(int) != sizeof(PetscInt)) { 3371 ierr = PetscFree(pintCells);CHKERRQ(ierr); 3372 } 3373 if (sizeof(double) != sizeof(PetscReal)) { 3374 ierr = PetscFree(prealVC);CHKERRQ(ierr); 3375 } 3376 PetscFunctionReturn(0); 3377 } 3378 3379 /*@ 3380 DMPlexCreateFromDAG - This takes as input the adjacency-list representation of the Directed Acyclic Graph (Hasse Diagram) encoding a mesh, and produces a DM 3381 3382 Input Parameters: 3383 + dm - The empty DM object, usually from DMCreate() and DMSetDimension() 3384 . depth - The depth of the DAG 3385 . numPoints - Array of size depth + 1 containing the number of points at each depth 3386 . coneSize - The cone size of each point 3387 . cones - The concatenation of the cone points for each point, the cone list must be oriented correctly for each point 3388 . coneOrientations - The orientation of each cone point 3389 - vertexCoords - An array of numPoints[0]*spacedim numbers representing the coordinates of each vertex, with spacedim the value set via DMSetCoordinateDim() 3390 3391 Output Parameter: 3392 . dm - The DM 3393 3394 Note: Two triangles sharing a face would have input 3395 $ depth = 1, numPoints = [4 2], coneSize = [3 3 0 0 0 0] 3396 $ cones = [2 3 4 3 5 4], coneOrientations = [0 0 0 0 0 0] 3397 $ vertexCoords = [-1.0 0.0 0.0 -1.0 0.0 1.0 1.0 0.0] 3398 $ 3399 which would result in the DMPlex 3400 $ 3401 $ 4 3402 $ / | \ 3403 $ / | \ 3404 $ / | \ 3405 $ 2 0 | 1 5 3406 $ \ | / 3407 $ \ | / 3408 $ \ | / 3409 $ 3 3410 $ 3411 $ Notice that all points are numbered consecutively, unlike DMPlexCreateFromCellListPetsc() 3412 3413 Level: advanced 3414 3415 .seealso: DMPlexCreateFromCellListPetsc(), DMPlexCreate() 3416 @*/ 3417 PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[]) 3418 { 3419 Vec coordinates; 3420 PetscSection coordSection; 3421 PetscScalar *coords; 3422 PetscInt coordSize, firstVertex = -1, pStart = 0, pEnd = 0, p, v, dim, dimEmbed, d, off; 3423 PetscErrorCode ierr; 3424 3425 PetscFunctionBegin; 3426 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3427 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 3428 if (dimEmbed < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Embedding dimension %D cannot be less than intrinsic dimension %d",dimEmbed,dim); 3429 for (d = 0; d <= depth; ++d) pEnd += numPoints[d]; 3430 ierr = DMPlexSetChart(dm, pStart, pEnd);CHKERRQ(ierr); 3431 for (p = pStart; p < pEnd; ++p) { 3432 ierr = DMPlexSetConeSize(dm, p, coneSize[p-pStart]);CHKERRQ(ierr); 3433 if (firstVertex < 0 && !coneSize[p - pStart]) { 3434 firstVertex = p - pStart; 3435 } 3436 } 3437 if (firstVertex < 0 && numPoints[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected %D vertices but could not find any", numPoints[0]); 3438 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 3439 for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) { 3440 ierr = DMPlexSetCone(dm, p, &cones[off]);CHKERRQ(ierr); 3441 ierr = DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);CHKERRQ(ierr); 3442 } 3443 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 3444 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 3445 /* Build coordinates */ 3446 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 3447 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 3448 ierr = PetscSectionSetFieldComponents(coordSection, 0, dimEmbed);CHKERRQ(ierr); 3449 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0]);CHKERRQ(ierr); 3450 for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) { 3451 ierr = PetscSectionSetDof(coordSection, v, dimEmbed);CHKERRQ(ierr); 3452 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed);CHKERRQ(ierr); 3453 } 3454 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 3455 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 3456 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 3457 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 3458 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 3459 ierr = VecSetBlockSize(coordinates, dimEmbed);CHKERRQ(ierr); 3460 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 3461 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 3462 for (v = 0; v < numPoints[0]; ++v) { 3463 PetscInt off; 3464 3465 ierr = PetscSectionGetOffset(coordSection, v+firstVertex, &off);CHKERRQ(ierr); 3466 for (d = 0; d < dimEmbed; ++d) { 3467 coords[off+d] = vertexCoords[v*dimEmbed+d]; 3468 } 3469 } 3470 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 3471 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 3472 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 3473 PetscFunctionReturn(0); 3474 } 3475 3476 /*@C 3477 DMPlexCreateCellVertexFromFile - Create a DMPlex mesh from a simple cell-vertex file. 3478 3479 + comm - The MPI communicator 3480 . filename - Name of the .dat file 3481 - interpolate - Create faces and edges in the mesh 3482 3483 Output Parameter: 3484 . dm - The DM object representing the mesh 3485 3486 Note: The format is the simplest possible: 3487 $ Ne 3488 $ v0 v1 ... vk 3489 $ Nv 3490 $ x y z marker 3491 3492 Level: beginner 3493 3494 .seealso: DMPlexCreateFromFile(), DMPlexCreateMedFromFile(), DMPlexCreateGmsh(), DMPlexCreate() 3495 @*/ 3496 PetscErrorCode DMPlexCreateCellVertexFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 3497 { 3498 DMLabel marker; 3499 PetscViewer viewer; 3500 Vec coordinates; 3501 PetscSection coordSection; 3502 PetscScalar *coords; 3503 char line[PETSC_MAX_PATH_LEN]; 3504 PetscInt dim = 3, cdim = 3, coordSize, v, c, d; 3505 PetscMPIInt rank; 3506 int snum, Nv, Nc; 3507 PetscErrorCode ierr; 3508 3509 PetscFunctionBegin; 3510 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 3511 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 3512 ierr = PetscViewerSetType(viewer, PETSCVIEWERASCII);CHKERRQ(ierr); 3513 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 3514 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 3515 if (!rank) { 3516 ierr = PetscViewerRead(viewer, line, 2, NULL, PETSC_STRING);CHKERRQ(ierr); 3517 snum = sscanf(line, "%d %d", &Nc, &Nv); 3518 if (snum != 2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 3519 } else { 3520 Nc = Nv = 0; 3521 } 3522 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 3523 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 3524 ierr = DMPlexSetChart(*dm, 0, Nc+Nv);CHKERRQ(ierr); 3525 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 3526 ierr = DMSetCoordinateDim(*dm, cdim);CHKERRQ(ierr); 3527 /* Read topology */ 3528 if (!rank) { 3529 PetscInt cone[8], corners = 8; 3530 int vbuf[8], v; 3531 3532 for (c = 0; c < Nc; ++c) {ierr = DMPlexSetConeSize(*dm, c, corners);CHKERRQ(ierr);} 3533 ierr = DMSetUp(*dm);CHKERRQ(ierr); 3534 for (c = 0; c < Nc; ++c) { 3535 ierr = PetscViewerRead(viewer, line, corners, NULL, PETSC_STRING);CHKERRQ(ierr); 3536 snum = sscanf(line, "%d %d %d %d %d %d %d %d", &vbuf[0], &vbuf[1], &vbuf[2], &vbuf[3], &vbuf[4], &vbuf[5], &vbuf[6], &vbuf[7]); 3537 if (snum != corners) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 3538 for (v = 0; v < corners; ++v) cone[v] = vbuf[v] + Nc; 3539 /* Hexahedra are inverted */ 3540 { 3541 PetscInt tmp = cone[1]; 3542 cone[1] = cone[3]; 3543 cone[3] = tmp; 3544 } 3545 ierr = DMPlexSetCone(*dm, c, cone);CHKERRQ(ierr); 3546 } 3547 } 3548 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 3549 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 3550 /* Read coordinates */ 3551 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 3552 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 3553 ierr = PetscSectionSetFieldComponents(coordSection, 0, cdim);CHKERRQ(ierr); 3554 ierr = PetscSectionSetChart(coordSection, Nc, Nc + Nv);CHKERRQ(ierr); 3555 for (v = Nc; v < Nc+Nv; ++v) { 3556 ierr = PetscSectionSetDof(coordSection, v, cdim);CHKERRQ(ierr); 3557 ierr = PetscSectionSetFieldDof(coordSection, v, 0, cdim);CHKERRQ(ierr); 3558 } 3559 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 3560 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 3561 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 3562 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 3563 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 3564 ierr = VecSetBlockSize(coordinates, cdim);CHKERRQ(ierr); 3565 ierr = VecSetType(coordinates, VECSTANDARD);CHKERRQ(ierr); 3566 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 3567 if (!rank) { 3568 double x[3]; 3569 int val; 3570 3571 ierr = DMCreateLabel(*dm, "marker");CHKERRQ(ierr); 3572 ierr = DMGetLabel(*dm, "marker", &marker);CHKERRQ(ierr); 3573 for (v = 0; v < Nv; ++v) { 3574 ierr = PetscViewerRead(viewer, line, 4, NULL, PETSC_STRING);CHKERRQ(ierr); 3575 snum = sscanf(line, "%lg %lg %lg %d", &x[0], &x[1], &x[2], &val); 3576 if (snum != 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 3577 for (d = 0; d < cdim; ++d) coords[v*cdim+d] = x[d]; 3578 if (val) {ierr = DMLabelSetValue(marker, v+Nc, val);CHKERRQ(ierr);} 3579 } 3580 } 3581 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 3582 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 3583 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 3584 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3585 if (interpolate) { 3586 DM idm; 3587 DMLabel bdlabel; 3588 3589 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 3590 ierr = DMDestroy(dm);CHKERRQ(ierr); 3591 *dm = idm; 3592 3593 ierr = DMGetLabel(*dm, "marker", &bdlabel);CHKERRQ(ierr); 3594 ierr = DMPlexMarkBoundaryFaces(*dm, PETSC_DETERMINE, bdlabel);CHKERRQ(ierr); 3595 ierr = DMPlexLabelComplete(*dm, bdlabel);CHKERRQ(ierr); 3596 } 3597 PetscFunctionReturn(0); 3598 } 3599 3600 /*@C 3601 DMPlexCreateFromFile - This takes a filename and produces a DM 3602 3603 Input Parameters: 3604 + comm - The communicator 3605 . filename - A file name 3606 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 3607 3608 Output Parameter: 3609 . dm - The DM 3610 3611 Options Database Keys: 3612 . -dm_plex_create_from_hdf5_xdmf - use the PETSC_VIEWER_HDF5_XDMF format for reading HDF5 3613 3614 Use -dm_plex_create_ prefix to pass options to the internal PetscViewer, e.g. 3615 $ -dm_plex_create_viewer_hdf5_collective 3616 3617 Level: beginner 3618 3619 .seealso: DMPlexCreateFromDAG(), DMPlexCreateFromCellListPetsc(), DMPlexCreate() 3620 @*/ 3621 PetscErrorCode DMPlexCreateFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 3622 { 3623 const char *extGmsh = ".msh"; 3624 const char *extGmsh2 = ".msh2"; 3625 const char *extGmsh4 = ".msh4"; 3626 const char *extCGNS = ".cgns"; 3627 const char *extExodus = ".exo"; 3628 const char *extGenesis = ".gen"; 3629 const char *extFluent = ".cas"; 3630 const char *extHDF5 = ".h5"; 3631 const char *extMed = ".med"; 3632 const char *extPLY = ".ply"; 3633 const char *extEGADS = ".egadslite"; 3634 const char *extCV = ".dat"; 3635 size_t len; 3636 PetscBool isGmsh, isGmsh2, isGmsh4, isCGNS, isExodus, isGenesis, isFluent, isHDF5, isMed, isPLY, isEGADS, isCV; 3637 PetscMPIInt rank; 3638 PetscErrorCode ierr; 3639 3640 PetscFunctionBegin; 3641 PetscValidCharPointer(filename, 2); 3642 PetscValidPointer(dm, 4); 3643 ierr = DMInitializePackage();CHKERRQ(ierr); 3644 ierr = PetscLogEventBegin(DMPLEX_CreateFromFile,0,0,0,0);CHKERRQ(ierr); 3645 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 3646 ierr = PetscStrlen(filename, &len);CHKERRQ(ierr); 3647 if (!len) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Filename must be a valid path"); 3648 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGmsh, 4, &isGmsh);CHKERRQ(ierr); 3649 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extGmsh2, 5, &isGmsh2);CHKERRQ(ierr); 3650 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extGmsh4, 5, &isGmsh4);CHKERRQ(ierr); 3651 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extCGNS, 5, &isCGNS);CHKERRQ(ierr); 3652 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extExodus, 4, &isExodus);CHKERRQ(ierr); 3653 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGenesis, 4, &isGenesis);CHKERRQ(ierr); 3654 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extFluent, 4, &isFluent);CHKERRQ(ierr); 3655 ierr = PetscStrncmp(&filename[PetscMax(0,len-3)], extHDF5, 3, &isHDF5);CHKERRQ(ierr); 3656 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extMed, 4, &isMed);CHKERRQ(ierr); 3657 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extPLY, 4, &isPLY);CHKERRQ(ierr); 3658 ierr = PetscStrncmp(&filename[PetscMax(0,len-10)], extEGADS, 9, &isEGADS);CHKERRQ(ierr); 3659 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extCV, 4, &isCV);CHKERRQ(ierr); 3660 if (isGmsh || isGmsh2 || isGmsh4) { 3661 ierr = DMPlexCreateGmshFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3662 } else if (isCGNS) { 3663 ierr = DMPlexCreateCGNSFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3664 } else if (isExodus || isGenesis) { 3665 ierr = DMPlexCreateExodusFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3666 } else if (isFluent) { 3667 ierr = DMPlexCreateFluentFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3668 } else if (isHDF5) { 3669 PetscBool load_hdf5_xdmf = PETSC_FALSE; 3670 PetscViewer viewer; 3671 3672 /* PETSC_VIEWER_HDF5_XDMF is used if the filename ends with .xdmf.h5, or if -dm_plex_create_from_hdf5_xdmf option is present */ 3673 ierr = PetscStrncmp(&filename[PetscMax(0,len-8)], ".xdmf", 5, &load_hdf5_xdmf);CHKERRQ(ierr); 3674 ierr = PetscOptionsGetBool(NULL, NULL, "-dm_plex_create_from_hdf5_xdmf", &load_hdf5_xdmf, NULL);CHKERRQ(ierr); 3675 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 3676 ierr = PetscViewerSetType(viewer, PETSCVIEWERHDF5);CHKERRQ(ierr); 3677 ierr = PetscViewerSetOptionsPrefix(viewer, "dm_plex_create_");CHKERRQ(ierr); 3678 ierr = PetscViewerSetFromOptions(viewer);CHKERRQ(ierr); 3679 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 3680 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 3681 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 3682 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 3683 if (load_hdf5_xdmf) {ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_XDMF);CHKERRQ(ierr);} 3684 ierr = DMLoad(*dm, viewer);CHKERRQ(ierr); 3685 if (load_hdf5_xdmf) {ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);} 3686 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3687 3688 if (interpolate) { 3689 DM idm; 3690 3691 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 3692 ierr = DMDestroy(dm);CHKERRQ(ierr); 3693 *dm = idm; 3694 } 3695 } else if (isMed) { 3696 ierr = DMPlexCreateMedFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3697 } else if (isPLY) { 3698 ierr = DMPlexCreatePLYFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3699 } else if (isEGADS) { 3700 ierr = DMPlexCreateEGADSFromFile(comm, filename, dm);CHKERRQ(ierr); 3701 if (!interpolate) { 3702 DM udm; 3703 3704 ierr = DMPlexUninterpolate(*dm, &udm);CHKERRQ(ierr); 3705 ierr = DMDestroy(dm);CHKERRQ(ierr); 3706 *dm = udm; 3707 } 3708 } else if (isCV) { 3709 ierr = DMPlexCreateCellVertexFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3710 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot load file %s: unrecognized extension", filename); 3711 ierr = PetscLogEventEnd(DMPLEX_CreateFromFile,0,0,0,0);CHKERRQ(ierr); 3712 PetscFunctionReturn(0); 3713 } 3714 3715 /*@ 3716 DMPlexCreateReferenceCellByType - Create a DMPLEX with the appropriate FEM reference cell 3717 3718 Collective 3719 3720 Input Parameters: 3721 + comm - The communicator 3722 - ct - The cell type of the reference cell 3723 3724 Output Parameter: 3725 . refdm - The reference cell 3726 3727 Options Database Keys: 3728 . -dm_plex_ref_type <ct> - Specify the celltyoe for the reference cell 3729 3730 Level: intermediate 3731 3732 .seealso: DMPlexCreateReferenceCell(), DMPlexCreateBoxMesh() 3733 @*/ 3734 PetscErrorCode DMPlexCreateReferenceCellByType(MPI_Comm comm, DMPolytopeType ct, DM *refdm) 3735 { 3736 DM rdm; 3737 Vec coords; 3738 PetscErrorCode ierr; 3739 3740 PetscFunctionBegin; 3741 ierr = PetscOptionsGetEnum(NULL, NULL, "-dm_plex_ref_type", DMPolytopeTypes, (PetscEnum *) &ct, NULL);CHKERRQ(ierr); 3742 ierr = DMCreate(comm, &rdm);CHKERRQ(ierr); 3743 ierr = DMSetType(rdm, DMPLEX);CHKERRQ(ierr); 3744 switch (ct) { 3745 case DM_POLYTOPE_POINT: 3746 { 3747 PetscInt numPoints[1] = {1}; 3748 PetscInt coneSize[1] = {0}; 3749 PetscInt cones[1] = {0}; 3750 PetscInt coneOrientations[1] = {0}; 3751 PetscScalar vertexCoords[1] = {0.0}; 3752 3753 ierr = DMSetDimension(rdm, 0);CHKERRQ(ierr); 3754 ierr = DMPlexCreateFromDAG(rdm, 0, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3755 } 3756 break; 3757 case DM_POLYTOPE_SEGMENT: 3758 { 3759 PetscInt numPoints[2] = {2, 1}; 3760 PetscInt coneSize[3] = {2, 0, 0}; 3761 PetscInt cones[2] = {1, 2}; 3762 PetscInt coneOrientations[2] = {0, 0}; 3763 PetscScalar vertexCoords[2] = {-1.0, 1.0}; 3764 3765 ierr = DMSetDimension(rdm, 1);CHKERRQ(ierr); 3766 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3767 } 3768 break; 3769 case DM_POLYTOPE_TRIANGLE: 3770 { 3771 PetscInt numPoints[2] = {3, 1}; 3772 PetscInt coneSize[4] = {3, 0, 0, 0}; 3773 PetscInt cones[3] = {1, 2, 3}; 3774 PetscInt coneOrientations[3] = {0, 0, 0}; 3775 PetscScalar vertexCoords[6] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0}; 3776 3777 ierr = DMSetDimension(rdm, 2);CHKERRQ(ierr); 3778 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3779 } 3780 break; 3781 case DM_POLYTOPE_QUADRILATERAL: 3782 { 3783 PetscInt numPoints[2] = {4, 1}; 3784 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 3785 PetscInt cones[4] = {1, 2, 3, 4}; 3786 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 3787 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0}; 3788 3789 ierr = DMSetDimension(rdm, 2);CHKERRQ(ierr); 3790 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3791 } 3792 break; 3793 case DM_POLYTOPE_SEG_PRISM_TENSOR: 3794 { 3795 PetscInt numPoints[2] = {4, 1}; 3796 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 3797 PetscInt cones[4] = {1, 2, 3, 4}; 3798 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 3799 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0}; 3800 3801 ierr = DMSetDimension(rdm, 2);CHKERRQ(ierr); 3802 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3803 } 3804 break; 3805 case DM_POLYTOPE_TETRAHEDRON: 3806 { 3807 PetscInt numPoints[2] = {4, 1}; 3808 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 3809 PetscInt cones[4] = {1, 3, 2, 4}; 3810 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 3811 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}; 3812 3813 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 3814 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3815 } 3816 break; 3817 case DM_POLYTOPE_HEXAHEDRON: 3818 { 3819 PetscInt numPoints[2] = {8, 1}; 3820 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0}; 3821 PetscInt cones[8] = {1, 4, 3, 2, 5, 6, 7, 8}; 3822 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 3823 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, 3824 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0}; 3825 3826 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 3827 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3828 } 3829 break; 3830 case DM_POLYTOPE_TRI_PRISM: 3831 { 3832 PetscInt numPoints[2] = {6, 1}; 3833 PetscInt coneSize[7] = {6, 0, 0, 0, 0, 0, 0}; 3834 PetscInt cones[6] = {1, 3, 2, 4, 5, 6}; 3835 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 3836 PetscScalar vertexCoords[18] = {-1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 3837 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0}; 3838 3839 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 3840 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3841 } 3842 break; 3843 case DM_POLYTOPE_TRI_PRISM_TENSOR: 3844 { 3845 PetscInt numPoints[2] = {6, 1}; 3846 PetscInt coneSize[7] = {6, 0, 0, 0, 0, 0, 0}; 3847 PetscInt cones[6] = {1, 2, 3, 4, 5, 6}; 3848 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 3849 PetscScalar vertexCoords[18] = {-1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 3850 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0}; 3851 3852 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 3853 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3854 } 3855 break; 3856 case DM_POLYTOPE_QUAD_PRISM_TENSOR: 3857 { 3858 PetscInt numPoints[2] = {8, 1}; 3859 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0}; 3860 PetscInt cones[8] = {1, 2, 3, 4, 5, 6, 7, 8}; 3861 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 3862 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, 3863 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0}; 3864 3865 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 3866 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3867 } 3868 break; 3869 case DM_POLYTOPE_PYRAMID: 3870 { 3871 PetscInt numPoints[2] = {5, 1}; 3872 PetscInt coneSize[6] = {5, 0, 0, 0, 0, 0}; 3873 PetscInt cones[5] = {1, 4, 3, 2, 5}; 3874 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 3875 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, 3876 0.0, 0.0, 1.0}; 3877 3878 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 3879 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3880 } 3881 break; 3882 default: SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Cannot create reference cell for cell type %s", DMPolytopeTypes[ct]); 3883 } 3884 { 3885 PetscInt Nv, v; 3886 3887 /* Must create the celltype label here so that we do not automatically try to compute the types */ 3888 ierr = DMCreateLabel(rdm, "celltype");CHKERRQ(ierr); 3889 ierr = DMPlexSetCellType(rdm, 0, ct);CHKERRQ(ierr); 3890 ierr = DMPlexGetChart(rdm, NULL, &Nv);CHKERRQ(ierr); 3891 for (v = 1; v < Nv; ++v) {ierr = DMPlexSetCellType(rdm, v, DM_POLYTOPE_POINT);CHKERRQ(ierr);} 3892 } 3893 ierr = DMPlexInterpolate(rdm, refdm);CHKERRQ(ierr); 3894 if (rdm->coordinateDM) { 3895 DM ncdm; 3896 PetscSection cs; 3897 PetscInt pEnd = -1; 3898 3899 ierr = DMGetLocalSection(rdm->coordinateDM, &cs);CHKERRQ(ierr); 3900 if (cs) {ierr = PetscSectionGetChart(cs, NULL, &pEnd);CHKERRQ(ierr);} 3901 if (pEnd >= 0) { 3902 ierr = DMClone(*refdm, &ncdm);CHKERRQ(ierr); 3903 ierr = DMCopyDisc(rdm->coordinateDM, ncdm);CHKERRQ(ierr); 3904 ierr = DMSetLocalSection(ncdm, cs);CHKERRQ(ierr); 3905 ierr = DMSetCoordinateDM(*refdm, ncdm);CHKERRQ(ierr); 3906 ierr = DMDestroy(&ncdm);CHKERRQ(ierr); 3907 } 3908 } 3909 ierr = DMGetCoordinatesLocal(rdm, &coords);CHKERRQ(ierr); 3910 if (coords) { 3911 ierr = DMSetCoordinatesLocal(*refdm, coords);CHKERRQ(ierr); 3912 } else { 3913 ierr = DMGetCoordinates(rdm, &coords);CHKERRQ(ierr); 3914 if (coords) {ierr = DMSetCoordinates(*refdm, coords);CHKERRQ(ierr);} 3915 } 3916 ierr = DMDestroy(&rdm);CHKERRQ(ierr); 3917 PetscFunctionReturn(0); 3918 } 3919 3920 /*@ 3921 DMPlexCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell 3922 3923 Collective 3924 3925 Input Parameters: 3926 + comm - The communicator 3927 . dim - The spatial dimension 3928 - simplex - Flag for simplex, otherwise use a tensor-product cell 3929 3930 Output Parameter: 3931 . refdm - The reference cell 3932 3933 Level: intermediate 3934 3935 .seealso: DMPlexCreateReferenceCellByType(), DMPlexCreateBoxMesh() 3936 @*/ 3937 PetscErrorCode DMPlexCreateReferenceCell(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *refdm) 3938 { 3939 PetscErrorCode ierr; 3940 3941 PetscFunctionBeginUser; 3942 switch (dim) { 3943 case 0: ierr = DMPlexCreateReferenceCellByType(comm, DM_POLYTOPE_POINT, refdm);CHKERRQ(ierr);break; 3944 case 1: ierr = DMPlexCreateReferenceCellByType(comm, DM_POLYTOPE_SEGMENT, refdm);CHKERRQ(ierr);break; 3945 case 2: 3946 if (simplex) {ierr = DMPlexCreateReferenceCellByType(comm, DM_POLYTOPE_TRIANGLE, refdm);CHKERRQ(ierr);} 3947 else {ierr = DMPlexCreateReferenceCellByType(comm, DM_POLYTOPE_QUADRILATERAL, refdm);CHKERRQ(ierr);} 3948 break; 3949 case 3: 3950 if (simplex) {ierr = DMPlexCreateReferenceCellByType(comm, DM_POLYTOPE_TETRAHEDRON, refdm);CHKERRQ(ierr);} 3951 else {ierr = DMPlexCreateReferenceCellByType(comm, DM_POLYTOPE_HEXAHEDRON, refdm);CHKERRQ(ierr);} 3952 break; 3953 default: 3954 SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Cannot create reference cell for dimension %D", dim); 3955 } 3956 PetscFunctionReturn(0); 3957 } 3958