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 = PetscSFGetRootRanks(sf, &njranks, &jranks, NULL, NULL, NULL);CHKERRQ(ierr); 2604 ierr = PetscSFGetLeafRanks(sf, &niranks, &iranks, NULL, NULL);CHKERRQ(ierr); 2605 ierr = PetscMalloc1(njranks + niranks + 1, &data->neighbors);CHKERRQ(ierr); 2606 ierr = PetscArraycpy(data->neighbors + 1, jranks, njranks);CHKERRQ(ierr); 2607 ierr = PetscArraycpy(data->neighbors + njranks + 1, iranks, niranks);CHKERRQ(ierr); 2608 n = njranks + niranks; 2609 ierr = PetscSortRemoveDupsMPIInt(&n, data->neighbors + 1);CHKERRQ(ierr); 2610 /* The following cast should never fail: can't have more neighbors than PETSC_MPI_INT_MAX */ 2611 ierr = PetscMPIIntCast(n, data->neighbors);CHKERRQ(ierr); 2612 } 2613 if (nranks) *nranks = data->neighbors[0]; 2614 if (ranks) { 2615 if (data->neighbors[0]) *ranks = data->neighbors + 1; 2616 else *ranks = NULL; 2617 } 2618 PetscFunctionReturn(0); 2619 } 2620 2621 static PetscErrorCode DMInitialize_Plex(DM dm) 2622 { 2623 PetscErrorCode ierr; 2624 2625 PetscFunctionBegin; 2626 dm->ops->view = DMView_Plex; 2627 dm->ops->load = DMLoad_Plex; 2628 dm->ops->setfromoptions = DMSetFromOptions_Plex; 2629 dm->ops->clone = DMClone_Plex; 2630 dm->ops->setup = DMSetUp_Plex; 2631 dm->ops->createlocalsection = DMCreateLocalSection_Plex; 2632 dm->ops->createdefaultconstraints = DMCreateDefaultConstraints_Plex; 2633 dm->ops->createglobalvector = DMCreateGlobalVector_Plex; 2634 dm->ops->createlocalvector = DMCreateLocalVector_Plex; 2635 dm->ops->getlocaltoglobalmapping = NULL; 2636 dm->ops->createfieldis = NULL; 2637 dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex; 2638 dm->ops->createcoordinatefield = DMCreateCoordinateField_Plex; 2639 dm->ops->getcoloring = NULL; 2640 dm->ops->creatematrix = DMCreateMatrix_Plex; 2641 dm->ops->createinterpolation = DMCreateInterpolation_Plex; 2642 dm->ops->createmassmatrix = DMCreateMassMatrix_Plex; 2643 dm->ops->createinjection = DMCreateInjection_Plex; 2644 dm->ops->refine = DMRefine_Plex; 2645 dm->ops->coarsen = DMCoarsen_Plex; 2646 dm->ops->refinehierarchy = DMRefineHierarchy_Plex; 2647 dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Plex; 2648 dm->ops->adaptlabel = DMAdaptLabel_Plex; 2649 dm->ops->adaptmetric = DMAdaptMetric_Plex; 2650 dm->ops->globaltolocalbegin = NULL; 2651 dm->ops->globaltolocalend = NULL; 2652 dm->ops->localtoglobalbegin = NULL; 2653 dm->ops->localtoglobalend = NULL; 2654 dm->ops->destroy = DMDestroy_Plex; 2655 dm->ops->createsubdm = DMCreateSubDM_Plex; 2656 dm->ops->createsuperdm = DMCreateSuperDM_Plex; 2657 dm->ops->getdimpoints = DMGetDimPoints_Plex; 2658 dm->ops->locatepoints = DMLocatePoints_Plex; 2659 dm->ops->projectfunctionlocal = DMProjectFunctionLocal_Plex; 2660 dm->ops->projectfunctionlabellocal = DMProjectFunctionLabelLocal_Plex; 2661 dm->ops->projectfieldlocal = DMProjectFieldLocal_Plex; 2662 dm->ops->projectfieldlabellocal = DMProjectFieldLabelLocal_Plex; 2663 dm->ops->projectbdfieldlabellocal = DMProjectBdFieldLabelLocal_Plex; 2664 dm->ops->computel2diff = DMComputeL2Diff_Plex; 2665 dm->ops->computel2gradientdiff = DMComputeL2GradientDiff_Plex; 2666 dm->ops->computel2fielddiff = DMComputeL2FieldDiff_Plex; 2667 dm->ops->getneighbors = DMGetNeighbors_Plex; 2668 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMPlexInsertBoundaryValues_C",DMPlexInsertBoundaryValues_Plex);CHKERRQ(ierr); 2669 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMPlexInsertTimeDerviativeBoundaryValues_C",DMPlexInsertTimeDerivativeBoundaryValues_Plex);CHKERRQ(ierr); 2670 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMSetUpGLVisViewer_C",DMSetUpGLVisViewer_Plex);CHKERRQ(ierr); 2671 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMCreateNeumannOverlap_C",DMCreateNeumannOverlap_Plex);CHKERRQ(ierr); 2672 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMPlexGetOverlap_C",DMPlexGetOverlap_Plex);CHKERRQ(ierr); 2673 PetscFunctionReturn(0); 2674 } 2675 2676 PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm) 2677 { 2678 DM_Plex *mesh = (DM_Plex *) dm->data; 2679 PetscErrorCode ierr; 2680 2681 PetscFunctionBegin; 2682 mesh->refct++; 2683 (*newdm)->data = mesh; 2684 ierr = PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);CHKERRQ(ierr); 2685 ierr = DMInitialize_Plex(*newdm);CHKERRQ(ierr); 2686 PetscFunctionReturn(0); 2687 } 2688 2689 /*MC 2690 DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram. 2691 In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is 2692 specified by a PetscSection object. Ownership in the global representation is determined by 2693 ownership of the underlying DMPlex points. This is specified by another PetscSection object. 2694 2695 Options Database Keys: 2696 + -dm_refine_pre - Refine mesh before distribution 2697 + -dm_refine_uniform_pre - Choose uniform or generator-based refinement 2698 + -dm_refine_volume_limit_pre - Cell volume limit after pre-refinement using generator 2699 . -dm_distribute - Distribute mesh across processes 2700 . -dm_distribute_overlap - Number of cells to overlap for distribution 2701 . -dm_refine - Refine mesh after distribution 2702 . -dm_plex_hash_location - Use grid hashing for point location 2703 . -dm_plex_partition_balance - Attempt to evenly divide points on partition boundary between processes 2704 . -dm_plex_remesh_bd - Allow changes to the boundary on remeshing 2705 . -dm_plex_max_projection_height - Maxmimum mesh point height used to project locally 2706 . -dm_plex_regular_refinement - Use special nested projection algorithm for regular refinement 2707 . -dm_plex_check_all - Perform all shecks below 2708 . -dm_plex_check_symmetry - Check that the adjacency information in the mesh is symmetric 2709 . -dm_plex_check_skeleton <celltype> - Check that each cell has the correct number of vertices 2710 . -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 2711 . -dm_plex_check_geometry - Check that cells have positive volume 2712 . -dm_view :mesh.tex:ascii_latex - View the mesh in LaTeX/TikZ 2713 . -dm_plex_view_scale <num> - Scale the TikZ 2714 - -dm_plex_print_fem <num> - View FEM assembly information, such as element vectors and matrices 2715 2716 2717 Level: intermediate 2718 2719 .seealso: DMType, DMPlexCreate(), DMCreate(), DMSetType() 2720 M*/ 2721 2722 PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm) 2723 { 2724 DM_Plex *mesh; 2725 PetscInt unit; 2726 PetscErrorCode ierr; 2727 2728 PetscFunctionBegin; 2729 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2730 ierr = PetscNewLog(dm,&mesh);CHKERRQ(ierr); 2731 dm->dim = 0; 2732 dm->data = mesh; 2733 2734 mesh->refct = 1; 2735 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);CHKERRQ(ierr); 2736 mesh->maxConeSize = 0; 2737 mesh->cones = NULL; 2738 mesh->coneOrientations = NULL; 2739 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);CHKERRQ(ierr); 2740 mesh->maxSupportSize = 0; 2741 mesh->supports = NULL; 2742 mesh->refinementUniform = PETSC_TRUE; 2743 mesh->refinementLimit = -1.0; 2744 mesh->interpolated = DMPLEX_INTERPOLATED_INVALID; 2745 mesh->interpolatedCollective = DMPLEX_INTERPOLATED_INVALID; 2746 2747 mesh->facesTmp = NULL; 2748 2749 mesh->tetgenOpts = NULL; 2750 mesh->triangleOpts = NULL; 2751 ierr = PetscPartitionerCreate(PetscObjectComm((PetscObject)dm), &mesh->partitioner);CHKERRQ(ierr); 2752 mesh->remeshBd = PETSC_FALSE; 2753 2754 mesh->subpointMap = NULL; 2755 2756 for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0; 2757 2758 mesh->regularRefinement = PETSC_FALSE; 2759 mesh->depthState = -1; 2760 mesh->celltypeState = -1; 2761 mesh->globalVertexNumbers = NULL; 2762 mesh->globalCellNumbers = NULL; 2763 mesh->anchorSection = NULL; 2764 mesh->anchorIS = NULL; 2765 mesh->createanchors = NULL; 2766 mesh->computeanchormatrix = NULL; 2767 mesh->parentSection = NULL; 2768 mesh->parents = NULL; 2769 mesh->childIDs = NULL; 2770 mesh->childSection = NULL; 2771 mesh->children = NULL; 2772 mesh->referenceTree = NULL; 2773 mesh->getchildsymmetry = NULL; 2774 mesh->vtkCellHeight = 0; 2775 mesh->useAnchors = PETSC_FALSE; 2776 2777 mesh->maxProjectionHeight = 0; 2778 2779 mesh->neighbors = NULL; 2780 2781 mesh->printSetValues = PETSC_FALSE; 2782 mesh->printFEM = 0; 2783 mesh->printTol = 1.0e-10; 2784 2785 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 2786 PetscFunctionReturn(0); 2787 } 2788 2789 /*@ 2790 DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram. 2791 2792 Collective 2793 2794 Input Parameter: 2795 . comm - The communicator for the DMPlex object 2796 2797 Output Parameter: 2798 . mesh - The DMPlex object 2799 2800 Level: beginner 2801 2802 @*/ 2803 PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh) 2804 { 2805 PetscErrorCode ierr; 2806 2807 PetscFunctionBegin; 2808 PetscValidPointer(mesh,2); 2809 ierr = DMCreate(comm, mesh);CHKERRQ(ierr); 2810 ierr = DMSetType(*mesh, DMPLEX);CHKERRQ(ierr); 2811 PetscFunctionReturn(0); 2812 } 2813 2814 /*@C 2815 DMPlexBuildFromCellListParallel - Build distributed DMPLEX topology from a list of vertices for each cell (common mesh generator output) 2816 2817 Input Parameters: 2818 + dm - The DM 2819 . numCells - The number of cells owned by this process 2820 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE 2821 . NVertices - The global number of vertices, or PETSC_DECIDE 2822 . numCorners - The number of vertices for each cell 2823 - cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 2824 2825 Output Parameter: 2826 . vertexSF - (Optional) SF describing complete vertex ownership 2827 2828 Notes: 2829 Two triangles sharing a face 2830 $ 2831 $ 2 2832 $ / | \ 2833 $ / | \ 2834 $ / | \ 2835 $ 0 0 | 1 3 2836 $ \ | / 2837 $ \ | / 2838 $ \ | / 2839 $ 1 2840 would have input 2841 $ numCells = 2, numVertices = 4 2842 $ cells = [0 1 2 1 3 2] 2843 $ 2844 which would result in the DMPlex 2845 $ 2846 $ 4 2847 $ / | \ 2848 $ / | \ 2849 $ / | \ 2850 $ 2 0 | 1 5 2851 $ \ | / 2852 $ \ | / 2853 $ \ | / 2854 $ 3 2855 2856 Vertices are implicitly numbered consecutively 0,...,NVertices. 2857 Each rank owns a chunk of numVertices consecutive vertices. 2858 If numVertices is PETSC_DECIDE, PETSc will distribute them as evenly as possible using PetscLayout. 2859 If both NVertices and numVertices are PETSC_DECIDE, NVertices is computed by PETSc as the maximum vertex index in cells + 1. 2860 If only NVertices is PETSC_DECIDE, it is computed as the sum of numVertices over all ranks. 2861 2862 The cell distribution is arbitrary non-overlapping, independent of the vertex distribution. 2863 2864 Not currently supported in Fortran. 2865 2866 Level: advanced 2867 2868 .seealso: DMPlexBuildFromCellList(), DMPlexCreateFromCellListParallelPetsc(), DMPlexBuildCoordinatesFromCellListParallel() 2869 @*/ 2870 PetscErrorCode DMPlexBuildFromCellListParallel(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt NVertices, PetscInt numCorners, const PetscInt cells[], PetscSF *vertexSF) 2871 { 2872 PetscSF sfPoint, sfVert; 2873 PetscLayout vLayout; 2874 PetscSFNode *vertexLocal, *vertexOwner, *remoteVertex; 2875 PetscInt numVerticesAdj, *verticesAdj, numVerticesGhost = 0, *localVertex, *cones, c, p, v, g, dim; 2876 PetscMPIInt rank, size; 2877 PetscErrorCode ierr; 2878 2879 PetscFunctionBegin; 2880 PetscValidLogicalCollectiveInt(dm,NVertices,4); 2881 ierr = PetscLogEventBegin(DMPLEX_BuildFromCellList,dm,0,0,0);CHKERRQ(ierr); 2882 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRMPI(ierr); 2883 ierr = MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);CHKERRMPI(ierr); 2884 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2885 /* Get/check global number of vertices */ 2886 { 2887 PetscInt NVerticesInCells, i; 2888 const PetscInt len = numCells * numCorners; 2889 2890 /* NVerticesInCells = max(cells) + 1 */ 2891 NVerticesInCells = PETSC_MIN_INT; 2892 for (i=0; i<len; i++) if (cells[i] > NVerticesInCells) NVerticesInCells = cells[i]; 2893 ++NVerticesInCells; 2894 ierr = MPI_Allreduce(MPI_IN_PLACE, &NVerticesInCells, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject) dm));CHKERRMPI(ierr); 2895 2896 if (numVertices == PETSC_DECIDE && NVertices == PETSC_DECIDE) NVertices = NVerticesInCells; 2897 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); 2898 } 2899 /* Partition vertices */ 2900 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) dm), &vLayout);CHKERRQ(ierr); 2901 ierr = PetscLayoutSetLocalSize(vLayout, numVertices);CHKERRQ(ierr); 2902 ierr = PetscLayoutSetSize(vLayout, NVertices);CHKERRQ(ierr); 2903 ierr = PetscLayoutSetBlockSize(vLayout, 1);CHKERRQ(ierr); 2904 ierr = PetscLayoutSetUp(vLayout);CHKERRQ(ierr); 2905 /* Count locally unique vertices */ 2906 { 2907 PetscHSetI vhash; 2908 PetscInt off = 0; 2909 2910 ierr = PetscHSetICreate(&vhash);CHKERRQ(ierr); 2911 for (c = 0; c < numCells; ++c) { 2912 for (p = 0; p < numCorners; ++p) { 2913 ierr = PetscHSetIAdd(vhash, cells[c*numCorners+p]);CHKERRQ(ierr); 2914 } 2915 } 2916 ierr = PetscHSetIGetSize(vhash, &numVerticesAdj);CHKERRQ(ierr); 2917 ierr = PetscMalloc1(numVerticesAdj, &verticesAdj);CHKERRQ(ierr); 2918 ierr = PetscHSetIGetElems(vhash, &off, verticesAdj);CHKERRQ(ierr); 2919 ierr = PetscHSetIDestroy(&vhash);CHKERRQ(ierr); 2920 if (off != numVerticesAdj) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid number of local vertices %D should be %D", off, numVerticesAdj); 2921 } 2922 ierr = PetscSortInt(numVerticesAdj, verticesAdj);CHKERRQ(ierr); 2923 /* Create cones */ 2924 ierr = DMPlexSetChart(dm, 0, numCells+numVerticesAdj);CHKERRQ(ierr); 2925 for (c = 0; c < numCells; ++c) {ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr);} 2926 ierr = DMSetUp(dm);CHKERRQ(ierr); 2927 ierr = DMPlexGetCones(dm,&cones);CHKERRQ(ierr); 2928 for (c = 0; c < numCells; ++c) { 2929 for (p = 0; p < numCorners; ++p) { 2930 const PetscInt gv = cells[c*numCorners+p]; 2931 PetscInt lv; 2932 2933 /* Positions within verticesAdj form 0-based local vertex numbering; 2934 we need to shift it by numCells to get correct DAG points (cells go first) */ 2935 ierr = PetscFindInt(gv, numVerticesAdj, verticesAdj, &lv);CHKERRQ(ierr); 2936 if (lv < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not find global vertex %D in local connectivity", gv); 2937 cones[c*numCorners+p] = lv+numCells; 2938 } 2939 } 2940 /* Create SF for vertices */ 2941 ierr = PetscSFCreate(PetscObjectComm((PetscObject)dm), &sfVert);CHKERRQ(ierr); 2942 ierr = PetscObjectSetName((PetscObject) sfVert, "Vertex Ownership SF");CHKERRQ(ierr); 2943 ierr = PetscSFSetFromOptions(sfVert);CHKERRQ(ierr); 2944 ierr = PetscSFSetGraphLayout(sfVert, vLayout, numVerticesAdj, NULL, PETSC_OWN_POINTER, verticesAdj);CHKERRQ(ierr); 2945 ierr = PetscFree(verticesAdj);CHKERRQ(ierr); 2946 /* Build pointSF */ 2947 ierr = PetscMalloc2(numVerticesAdj, &vertexLocal, numVertices, &vertexOwner);CHKERRQ(ierr); 2948 for (v = 0; v < numVerticesAdj; ++v) {vertexLocal[v].index = v+numCells; vertexLocal[v].rank = rank;} 2949 for (v = 0; v < numVertices; ++v) {vertexOwner[v].index = -1; vertexOwner[v].rank = -1;} 2950 ierr = PetscSFReduceBegin(sfVert, MPIU_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);CHKERRQ(ierr); 2951 ierr = PetscSFReduceEnd(sfVert, MPIU_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);CHKERRQ(ierr); 2952 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); 2953 ierr = PetscSFBcastBegin(sfVert, MPIU_2INT, vertexOwner, vertexLocal);CHKERRQ(ierr); 2954 ierr = PetscSFBcastEnd(sfVert, MPIU_2INT, vertexOwner, vertexLocal);CHKERRQ(ierr); 2955 for (v = 0; v < numVerticesAdj; ++v) if (vertexLocal[v].rank != rank) ++numVerticesGhost; 2956 ierr = PetscMalloc1(numVerticesGhost, &localVertex);CHKERRQ(ierr); 2957 ierr = PetscMalloc1(numVerticesGhost, &remoteVertex);CHKERRQ(ierr); 2958 for (v = 0, g = 0; v < numVerticesAdj; ++v) { 2959 if (vertexLocal[v].rank != rank) { 2960 localVertex[g] = v+numCells; 2961 remoteVertex[g].index = vertexLocal[v].index; 2962 remoteVertex[g].rank = vertexLocal[v].rank; 2963 ++g; 2964 } 2965 } 2966 ierr = PetscFree2(vertexLocal, vertexOwner);CHKERRQ(ierr); 2967 if (g != numVerticesGhost) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of vertex ghosts %D should be %D", g, numVerticesGhost); 2968 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 2969 ierr = PetscObjectSetName((PetscObject) sfPoint, "point SF");CHKERRQ(ierr); 2970 ierr = PetscSFSetGraph(sfPoint, numCells+numVerticesAdj, numVerticesGhost, localVertex, PETSC_OWN_POINTER, remoteVertex, PETSC_OWN_POINTER);CHKERRQ(ierr); 2971 ierr = PetscLayoutDestroy(&vLayout);CHKERRQ(ierr); 2972 /* Fill in the rest of the topology structure */ 2973 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 2974 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 2975 ierr = PetscLogEventEnd(DMPLEX_BuildFromCellList,dm,0,0,0);CHKERRQ(ierr); 2976 if (vertexSF) *vertexSF = sfVert; 2977 else {ierr = PetscSFDestroy(&sfVert);CHKERRQ(ierr);} 2978 PetscFunctionReturn(0); 2979 } 2980 2981 /*@C 2982 DMPlexBuildCoordinatesFromCellListParallel - Build DM coordinates from a list of coordinates for each owned vertex (common mesh generator output) 2983 2984 Input Parameters: 2985 + dm - The DM 2986 . spaceDim - The spatial dimension used for coordinates 2987 . sfVert - SF describing complete vertex ownership 2988 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 2989 2990 Level: advanced 2991 2992 Notes: 2993 Not currently supported in Fortran. 2994 2995 .seealso: DMPlexBuildCoordinatesFromCellList(), DMPlexCreateFromCellListParallelPetsc(), DMPlexBuildFromCellListParallel() 2996 @*/ 2997 PetscErrorCode DMPlexBuildCoordinatesFromCellListParallel(DM dm, PetscInt spaceDim, PetscSF sfVert, const PetscReal vertexCoords[]) 2998 { 2999 PetscSection coordSection; 3000 Vec coordinates; 3001 PetscScalar *coords; 3002 PetscInt numVertices, numVerticesAdj, coordSize, v, vStart, vEnd; 3003 PetscErrorCode ierr; 3004 3005 PetscFunctionBegin; 3006 ierr = PetscLogEventBegin(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0);CHKERRQ(ierr); 3007 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 3008 if (vStart < 0 || vEnd < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM is not set up properly. DMPlexBuildFromCellList() should be called first."); 3009 ierr = DMSetCoordinateDim(dm, spaceDim);CHKERRQ(ierr); 3010 ierr = PetscSFGetGraph(sfVert, &numVertices, &numVerticesAdj, NULL, NULL);CHKERRQ(ierr); 3011 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); 3012 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 3013 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 3014 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 3015 ierr = PetscSectionSetChart(coordSection, vStart, vEnd);CHKERRQ(ierr); 3016 for (v = vStart; v < vEnd; ++v) { 3017 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 3018 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 3019 } 3020 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 3021 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 3022 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 3023 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 3024 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 3025 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 3026 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 3027 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 3028 { 3029 MPI_Datatype coordtype; 3030 3031 /* Need a temp buffer for coords if we have complex/single */ 3032 ierr = MPI_Type_contiguous(spaceDim, MPIU_SCALAR, &coordtype);CHKERRMPI(ierr); 3033 ierr = MPI_Type_commit(&coordtype);CHKERRMPI(ierr); 3034 #if defined(PETSC_USE_COMPLEX) 3035 { 3036 PetscScalar *svertexCoords; 3037 PetscInt i; 3038 ierr = PetscMalloc1(numVertices*spaceDim,&svertexCoords);CHKERRQ(ierr); 3039 for (i=0; i<numVertices*spaceDim; i++) svertexCoords[i] = vertexCoords[i]; 3040 ierr = PetscSFBcastBegin(sfVert, coordtype, svertexCoords, coords);CHKERRQ(ierr); 3041 ierr = PetscSFBcastEnd(sfVert, coordtype, svertexCoords, coords);CHKERRQ(ierr); 3042 ierr = PetscFree(svertexCoords);CHKERRQ(ierr); 3043 } 3044 #else 3045 ierr = PetscSFBcastBegin(sfVert, coordtype, vertexCoords, coords);CHKERRQ(ierr); 3046 ierr = PetscSFBcastEnd(sfVert, coordtype, vertexCoords, coords);CHKERRQ(ierr); 3047 #endif 3048 ierr = MPI_Type_free(&coordtype);CHKERRMPI(ierr); 3049 } 3050 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 3051 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 3052 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 3053 ierr = PetscLogEventEnd(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0);CHKERRQ(ierr); 3054 PetscFunctionReturn(0); 3055 } 3056 3057 /*@ 3058 DMPlexCreateFromCellListParallelPetsc - Create distributed DMPLEX from a list of vertices for each cell (common mesh generator output) 3059 3060 Input Parameters: 3061 + comm - The communicator 3062 . dim - The topological dimension of the mesh 3063 . numCells - The number of cells owned by this process 3064 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE 3065 . NVertices - The global number of vertices, or PETSC_DECIDE 3066 . numCorners - The number of vertices for each cell 3067 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 3068 . cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 3069 . spaceDim - The spatial dimension used for coordinates 3070 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 3071 3072 Output Parameter: 3073 + dm - The DM 3074 - vertexSF - (Optional) SF describing complete vertex ownership 3075 3076 Notes: 3077 This function is just a convenient sequence of DMCreate(), DMSetType(), DMSetDimension(), 3078 DMPlexBuildFromCellListParallel(), DMPlexInterpolate(), DMPlexBuildCoordinatesFromCellListParallel() 3079 3080 See DMPlexBuildFromCellListParallel() for an example and details about the topology-related parameters. 3081 See DMPlexBuildCoordinatesFromCellListParallel() for details about the geometry-related parameters. 3082 3083 Level: intermediate 3084 3085 .seealso: DMPlexCreateFromCellListPetsc(), DMPlexBuildFromCellListParallel(), DMPlexBuildCoordinatesFromCellListParallel(), DMPlexCreateFromDAG(), DMPlexCreate() 3086 @*/ 3087 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) 3088 { 3089 PetscSF sfVert; 3090 PetscErrorCode ierr; 3091 3092 PetscFunctionBegin; 3093 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 3094 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 3095 PetscValidLogicalCollectiveInt(*dm, dim, 2); 3096 PetscValidLogicalCollectiveInt(*dm, spaceDim, 8); 3097 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 3098 ierr = DMPlexBuildFromCellListParallel(*dm, numCells, numVertices, NVertices, numCorners, cells, &sfVert);CHKERRQ(ierr); 3099 if (interpolate) { 3100 DM idm; 3101 3102 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 3103 ierr = DMDestroy(dm);CHKERRQ(ierr); 3104 *dm = idm; 3105 } 3106 ierr = DMPlexBuildCoordinatesFromCellListParallel(*dm, spaceDim, sfVert, vertexCoords);CHKERRQ(ierr); 3107 if (vertexSF) *vertexSF = sfVert; 3108 else {ierr = PetscSFDestroy(&sfVert);CHKERRQ(ierr);} 3109 PetscFunctionReturn(0); 3110 } 3111 3112 3113 /*@ 3114 DMPlexCreateFromCellListParallel - Deprecated, use DMPlexCreateFromCellListParallelPetsc() 3115 3116 Level: deprecated 3117 3118 .seealso: DMPlexCreateFromCellListParallelPetsc() 3119 @*/ 3120 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) 3121 { 3122 PetscErrorCode ierr; 3123 PetscInt i; 3124 PetscInt *pintCells; 3125 3126 PetscFunctionBegin; 3127 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)); 3128 if (sizeof(int) == sizeof(PetscInt)) { 3129 pintCells = (PetscInt *) cells; 3130 } else { 3131 ierr = PetscMalloc1(numCells*numCorners, &pintCells);CHKERRQ(ierr); 3132 for (i = 0; i < numCells*numCorners; i++) { 3133 pintCells[i] = (PetscInt) cells[i]; 3134 } 3135 } 3136 ierr = DMPlexCreateFromCellListParallelPetsc(comm, dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, pintCells, spaceDim, vertexCoords, vertexSF, dm);CHKERRQ(ierr); 3137 if (sizeof(int) != sizeof(PetscInt)) { 3138 ierr = PetscFree(pintCells);CHKERRQ(ierr); 3139 } 3140 PetscFunctionReturn(0); 3141 } 3142 3143 /*@C 3144 DMPlexBuildFromCellList - Build DMPLEX topology from a list of vertices for each cell (common mesh generator output) 3145 3146 Input Parameters: 3147 + dm - The DM 3148 . numCells - The number of cells owned by this process 3149 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE 3150 . numCorners - The number of vertices for each cell 3151 - cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 3152 3153 Level: advanced 3154 3155 Notes: 3156 Two triangles sharing a face 3157 $ 3158 $ 2 3159 $ / | \ 3160 $ / | \ 3161 $ / | \ 3162 $ 0 0 | 1 3 3163 $ \ | / 3164 $ \ | / 3165 $ \ | / 3166 $ 1 3167 would have input 3168 $ numCells = 2, numVertices = 4 3169 $ cells = [0 1 2 1 3 2] 3170 $ 3171 which would result in the DMPlex 3172 $ 3173 $ 4 3174 $ / | \ 3175 $ / | \ 3176 $ / | \ 3177 $ 2 0 | 1 5 3178 $ \ | / 3179 $ \ | / 3180 $ \ | / 3181 $ 3 3182 3183 If numVertices is PETSC_DECIDE, it is computed by PETSc as the maximum vertex index in cells + 1. 3184 3185 Not currently supported in Fortran. 3186 3187 .seealso: DMPlexBuildFromCellListParallel(), DMPlexBuildCoordinatesFromCellList(), DMPlexCreateFromCellListPetsc() 3188 @*/ 3189 PetscErrorCode DMPlexBuildFromCellList(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const PetscInt cells[]) 3190 { 3191 PetscInt *cones, c, p, dim; 3192 PetscErrorCode ierr; 3193 3194 PetscFunctionBegin; 3195 ierr = PetscLogEventBegin(DMPLEX_BuildFromCellList,dm,0,0,0);CHKERRQ(ierr); 3196 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3197 /* Get/check global number of vertices */ 3198 { 3199 PetscInt NVerticesInCells, i; 3200 const PetscInt len = numCells * numCorners; 3201 3202 /* NVerticesInCells = max(cells) + 1 */ 3203 NVerticesInCells = PETSC_MIN_INT; 3204 for (i=0; i<len; i++) if (cells[i] > NVerticesInCells) NVerticesInCells = cells[i]; 3205 ++NVerticesInCells; 3206 3207 if (numVertices == PETSC_DECIDE) numVertices = NVerticesInCells; 3208 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); 3209 } 3210 ierr = DMPlexSetChart(dm, 0, numCells+numVertices);CHKERRQ(ierr); 3211 for (c = 0; c < numCells; ++c) { 3212 ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr); 3213 } 3214 ierr = DMSetUp(dm);CHKERRQ(ierr); 3215 ierr = DMPlexGetCones(dm,&cones);CHKERRQ(ierr); 3216 for (c = 0; c < numCells; ++c) { 3217 for (p = 0; p < numCorners; ++p) { 3218 cones[c*numCorners+p] = cells[c*numCorners+p]+numCells; 3219 } 3220 } 3221 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 3222 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 3223 ierr = PetscLogEventEnd(DMPLEX_BuildFromCellList,dm,0,0,0);CHKERRQ(ierr); 3224 PetscFunctionReturn(0); 3225 } 3226 3227 /*@C 3228 DMPlexBuildCoordinatesFromCellList - Build DM coordinates from a list of coordinates for each owned vertex (common mesh generator output) 3229 3230 Input Parameters: 3231 + dm - The DM 3232 . spaceDim - The spatial dimension used for coordinates 3233 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 3234 3235 Level: advanced 3236 3237 Notes: 3238 Not currently supported in Fortran. 3239 3240 .seealso: DMPlexBuildCoordinatesFromCellListParallel(), DMPlexCreateFromCellListPetsc(), DMPlexBuildFromCellList() 3241 @*/ 3242 PetscErrorCode DMPlexBuildCoordinatesFromCellList(DM dm, PetscInt spaceDim, const PetscReal vertexCoords[]) 3243 { 3244 PetscSection coordSection; 3245 Vec coordinates; 3246 DM cdm; 3247 PetscScalar *coords; 3248 PetscInt v, vStart, vEnd, d; 3249 PetscErrorCode ierr; 3250 3251 PetscFunctionBegin; 3252 ierr = PetscLogEventBegin(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0);CHKERRQ(ierr); 3253 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 3254 if (vStart < 0 || vEnd < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM is not set up properly. DMPlexBuildFromCellList() should be called first."); 3255 ierr = DMSetCoordinateDim(dm, spaceDim);CHKERRQ(ierr); 3256 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 3257 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 3258 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 3259 ierr = PetscSectionSetChart(coordSection, vStart, vEnd);CHKERRQ(ierr); 3260 for (v = vStart; v < vEnd; ++v) { 3261 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 3262 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 3263 } 3264 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 3265 3266 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 3267 ierr = DMCreateLocalVector(cdm, &coordinates);CHKERRQ(ierr); 3268 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 3269 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 3270 ierr = VecGetArrayWrite(coordinates, &coords);CHKERRQ(ierr); 3271 for (v = 0; v < vEnd-vStart; ++v) { 3272 for (d = 0; d < spaceDim; ++d) { 3273 coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d]; 3274 } 3275 } 3276 ierr = VecRestoreArrayWrite(coordinates, &coords);CHKERRQ(ierr); 3277 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 3278 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 3279 ierr = PetscLogEventEnd(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0);CHKERRQ(ierr); 3280 PetscFunctionReturn(0); 3281 } 3282 3283 /*@ 3284 DMPlexCreateFromCellListPetsc - Create DMPLEX from a list of vertices for each cell (common mesh generator output) 3285 3286 Input Parameters: 3287 + comm - The communicator 3288 . dim - The topological dimension of the mesh 3289 . numCells - The number of cells 3290 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE 3291 . numCorners - The number of vertices for each cell 3292 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 3293 . cells - An array of numCells*numCorners numbers, the vertices for each cell 3294 . spaceDim - The spatial dimension used for coordinates 3295 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 3296 3297 Output Parameter: 3298 . dm - The DM 3299 3300 Notes: 3301 This function is just a convenient sequence of DMCreate(), DMSetType(), DMSetDimension(), DMPlexBuildFromCellList(), 3302 DMPlexInterpolate(), DMPlexBuildCoordinatesFromCellList() 3303 3304 See DMPlexBuildFromCellList() for an example and details about the topology-related parameters. 3305 See DMPlexBuildCoordinatesFromCellList() for details about the geometry-related parameters. 3306 3307 Level: intermediate 3308 3309 .seealso: DMPlexCreateFromCellListParallelPetsc(), DMPlexBuildFromCellList(), DMPlexBuildCoordinatesFromCellList(), DMPlexCreateFromDAG(), DMPlexCreate() 3310 @*/ 3311 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) 3312 { 3313 PetscErrorCode ierr; 3314 3315 PetscFunctionBegin; 3316 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."); 3317 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 3318 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 3319 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 3320 ierr = DMPlexBuildFromCellList(*dm, numCells, numVertices, numCorners, cells);CHKERRQ(ierr); 3321 if (interpolate) { 3322 DM idm; 3323 3324 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 3325 ierr = DMDestroy(dm);CHKERRQ(ierr); 3326 *dm = idm; 3327 } 3328 ierr = DMPlexBuildCoordinatesFromCellList(*dm, spaceDim, vertexCoords);CHKERRQ(ierr); 3329 PetscFunctionReturn(0); 3330 } 3331 3332 /*@ 3333 DMPlexCreateFromCellList - Deprecated, use DMPlexCreateFromCellListPetsc() 3334 3335 Level: deprecated 3336 3337 .seealso: DMPlexCreateFromCellListPetsc() 3338 @*/ 3339 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) 3340 { 3341 PetscErrorCode ierr; 3342 PetscInt i; 3343 PetscInt *pintCells; 3344 PetscReal *prealVC; 3345 3346 PetscFunctionBegin; 3347 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)); 3348 if (sizeof(int) == sizeof(PetscInt)) { 3349 pintCells = (PetscInt *) cells; 3350 } else { 3351 ierr = PetscMalloc1(numCells*numCorners, &pintCells);CHKERRQ(ierr); 3352 for (i = 0; i < numCells*numCorners; i++) { 3353 pintCells[i] = (PetscInt) cells[i]; 3354 } 3355 } 3356 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)); 3357 if (sizeof(double) == sizeof(PetscReal)) { 3358 prealVC = (PetscReal *) vertexCoords; 3359 } else { 3360 ierr = PetscMalloc1(numVertices*spaceDim, &prealVC);CHKERRQ(ierr); 3361 for (i = 0; i < numVertices*spaceDim; i++) { 3362 prealVC[i] = (PetscReal) vertexCoords[i]; 3363 } 3364 } 3365 ierr = DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, pintCells, spaceDim, prealVC, dm);CHKERRQ(ierr); 3366 if (sizeof(int) != sizeof(PetscInt)) { 3367 ierr = PetscFree(pintCells);CHKERRQ(ierr); 3368 } 3369 if (sizeof(double) != sizeof(PetscReal)) { 3370 ierr = PetscFree(prealVC);CHKERRQ(ierr); 3371 } 3372 PetscFunctionReturn(0); 3373 } 3374 3375 /*@ 3376 DMPlexCreateFromDAG - This takes as input the adjacency-list representation of the Directed Acyclic Graph (Hasse Diagram) encoding a mesh, and produces a DM 3377 3378 Input Parameters: 3379 + dm - The empty DM object, usually from DMCreate() and DMSetDimension() 3380 . depth - The depth of the DAG 3381 . numPoints - Array of size depth + 1 containing the number of points at each depth 3382 . coneSize - The cone size of each point 3383 . cones - The concatenation of the cone points for each point, the cone list must be oriented correctly for each point 3384 . coneOrientations - The orientation of each cone point 3385 - vertexCoords - An array of numPoints[0]*spacedim numbers representing the coordinates of each vertex, with spacedim the value set via DMSetCoordinateDim() 3386 3387 Output Parameter: 3388 . dm - The DM 3389 3390 Note: Two triangles sharing a face would have input 3391 $ depth = 1, numPoints = [4 2], coneSize = [3 3 0 0 0 0] 3392 $ cones = [2 3 4 3 5 4], coneOrientations = [0 0 0 0 0 0] 3393 $ vertexCoords = [-1.0 0.0 0.0 -1.0 0.0 1.0 1.0 0.0] 3394 $ 3395 which would result in the DMPlex 3396 $ 3397 $ 4 3398 $ / | \ 3399 $ / | \ 3400 $ / | \ 3401 $ 2 0 | 1 5 3402 $ \ | / 3403 $ \ | / 3404 $ \ | / 3405 $ 3 3406 $ 3407 $ Notice that all points are numbered consecutively, unlike DMPlexCreateFromCellListPetsc() 3408 3409 Level: advanced 3410 3411 .seealso: DMPlexCreateFromCellListPetsc(), DMPlexCreate() 3412 @*/ 3413 PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[]) 3414 { 3415 Vec coordinates; 3416 PetscSection coordSection; 3417 PetscScalar *coords; 3418 PetscInt coordSize, firstVertex = -1, pStart = 0, pEnd = 0, p, v, dim, dimEmbed, d, off; 3419 PetscErrorCode ierr; 3420 3421 PetscFunctionBegin; 3422 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3423 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 3424 if (dimEmbed < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Embedding dimension %D cannot be less than intrinsic dimension %d",dimEmbed,dim); 3425 for (d = 0; d <= depth; ++d) pEnd += numPoints[d]; 3426 ierr = DMPlexSetChart(dm, pStart, pEnd);CHKERRQ(ierr); 3427 for (p = pStart; p < pEnd; ++p) { 3428 ierr = DMPlexSetConeSize(dm, p, coneSize[p-pStart]);CHKERRQ(ierr); 3429 if (firstVertex < 0 && !coneSize[p - pStart]) { 3430 firstVertex = p - pStart; 3431 } 3432 } 3433 if (firstVertex < 0 && numPoints[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected %D vertices but could not find any", numPoints[0]); 3434 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 3435 for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) { 3436 ierr = DMPlexSetCone(dm, p, &cones[off]);CHKERRQ(ierr); 3437 ierr = DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);CHKERRQ(ierr); 3438 } 3439 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 3440 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 3441 /* Build coordinates */ 3442 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 3443 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 3444 ierr = PetscSectionSetFieldComponents(coordSection, 0, dimEmbed);CHKERRQ(ierr); 3445 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0]);CHKERRQ(ierr); 3446 for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) { 3447 ierr = PetscSectionSetDof(coordSection, v, dimEmbed);CHKERRQ(ierr); 3448 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed);CHKERRQ(ierr); 3449 } 3450 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 3451 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 3452 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 3453 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 3454 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 3455 ierr = VecSetBlockSize(coordinates, dimEmbed);CHKERRQ(ierr); 3456 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 3457 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 3458 for (v = 0; v < numPoints[0]; ++v) { 3459 PetscInt off; 3460 3461 ierr = PetscSectionGetOffset(coordSection, v+firstVertex, &off);CHKERRQ(ierr); 3462 for (d = 0; d < dimEmbed; ++d) { 3463 coords[off+d] = vertexCoords[v*dimEmbed+d]; 3464 } 3465 } 3466 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 3467 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 3468 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 3469 PetscFunctionReturn(0); 3470 } 3471 3472 /*@C 3473 DMPlexCreateCellVertexFromFile - Create a DMPlex mesh from a simple cell-vertex file. 3474 3475 + comm - The MPI communicator 3476 . filename - Name of the .dat file 3477 - interpolate - Create faces and edges in the mesh 3478 3479 Output Parameter: 3480 . dm - The DM object representing the mesh 3481 3482 Note: The format is the simplest possible: 3483 $ Ne 3484 $ v0 v1 ... vk 3485 $ Nv 3486 $ x y z marker 3487 3488 Level: beginner 3489 3490 .seealso: DMPlexCreateFromFile(), DMPlexCreateMedFromFile(), DMPlexCreateGmsh(), DMPlexCreate() 3491 @*/ 3492 PetscErrorCode DMPlexCreateCellVertexFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 3493 { 3494 DMLabel marker; 3495 PetscViewer viewer; 3496 Vec coordinates; 3497 PetscSection coordSection; 3498 PetscScalar *coords; 3499 char line[PETSC_MAX_PATH_LEN]; 3500 PetscInt dim = 3, cdim = 3, coordSize, v, c, d; 3501 PetscMPIInt rank; 3502 int snum, Nv, Nc; 3503 PetscErrorCode ierr; 3504 3505 PetscFunctionBegin; 3506 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 3507 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 3508 ierr = PetscViewerSetType(viewer, PETSCVIEWERASCII);CHKERRQ(ierr); 3509 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 3510 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 3511 if (!rank) { 3512 ierr = PetscViewerRead(viewer, line, 2, NULL, PETSC_STRING);CHKERRQ(ierr); 3513 snum = sscanf(line, "%d %d", &Nc, &Nv); 3514 if (snum != 2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 3515 } else { 3516 Nc = Nv = 0; 3517 } 3518 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 3519 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 3520 ierr = DMPlexSetChart(*dm, 0, Nc+Nv);CHKERRQ(ierr); 3521 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 3522 ierr = DMSetCoordinateDim(*dm, cdim);CHKERRQ(ierr); 3523 /* Read topology */ 3524 if (!rank) { 3525 PetscInt cone[8], corners = 8; 3526 int vbuf[8], v; 3527 3528 for (c = 0; c < Nc; ++c) {ierr = DMPlexSetConeSize(*dm, c, corners);CHKERRQ(ierr);} 3529 ierr = DMSetUp(*dm);CHKERRQ(ierr); 3530 for (c = 0; c < Nc; ++c) { 3531 ierr = PetscViewerRead(viewer, line, corners, NULL, PETSC_STRING);CHKERRQ(ierr); 3532 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]); 3533 if (snum != corners) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 3534 for (v = 0; v < corners; ++v) cone[v] = vbuf[v] + Nc; 3535 /* Hexahedra are inverted */ 3536 { 3537 PetscInt tmp = cone[1]; 3538 cone[1] = cone[3]; 3539 cone[3] = tmp; 3540 } 3541 ierr = DMPlexSetCone(*dm, c, cone);CHKERRQ(ierr); 3542 } 3543 } 3544 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 3545 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 3546 /* Read coordinates */ 3547 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 3548 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 3549 ierr = PetscSectionSetFieldComponents(coordSection, 0, cdim);CHKERRQ(ierr); 3550 ierr = PetscSectionSetChart(coordSection, Nc, Nc + Nv);CHKERRQ(ierr); 3551 for (v = Nc; v < Nc+Nv; ++v) { 3552 ierr = PetscSectionSetDof(coordSection, v, cdim);CHKERRQ(ierr); 3553 ierr = PetscSectionSetFieldDof(coordSection, v, 0, cdim);CHKERRQ(ierr); 3554 } 3555 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 3556 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 3557 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 3558 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 3559 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 3560 ierr = VecSetBlockSize(coordinates, cdim);CHKERRQ(ierr); 3561 ierr = VecSetType(coordinates, VECSTANDARD);CHKERRQ(ierr); 3562 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 3563 if (!rank) { 3564 double x[3]; 3565 int val; 3566 3567 ierr = DMCreateLabel(*dm, "marker");CHKERRQ(ierr); 3568 ierr = DMGetLabel(*dm, "marker", &marker);CHKERRQ(ierr); 3569 for (v = 0; v < Nv; ++v) { 3570 ierr = PetscViewerRead(viewer, line, 4, NULL, PETSC_STRING);CHKERRQ(ierr); 3571 snum = sscanf(line, "%lg %lg %lg %d", &x[0], &x[1], &x[2], &val); 3572 if (snum != 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 3573 for (d = 0; d < cdim; ++d) coords[v*cdim+d] = x[d]; 3574 if (val) {ierr = DMLabelSetValue(marker, v+Nc, val);CHKERRQ(ierr);} 3575 } 3576 } 3577 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 3578 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 3579 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 3580 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3581 if (interpolate) { 3582 DM idm; 3583 DMLabel bdlabel; 3584 3585 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 3586 ierr = DMDestroy(dm);CHKERRQ(ierr); 3587 *dm = idm; 3588 3589 ierr = DMGetLabel(*dm, "marker", &bdlabel);CHKERRQ(ierr); 3590 ierr = DMPlexMarkBoundaryFaces(*dm, PETSC_DETERMINE, bdlabel);CHKERRQ(ierr); 3591 ierr = DMPlexLabelComplete(*dm, bdlabel);CHKERRQ(ierr); 3592 } 3593 PetscFunctionReturn(0); 3594 } 3595 3596 /*@C 3597 DMPlexCreateFromFile - This takes a filename and produces a DM 3598 3599 Input Parameters: 3600 + comm - The communicator 3601 . filename - A file name 3602 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 3603 3604 Output Parameter: 3605 . dm - The DM 3606 3607 Options Database Keys: 3608 . -dm_plex_create_from_hdf5_xdmf - use the PETSC_VIEWER_HDF5_XDMF format for reading HDF5 3609 3610 Use -dm_plex_create_ prefix to pass options to the internal PetscViewer, e.g. 3611 $ -dm_plex_create_viewer_hdf5_collective 3612 3613 Level: beginner 3614 3615 .seealso: DMPlexCreateFromDAG(), DMPlexCreateFromCellListPetsc(), DMPlexCreate() 3616 @*/ 3617 PetscErrorCode DMPlexCreateFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 3618 { 3619 const char *extGmsh = ".msh"; 3620 const char *extGmsh2 = ".msh2"; 3621 const char *extGmsh4 = ".msh4"; 3622 const char *extCGNS = ".cgns"; 3623 const char *extExodus = ".exo"; 3624 const char *extGenesis = ".gen"; 3625 const char *extFluent = ".cas"; 3626 const char *extHDF5 = ".h5"; 3627 const char *extMed = ".med"; 3628 const char *extPLY = ".ply"; 3629 const char *extEGADS = ".egadslite"; 3630 const char *extCV = ".dat"; 3631 size_t len; 3632 PetscBool isGmsh, isGmsh2, isGmsh4, isCGNS, isExodus, isGenesis, isFluent, isHDF5, isMed, isPLY, isEGADS, isCV; 3633 PetscMPIInt rank; 3634 PetscErrorCode ierr; 3635 3636 PetscFunctionBegin; 3637 PetscValidCharPointer(filename, 2); 3638 PetscValidPointer(dm, 4); 3639 ierr = DMInitializePackage();CHKERRQ(ierr); 3640 ierr = PetscLogEventBegin(DMPLEX_CreateFromFile,0,0,0,0);CHKERRQ(ierr); 3641 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 3642 ierr = PetscStrlen(filename, &len);CHKERRQ(ierr); 3643 if (!len) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Filename must be a valid path"); 3644 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGmsh, 4, &isGmsh);CHKERRQ(ierr); 3645 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extGmsh2, 5, &isGmsh2);CHKERRQ(ierr); 3646 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extGmsh4, 5, &isGmsh4);CHKERRQ(ierr); 3647 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extCGNS, 5, &isCGNS);CHKERRQ(ierr); 3648 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extExodus, 4, &isExodus);CHKERRQ(ierr); 3649 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGenesis, 4, &isGenesis);CHKERRQ(ierr); 3650 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extFluent, 4, &isFluent);CHKERRQ(ierr); 3651 ierr = PetscStrncmp(&filename[PetscMax(0,len-3)], extHDF5, 3, &isHDF5);CHKERRQ(ierr); 3652 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extMed, 4, &isMed);CHKERRQ(ierr); 3653 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extPLY, 4, &isPLY);CHKERRQ(ierr); 3654 ierr = PetscStrncmp(&filename[PetscMax(0,len-10)], extEGADS, 9, &isEGADS);CHKERRQ(ierr); 3655 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extCV, 4, &isCV);CHKERRQ(ierr); 3656 if (isGmsh || isGmsh2 || isGmsh4) { 3657 ierr = DMPlexCreateGmshFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3658 } else if (isCGNS) { 3659 ierr = DMPlexCreateCGNSFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3660 } else if (isExodus || isGenesis) { 3661 ierr = DMPlexCreateExodusFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3662 } else if (isFluent) { 3663 ierr = DMPlexCreateFluentFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3664 } else if (isHDF5) { 3665 PetscBool load_hdf5_xdmf = PETSC_FALSE; 3666 PetscViewer viewer; 3667 3668 /* PETSC_VIEWER_HDF5_XDMF is used if the filename ends with .xdmf.h5, or if -dm_plex_create_from_hdf5_xdmf option is present */ 3669 ierr = PetscStrncmp(&filename[PetscMax(0,len-8)], ".xdmf", 5, &load_hdf5_xdmf);CHKERRQ(ierr); 3670 ierr = PetscOptionsGetBool(NULL, NULL, "-dm_plex_create_from_hdf5_xdmf", &load_hdf5_xdmf, NULL);CHKERRQ(ierr); 3671 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 3672 ierr = PetscViewerSetType(viewer, PETSCVIEWERHDF5);CHKERRQ(ierr); 3673 ierr = PetscViewerSetOptionsPrefix(viewer, "dm_plex_create_");CHKERRQ(ierr); 3674 ierr = PetscViewerSetFromOptions(viewer);CHKERRQ(ierr); 3675 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 3676 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 3677 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 3678 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 3679 if (load_hdf5_xdmf) {ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_XDMF);CHKERRQ(ierr);} 3680 ierr = DMLoad(*dm, viewer);CHKERRQ(ierr); 3681 if (load_hdf5_xdmf) {ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);} 3682 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3683 3684 if (interpolate) { 3685 DM idm; 3686 3687 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 3688 ierr = DMDestroy(dm);CHKERRQ(ierr); 3689 *dm = idm; 3690 } 3691 } else if (isMed) { 3692 ierr = DMPlexCreateMedFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3693 } else if (isPLY) { 3694 ierr = DMPlexCreatePLYFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3695 } else if (isEGADS) { 3696 ierr = DMPlexCreateEGADSFromFile(comm, filename, dm);CHKERRQ(ierr); 3697 if (!interpolate) { 3698 DM udm; 3699 3700 ierr = DMPlexUninterpolate(*dm, &udm);CHKERRQ(ierr); 3701 ierr = DMDestroy(dm);CHKERRQ(ierr); 3702 *dm = udm; 3703 } 3704 } else if (isCV) { 3705 ierr = DMPlexCreateCellVertexFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3706 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot load file %s: unrecognized extension", filename); 3707 ierr = PetscLogEventEnd(DMPLEX_CreateFromFile,0,0,0,0);CHKERRQ(ierr); 3708 PetscFunctionReturn(0); 3709 } 3710 3711 /*@ 3712 DMPlexCreateReferenceCellByType - Create a DMPLEX with the appropriate FEM reference cell 3713 3714 Collective 3715 3716 Input Parameters: 3717 + comm - The communicator 3718 - ct - The cell type of the reference cell 3719 3720 Output Parameter: 3721 . refdm - The reference cell 3722 3723 Options Database Keys: 3724 . -dm_plex_ref_type <ct> - Specify the celltyoe for the reference cell 3725 3726 Level: intermediate 3727 3728 .seealso: DMPlexCreateReferenceCell(), DMPlexCreateBoxMesh() 3729 @*/ 3730 PetscErrorCode DMPlexCreateReferenceCellByType(MPI_Comm comm, DMPolytopeType ct, DM *refdm) 3731 { 3732 DM rdm; 3733 Vec coords; 3734 PetscErrorCode ierr; 3735 3736 PetscFunctionBegin; 3737 ierr = PetscOptionsGetEnum(NULL, NULL, "-dm_plex_ref_type", DMPolytopeTypes, (PetscEnum *) &ct, NULL);CHKERRQ(ierr); 3738 ierr = DMCreate(comm, &rdm);CHKERRQ(ierr); 3739 ierr = DMSetType(rdm, DMPLEX);CHKERRQ(ierr); 3740 switch (ct) { 3741 case DM_POLYTOPE_POINT: 3742 { 3743 PetscInt numPoints[1] = {1}; 3744 PetscInt coneSize[1] = {0}; 3745 PetscInt cones[1] = {0}; 3746 PetscInt coneOrientations[1] = {0}; 3747 PetscScalar vertexCoords[1] = {0.0}; 3748 3749 ierr = DMSetDimension(rdm, 0);CHKERRQ(ierr); 3750 ierr = DMPlexCreateFromDAG(rdm, 0, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3751 } 3752 break; 3753 case DM_POLYTOPE_SEGMENT: 3754 { 3755 PetscInt numPoints[2] = {2, 1}; 3756 PetscInt coneSize[3] = {2, 0, 0}; 3757 PetscInt cones[2] = {1, 2}; 3758 PetscInt coneOrientations[2] = {0, 0}; 3759 PetscScalar vertexCoords[2] = {-1.0, 1.0}; 3760 3761 ierr = DMSetDimension(rdm, 1);CHKERRQ(ierr); 3762 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3763 } 3764 break; 3765 case DM_POLYTOPE_TRIANGLE: 3766 { 3767 PetscInt numPoints[2] = {3, 1}; 3768 PetscInt coneSize[4] = {3, 0, 0, 0}; 3769 PetscInt cones[3] = {1, 2, 3}; 3770 PetscInt coneOrientations[3] = {0, 0, 0}; 3771 PetscScalar vertexCoords[6] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0}; 3772 3773 ierr = DMSetDimension(rdm, 2);CHKERRQ(ierr); 3774 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3775 } 3776 break; 3777 case DM_POLYTOPE_QUADRILATERAL: 3778 { 3779 PetscInt numPoints[2] = {4, 1}; 3780 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 3781 PetscInt cones[4] = {1, 2, 3, 4}; 3782 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 3783 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0}; 3784 3785 ierr = DMSetDimension(rdm, 2);CHKERRQ(ierr); 3786 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3787 } 3788 break; 3789 case DM_POLYTOPE_SEG_PRISM_TENSOR: 3790 { 3791 PetscInt numPoints[2] = {4, 1}; 3792 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 3793 PetscInt cones[4] = {1, 2, 3, 4}; 3794 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 3795 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0}; 3796 3797 ierr = DMSetDimension(rdm, 2);CHKERRQ(ierr); 3798 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3799 } 3800 break; 3801 case DM_POLYTOPE_TETRAHEDRON: 3802 { 3803 PetscInt numPoints[2] = {4, 1}; 3804 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 3805 PetscInt cones[4] = {1, 3, 2, 4}; 3806 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 3807 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}; 3808 3809 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 3810 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3811 } 3812 break; 3813 case DM_POLYTOPE_HEXAHEDRON: 3814 { 3815 PetscInt numPoints[2] = {8, 1}; 3816 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0}; 3817 PetscInt cones[8] = {1, 4, 3, 2, 5, 6, 7, 8}; 3818 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 3819 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, 3820 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0}; 3821 3822 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 3823 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3824 } 3825 break; 3826 case DM_POLYTOPE_TRI_PRISM: 3827 { 3828 PetscInt numPoints[2] = {6, 1}; 3829 PetscInt coneSize[7] = {6, 0, 0, 0, 0, 0, 0}; 3830 PetscInt cones[6] = {1, 3, 2, 4, 5, 6}; 3831 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 3832 PetscScalar vertexCoords[18] = {-1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 3833 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0}; 3834 3835 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 3836 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3837 } 3838 break; 3839 case DM_POLYTOPE_TRI_PRISM_TENSOR: 3840 { 3841 PetscInt numPoints[2] = {6, 1}; 3842 PetscInt coneSize[7] = {6, 0, 0, 0, 0, 0, 0}; 3843 PetscInt cones[6] = {1, 2, 3, 4, 5, 6}; 3844 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 3845 PetscScalar vertexCoords[18] = {-1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 3846 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0}; 3847 3848 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 3849 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3850 } 3851 break; 3852 case DM_POLYTOPE_QUAD_PRISM_TENSOR: 3853 { 3854 PetscInt numPoints[2] = {8, 1}; 3855 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0}; 3856 PetscInt cones[8] = {1, 2, 3, 4, 5, 6, 7, 8}; 3857 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 3858 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, 3859 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0}; 3860 3861 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 3862 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3863 } 3864 break; 3865 case DM_POLYTOPE_PYRAMID: 3866 { 3867 PetscInt numPoints[2] = {5, 1}; 3868 PetscInt coneSize[6] = {5, 0, 0, 0, 0, 0}; 3869 PetscInt cones[5] = {1, 4, 3, 2, 5}; 3870 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 3871 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, 3872 0.0, 0.0, 1.0}; 3873 3874 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 3875 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3876 } 3877 break; 3878 default: SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Cannot create reference cell for cell type %s", DMPolytopeTypes[ct]); 3879 } 3880 { 3881 PetscInt Nv, v; 3882 3883 /* Must create the celltype label here so that we do not automatically try to compute the types */ 3884 ierr = DMCreateLabel(rdm, "celltype");CHKERRQ(ierr); 3885 ierr = DMPlexSetCellType(rdm, 0, ct);CHKERRQ(ierr); 3886 ierr = DMPlexGetChart(rdm, NULL, &Nv);CHKERRQ(ierr); 3887 for (v = 1; v < Nv; ++v) {ierr = DMPlexSetCellType(rdm, v, DM_POLYTOPE_POINT);CHKERRQ(ierr);} 3888 } 3889 ierr = DMPlexInterpolate(rdm, refdm);CHKERRQ(ierr); 3890 if (rdm->coordinateDM) { 3891 DM ncdm; 3892 PetscSection cs; 3893 PetscInt pEnd = -1; 3894 3895 ierr = DMGetLocalSection(rdm->coordinateDM, &cs);CHKERRQ(ierr); 3896 if (cs) {ierr = PetscSectionGetChart(cs, NULL, &pEnd);CHKERRQ(ierr);} 3897 if (pEnd >= 0) { 3898 ierr = DMClone(*refdm, &ncdm);CHKERRQ(ierr); 3899 ierr = DMCopyDisc(rdm->coordinateDM, ncdm);CHKERRQ(ierr); 3900 ierr = DMSetLocalSection(ncdm, cs);CHKERRQ(ierr); 3901 ierr = DMSetCoordinateDM(*refdm, ncdm);CHKERRQ(ierr); 3902 ierr = DMDestroy(&ncdm);CHKERRQ(ierr); 3903 } 3904 } 3905 ierr = DMGetCoordinatesLocal(rdm, &coords);CHKERRQ(ierr); 3906 if (coords) { 3907 ierr = DMSetCoordinatesLocal(*refdm, coords);CHKERRQ(ierr); 3908 } else { 3909 ierr = DMGetCoordinates(rdm, &coords);CHKERRQ(ierr); 3910 if (coords) {ierr = DMSetCoordinates(*refdm, coords);CHKERRQ(ierr);} 3911 } 3912 ierr = DMDestroy(&rdm);CHKERRQ(ierr); 3913 PetscFunctionReturn(0); 3914 } 3915 3916 /*@ 3917 DMPlexCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell 3918 3919 Collective 3920 3921 Input Parameters: 3922 + comm - The communicator 3923 . dim - The spatial dimension 3924 - simplex - Flag for simplex, otherwise use a tensor-product cell 3925 3926 Output Parameter: 3927 . refdm - The reference cell 3928 3929 Level: intermediate 3930 3931 .seealso: DMPlexCreateReferenceCellByType(), DMPlexCreateBoxMesh() 3932 @*/ 3933 PetscErrorCode DMPlexCreateReferenceCell(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *refdm) 3934 { 3935 PetscErrorCode ierr; 3936 3937 PetscFunctionBeginUser; 3938 switch (dim) { 3939 case 0: ierr = DMPlexCreateReferenceCellByType(comm, DM_POLYTOPE_POINT, refdm);CHKERRQ(ierr);break; 3940 case 1: ierr = DMPlexCreateReferenceCellByType(comm, DM_POLYTOPE_SEGMENT, refdm);CHKERRQ(ierr);break; 3941 case 2: 3942 if (simplex) {ierr = DMPlexCreateReferenceCellByType(comm, DM_POLYTOPE_TRIANGLE, refdm);CHKERRQ(ierr);} 3943 else {ierr = DMPlexCreateReferenceCellByType(comm, DM_POLYTOPE_QUADRILATERAL, refdm);CHKERRQ(ierr);} 3944 break; 3945 case 3: 3946 if (simplex) {ierr = DMPlexCreateReferenceCellByType(comm, DM_POLYTOPE_TETRAHEDRON, refdm);CHKERRQ(ierr);} 3947 else {ierr = DMPlexCreateReferenceCellByType(comm, DM_POLYTOPE_HEXAHEDRON, refdm);CHKERRQ(ierr);} 3948 break; 3949 default: 3950 SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Cannot create reference cell for dimension %D", dim); 3951 } 3952 PetscFunctionReturn(0); 3953 } 3954