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