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