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 - -dm_plex_box_interpolate <bool> - PETSC_TRUE turns on topological interpolation (creating edges and faces) 1013 1014 Notes: 1015 The options database keys above take lists of length d in d dimensions. 1016 1017 Here is the numbering returned for 2 faces in each direction for tensor cells: 1018 $ 10---17---11---18----12 1019 $ | | | 1020 $ | | | 1021 $ 20 2 22 3 24 1022 $ | | | 1023 $ | | | 1024 $ 7---15----8---16----9 1025 $ | | | 1026 $ | | | 1027 $ 19 0 21 1 23 1028 $ | | | 1029 $ | | | 1030 $ 4---13----5---14----6 1031 1032 and for simplicial cells 1033 1034 $ 14----8---15----9----16 1035 $ |\ 5 |\ 7 | 1036 $ | \ | \ | 1037 $ 13 2 14 3 15 1038 $ | 4 \ | 6 \ | 1039 $ | \ | \ | 1040 $ 11----6---12----7----13 1041 $ |\ |\ | 1042 $ | \ 1 | \ 3 | 1043 $ 10 0 11 1 12 1044 $ | 0 \ | 2 \ | 1045 $ | \ | \ | 1046 $ 8----4----9----5----10 1047 1048 Level: beginner 1049 1050 .seealso: DMPlexCreateFromFile(), DMPlexCreateHexCylinderMesh(), DMSetType(), DMCreate() 1051 @*/ 1052 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) 1053 { 1054 PetscInt fac[3] = {0, 0, 0}; 1055 PetscReal low[3] = {0, 0, 0}; 1056 PetscReal upp[3] = {1, 1, 1}; 1057 DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE}; 1058 PetscInt i, n; 1059 PetscBool flg; 1060 PetscErrorCode ierr; 1061 1062 PetscFunctionBegin; 1063 ierr = PetscOptionsGetInt(NULL, NULL, "-dm_plex_box_dim", &dim, &flg);CHKERRQ(ierr); 1064 if ((dim < 0) || (dim > 3)) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %D should be in [1, 3]", dim); 1065 ierr = PetscOptionsGetBool(NULL, NULL, "-dm_plex_box_simplex", &simplex, &flg);CHKERRQ(ierr); 1066 n = 3; 1067 ierr = PetscOptionsGetIntArray(NULL, NULL, "-dm_plex_box_faces", fac, &n, &flg);CHKERRQ(ierr); 1068 for (i = 0; i < dim; ++i) fac[i] = faces ? faces[i] : (flg && i < n ? fac[i] : (dim == 1 ? 1 : 4-dim)); 1069 if (lower) for (i = 0; i < dim; ++i) low[i] = lower[i]; 1070 if (upper) for (i = 0; i < dim; ++i) upp[i] = upper[i]; 1071 if (periodicity) for (i = 0; i < dim; ++i) bdt[i] = periodicity[i]; 1072 /* Allow bounds to be specified from the command line */ 1073 n = 3; 1074 ierr = PetscOptionsGetRealArray(NULL, NULL, "-dm_plex_box_lower", low, &n, &flg);CHKERRQ(ierr); 1075 if (flg && (n != dim)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Lower box point had %D values, should have been %D", n, dim); 1076 n = 3; 1077 ierr = PetscOptionsGetRealArray(NULL, NULL, "-dm_plex_box_upper", upp, &n, &flg);CHKERRQ(ierr); 1078 if (flg && (n != dim)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Upper box point had %D values, should have been %D", n, dim); 1079 n = 3; 1080 ierr = PetscOptionsGetEnumArray(NULL, NULL, "-dm_plex_box_bd", DMBoundaryTypes, (PetscEnum *) bdt, &n, &flg);CHKERRQ(ierr); 1081 if (flg && (n != dim)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Box boundary types had %D values, should have been %D", n, dim); 1082 ierr = PetscOptionsGetBool(NULL, NULL, "-dm_plex_box_interpolate", &interpolate, &flg);CHKERRQ(ierr); 1083 1084 if (dim == 1) {ierr = DMPlexCreateLineMesh_Internal(comm, fac[0], low[0], upp[0], bdt[0], dm);CHKERRQ(ierr);} 1085 else if (simplex) {ierr = DMPlexCreateBoxMesh_Simplex_Internal(comm, dim, fac, low, upp, bdt, interpolate, dm);CHKERRQ(ierr);} 1086 else {ierr = DMPlexCreateBoxMesh_Tensor_Internal(comm, dim, fac, low, upp, bdt, interpolate, dm);CHKERRQ(ierr);} 1087 PetscFunctionReturn(0); 1088 } 1089 1090 /*@ 1091 DMPlexCreateWedgeBoxMesh - Creates a 3-D mesh tesselating the (x,y) plane and extruding in the third direction using wedge cells. 1092 1093 Collective 1094 1095 Input Parameters: 1096 + comm - The communicator for the DM object 1097 . faces - Number of faces per dimension, or NULL for (1, 1, 1) 1098 . lower - The lower left corner, or NULL for (0, 0, 0) 1099 . upper - The upper right corner, or NULL for (1, 1, 1) 1100 . periodicity - The boundary type for the X,Y,Z direction, or NULL for DM_BOUNDARY_NONE 1101 . orderHeight - If PETSC_TRUE, orders the extruded cells in the height first. Otherwise, orders the cell on the layers first 1102 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 1103 1104 Output Parameter: 1105 . dm - The DM object 1106 1107 Level: beginner 1108 1109 .seealso: DMPlexCreateHexCylinderMesh(), DMPlexCreateWedgeCylinderMesh(), DMPlexExtrude(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 1110 @*/ 1111 PetscErrorCode DMPlexCreateWedgeBoxMesh(MPI_Comm comm, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool orderHeight, PetscBool interpolate, DM *dm) 1112 { 1113 DM bdm, botdm; 1114 PetscInt i; 1115 PetscInt fac[3] = {0, 0, 0}; 1116 PetscReal low[3] = {0, 0, 0}; 1117 PetscReal upp[3] = {1, 1, 1}; 1118 DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE}; 1119 PetscErrorCode ierr; 1120 1121 PetscFunctionBegin; 1122 for (i = 0; i < 3; ++i) fac[i] = faces ? (faces[i] > 0 ? faces[i] : 1) : 1; 1123 if (lower) for (i = 0; i < 3; ++i) low[i] = lower[i]; 1124 if (upper) for (i = 0; i < 3; ++i) upp[i] = upper[i]; 1125 if (periodicity) for (i = 0; i < 3; ++i) bdt[i] = periodicity[i]; 1126 for (i = 0; i < 3; ++i) if (bdt[i] != DM_BOUNDARY_NONE) SETERRQ(comm, PETSC_ERR_SUP, "Periodicity not yet supported"); 1127 1128 ierr = DMCreate(comm, &bdm);CHKERRQ(ierr); 1129 ierr = DMSetType(bdm, DMPLEX);CHKERRQ(ierr); 1130 ierr = DMSetDimension(bdm, 1);CHKERRQ(ierr); 1131 ierr = DMSetCoordinateDim(bdm, 2);CHKERRQ(ierr); 1132 ierr = DMPlexCreateSquareBoundary(bdm, low, upp, fac);CHKERRQ(ierr); 1133 ierr = DMPlexGenerate(bdm, NULL, PETSC_FALSE, &botdm);CHKERRQ(ierr); 1134 ierr = DMDestroy(&bdm);CHKERRQ(ierr); 1135 ierr = DMPlexExtrude(botdm, fac[2], upp[2] - low[2], orderHeight, NULL, interpolate, dm);CHKERRQ(ierr); 1136 if (low[2] != 0.0) { 1137 Vec v; 1138 PetscScalar *x; 1139 PetscInt cDim, n; 1140 1141 ierr = DMGetCoordinatesLocal(*dm, &v);CHKERRQ(ierr); 1142 ierr = VecGetBlockSize(v, &cDim);CHKERRQ(ierr); 1143 ierr = VecGetLocalSize(v, &n);CHKERRQ(ierr); 1144 ierr = VecGetArray(v, &x);CHKERRQ(ierr); 1145 x += cDim; 1146 for (i=0; i<n; i+=cDim) x[i] += low[2]; 1147 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1148 ierr = DMSetCoordinatesLocal(*dm, v);CHKERRQ(ierr); 1149 } 1150 ierr = DMDestroy(&botdm);CHKERRQ(ierr); 1151 PetscFunctionReturn(0); 1152 } 1153 1154 /*@C 1155 DMPlexExtrude - Creates a (d+1)-D mesh by extruding a d-D mesh in the normal direction using prismatic cells. 1156 1157 Collective on idm 1158 1159 Input Parameters: 1160 + idm - The mesh to be extruded 1161 . layers - The number of layers, or PETSC_DETERMINE to use the default 1162 . height - The total height of the extrusion, or PETSC_DETERMINE to use the default 1163 . orderHeight - If PETSC_TRUE, orders the extruded cells in the height first. Otherwise, orders the cell on the layers first 1164 . extNormal - The normal direction in which the mesh should be extruded, or NULL to extrude using the surface normal 1165 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 1166 1167 Output Parameter: 1168 . dm - The DM object 1169 1170 Notes: 1171 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. 1172 1173 Options Database Keys: 1174 + -dm_plex_extrude_layers <k> - Sets the number of layers k 1175 . -dm_plex_extrude_height <h> - Sets the total height of the extrusion 1176 . -dm_plex_extrude_heights <h0,h1,...> - Sets the height of each layer 1177 . -dm_plex_extrude_order_height - If true, order cells by height first 1178 - -dm_plex_extrude_normal <n0,...,nd> - Sets the normal vector along which to extrude 1179 1180 Level: advanced 1181 1182 .seealso: DMPlexCreateWedgeCylinderMesh(), DMPlexCreateWedgeBoxMesh(), DMSetType(), DMCreate() 1183 @*/ 1184 PetscErrorCode DMPlexExtrude(DM idm, PetscInt layers, PetscReal height, PetscBool orderHeight, const PetscReal extNormal[], PetscBool interpolate, DM* dm) 1185 { 1186 PetscScalar *coordsB; 1187 const PetscScalar *coordsA; 1188 PetscReal *normals = NULL, *heights = NULL; 1189 PetscReal clNormal[3]; 1190 Vec coordinatesA, coordinatesB; 1191 PetscSection coordSectionA, coordSectionB; 1192 PetscInt dim, cDim, cDimB, c, l, v, coordSize, *newCone, nl; 1193 PetscInt cStart, cEnd, vStart, vEnd, cellV, numCells, numVertices; 1194 const char *prefix; 1195 PetscBool haveCLNormal, flg; 1196 PetscErrorCode ierr; 1197 1198 PetscFunctionBegin; 1199 PetscValidHeaderSpecific(idm, DM_CLASSID, 1); 1200 PetscValidLogicalCollectiveInt(idm, layers, 2); 1201 PetscValidLogicalCollectiveReal(idm, height, 3); 1202 PetscValidLogicalCollectiveBool(idm, interpolate, 4); 1203 ierr = DMGetDimension(idm, &dim);CHKERRQ(ierr); 1204 ierr = DMGetCoordinateDim(idm, &cDim);CHKERRQ(ierr); 1205 cDimB = cDim == dim ? cDim+1 : cDim; 1206 if (dim < 1 || dim > 3) SETERRQ1(PetscObjectComm((PetscObject)idm), PETSC_ERR_SUP, "Support for dimension %D not coded", dim); 1207 1208 ierr = PetscObjectGetOptionsPrefix((PetscObject) idm, &prefix);CHKERRQ(ierr); 1209 if (layers < 0) layers = 1; 1210 ierr = PetscOptionsGetInt(NULL, prefix, "-dm_plex_extrude_layers", &layers, NULL);CHKERRQ(ierr); 1211 if (layers <= 0) SETERRQ1(PetscObjectComm((PetscObject) idm), PETSC_ERR_ARG_OUTOFRANGE, "Number of layers %D must be positive", layers); 1212 if (height < 0.) height = 1.; 1213 ierr = PetscOptionsGetReal(NULL, prefix, "-dm_plex_extrude_height", &height, NULL);CHKERRQ(ierr); 1214 if (height <= 0.) SETERRQ1(PetscObjectComm((PetscObject) idm), PETSC_ERR_ARG_OUTOFRANGE, "Height of layers %g must be positive", (double) height); 1215 ierr = PetscMalloc1(layers, &heights);CHKERRQ(ierr); 1216 nl = layers; 1217 ierr = PetscOptionsGetRealArray(NULL, prefix, "-dm_plex_extrude_heights", heights, &nl, &flg);CHKERRQ(ierr); 1218 if (flg) { 1219 if (!nl) SETERRQ(PetscObjectComm((PetscObject) idm), PETSC_ERR_ARG_OUTOFRANGE, "Must give at least one height for -dm_plex_extrude_heights"); 1220 for (l = nl; l < layers; ++l) heights[l] = heights[l-1]; 1221 for (l = 0; l < layers; ++l) if (heights[l] <= 0.) SETERRQ2(PetscObjectComm((PetscObject) idm), PETSC_ERR_ARG_OUTOFRANGE, "Height %g of layers %D must be positive", (double) heights[l], l); 1222 } else { 1223 for (l = 0; l < layers; ++l) heights[l] = height/layers; 1224 } 1225 ierr = PetscOptionsGetBool(NULL, prefix, "-dm_plex_extrude_order_height", &orderHeight, NULL);CHKERRQ(ierr); 1226 c = 3; 1227 ierr = PetscOptionsGetRealArray(NULL, prefix, "-dm_plex_extrude_normal", clNormal, &c, &haveCLNormal);CHKERRQ(ierr); 1228 if (haveCLNormal && c != cDimB) SETERRQ2(PetscObjectComm((PetscObject)idm), PETSC_ERR_ARG_SIZ, "Input normal has size %D != %D extruded coordinate dimension", c, cDimB); 1229 1230 ierr = DMPlexGetHeightStratum(idm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1231 ierr = DMPlexGetDepthStratum(idm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1232 numCells = (cEnd - cStart)*layers; 1233 numVertices = (vEnd - vStart)*(layers+1); 1234 ierr = DMCreate(PetscObjectComm((PetscObject)idm), dm);CHKERRQ(ierr); 1235 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1236 ierr = DMSetDimension(*dm, dim+1);CHKERRQ(ierr); 1237 ierr = DMPlexSetChart(*dm, 0, numCells+numVertices);CHKERRQ(ierr); 1238 /* Must create the celltype label here so that we do not automatically try to compute the types */ 1239 ierr = DMCreateLabel(*dm, "celltype");CHKERRQ(ierr); 1240 for (c = cStart, cellV = 0; c < cEnd; ++c) { 1241 DMPolytopeType ct, nct; 1242 PetscInt *closure = NULL; 1243 PetscInt closureSize, numCorners = 0; 1244 1245 ierr = DMPlexGetCellType(idm, c, &ct);CHKERRQ(ierr); 1246 switch (ct) { 1247 case DM_POLYTOPE_SEGMENT: nct = DM_POLYTOPE_SEG_PRISM_TENSOR;break; 1248 case DM_POLYTOPE_TRIANGLE: nct = DM_POLYTOPE_TRI_PRISM_TENSOR;break; 1249 case DM_POLYTOPE_QUADRILATERAL: nct = DM_POLYTOPE_QUAD_PRISM_TENSOR;break; 1250 default: nct = DM_POLYTOPE_UNKNOWN; 1251 } 1252 ierr = DMPlexGetTransitiveClosure(idm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1253 for (v = 0; v < closureSize*2; v += 2) if ((closure[v] >= vStart) && (closure[v] < vEnd)) numCorners++; 1254 ierr = DMPlexRestoreTransitiveClosure(idm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1255 for (l = 0; l < layers; ++l) { 1256 const PetscInt cell = orderHeight ? layers*(c - cStart) + l : l*(cEnd - cStart) + c - cStart; 1257 1258 ierr = DMPlexSetConeSize(*dm, cell, 2*numCorners);CHKERRQ(ierr); 1259 ierr = DMPlexSetCellType(*dm, cell, nct);CHKERRQ(ierr); 1260 } 1261 cellV = PetscMax(numCorners,cellV); 1262 } 1263 ierr = DMSetUp(*dm);CHKERRQ(ierr); 1264 1265 if (dim != cDim && !(extNormal || haveCLNormal)) {ierr = PetscCalloc1(cDim*(vEnd - vStart), &normals);CHKERRQ(ierr);} 1266 ierr = PetscMalloc1(3*cellV,&newCone);CHKERRQ(ierr); 1267 for (c = cStart; c < cEnd; ++c) { 1268 PetscInt *closure = NULL; 1269 PetscInt closureSize, numCorners = 0, l; 1270 PetscReal normal[3] = {0, 0, 0}; 1271 1272 if (normals) {ierr = DMPlexComputeCellGeometryFVM(idm, c, NULL, NULL, normal);CHKERRQ(ierr);} 1273 ierr = DMPlexGetTransitiveClosure(idm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1274 for (v = 0; v < closureSize*2; v += 2) { 1275 if ((closure[v] >= vStart) && (closure[v] < vEnd)) { 1276 PetscInt d; 1277 1278 newCone[numCorners++] = closure[v] - vStart; 1279 if (normals) {for (d = 0; d < cDim; ++d) normals[cDim*(closure[v]-vStart)+d] += normal[d];} 1280 } 1281 } 1282 ierr = DMPlexRestoreTransitiveClosure(idm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1283 for (l = 0; l < layers; ++l) { 1284 PetscInt i; 1285 1286 for (i = 0; i < numCorners; ++i) { 1287 newCone[ numCorners + i] = orderHeight ? (layers+1)*newCone[i] + l + numCells : l*(vEnd - vStart) + newCone[i] + numCells; 1288 newCone[2*numCorners + i] = orderHeight ? (layers+1)*newCone[i] + l + 1 + numCells : (l+1)*(vEnd - vStart) + newCone[i] + numCells; 1289 } 1290 ierr = DMPlexSetCone(*dm, orderHeight ? layers*(c - cStart) + l : l*(cEnd - cStart) + c - cStart, newCone + numCorners);CHKERRQ(ierr); 1291 } 1292 } 1293 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1294 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1295 ierr = PetscFree(newCone);CHKERRQ(ierr); 1296 1297 ierr = DMGetCoordinateSection(*dm, &coordSectionB);CHKERRQ(ierr); 1298 ierr = PetscSectionSetNumFields(coordSectionB, 1);CHKERRQ(ierr); 1299 ierr = PetscSectionSetFieldComponents(coordSectionB, 0, cDimB);CHKERRQ(ierr); 1300 ierr = PetscSectionSetChart(coordSectionB, numCells, numCells+numVertices);CHKERRQ(ierr); 1301 for (v = numCells; v < numCells+numVertices; ++v) { 1302 ierr = PetscSectionSetDof(coordSectionB, v, cDimB);CHKERRQ(ierr); 1303 ierr = PetscSectionSetFieldDof(coordSectionB, v, 0, cDimB);CHKERRQ(ierr); 1304 ierr = DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT);CHKERRQ(ierr); 1305 } 1306 ierr = PetscSectionSetUp(coordSectionB);CHKERRQ(ierr); 1307 ierr = PetscSectionGetStorageSize(coordSectionB, &coordSize);CHKERRQ(ierr); 1308 ierr = VecCreate(PETSC_COMM_SELF, &coordinatesB);CHKERRQ(ierr); 1309 ierr = PetscObjectSetName((PetscObject) coordinatesB, "coordinates");CHKERRQ(ierr); 1310 ierr = VecSetSizes(coordinatesB, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1311 ierr = VecSetBlockSize(coordinatesB, cDimB);CHKERRQ(ierr); 1312 ierr = VecSetType(coordinatesB,VECSTANDARD);CHKERRQ(ierr); 1313 1314 ierr = DMGetCoordinateSection(idm, &coordSectionA);CHKERRQ(ierr); 1315 ierr = DMGetCoordinatesLocal(idm, &coordinatesA);CHKERRQ(ierr); 1316 ierr = VecGetArray(coordinatesB, &coordsB);CHKERRQ(ierr); 1317 ierr = VecGetArrayRead(coordinatesA, &coordsA);CHKERRQ(ierr); 1318 for (v = vStart; v < vEnd; ++v) { 1319 const PetscScalar *cptr; 1320 PetscReal ones2[2] = { 0., 1.}, ones3[3] = { 0., 0., 1.}; 1321 PetscReal normal[3]; 1322 PetscReal norm; 1323 PetscInt offA, d, cDimA = cDim; 1324 1325 if (normals) {for (d = 0; d < cDimB; ++d) normal[d] = normals[cDimB*(v - vStart)+d];} 1326 else if (haveCLNormal) {for (d = 0; d < cDimB; ++d) normal[d] = clNormal[d];} 1327 else if (extNormal) {for (d = 0; d < cDimB; ++d) normal[d] = extNormal[d];} 1328 else if (cDimB == 2) {for (d = 0; d < cDimB; ++d) normal[d] = ones2[d];} 1329 else if (cDimB == 3) {for (d = 0; d < cDimB; ++d) normal[d] = ones3[d];} 1330 else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unable to determine normal for extrusion"); 1331 for (d = 0, norm = 0.0; d < cDimB; ++d) norm += normal[d]*normal[d]; 1332 for (d = 0; d < cDimB; ++d) normal[d] *= 1./PetscSqrtReal(norm); 1333 1334 ierr = PetscSectionGetOffset(coordSectionA, v, &offA);CHKERRQ(ierr); 1335 cptr = coordsA + offA; 1336 for (l = 0; l <= layers; ++l) { 1337 PetscInt offB, d, newV; 1338 1339 newV = orderHeight ? (layers+1)*(v -vStart) + l + numCells : (vEnd -vStart)*l + (v -vStart) + numCells; 1340 ierr = PetscSectionGetOffset(coordSectionB, newV, &offB);CHKERRQ(ierr); 1341 for (d = 0; d < cDimA; ++d) { coordsB[offB+d] = cptr[d]; } 1342 for (d = 0; d < cDimB; ++d) { coordsB[offB+d] += l ? normal[d]*heights[l-1] : 0.0; } 1343 cptr = coordsB + offB; 1344 cDimA = cDimB; 1345 } 1346 } 1347 ierr = VecRestoreArrayRead(coordinatesA, &coordsA);CHKERRQ(ierr); 1348 ierr = VecRestoreArray(coordinatesB, &coordsB);CHKERRQ(ierr); 1349 ierr = DMSetCoordinatesLocal(*dm, coordinatesB);CHKERRQ(ierr); 1350 ierr = VecDestroy(&coordinatesB);CHKERRQ(ierr); 1351 ierr = PetscFree(normals);CHKERRQ(ierr); 1352 ierr = PetscFree(heights);CHKERRQ(ierr); 1353 if (interpolate) { 1354 DM idm; 1355 1356 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1357 ierr = DMPlexCopyCoordinates(*dm, idm);CHKERRQ(ierr); 1358 ierr = DMDestroy(dm);CHKERRQ(ierr); 1359 *dm = idm; 1360 } 1361 PetscFunctionReturn(0); 1362 } 1363 1364 /*@C 1365 DMPlexSetOptionsPrefix - Sets the prefix used for searching for all DM options in the database. 1366 1367 Logically Collective on dm 1368 1369 Input Parameters: 1370 + dm - the DM context 1371 - prefix - the prefix to prepend to all option names 1372 1373 Notes: 1374 A hyphen (-) must NOT be given at the beginning of the prefix name. 1375 The first character of all runtime options is AUTOMATICALLY the hyphen. 1376 1377 Level: advanced 1378 1379 .seealso: SNESSetFromOptions() 1380 @*/ 1381 PetscErrorCode DMPlexSetOptionsPrefix(DM dm, const char prefix[]) 1382 { 1383 DM_Plex *mesh = (DM_Plex *) dm->data; 1384 PetscErrorCode ierr; 1385 1386 PetscFunctionBegin; 1387 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1388 ierr = PetscObjectSetOptionsPrefix((PetscObject) dm, prefix);CHKERRQ(ierr); 1389 ierr = PetscObjectSetOptionsPrefix((PetscObject) mesh->partitioner, prefix);CHKERRQ(ierr); 1390 PetscFunctionReturn(0); 1391 } 1392 1393 /*@ 1394 DMPlexCreateHexCylinderMesh - Creates a mesh on the tensor product of the unit interval with the circle (cylinder) using hexahedra. 1395 1396 Collective 1397 1398 Input Parameters: 1399 + comm - The communicator for the DM object 1400 . numRefine - The number of regular refinements to the basic 5 cell structure 1401 - periodicZ - The boundary type for the Z direction 1402 1403 Output Parameter: 1404 . dm - The DM object 1405 1406 Note: Here is the output numbering looking from the bottom of the cylinder: 1407 $ 17-----14 1408 $ | | 1409 $ | 2 | 1410 $ | | 1411 $ 17-----8-----7-----14 1412 $ | | | | 1413 $ | 3 | 0 | 1 | 1414 $ | | | | 1415 $ 19-----5-----6-----13 1416 $ | | 1417 $ | 4 | 1418 $ | | 1419 $ 19-----13 1420 $ 1421 $ and up through the top 1422 $ 1423 $ 18-----16 1424 $ | | 1425 $ | 2 | 1426 $ | | 1427 $ 18----10----11-----16 1428 $ | | | | 1429 $ | 3 | 0 | 1 | 1430 $ | | | | 1431 $ 20-----9----12-----15 1432 $ | | 1433 $ | 4 | 1434 $ | | 1435 $ 20-----15 1436 1437 Level: beginner 1438 1439 .seealso: DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 1440 @*/ 1441 PetscErrorCode DMPlexCreateHexCylinderMesh(MPI_Comm comm, PetscInt numRefine, DMBoundaryType periodicZ, DM *dm) 1442 { 1443 const PetscInt dim = 3; 1444 PetscInt numCells, numVertices, r; 1445 PetscMPIInt rank; 1446 PetscErrorCode ierr; 1447 1448 PetscFunctionBegin; 1449 PetscValidPointer(dm, 4); 1450 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 1451 if (numRefine < 0) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Number of refinements %D cannot be negative", numRefine); 1452 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1453 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1454 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1455 /* Create topology */ 1456 { 1457 PetscInt cone[8], c; 1458 1459 numCells = !rank ? 5 : 0; 1460 numVertices = !rank ? 16 : 0; 1461 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1462 numCells *= 3; 1463 numVertices = !rank ? 24 : 0; 1464 } 1465 ierr = DMPlexSetChart(*dm, 0, numCells+numVertices);CHKERRQ(ierr); 1466 for (c = 0; c < numCells; c++) {ierr = DMPlexSetConeSize(*dm, c, 8);CHKERRQ(ierr);} 1467 ierr = DMSetUp(*dm);CHKERRQ(ierr); 1468 if (!rank) { 1469 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1470 cone[0] = 15; cone[1] = 18; cone[2] = 17; cone[3] = 16; 1471 cone[4] = 31; cone[5] = 32; cone[6] = 33; cone[7] = 34; 1472 ierr = DMPlexSetCone(*dm, 0, cone);CHKERRQ(ierr); 1473 cone[0] = 16; cone[1] = 17; cone[2] = 24; cone[3] = 23; 1474 cone[4] = 32; cone[5] = 36; cone[6] = 37; cone[7] = 33; /* 22 25 26 21 */ 1475 ierr = DMPlexSetCone(*dm, 1, cone);CHKERRQ(ierr); 1476 cone[0] = 18; cone[1] = 27; cone[2] = 24; cone[3] = 17; 1477 cone[4] = 34; cone[5] = 33; cone[6] = 37; cone[7] = 38; 1478 ierr = DMPlexSetCone(*dm, 2, cone);CHKERRQ(ierr); 1479 cone[0] = 29; cone[1] = 27; cone[2] = 18; cone[3] = 15; 1480 cone[4] = 35; cone[5] = 31; cone[6] = 34; cone[7] = 38; 1481 ierr = DMPlexSetCone(*dm, 3, cone);CHKERRQ(ierr); 1482 cone[0] = 29; cone[1] = 15; cone[2] = 16; cone[3] = 23; 1483 cone[4] = 35; cone[5] = 36; cone[6] = 32; cone[7] = 31; 1484 ierr = DMPlexSetCone(*dm, 4, cone);CHKERRQ(ierr); 1485 1486 cone[0] = 31; cone[1] = 34; cone[2] = 33; cone[3] = 32; 1487 cone[4] = 19; cone[5] = 22; cone[6] = 21; cone[7] = 20; 1488 ierr = DMPlexSetCone(*dm, 5, cone);CHKERRQ(ierr); 1489 cone[0] = 32; cone[1] = 33; cone[2] = 37; cone[3] = 36; 1490 cone[4] = 22; cone[5] = 25; cone[6] = 26; cone[7] = 21; 1491 ierr = DMPlexSetCone(*dm, 6, cone);CHKERRQ(ierr); 1492 cone[0] = 34; cone[1] = 38; cone[2] = 37; cone[3] = 33; 1493 cone[4] = 20; cone[5] = 21; cone[6] = 26; cone[7] = 28; 1494 ierr = DMPlexSetCone(*dm, 7, cone);CHKERRQ(ierr); 1495 cone[0] = 35; cone[1] = 38; cone[2] = 34; cone[3] = 31; 1496 cone[4] = 30; cone[5] = 19; cone[6] = 20; cone[7] = 28; 1497 ierr = DMPlexSetCone(*dm, 8, cone);CHKERRQ(ierr); 1498 cone[0] = 35; cone[1] = 31; cone[2] = 32; cone[3] = 36; 1499 cone[4] = 30; cone[5] = 25; cone[6] = 22; cone[7] = 19; 1500 ierr = DMPlexSetCone(*dm, 9, cone);CHKERRQ(ierr); 1501 1502 cone[0] = 19; cone[1] = 20; cone[2] = 21; cone[3] = 22; 1503 cone[4] = 15; cone[5] = 16; cone[6] = 17; cone[7] = 18; 1504 ierr = DMPlexSetCone(*dm, 10, cone);CHKERRQ(ierr); 1505 cone[0] = 22; cone[1] = 21; cone[2] = 26; cone[3] = 25; 1506 cone[4] = 16; cone[5] = 23; cone[6] = 24; cone[7] = 17; 1507 ierr = DMPlexSetCone(*dm, 11, cone);CHKERRQ(ierr); 1508 cone[0] = 20; cone[1] = 28; cone[2] = 26; cone[3] = 21; 1509 cone[4] = 18; cone[5] = 17; cone[6] = 24; cone[7] = 27; 1510 ierr = DMPlexSetCone(*dm, 12, cone);CHKERRQ(ierr); 1511 cone[0] = 30; cone[1] = 28; cone[2] = 20; cone[3] = 19; 1512 cone[4] = 29; cone[5] = 15; cone[6] = 18; cone[7] = 27; 1513 ierr = DMPlexSetCone(*dm, 13, cone);CHKERRQ(ierr); 1514 cone[0] = 30; cone[1] = 19; cone[2] = 22; cone[3] = 25; 1515 cone[4] = 29; cone[5] = 23; cone[6] = 16; cone[7] = 15; 1516 ierr = DMPlexSetCone(*dm, 14, cone);CHKERRQ(ierr); 1517 } else { 1518 cone[0] = 5; cone[1] = 8; cone[2] = 7; cone[3] = 6; 1519 cone[4] = 9; cone[5] = 12; cone[6] = 11; cone[7] = 10; 1520 ierr = DMPlexSetCone(*dm, 0, cone);CHKERRQ(ierr); 1521 cone[0] = 6; cone[1] = 7; cone[2] = 14; cone[3] = 13; 1522 cone[4] = 12; cone[5] = 15; cone[6] = 16; cone[7] = 11; 1523 ierr = DMPlexSetCone(*dm, 1, cone);CHKERRQ(ierr); 1524 cone[0] = 8; cone[1] = 17; cone[2] = 14; cone[3] = 7; 1525 cone[4] = 10; cone[5] = 11; cone[6] = 16; cone[7] = 18; 1526 ierr = DMPlexSetCone(*dm, 2, cone);CHKERRQ(ierr); 1527 cone[0] = 19; cone[1] = 17; cone[2] = 8; cone[3] = 5; 1528 cone[4] = 20; cone[5] = 9; cone[6] = 10; cone[7] = 18; 1529 ierr = DMPlexSetCone(*dm, 3, cone);CHKERRQ(ierr); 1530 cone[0] = 19; cone[1] = 5; cone[2] = 6; cone[3] = 13; 1531 cone[4] = 20; cone[5] = 15; cone[6] = 12; cone[7] = 9; 1532 ierr = DMPlexSetCone(*dm, 4, cone);CHKERRQ(ierr); 1533 } 1534 } 1535 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1536 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1537 } 1538 /* Interpolate */ 1539 { 1540 DM idm; 1541 1542 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1543 ierr = DMDestroy(dm);CHKERRQ(ierr); 1544 *dm = idm; 1545 } 1546 /* Create cube geometry */ 1547 { 1548 Vec coordinates; 1549 PetscSection coordSection; 1550 PetscScalar *coords; 1551 PetscInt coordSize, v; 1552 const PetscReal dis = 1.0/PetscSqrtReal(2.0); 1553 const PetscReal ds2 = dis/2.0; 1554 1555 /* Build coordinates */ 1556 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 1557 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1558 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 1559 ierr = PetscSectionSetChart(coordSection, numCells, numCells+numVertices);CHKERRQ(ierr); 1560 for (v = numCells; v < numCells+numVertices; ++v) { 1561 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 1562 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 1563 } 1564 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1565 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1566 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 1567 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1568 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1569 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 1570 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1571 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1572 if (!rank) { 1573 coords[0*dim+0] = -ds2; coords[0*dim+1] = -ds2; coords[0*dim+2] = 0.0; 1574 coords[1*dim+0] = ds2; coords[1*dim+1] = -ds2; coords[1*dim+2] = 0.0; 1575 coords[2*dim+0] = ds2; coords[2*dim+1] = ds2; coords[2*dim+2] = 0.0; 1576 coords[3*dim+0] = -ds2; coords[3*dim+1] = ds2; coords[3*dim+2] = 0.0; 1577 coords[4*dim+0] = -ds2; coords[4*dim+1] = -ds2; coords[4*dim+2] = 1.0; 1578 coords[5*dim+0] = -ds2; coords[5*dim+1] = ds2; coords[5*dim+2] = 1.0; 1579 coords[6*dim+0] = ds2; coords[6*dim+1] = ds2; coords[6*dim+2] = 1.0; 1580 coords[7*dim+0] = ds2; coords[7*dim+1] = -ds2; coords[7*dim+2] = 1.0; 1581 coords[ 8*dim+0] = dis; coords[ 8*dim+1] = -dis; coords[ 8*dim+2] = 0.0; 1582 coords[ 9*dim+0] = dis; coords[ 9*dim+1] = dis; coords[ 9*dim+2] = 0.0; 1583 coords[10*dim+0] = dis; coords[10*dim+1] = -dis; coords[10*dim+2] = 1.0; 1584 coords[11*dim+0] = dis; coords[11*dim+1] = dis; coords[11*dim+2] = 1.0; 1585 coords[12*dim+0] = -dis; coords[12*dim+1] = dis; coords[12*dim+2] = 0.0; 1586 coords[13*dim+0] = -dis; coords[13*dim+1] = dis; coords[13*dim+2] = 1.0; 1587 coords[14*dim+0] = -dis; coords[14*dim+1] = -dis; coords[14*dim+2] = 0.0; 1588 coords[15*dim+0] = -dis; coords[15*dim+1] = -dis; coords[15*dim+2] = 1.0; 1589 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1590 /* 15 31 19 */ coords[16*dim+0] = -ds2; coords[16*dim+1] = -ds2; coords[16*dim+2] = 0.5; 1591 /* 16 32 22 */ coords[17*dim+0] = ds2; coords[17*dim+1] = -ds2; coords[17*dim+2] = 0.5; 1592 /* 17 33 21 */ coords[18*dim+0] = ds2; coords[18*dim+1] = ds2; coords[18*dim+2] = 0.5; 1593 /* 18 34 20 */ coords[19*dim+0] = -ds2; coords[19*dim+1] = ds2; coords[19*dim+2] = 0.5; 1594 /* 29 35 30 */ coords[20*dim+0] = -dis; coords[20*dim+1] = -dis; coords[20*dim+2] = 0.5; 1595 /* 23 36 25 */ coords[21*dim+0] = dis; coords[21*dim+1] = -dis; coords[21*dim+2] = 0.5; 1596 /* 24 37 26 */ coords[22*dim+0] = dis; coords[22*dim+1] = dis; coords[22*dim+2] = 0.5; 1597 /* 27 38 28 */ coords[23*dim+0] = -dis; coords[23*dim+1] = dis; coords[23*dim+2] = 0.5; 1598 } 1599 } 1600 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1601 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 1602 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1603 } 1604 /* Create periodicity */ 1605 if (periodicZ == DM_BOUNDARY_PERIODIC || periodicZ == DM_BOUNDARY_TWIST) { 1606 PetscReal L[3]; 1607 PetscReal maxCell[3]; 1608 DMBoundaryType bdType[3]; 1609 PetscReal lower[3] = {0.0, 0.0, 0.0}; 1610 PetscReal upper[3] = {1.0, 1.0, 1.5}; 1611 PetscInt i, numZCells = 3; 1612 1613 bdType[0] = DM_BOUNDARY_NONE; 1614 bdType[1] = DM_BOUNDARY_NONE; 1615 bdType[2] = periodicZ; 1616 for (i = 0; i < dim; i++) { 1617 L[i] = upper[i] - lower[i]; 1618 maxCell[i] = 1.1 * (L[i] / numZCells); 1619 } 1620 ierr = DMSetPeriodicity(*dm, PETSC_TRUE, maxCell, L, bdType);CHKERRQ(ierr); 1621 } 1622 /* Refine topology */ 1623 for (r = 0; r < numRefine; ++r) { 1624 DM rdm = NULL; 1625 1626 ierr = DMRefine(*dm, comm, &rdm);CHKERRQ(ierr); 1627 ierr = DMDestroy(dm);CHKERRQ(ierr); 1628 *dm = rdm; 1629 } 1630 /* Remap geometry to cylinder 1631 Interior square: Linear interpolation is correct 1632 The other cells all have vertices on rays from the origin. We want to uniformly expand the spacing 1633 such that the last vertex is on the unit circle. So the closest and farthest vertices are at distance 1634 1635 phi = arctan(y/x) 1636 d_close = sqrt(1/8 + 1/4 sin^2(phi)) 1637 d_far = sqrt(1/2 + sin^2(phi)) 1638 1639 so we remap them using 1640 1641 x_new = x_close + (x - x_close) (1 - d_close) / (d_far - d_close) 1642 y_new = y_close + (y - y_close) (1 - d_close) / (d_far - d_close) 1643 1644 If pi/4 < phi < 3pi/4 or -3pi/4 < phi < -pi/4, then we switch x and y. 1645 */ 1646 { 1647 Vec coordinates; 1648 PetscSection coordSection; 1649 PetscScalar *coords; 1650 PetscInt vStart, vEnd, v; 1651 const PetscReal dis = 1.0/PetscSqrtReal(2.0); 1652 const PetscReal ds2 = 0.5*dis; 1653 1654 ierr = DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1655 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 1656 ierr = DMGetCoordinatesLocal(*dm, &coordinates);CHKERRQ(ierr); 1657 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1658 for (v = vStart; v < vEnd; ++v) { 1659 PetscReal phi, sinp, cosp, dc, df, x, y, xc, yc; 1660 PetscInt off; 1661 1662 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr); 1663 if ((PetscAbsScalar(coords[off+0]) <= ds2) && (PetscAbsScalar(coords[off+1]) <= ds2)) continue; 1664 x = PetscRealPart(coords[off]); 1665 y = PetscRealPart(coords[off+1]); 1666 phi = PetscAtan2Real(y, x); 1667 sinp = PetscSinReal(phi); 1668 cosp = PetscCosReal(phi); 1669 if ((PetscAbsReal(phi) > PETSC_PI/4.0) && (PetscAbsReal(phi) < 3.0*PETSC_PI/4.0)) { 1670 dc = PetscAbsReal(ds2/sinp); 1671 df = PetscAbsReal(dis/sinp); 1672 xc = ds2*x/PetscAbsReal(y); 1673 yc = ds2*PetscSignReal(y); 1674 } else { 1675 dc = PetscAbsReal(ds2/cosp); 1676 df = PetscAbsReal(dis/cosp); 1677 xc = ds2*PetscSignReal(x); 1678 yc = ds2*y/PetscAbsReal(x); 1679 } 1680 coords[off+0] = xc + (coords[off+0] - xc)*(1.0 - dc)/(df - dc); 1681 coords[off+1] = yc + (coords[off+1] - yc)*(1.0 - dc)/(df - dc); 1682 } 1683 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1684 if (periodicZ == DM_BOUNDARY_PERIODIC || periodicZ == DM_BOUNDARY_TWIST) { 1685 ierr = DMLocalizeCoordinates(*dm);CHKERRQ(ierr); 1686 } 1687 } 1688 PetscFunctionReturn(0); 1689 } 1690 1691 /*@ 1692 DMPlexCreateWedgeCylinderMesh - Creates a mesh on the tensor product of the unit interval with the circle (cylinder) using wedges. 1693 1694 Collective 1695 1696 Input Parameters: 1697 + comm - The communicator for the DM object 1698 . n - The number of wedges around the origin 1699 - interpolate - Create edges and faces 1700 1701 Output Parameter: 1702 . dm - The DM object 1703 1704 Level: beginner 1705 1706 .seealso: DMPlexCreateHexCylinderMesh(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 1707 @*/ 1708 PetscErrorCode DMPlexCreateWedgeCylinderMesh(MPI_Comm comm, PetscInt n, PetscBool interpolate, DM *dm) 1709 { 1710 const PetscInt dim = 3; 1711 PetscInt numCells, numVertices, v; 1712 PetscMPIInt rank; 1713 PetscErrorCode ierr; 1714 1715 PetscFunctionBegin; 1716 PetscValidPointer(dm, 3); 1717 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 1718 if (n < 0) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Number of wedges %D cannot be negative", n); 1719 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1720 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1721 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1722 /* Must create the celltype label here so that we do not automatically try to compute the types */ 1723 ierr = DMCreateLabel(*dm, "celltype");CHKERRQ(ierr); 1724 /* Create topology */ 1725 { 1726 PetscInt cone[6], c; 1727 1728 numCells = !rank ? n : 0; 1729 numVertices = !rank ? 2*(n+1) : 0; 1730 ierr = DMPlexSetChart(*dm, 0, numCells+numVertices);CHKERRQ(ierr); 1731 for (c = 0; c < numCells; c++) {ierr = DMPlexSetConeSize(*dm, c, 6);CHKERRQ(ierr);} 1732 ierr = DMSetUp(*dm);CHKERRQ(ierr); 1733 for (c = 0; c < numCells; c++) { 1734 cone[0] = c+n*1; cone[1] = (c+1)%n+n*1; cone[2] = 0+3*n; 1735 cone[3] = c+n*2; cone[4] = (c+1)%n+n*2; cone[5] = 1+3*n; 1736 ierr = DMPlexSetCone(*dm, c, cone);CHKERRQ(ierr); 1737 ierr = DMPlexSetCellType(*dm, c, DM_POLYTOPE_TRI_PRISM_TENSOR);CHKERRQ(ierr); 1738 } 1739 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1740 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1741 } 1742 for (v = numCells; v < numCells+numVertices; ++v) { 1743 ierr = DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT);CHKERRQ(ierr); 1744 } 1745 /* Interpolate */ 1746 if (interpolate) { 1747 DM idm; 1748 1749 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1750 ierr = DMDestroy(dm);CHKERRQ(ierr); 1751 *dm = idm; 1752 } 1753 /* Create cylinder geometry */ 1754 { 1755 Vec coordinates; 1756 PetscSection coordSection; 1757 PetscScalar *coords; 1758 PetscInt coordSize, c; 1759 1760 /* Build coordinates */ 1761 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 1762 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1763 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 1764 ierr = PetscSectionSetChart(coordSection, numCells, numCells+numVertices);CHKERRQ(ierr); 1765 for (v = numCells; v < numCells+numVertices; ++v) { 1766 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 1767 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 1768 } 1769 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1770 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1771 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 1772 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1773 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1774 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 1775 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1776 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1777 for (c = 0; c < numCells; c++) { 1778 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; 1779 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; 1780 } 1781 if (!rank) { 1782 coords[(2*n+0)*dim+0] = 0.0; coords[(2*n+0)*dim+1] = 0.0; coords[(2*n+0)*dim+2] = 1.0; 1783 coords[(2*n+1)*dim+0] = 0.0; coords[(2*n+1)*dim+1] = 0.0; coords[(2*n+1)*dim+2] = 0.0; 1784 } 1785 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1786 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 1787 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1788 } 1789 PetscFunctionReturn(0); 1790 } 1791 1792 PETSC_STATIC_INLINE PetscReal DiffNormReal(PetscInt dim, const PetscReal x[], const PetscReal y[]) 1793 { 1794 PetscReal prod = 0.0; 1795 PetscInt i; 1796 for (i = 0; i < dim; ++i) prod += PetscSqr(x[i] - y[i]); 1797 return PetscSqrtReal(prod); 1798 } 1799 PETSC_STATIC_INLINE PetscReal DotReal(PetscInt dim, const PetscReal x[], const PetscReal y[]) 1800 { 1801 PetscReal prod = 0.0; 1802 PetscInt i; 1803 for (i = 0; i < dim; ++i) prod += x[i]*y[i]; 1804 return prod; 1805 } 1806 1807 /* The first constant is the sphere radius */ 1808 static void snapToSphere(PetscInt dim, PetscInt Nf, PetscInt NfAux, 1809 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 1810 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 1811 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 1812 { 1813 PetscReal r = PetscRealPart(constants[0]); 1814 PetscReal norm2 = 0.0, fac; 1815 PetscInt n = uOff[1] - uOff[0], d; 1816 1817 for (d = 0; d < n; ++d) norm2 += PetscSqr(PetscRealPart(u[d])); 1818 fac = r/PetscSqrtReal(norm2); 1819 for (d = 0; d < n; ++d) f0[d] = u[d]*fac; 1820 } 1821 1822 /*@ 1823 DMPlexCreateSphereMesh - Creates a mesh on the d-dimensional sphere, S^d. 1824 1825 Collective 1826 1827 Input Parameters: 1828 + comm - The communicator for the DM object 1829 . dim - The dimension 1830 . simplex - Use simplices, or tensor product cells 1831 - R - The radius 1832 1833 Output Parameter: 1834 . dm - The DM object 1835 1836 Level: beginner 1837 1838 .seealso: DMPlexCreateBallMesh(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 1839 @*/ 1840 PetscErrorCode DMPlexCreateSphereMesh(MPI_Comm comm, PetscInt dim, PetscBool simplex, PetscReal R, DM *dm) 1841 { 1842 const PetscInt embedDim = dim+1; 1843 PetscSection coordSection; 1844 Vec coordinates; 1845 PetscScalar *coords; 1846 PetscReal *coordsIn; 1847 PetscInt numCells, numEdges, numVerts, firstVertex, v, firstEdge, coordSize, d, c, e; 1848 PetscMPIInt rank; 1849 PetscErrorCode ierr; 1850 1851 PetscFunctionBegin; 1852 PetscValidPointer(dm, 4); 1853 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1854 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1855 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1856 ierr = DMSetCoordinateDim(*dm, dim+1);CHKERRQ(ierr); 1857 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) *dm), &rank);CHKERRQ(ierr); 1858 switch (dim) { 1859 case 2: 1860 if (simplex) { 1861 DM idm; 1862 const PetscReal radius = PetscSqrtReal(1 + PETSC_PHI*PETSC_PHI)/(1.0 + PETSC_PHI); 1863 const PetscReal edgeLen = 2.0/(1.0 + PETSC_PHI) * (R/radius); 1864 const PetscInt degree = 5; 1865 PetscReal vertex[3] = {0.0, 1.0/(1.0 + PETSC_PHI), PETSC_PHI/(1.0 + PETSC_PHI)}; 1866 PetscInt s[3] = {1, 1, 1}; 1867 PetscInt cone[3]; 1868 PetscInt *graph, p, i, j, k; 1869 1870 vertex[0] *= R/radius; vertex[1] *= R/radius; vertex[2] *= R/radius; 1871 numCells = !rank ? 20 : 0; 1872 numVerts = !rank ? 12 : 0; 1873 firstVertex = numCells; 1874 /* Use icosahedron, which for a R-sphere has coordinates which are all cyclic permutations of 1875 1876 (0, \pm 1/\phi+1, \pm \phi/\phi+1) 1877 1878 where \phi^2 - \phi - 1 = 0, meaning \phi is the golden ratio \frac{1 + \sqrt{5}}{2}. The edge 1879 length is then given by 2/(1+\phi) = 2 * 0.38197 = 0.76393. 1880 */ 1881 /* Construct vertices */ 1882 ierr = PetscCalloc1(numVerts * embedDim, &coordsIn);CHKERRQ(ierr); 1883 if (!rank) { 1884 for (p = 0, i = 0; p < embedDim; ++p) { 1885 for (s[1] = -1; s[1] < 2; s[1] += 2) { 1886 for (s[2] = -1; s[2] < 2; s[2] += 2) { 1887 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[(d+p)%embedDim]*vertex[(d+p)%embedDim]; 1888 ++i; 1889 } 1890 } 1891 } 1892 } 1893 /* Construct graph */ 1894 ierr = PetscCalloc1(numVerts * numVerts, &graph);CHKERRQ(ierr); 1895 for (i = 0; i < numVerts; ++i) { 1896 for (j = 0, k = 0; j < numVerts; ++j) { 1897 if (PetscAbsReal(DiffNormReal(embedDim, &coordsIn[i*embedDim], &coordsIn[j*embedDim]) - edgeLen) < PETSC_SMALL) {graph[i*numVerts+j] = 1; ++k;} 1898 } 1899 if (k != degree) SETERRQ3(comm, PETSC_ERR_PLIB, "Invalid icosahedron, vertex %D degree %D != %D", i, k, degree); 1900 } 1901 /* Build Topology */ 1902 ierr = DMPlexSetChart(*dm, 0, numCells+numVerts);CHKERRQ(ierr); 1903 for (c = 0; c < numCells; c++) { 1904 ierr = DMPlexSetConeSize(*dm, c, embedDim);CHKERRQ(ierr); 1905 } 1906 ierr = DMSetUp(*dm);CHKERRQ(ierr); /* Allocate space for cones */ 1907 /* Cells */ 1908 for (i = 0, c = 0; i < numVerts; ++i) { 1909 for (j = 0; j < i; ++j) { 1910 for (k = 0; k < j; ++k) { 1911 if (graph[i*numVerts+j] && graph[j*numVerts+k] && graph[k*numVerts+i]) { 1912 cone[0] = firstVertex+i; cone[1] = firstVertex+j; cone[2] = firstVertex+k; 1913 /* Check orientation */ 1914 { 1915 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}}}; 1916 PetscReal normal[3]; 1917 PetscInt e, f; 1918 1919 for (d = 0; d < embedDim; ++d) { 1920 normal[d] = 0.0; 1921 for (e = 0; e < embedDim; ++e) { 1922 for (f = 0; f < embedDim; ++f) { 1923 normal[d] += epsilon[d][e][f]*(coordsIn[j*embedDim+e] - coordsIn[i*embedDim+e])*(coordsIn[k*embedDim+f] - coordsIn[i*embedDim+f]); 1924 } 1925 } 1926 } 1927 if (DotReal(embedDim, normal, &coordsIn[i*embedDim]) < 0) {PetscInt tmp = cone[1]; cone[1] = cone[2]; cone[2] = tmp;} 1928 } 1929 ierr = DMPlexSetCone(*dm, c++, cone);CHKERRQ(ierr); 1930 } 1931 } 1932 } 1933 } 1934 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1935 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1936 ierr = PetscFree(graph);CHKERRQ(ierr); 1937 /* Interpolate mesh */ 1938 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1939 ierr = DMDestroy(dm);CHKERRQ(ierr); 1940 *dm = idm; 1941 } else { 1942 /* 1943 12-21--13 1944 | | 1945 25 4 24 1946 | | 1947 12-25--9-16--8-24--13 1948 | | | | 1949 23 5 17 0 15 3 22 1950 | | | | 1951 10-20--6-14--7-19--11 1952 | | 1953 20 1 19 1954 | | 1955 10-18--11 1956 | | 1957 23 2 22 1958 | | 1959 12-21--13 1960 */ 1961 PetscInt cone[4], ornt[4]; 1962 1963 numCells = !rank ? 6 : 0; 1964 numEdges = !rank ? 12 : 0; 1965 numVerts = !rank ? 8 : 0; 1966 firstVertex = numCells; 1967 firstEdge = numCells + numVerts; 1968 /* Build Topology */ 1969 ierr = DMPlexSetChart(*dm, 0, numCells+numEdges+numVerts);CHKERRQ(ierr); 1970 for (c = 0; c < numCells; c++) { 1971 ierr = DMPlexSetConeSize(*dm, c, 4);CHKERRQ(ierr); 1972 } 1973 for (e = firstEdge; e < firstEdge+numEdges; ++e) { 1974 ierr = DMPlexSetConeSize(*dm, e, 2);CHKERRQ(ierr); 1975 } 1976 ierr = DMSetUp(*dm);CHKERRQ(ierr); /* Allocate space for cones */ 1977 if (!rank) { 1978 /* Cell 0 */ 1979 cone[0] = 14; cone[1] = 15; cone[2] = 16; cone[3] = 17; 1980 ierr = DMPlexSetCone(*dm, 0, cone);CHKERRQ(ierr); 1981 ornt[0] = 0; ornt[1] = 0; ornt[2] = 0; ornt[3] = 0; 1982 ierr = DMPlexSetConeOrientation(*dm, 0, ornt);CHKERRQ(ierr); 1983 /* Cell 1 */ 1984 cone[0] = 18; cone[1] = 19; cone[2] = 14; cone[3] = 20; 1985 ierr = DMPlexSetCone(*dm, 1, cone);CHKERRQ(ierr); 1986 ornt[0] = 0; ornt[1] = 0; ornt[2] = -2; ornt[3] = 0; 1987 ierr = DMPlexSetConeOrientation(*dm, 1, ornt);CHKERRQ(ierr); 1988 /* Cell 2 */ 1989 cone[0] = 21; cone[1] = 22; cone[2] = 18; cone[3] = 23; 1990 ierr = DMPlexSetCone(*dm, 2, cone);CHKERRQ(ierr); 1991 ornt[0] = 0; ornt[1] = 0; ornt[2] = -2; ornt[3] = 0; 1992 ierr = DMPlexSetConeOrientation(*dm, 2, ornt);CHKERRQ(ierr); 1993 /* Cell 3 */ 1994 cone[0] = 19; cone[1] = 22; cone[2] = 24; cone[3] = 15; 1995 ierr = DMPlexSetCone(*dm, 3, cone);CHKERRQ(ierr); 1996 ornt[0] = -2; ornt[1] = -2; ornt[2] = 0; ornt[3] = -2; 1997 ierr = DMPlexSetConeOrientation(*dm, 3, ornt);CHKERRQ(ierr); 1998 /* Cell 4 */ 1999 cone[0] = 16; cone[1] = 24; cone[2] = 21; cone[3] = 25; 2000 ierr = DMPlexSetCone(*dm, 4, cone);CHKERRQ(ierr); 2001 ornt[0] = -2; ornt[1] = -2; ornt[2] = -2; ornt[3] = 0; 2002 ierr = DMPlexSetConeOrientation(*dm, 4, ornt);CHKERRQ(ierr); 2003 /* Cell 5 */ 2004 cone[0] = 20; cone[1] = 17; cone[2] = 25; cone[3] = 23; 2005 ierr = DMPlexSetCone(*dm, 5, cone);CHKERRQ(ierr); 2006 ornt[0] = -2; ornt[1] = -2; ornt[2] = -2; ornt[3] = -2; 2007 ierr = DMPlexSetConeOrientation(*dm, 5, ornt);CHKERRQ(ierr); 2008 /* Edges */ 2009 cone[0] = 6; cone[1] = 7; 2010 ierr = DMPlexSetCone(*dm, 14, cone);CHKERRQ(ierr); 2011 cone[0] = 7; cone[1] = 8; 2012 ierr = DMPlexSetCone(*dm, 15, cone);CHKERRQ(ierr); 2013 cone[0] = 8; cone[1] = 9; 2014 ierr = DMPlexSetCone(*dm, 16, cone);CHKERRQ(ierr); 2015 cone[0] = 9; cone[1] = 6; 2016 ierr = DMPlexSetCone(*dm, 17, cone);CHKERRQ(ierr); 2017 cone[0] = 10; cone[1] = 11; 2018 ierr = DMPlexSetCone(*dm, 18, cone);CHKERRQ(ierr); 2019 cone[0] = 11; cone[1] = 7; 2020 ierr = DMPlexSetCone(*dm, 19, cone);CHKERRQ(ierr); 2021 cone[0] = 6; cone[1] = 10; 2022 ierr = DMPlexSetCone(*dm, 20, cone);CHKERRQ(ierr); 2023 cone[0] = 12; cone[1] = 13; 2024 ierr = DMPlexSetCone(*dm, 21, cone);CHKERRQ(ierr); 2025 cone[0] = 13; cone[1] = 11; 2026 ierr = DMPlexSetCone(*dm, 22, cone);CHKERRQ(ierr); 2027 cone[0] = 10; cone[1] = 12; 2028 ierr = DMPlexSetCone(*dm, 23, cone);CHKERRQ(ierr); 2029 cone[0] = 13; cone[1] = 8; 2030 ierr = DMPlexSetCone(*dm, 24, cone);CHKERRQ(ierr); 2031 cone[0] = 12; cone[1] = 9; 2032 ierr = DMPlexSetCone(*dm, 25, cone);CHKERRQ(ierr); 2033 } 2034 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 2035 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 2036 /* Build coordinates */ 2037 ierr = PetscCalloc1(numVerts * embedDim, &coordsIn);CHKERRQ(ierr); 2038 if (!rank) { 2039 coordsIn[0*embedDim+0] = -R; coordsIn[0*embedDim+1] = R; coordsIn[0*embedDim+2] = -R; 2040 coordsIn[1*embedDim+0] = R; coordsIn[1*embedDim+1] = R; coordsIn[1*embedDim+2] = -R; 2041 coordsIn[2*embedDim+0] = R; coordsIn[2*embedDim+1] = -R; coordsIn[2*embedDim+2] = -R; 2042 coordsIn[3*embedDim+0] = -R; coordsIn[3*embedDim+1] = -R; coordsIn[3*embedDim+2] = -R; 2043 coordsIn[4*embedDim+0] = -R; coordsIn[4*embedDim+1] = R; coordsIn[4*embedDim+2] = R; 2044 coordsIn[5*embedDim+0] = R; coordsIn[5*embedDim+1] = R; coordsIn[5*embedDim+2] = R; 2045 coordsIn[6*embedDim+0] = -R; coordsIn[6*embedDim+1] = -R; coordsIn[6*embedDim+2] = R; 2046 coordsIn[7*embedDim+0] = R; coordsIn[7*embedDim+1] = -R; coordsIn[7*embedDim+2] = R; 2047 } 2048 } 2049 break; 2050 case 3: 2051 if (simplex) { 2052 DM idm; 2053 const PetscReal edgeLen = 1.0/PETSC_PHI; 2054 PetscReal vertexA[4] = {0.5, 0.5, 0.5, 0.5}; 2055 PetscReal vertexB[4] = {1.0, 0.0, 0.0, 0.0}; 2056 PetscReal vertexC[4] = {0.5, 0.5*PETSC_PHI, 0.5/PETSC_PHI, 0.0}; 2057 const PetscInt degree = 12; 2058 PetscInt s[4] = {1, 1, 1}; 2059 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}, 2060 {2, 0, 1, 3}, {2, 1, 3, 0}, {2, 3, 0, 1}, {3, 0, 2, 1}, {3, 1, 0, 2}, {3, 2, 1, 0}}; 2061 PetscInt cone[4]; 2062 PetscInt *graph, p, i, j, k, l; 2063 2064 vertexA[0] *= R; vertexA[1] *= R; vertexA[2] *= R; vertexA[3] *= R; 2065 vertexB[0] *= R; vertexB[1] *= R; vertexB[2] *= R; vertexB[3] *= R; 2066 vertexC[0] *= R; vertexC[1] *= R; vertexC[2] *= R; vertexC[3] *= R; 2067 numCells = !rank ? 600 : 0; 2068 numVerts = !rank ? 120 : 0; 2069 firstVertex = numCells; 2070 /* Use the 600-cell, which for a unit sphere has coordinates which are 2071 2072 1/2 (\pm 1, \pm 1, \pm 1, \pm 1) 16 2073 (\pm 1, 0, 0, 0) all cyclic permutations 8 2074 1/2 (\pm 1, \pm phi, \pm 1/phi, 0) all even permutations 96 2075 2076 where \phi^2 - \phi - 1 = 0, meaning \phi is the golden ratio \frac{1 + \sqrt{5}}{2}. The edge 2077 length is then given by 1/\phi = 0.61803. 2078 2079 http://buzzard.pugetsound.edu/sage-practice/ch03s03.html 2080 http://mathworld.wolfram.com/600-Cell.html 2081 */ 2082 /* Construct vertices */ 2083 ierr = PetscCalloc1(numVerts * embedDim, &coordsIn);CHKERRQ(ierr); 2084 i = 0; 2085 if (!rank) { 2086 for (s[0] = -1; s[0] < 2; s[0] += 2) { 2087 for (s[1] = -1; s[1] < 2; s[1] += 2) { 2088 for (s[2] = -1; s[2] < 2; s[2] += 2) { 2089 for (s[3] = -1; s[3] < 2; s[3] += 2) { 2090 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[d]*vertexA[d]; 2091 ++i; 2092 } 2093 } 2094 } 2095 } 2096 for (p = 0; p < embedDim; ++p) { 2097 s[1] = s[2] = s[3] = 1; 2098 for (s[0] = -1; s[0] < 2; s[0] += 2) { 2099 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[(d+p)%embedDim]*vertexB[(d+p)%embedDim]; 2100 ++i; 2101 } 2102 } 2103 for (p = 0; p < 12; ++p) { 2104 s[3] = 1; 2105 for (s[0] = -1; s[0] < 2; s[0] += 2) { 2106 for (s[1] = -1; s[1] < 2; s[1] += 2) { 2107 for (s[2] = -1; s[2] < 2; s[2] += 2) { 2108 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[evenPerm[p][d]]*vertexC[evenPerm[p][d]]; 2109 ++i; 2110 } 2111 } 2112 } 2113 } 2114 } 2115 if (i != numVerts) SETERRQ2(comm, PETSC_ERR_PLIB, "Invalid 600-cell, vertices %D != %D", i, numVerts); 2116 /* Construct graph */ 2117 ierr = PetscCalloc1(numVerts * numVerts, &graph);CHKERRQ(ierr); 2118 for (i = 0; i < numVerts; ++i) { 2119 for (j = 0, k = 0; j < numVerts; ++j) { 2120 if (PetscAbsReal(DiffNormReal(embedDim, &coordsIn[i*embedDim], &coordsIn[j*embedDim]) - edgeLen) < PETSC_SMALL) {graph[i*numVerts+j] = 1; ++k;} 2121 } 2122 if (k != degree) SETERRQ3(comm, PETSC_ERR_PLIB, "Invalid 600-cell, vertex %D degree %D != %D", i, k, degree); 2123 } 2124 /* Build Topology */ 2125 ierr = DMPlexSetChart(*dm, 0, numCells+numVerts);CHKERRQ(ierr); 2126 for (c = 0; c < numCells; c++) { 2127 ierr = DMPlexSetConeSize(*dm, c, embedDim);CHKERRQ(ierr); 2128 } 2129 ierr = DMSetUp(*dm);CHKERRQ(ierr); /* Allocate space for cones */ 2130 /* Cells */ 2131 if (!rank) { 2132 for (i = 0, c = 0; i < numVerts; ++i) { 2133 for (j = 0; j < i; ++j) { 2134 for (k = 0; k < j; ++k) { 2135 for (l = 0; l < k; ++l) { 2136 if (graph[i*numVerts+j] && graph[j*numVerts+k] && graph[k*numVerts+i] && 2137 graph[l*numVerts+i] && graph[l*numVerts+j] && graph[l*numVerts+k]) { 2138 cone[0] = firstVertex+i; cone[1] = firstVertex+j; cone[2] = firstVertex+k; cone[3] = firstVertex+l; 2139 /* Check orientation: https://ef.gy/linear-algebra:normal-vectors-in-higher-dimensional-spaces */ 2140 { 2141 const PetscInt epsilon[4][4][4][4] = {{{{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 2142 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 1}, { 0, 0, -1, 0}}, 2143 {{0, 0, 0, 0}, { 0, 0, 0, -1}, { 0, 0, 0, 0}, { 0, 1, 0, 0}}, 2144 {{0, 0, 0, 0}, { 0, 0, 1, 0}, { 0, -1, 0, 0}, { 0, 0, 0, 0}}}, 2145 2146 {{{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, -1}, { 0, 0, 1, 0}}, 2147 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 2148 {{0, 0, 0, 1}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, {-1, 0, 0, 0}}, 2149 {{0, 0, -1, 0}, { 0, 0, 0, 0}, { 1, 0, 0, 0}, { 0, 0, 0, 0}}}, 2150 2151 {{{0, 0, 0, 0}, { 0, 0, 0, 1}, { 0, 0, 0, 0}, { 0, -1, 0, 0}}, 2152 {{0, 0, 0, -1}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 1, 0, 0, 0}}, 2153 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 2154 {{0, 1, 0, 0}, {-1, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}}, 2155 2156 {{{0, 0, 0, 0}, { 0, 0, -1, 0}, { 0, 1, 0, 0}, { 0, 0, 0, 0}}, 2157 {{0, 0, 1, 0}, { 0, 0, 0, 0}, {-1, 0, 0, 0}, { 0, 0, 0, 0}}, 2158 {{0, -1, 0, 0}, { 1, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 2159 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}}}; 2160 PetscReal normal[4]; 2161 PetscInt e, f, g; 2162 2163 for (d = 0; d < embedDim; ++d) { 2164 normal[d] = 0.0; 2165 for (e = 0; e < embedDim; ++e) { 2166 for (f = 0; f < embedDim; ++f) { 2167 for (g = 0; g < embedDim; ++g) { 2168 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]); 2169 } 2170 } 2171 } 2172 } 2173 if (DotReal(embedDim, normal, &coordsIn[i*embedDim]) < 0) {PetscInt tmp = cone[1]; cone[1] = cone[2]; cone[2] = tmp;} 2174 } 2175 ierr = DMPlexSetCone(*dm, c++, cone);CHKERRQ(ierr); 2176 } 2177 } 2178 } 2179 } 2180 } 2181 } 2182 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 2183 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 2184 ierr = PetscFree(graph);CHKERRQ(ierr); 2185 /* Interpolate mesh */ 2186 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 2187 ierr = DMDestroy(dm);CHKERRQ(ierr); 2188 *dm = idm; 2189 break; 2190 } 2191 default: SETERRQ1(comm, PETSC_ERR_SUP, "Unsupported dimension for sphere: %D", dim); 2192 } 2193 /* Create coordinates */ 2194 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 2195 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 2196 ierr = PetscSectionSetFieldComponents(coordSection, 0, embedDim);CHKERRQ(ierr); 2197 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVerts);CHKERRQ(ierr); 2198 for (v = firstVertex; v < firstVertex+numVerts; ++v) { 2199 ierr = PetscSectionSetDof(coordSection, v, embedDim);CHKERRQ(ierr); 2200 ierr = PetscSectionSetFieldDof(coordSection, v, 0, embedDim);CHKERRQ(ierr); 2201 } 2202 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 2203 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 2204 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 2205 ierr = VecSetBlockSize(coordinates, embedDim);CHKERRQ(ierr); 2206 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 2207 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 2208 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 2209 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2210 for (v = 0; v < numVerts; ++v) for (d = 0; d < embedDim; ++d) {coords[v*embedDim+d] = coordsIn[v*embedDim+d];} 2211 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2212 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 2213 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 2214 ierr = PetscFree(coordsIn);CHKERRQ(ierr); 2215 /* Create coordinate function space */ 2216 { 2217 DM cdm; 2218 PetscDS cds; 2219 PetscFE fe; 2220 PetscScalar radius = R; 2221 PetscInt dT, dE; 2222 2223 ierr = DMGetCoordinateDM(*dm, &cdm);CHKERRQ(ierr); 2224 ierr = DMGetDimension(*dm, &dT);CHKERRQ(ierr); 2225 ierr = DMGetCoordinateDim(*dm, &dE);CHKERRQ(ierr); 2226 ierr = PetscFECreateLagrange(PETSC_COMM_SELF, dT, dE, simplex, 1, -1, &fe);CHKERRQ(ierr); 2227 ierr = DMSetField(cdm, 0, NULL, (PetscObject) fe);CHKERRQ(ierr); 2228 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 2229 ierr = DMCreateDS(cdm);CHKERRQ(ierr); 2230 2231 ierr = DMGetDS(cdm, &cds);CHKERRQ(ierr); 2232 ierr = PetscDSSetConstants(cds, 1, &radius);CHKERRQ(ierr); 2233 } 2234 ((DM_Plex *) (*dm)->data)->coordFunc = snapToSphere; 2235 PetscFunctionReturn(0); 2236 } 2237 2238 /*@ 2239 DMPlexCreateBallMesh - Creates a simplex mesh on the d-dimensional ball, B^d. 2240 2241 Collective 2242 2243 Input Parameters: 2244 + comm - The communicator for the DM object 2245 . dim - The dimension 2246 - R - The radius 2247 2248 Output Parameter: 2249 . dm - The DM object 2250 2251 Options Database Keys: 2252 - bd_dm_refine - This will refine the surface mesh preserving the sphere geometry 2253 2254 Level: beginner 2255 2256 .seealso: DMPlexCreateSphereMesh(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 2257 @*/ 2258 PetscErrorCode DMPlexCreateBallMesh(MPI_Comm comm, PetscInt dim, PetscReal R, DM *dm) 2259 { 2260 DM sdm; 2261 DMLabel bdlabel; 2262 PetscErrorCode ierr; 2263 2264 PetscFunctionBegin; 2265 ierr = DMPlexCreateSphereMesh(comm, dim-1, PETSC_TRUE, R, &sdm);CHKERRQ(ierr); 2266 ierr = PetscObjectSetOptionsPrefix((PetscObject) sdm, "bd_");CHKERRQ(ierr); 2267 ierr = DMSetFromOptions(sdm);CHKERRQ(ierr); 2268 ierr = DMPlexGenerate(sdm, NULL, PETSC_TRUE, dm);CHKERRQ(ierr); 2269 ierr = DMDestroy(&sdm);CHKERRQ(ierr); 2270 ierr = DMCreateLabel(*dm, "marker");CHKERRQ(ierr); 2271 ierr = DMGetLabel(*dm, "marker", &bdlabel);CHKERRQ(ierr); 2272 ierr = DMPlexMarkBoundaryFaces(*dm, PETSC_DETERMINE, bdlabel);CHKERRQ(ierr); 2273 ierr = DMPlexLabelComplete(*dm, bdlabel);CHKERRQ(ierr); 2274 PetscFunctionReturn(0); 2275 } 2276 2277 /* External function declarations here */ 2278 extern PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling); 2279 extern PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, Mat *mat); 2280 extern PetscErrorCode DMCreateMassMatrix_Plex(DM dmCoarse, DM dmFine, Mat *mat); 2281 extern PetscErrorCode DMCreateLocalSection_Plex(DM dm); 2282 extern PetscErrorCode DMCreateDefaultConstraints_Plex(DM dm); 2283 extern PetscErrorCode DMCreateMatrix_Plex(DM dm, Mat *J); 2284 extern PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm); 2285 extern PetscErrorCode DMCreateCoordinateField_Plex(DM dm, DMField *field); 2286 PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm); 2287 extern PetscErrorCode DMSetUp_Plex(DM dm); 2288 extern PetscErrorCode DMDestroy_Plex(DM dm); 2289 extern PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer); 2290 extern PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer); 2291 extern PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm); 2292 extern PetscErrorCode DMCreateSuperDM_Plex(DM dms[], PetscInt len, IS **is, DM *superdm); 2293 static PetscErrorCode DMInitialize_Plex(DM dm); 2294 2295 /* Replace dm with the contents of dmNew 2296 - Share the DM_Plex structure 2297 - Share the coordinates 2298 - Share the SF 2299 */ 2300 static PetscErrorCode DMPlexReplace_Static(DM dm, DM dmNew) 2301 { 2302 PetscSF sf; 2303 DM coordDM, coarseDM; 2304 DMField coordField; 2305 Vec coords; 2306 PetscBool isper; 2307 const PetscReal *maxCell, *L; 2308 const DMBoundaryType *bd; 2309 PetscErrorCode ierr; 2310 2311 PetscFunctionBegin; 2312 ierr = DMGetPointSF(dmNew, &sf);CHKERRQ(ierr); 2313 ierr = DMSetPointSF(dm, sf);CHKERRQ(ierr); 2314 ierr = DMGetCoordinateDM(dmNew, &coordDM);CHKERRQ(ierr); 2315 ierr = DMGetCoordinatesLocal(dmNew, &coords);CHKERRQ(ierr); 2316 ierr = DMGetCoordinateField(dmNew, &coordField);CHKERRQ(ierr); 2317 ierr = DMSetCoordinateDM(dm, coordDM);CHKERRQ(ierr); 2318 ierr = DMSetCoordinatesLocal(dm, coords);CHKERRQ(ierr); 2319 ierr = DMSetCoordinateField(dm, coordField);CHKERRQ(ierr); 2320 ierr = DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);CHKERRQ(ierr); 2321 ierr = DMSetPeriodicity(dmNew, isper, maxCell, L, bd);CHKERRQ(ierr); 2322 ierr = DMDestroy_Plex(dm);CHKERRQ(ierr); 2323 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 2324 dm->data = dmNew->data; 2325 ((DM_Plex *) dmNew->data)->refct++; 2326 ierr = DMDestroyLabelLinkList_Internal(dm);CHKERRQ(ierr); 2327 ierr = DMCopyLabels(dmNew, dm, PETSC_OWN_POINTER, PETSC_TRUE);CHKERRQ(ierr); 2328 ierr = DMGetCoarseDM(dmNew,&coarseDM);CHKERRQ(ierr); 2329 ierr = DMSetCoarseDM(dm,coarseDM);CHKERRQ(ierr); 2330 PetscFunctionReturn(0); 2331 } 2332 2333 /* Swap dm with the contents of dmNew 2334 - Swap the DM_Plex structure 2335 - Swap the coordinates 2336 - Swap the point PetscSF 2337 */ 2338 static PetscErrorCode DMPlexSwap_Static(DM dmA, DM dmB) 2339 { 2340 DM coordDMA, coordDMB; 2341 Vec coordsA, coordsB; 2342 PetscSF sfA, sfB; 2343 void *tmp; 2344 DMLabelLink listTmp; 2345 DMLabel depthTmp; 2346 PetscInt tmpI; 2347 PetscErrorCode ierr; 2348 2349 PetscFunctionBegin; 2350 ierr = DMGetPointSF(dmA, &sfA);CHKERRQ(ierr); 2351 ierr = DMGetPointSF(dmB, &sfB);CHKERRQ(ierr); 2352 ierr = PetscObjectReference((PetscObject) sfA);CHKERRQ(ierr); 2353 ierr = DMSetPointSF(dmA, sfB);CHKERRQ(ierr); 2354 ierr = DMSetPointSF(dmB, sfA);CHKERRQ(ierr); 2355 ierr = PetscObjectDereference((PetscObject) sfA);CHKERRQ(ierr); 2356 2357 ierr = DMGetCoordinateDM(dmA, &coordDMA);CHKERRQ(ierr); 2358 ierr = DMGetCoordinateDM(dmB, &coordDMB);CHKERRQ(ierr); 2359 ierr = PetscObjectReference((PetscObject) coordDMA);CHKERRQ(ierr); 2360 ierr = DMSetCoordinateDM(dmA, coordDMB);CHKERRQ(ierr); 2361 ierr = DMSetCoordinateDM(dmB, coordDMA);CHKERRQ(ierr); 2362 ierr = PetscObjectDereference((PetscObject) coordDMA);CHKERRQ(ierr); 2363 2364 ierr = DMGetCoordinatesLocal(dmA, &coordsA);CHKERRQ(ierr); 2365 ierr = DMGetCoordinatesLocal(dmB, &coordsB);CHKERRQ(ierr); 2366 ierr = PetscObjectReference((PetscObject) coordsA);CHKERRQ(ierr); 2367 ierr = DMSetCoordinatesLocal(dmA, coordsB);CHKERRQ(ierr); 2368 ierr = DMSetCoordinatesLocal(dmB, coordsA);CHKERRQ(ierr); 2369 ierr = PetscObjectDereference((PetscObject) coordsA);CHKERRQ(ierr); 2370 2371 tmp = dmA->data; 2372 dmA->data = dmB->data; 2373 dmB->data = tmp; 2374 listTmp = dmA->labels; 2375 dmA->labels = dmB->labels; 2376 dmB->labels = listTmp; 2377 depthTmp = dmA->depthLabel; 2378 dmA->depthLabel = dmB->depthLabel; 2379 dmB->depthLabel = depthTmp; 2380 depthTmp = dmA->celltypeLabel; 2381 dmA->celltypeLabel = dmB->celltypeLabel; 2382 dmB->celltypeLabel = depthTmp; 2383 tmpI = dmA->levelup; 2384 dmA->levelup = dmB->levelup; 2385 dmB->levelup = tmpI; 2386 PetscFunctionReturn(0); 2387 } 2388 2389 PetscErrorCode DMSetFromOptions_NonRefinement_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 2390 { 2391 DM_Plex *mesh = (DM_Plex*) dm->data; 2392 PetscBool flg; 2393 PetscErrorCode ierr; 2394 2395 PetscFunctionBegin; 2396 /* Handle viewing */ 2397 ierr = PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMPlexMatSetClosure", PETSC_FALSE, &mesh->printSetValues, NULL);CHKERRQ(ierr); 2398 ierr = PetscOptionsBoundedInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMPlexSNESComputeResidualFEM", 0, &mesh->printFEM, NULL,0);CHKERRQ(ierr); 2399 ierr = PetscOptionsReal("-dm_plex_print_tol", "Tolerance for FEM output", "DMPlexSNESComputeResidualFEM", mesh->printTol, &mesh->printTol, NULL);CHKERRQ(ierr); 2400 ierr = PetscOptionsBoundedInt("-dm_plex_print_l2", "Debug output level all L2 diff computations", "DMComputeL2Diff", 0, &mesh->printL2, NULL,0);CHKERRQ(ierr); 2401 ierr = DMMonitorSetFromOptions(dm, "-dm_plex_monitor_throughput", "Monitor the simulation throughput", "DMPlexMonitorThroughput", DMPlexMonitorThroughput, NULL, &flg);CHKERRQ(ierr); 2402 if (flg) {ierr = PetscLogDefaultBegin();CHKERRQ(ierr);} 2403 /* Point Location */ 2404 ierr = PetscOptionsBool("-dm_plex_hash_location", "Use grid hashing for point location", "DMInterpolate", PETSC_FALSE, &mesh->useHashLocation, NULL);CHKERRQ(ierr); 2405 /* Partitioning and distribution */ 2406 ierr = PetscOptionsBool("-dm_plex_partition_balance", "Attempt to evenly divide points on partition boundary between processes", "DMPlexSetPartitionBalance", PETSC_FALSE, &mesh->partitionBalance, NULL);CHKERRQ(ierr); 2407 /* Generation and remeshing */ 2408 ierr = PetscOptionsBool("-dm_plex_remesh_bd", "Allow changes to the boundary on remeshing", "DMAdapt", PETSC_FALSE, &mesh->remeshBd, NULL);CHKERRQ(ierr); 2409 /* Projection behavior */ 2410 ierr = PetscOptionsBoundedInt("-dm_plex_max_projection_height", "Maxmimum mesh point height used to project locally", "DMPlexSetMaxProjectionHeight", 0, &mesh->maxProjectionHeight, NULL,0);CHKERRQ(ierr); 2411 ierr = PetscOptionsBool("-dm_plex_regular_refinement", "Use special nested projection algorithm for regular refinement", "DMPlexSetRegularRefinement", mesh->regularRefinement, &mesh->regularRefinement, NULL);CHKERRQ(ierr); 2412 ierr = PetscOptionsEnum("-dm_plex_cell_refiner", "Strategy for cell refinment", "ex40.c", DMPlexCellRefinerTypes, (PetscEnum) mesh->cellRefiner, (PetscEnum *) &mesh->cellRefiner, NULL);CHKERRQ(ierr); 2413 /* Checking structure */ 2414 { 2415 PetscBool flg = PETSC_FALSE, flg2 = PETSC_FALSE, all = PETSC_FALSE; 2416 2417 ierr = PetscOptionsBool("-dm_plex_check_all", "Perform all checks", NULL, PETSC_FALSE, &all, &flg2);CHKERRQ(ierr); 2418 ierr = PetscOptionsBool("-dm_plex_check_symmetry", "Check that the adjacency information in the mesh is symmetric", "DMPlexCheckSymmetry", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr); 2419 if (all || (flg && flg2)) {ierr = DMPlexCheckSymmetry(dm);CHKERRQ(ierr);} 2420 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); 2421 if (all || (flg && flg2)) {ierr = DMPlexCheckSkeleton(dm, 0);CHKERRQ(ierr);} 2422 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); 2423 if (all || (flg && flg2)) {ierr = DMPlexCheckFaces(dm, 0);CHKERRQ(ierr);} 2424 ierr = PetscOptionsBool("-dm_plex_check_geometry", "Check that cells have positive volume", "DMPlexCheckGeometry", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr); 2425 if (all || (flg && flg2)) {ierr = DMPlexCheckGeometry(dm);CHKERRQ(ierr);} 2426 ierr = PetscOptionsBool("-dm_plex_check_pointsf", "Check some necessary conditions for PointSF", "DMPlexCheckPointSF", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr); 2427 if (all || (flg && flg2)) {ierr = DMPlexCheckPointSF(dm);CHKERRQ(ierr);} 2428 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); 2429 if (all || (flg && flg2)) {ierr = DMPlexCheckInterfaceCones(dm);CHKERRQ(ierr);} 2430 ierr = PetscOptionsBool("-dm_plex_check_cell_shape", "Check cell shape", "DMPlexCheckCellShape", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr); 2431 if (flg && flg2) {ierr = DMPlexCheckCellShape(dm, PETSC_TRUE, PETSC_DETERMINE);CHKERRQ(ierr);} 2432 } 2433 2434 ierr = PetscPartitionerSetFromOptions(mesh->partitioner);CHKERRQ(ierr); 2435 PetscFunctionReturn(0); 2436 } 2437 2438 static PetscErrorCode DMSetFromOptions_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 2439 { 2440 PetscReal volume = -1.0; 2441 PetscInt prerefine = 0, refine = 0, r, coarsen = 0, overlap = 0; 2442 PetscBool uniformOrig, uniform = PETSC_TRUE, distribute = PETSC_FALSE, isHierarchy, flg; 2443 PetscErrorCode ierr; 2444 2445 PetscFunctionBegin; 2446 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2447 ierr = PetscOptionsHead(PetscOptionsObject,"DMPlex Options");CHKERRQ(ierr); 2448 /* Handle DMPlex refinement before distribution */ 2449 ierr = DMPlexGetRefinementUniform(dm, &uniformOrig);CHKERRQ(ierr); 2450 ierr = PetscOptionsBool("-dm_refine_uniform_pre", "Flag for uniform refinement before distribution", "DMCreate", uniform, &uniform, &flg);CHKERRQ(ierr); 2451 if (flg) {ierr = DMPlexSetRefinementUniform(dm, uniform);CHKERRQ(ierr);} 2452 ierr = PetscOptionsReal("-dm_refine_volume_limit_pre", "The maximum cell volume after refinement before distribution", "DMCreate", volume, &volume, &flg);CHKERRQ(ierr); 2453 if (flg) {ierr = DMPlexSetRefinementLimit(dm, volume);CHKERRQ(ierr);} 2454 ierr = PetscOptionsBoundedInt("-dm_refine_pre", "The number of refinements before distribution", "DMCreate", prerefine, &prerefine, NULL,0);CHKERRQ(ierr); 2455 for (r = 0; r < prerefine; ++r) { 2456 DM rdm; 2457 PetscPointFunc coordFunc = ((DM_Plex*) dm->data)->coordFunc; 2458 2459 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2460 ierr = DMRefine(dm, PetscObjectComm((PetscObject) dm), &rdm);CHKERRQ(ierr); 2461 /* Total hack since we do not pass in a pointer */ 2462 ierr = DMPlexReplace_Static(dm, rdm);CHKERRQ(ierr); 2463 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2464 if (coordFunc) { 2465 ierr = DMPlexRemapGeometry(dm, 0.0, coordFunc);CHKERRQ(ierr); 2466 ((DM_Plex*) dm->data)->coordFunc = coordFunc; 2467 } 2468 ierr = DMDestroy(&rdm);CHKERRQ(ierr); 2469 } 2470 ierr = DMPlexSetRefinementUniform(dm, uniformOrig);CHKERRQ(ierr); 2471 /* Handle DMPlex distribution */ 2472 ierr = PetscOptionsBool("-dm_distribute", "Flag to redistribute a mesh among processes", "DMCreate", distribute, &distribute, NULL);CHKERRQ(ierr); 2473 ierr = PetscOptionsBoundedInt("-dm_distribute_overlap", "The size of the overlap halo", "DMCreate", overlap, &overlap, NULL, 0);CHKERRQ(ierr); 2474 if (distribute) { 2475 DM pdm = NULL; 2476 PetscPartitioner part; 2477 2478 ierr = DMPlexGetPartitioner(dm, &part);CHKERRQ(ierr); 2479 ierr = PetscPartitionerSetFromOptions(part);CHKERRQ(ierr); 2480 ierr = DMPlexDistribute(dm, overlap, NULL, &pdm);CHKERRQ(ierr); 2481 if (pdm) { 2482 ierr = DMPlexReplace_Static(dm, pdm);CHKERRQ(ierr); 2483 ierr = DMDestroy(&pdm);CHKERRQ(ierr); 2484 } 2485 } 2486 /* Handle DMPlex refinement */ 2487 ierr = PetscOptionsBoundedInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL,0);CHKERRQ(ierr); 2488 ierr = PetscOptionsBoundedInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy,0);CHKERRQ(ierr); 2489 if (refine) {ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE);CHKERRQ(ierr);} 2490 if (refine && isHierarchy) { 2491 DM *dms, coarseDM; 2492 2493 ierr = DMGetCoarseDM(dm, &coarseDM);CHKERRQ(ierr); 2494 ierr = PetscObjectReference((PetscObject)coarseDM);CHKERRQ(ierr); 2495 ierr = PetscMalloc1(refine,&dms);CHKERRQ(ierr); 2496 ierr = DMRefineHierarchy(dm, refine, dms);CHKERRQ(ierr); 2497 /* Total hack since we do not pass in a pointer */ 2498 ierr = DMPlexSwap_Static(dm, dms[refine-1]);CHKERRQ(ierr); 2499 if (refine == 1) { 2500 ierr = DMSetCoarseDM(dm, dms[0]);CHKERRQ(ierr); 2501 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr); 2502 } else { 2503 ierr = DMSetCoarseDM(dm, dms[refine-2]);CHKERRQ(ierr); 2504 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr); 2505 ierr = DMSetCoarseDM(dms[0], dms[refine-1]);CHKERRQ(ierr); 2506 ierr = DMPlexSetRegularRefinement(dms[0], PETSC_TRUE);CHKERRQ(ierr); 2507 } 2508 ierr = DMSetCoarseDM(dms[refine-1], coarseDM);CHKERRQ(ierr); 2509 ierr = PetscObjectDereference((PetscObject)coarseDM);CHKERRQ(ierr); 2510 /* Free DMs */ 2511 for (r = 0; r < refine; ++r) { 2512 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr); 2513 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 2514 } 2515 ierr = PetscFree(dms);CHKERRQ(ierr); 2516 } else { 2517 for (r = 0; r < refine; ++r) { 2518 DM refinedMesh; 2519 PetscPointFunc coordFunc = ((DM_Plex*) dm->data)->coordFunc; 2520 2521 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2522 ierr = DMRefine(dm, PetscObjectComm((PetscObject) dm), &refinedMesh);CHKERRQ(ierr); 2523 /* Total hack since we do not pass in a pointer */ 2524 ierr = DMPlexReplace_Static(dm, refinedMesh);CHKERRQ(ierr); 2525 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2526 if (coordFunc) { 2527 ierr = DMPlexRemapGeometry(dm, 0.0, coordFunc);CHKERRQ(ierr); 2528 ((DM_Plex*) dm->data)->coordFunc = coordFunc; 2529 } 2530 ierr = DMDestroy(&refinedMesh);CHKERRQ(ierr); 2531 } 2532 } 2533 /* Handle DMPlex coarsening */ 2534 ierr = PetscOptionsBoundedInt("-dm_coarsen", "Coarsen the mesh", "DMCreate", coarsen, &coarsen, NULL,0);CHKERRQ(ierr); 2535 ierr = PetscOptionsBoundedInt("-dm_coarsen_hierarchy", "The number of coarsenings", "DMCreate", coarsen, &coarsen, &isHierarchy,0);CHKERRQ(ierr); 2536 if (coarsen && isHierarchy) { 2537 DM *dms; 2538 2539 ierr = PetscMalloc1(coarsen, &dms);CHKERRQ(ierr); 2540 ierr = DMCoarsenHierarchy(dm, coarsen, dms);CHKERRQ(ierr); 2541 /* Free DMs */ 2542 for (r = 0; r < coarsen; ++r) { 2543 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr); 2544 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 2545 } 2546 ierr = PetscFree(dms);CHKERRQ(ierr); 2547 } else { 2548 for (r = 0; r < coarsen; ++r) { 2549 DM coarseMesh; 2550 2551 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2552 ierr = DMCoarsen(dm, PetscObjectComm((PetscObject) dm), &coarseMesh);CHKERRQ(ierr); 2553 /* Total hack since we do not pass in a pointer */ 2554 ierr = DMPlexReplace_Static(dm, coarseMesh);CHKERRQ(ierr); 2555 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2556 ierr = DMDestroy(&coarseMesh);CHKERRQ(ierr); 2557 } 2558 } 2559 /* Handle */ 2560 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2561 ierr = PetscOptionsTail();CHKERRQ(ierr); 2562 PetscFunctionReturn(0); 2563 } 2564 2565 static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec) 2566 { 2567 PetscErrorCode ierr; 2568 2569 PetscFunctionBegin; 2570 ierr = DMCreateGlobalVector_Section_Private(dm,vec);CHKERRQ(ierr); 2571 /* ierr = VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM);CHKERRQ(ierr); */ 2572 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex);CHKERRQ(ierr); 2573 ierr = VecSetOperation(*vec, VECOP_VIEWNATIVE, (void (*)(void)) VecView_Plex_Native);CHKERRQ(ierr); 2574 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex);CHKERRQ(ierr); 2575 ierr = VecSetOperation(*vec, VECOP_LOADNATIVE, (void (*)(void)) VecLoad_Plex_Native);CHKERRQ(ierr); 2576 PetscFunctionReturn(0); 2577 } 2578 2579 static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec) 2580 { 2581 PetscErrorCode ierr; 2582 2583 PetscFunctionBegin; 2584 ierr = DMCreateLocalVector_Section_Private(dm,vec);CHKERRQ(ierr); 2585 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex_Local);CHKERRQ(ierr); 2586 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex_Local);CHKERRQ(ierr); 2587 PetscFunctionReturn(0); 2588 } 2589 2590 static PetscErrorCode DMGetDimPoints_Plex(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd) 2591 { 2592 PetscInt depth, d; 2593 PetscErrorCode ierr; 2594 2595 PetscFunctionBegin; 2596 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 2597 if (depth == 1) { 2598 ierr = DMGetDimension(dm, &d);CHKERRQ(ierr); 2599 if (dim == 0) {ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr);} 2600 else if (dim == d) {ierr = DMPlexGetDepthStratum(dm, 1, pStart, pEnd);CHKERRQ(ierr);} 2601 else {*pStart = 0; *pEnd = 0;} 2602 } else { 2603 ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr); 2604 } 2605 PetscFunctionReturn(0); 2606 } 2607 2608 static PetscErrorCode DMGetNeighbors_Plex(DM dm, PetscInt *nranks, const PetscMPIInt *ranks[]) 2609 { 2610 PetscSF sf; 2611 PetscInt niranks, njranks, n; 2612 const PetscMPIInt *iranks, *jranks; 2613 DM_Plex *data = (DM_Plex*) dm->data; 2614 PetscErrorCode ierr; 2615 2616 PetscFunctionBegin; 2617 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 2618 if (!data->neighbors) { 2619 ierr = PetscSFGetRootRanks(sf, &njranks, &jranks, NULL, NULL, NULL);CHKERRQ(ierr); 2620 ierr = PetscSFGetLeafRanks(sf, &niranks, &iranks, NULL, NULL);CHKERRQ(ierr); 2621 ierr = PetscMalloc1(njranks + niranks + 1, &data->neighbors);CHKERRQ(ierr); 2622 ierr = PetscArraycpy(data->neighbors + 1, jranks, njranks);CHKERRQ(ierr); 2623 ierr = PetscArraycpy(data->neighbors + njranks + 1, iranks, niranks);CHKERRQ(ierr); 2624 n = njranks + niranks; 2625 ierr = PetscSortRemoveDupsMPIInt(&n, data->neighbors + 1);CHKERRQ(ierr); 2626 /* The following cast should never fail: can't have more neighbors than PETSC_MPI_INT_MAX */ 2627 ierr = PetscMPIIntCast(n, data->neighbors);CHKERRQ(ierr); 2628 } 2629 if (nranks) *nranks = data->neighbors[0]; 2630 if (ranks) { 2631 if (data->neighbors[0]) *ranks = data->neighbors + 1; 2632 else *ranks = NULL; 2633 } 2634 PetscFunctionReturn(0); 2635 } 2636 2637 static PetscErrorCode DMInitialize_Plex(DM dm) 2638 { 2639 PetscErrorCode ierr; 2640 2641 PetscFunctionBegin; 2642 dm->ops->view = DMView_Plex; 2643 dm->ops->load = DMLoad_Plex; 2644 dm->ops->setfromoptions = DMSetFromOptions_Plex; 2645 dm->ops->clone = DMClone_Plex; 2646 dm->ops->setup = DMSetUp_Plex; 2647 dm->ops->createlocalsection = DMCreateLocalSection_Plex; 2648 dm->ops->createdefaultconstraints = DMCreateDefaultConstraints_Plex; 2649 dm->ops->createglobalvector = DMCreateGlobalVector_Plex; 2650 dm->ops->createlocalvector = DMCreateLocalVector_Plex; 2651 dm->ops->getlocaltoglobalmapping = NULL; 2652 dm->ops->createfieldis = NULL; 2653 dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex; 2654 dm->ops->createcoordinatefield = DMCreateCoordinateField_Plex; 2655 dm->ops->getcoloring = NULL; 2656 dm->ops->creatematrix = DMCreateMatrix_Plex; 2657 dm->ops->createinterpolation = DMCreateInterpolation_Plex; 2658 dm->ops->createmassmatrix = DMCreateMassMatrix_Plex; 2659 dm->ops->createinjection = DMCreateInjection_Plex; 2660 dm->ops->refine = DMRefine_Plex; 2661 dm->ops->coarsen = DMCoarsen_Plex; 2662 dm->ops->refinehierarchy = DMRefineHierarchy_Plex; 2663 dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Plex; 2664 dm->ops->adaptlabel = DMAdaptLabel_Plex; 2665 dm->ops->adaptmetric = DMAdaptMetric_Plex; 2666 dm->ops->globaltolocalbegin = NULL; 2667 dm->ops->globaltolocalend = NULL; 2668 dm->ops->localtoglobalbegin = NULL; 2669 dm->ops->localtoglobalend = NULL; 2670 dm->ops->destroy = DMDestroy_Plex; 2671 dm->ops->createsubdm = DMCreateSubDM_Plex; 2672 dm->ops->createsuperdm = DMCreateSuperDM_Plex; 2673 dm->ops->getdimpoints = DMGetDimPoints_Plex; 2674 dm->ops->locatepoints = DMLocatePoints_Plex; 2675 dm->ops->projectfunctionlocal = DMProjectFunctionLocal_Plex; 2676 dm->ops->projectfunctionlabellocal = DMProjectFunctionLabelLocal_Plex; 2677 dm->ops->projectfieldlocal = DMProjectFieldLocal_Plex; 2678 dm->ops->projectfieldlabellocal = DMProjectFieldLabelLocal_Plex; 2679 dm->ops->projectbdfieldlabellocal = DMProjectBdFieldLabelLocal_Plex; 2680 dm->ops->computel2diff = DMComputeL2Diff_Plex; 2681 dm->ops->computel2gradientdiff = DMComputeL2GradientDiff_Plex; 2682 dm->ops->computel2fielddiff = DMComputeL2FieldDiff_Plex; 2683 dm->ops->getneighbors = DMGetNeighbors_Plex; 2684 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMPlexInsertBoundaryValues_C",DMPlexInsertBoundaryValues_Plex);CHKERRQ(ierr); 2685 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMPlexInsertTimeDerviativeBoundaryValues_C",DMPlexInsertTimeDerivativeBoundaryValues_Plex);CHKERRQ(ierr); 2686 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMSetUpGLVisViewer_C",DMSetUpGLVisViewer_Plex);CHKERRQ(ierr); 2687 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMCreateNeumannOverlap_C",DMCreateNeumannOverlap_Plex);CHKERRQ(ierr); 2688 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMPlexGetOverlap_C",DMPlexGetOverlap_Plex);CHKERRQ(ierr); 2689 PetscFunctionReturn(0); 2690 } 2691 2692 PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm) 2693 { 2694 DM_Plex *mesh = (DM_Plex *) dm->data; 2695 PetscErrorCode ierr; 2696 2697 PetscFunctionBegin; 2698 mesh->refct++; 2699 (*newdm)->data = mesh; 2700 ierr = PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);CHKERRQ(ierr); 2701 ierr = DMInitialize_Plex(*newdm);CHKERRQ(ierr); 2702 PetscFunctionReturn(0); 2703 } 2704 2705 /*MC 2706 DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram. 2707 In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is 2708 specified by a PetscSection object. Ownership in the global representation is determined by 2709 ownership of the underlying DMPlex points. This is specified by another PetscSection object. 2710 2711 Options Database Keys: 2712 + -dm_refine_pre - Refine mesh before distribution 2713 + -dm_refine_uniform_pre - Choose uniform or generator-based refinement 2714 + -dm_refine_volume_limit_pre - Cell volume limit after pre-refinement using generator 2715 . -dm_distribute - Distribute mesh across processes 2716 . -dm_distribute_overlap - Number of cells to overlap for distribution 2717 . -dm_refine - Refine mesh after distribution 2718 . -dm_plex_hash_location - Use grid hashing for point location 2719 . -dm_plex_partition_balance - Attempt to evenly divide points on partition boundary between processes 2720 . -dm_plex_remesh_bd - Allow changes to the boundary on remeshing 2721 . -dm_plex_max_projection_height - Maxmimum mesh point height used to project locally 2722 . -dm_plex_regular_refinement - Use special nested projection algorithm for regular refinement 2723 . -dm_plex_check_all - Perform all shecks below 2724 . -dm_plex_check_symmetry - Check that the adjacency information in the mesh is symmetric 2725 . -dm_plex_check_skeleton <celltype> - Check that each cell has the correct number of vertices 2726 . -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 2727 . -dm_plex_check_geometry - Check that cells have positive volume 2728 . -dm_view :mesh.tex:ascii_latex - View the mesh in LaTeX/TikZ 2729 . -dm_plex_view_scale <num> - Scale the TikZ 2730 - -dm_plex_print_fem <num> - View FEM assembly information, such as element vectors and matrices 2731 2732 2733 Level: intermediate 2734 2735 .seealso: DMType, DMPlexCreate(), DMCreate(), DMSetType() 2736 M*/ 2737 2738 PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm) 2739 { 2740 DM_Plex *mesh; 2741 PetscInt unit; 2742 PetscErrorCode ierr; 2743 2744 PetscFunctionBegin; 2745 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2746 ierr = PetscNewLog(dm,&mesh);CHKERRQ(ierr); 2747 dm->dim = 0; 2748 dm->data = mesh; 2749 2750 mesh->refct = 1; 2751 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);CHKERRQ(ierr); 2752 mesh->maxConeSize = 0; 2753 mesh->cones = NULL; 2754 mesh->coneOrientations = NULL; 2755 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);CHKERRQ(ierr); 2756 mesh->maxSupportSize = 0; 2757 mesh->supports = NULL; 2758 mesh->refinementUniform = PETSC_TRUE; 2759 mesh->refinementLimit = -1.0; 2760 mesh->interpolated = DMPLEX_INTERPOLATED_INVALID; 2761 mesh->interpolatedCollective = DMPLEX_INTERPOLATED_INVALID; 2762 2763 mesh->facesTmp = NULL; 2764 2765 mesh->tetgenOpts = NULL; 2766 mesh->triangleOpts = NULL; 2767 ierr = PetscPartitionerCreate(PetscObjectComm((PetscObject)dm), &mesh->partitioner);CHKERRQ(ierr); 2768 mesh->remeshBd = PETSC_FALSE; 2769 2770 mesh->subpointMap = NULL; 2771 2772 for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0; 2773 2774 mesh->regularRefinement = PETSC_FALSE; 2775 mesh->depthState = -1; 2776 mesh->celltypeState = -1; 2777 mesh->globalVertexNumbers = NULL; 2778 mesh->globalCellNumbers = NULL; 2779 mesh->anchorSection = NULL; 2780 mesh->anchorIS = NULL; 2781 mesh->createanchors = NULL; 2782 mesh->computeanchormatrix = NULL; 2783 mesh->parentSection = NULL; 2784 mesh->parents = NULL; 2785 mesh->childIDs = NULL; 2786 mesh->childSection = NULL; 2787 mesh->children = NULL; 2788 mesh->referenceTree = NULL; 2789 mesh->getchildsymmetry = NULL; 2790 mesh->vtkCellHeight = 0; 2791 mesh->useAnchors = PETSC_FALSE; 2792 2793 mesh->maxProjectionHeight = 0; 2794 2795 mesh->neighbors = NULL; 2796 2797 mesh->printSetValues = PETSC_FALSE; 2798 mesh->printFEM = 0; 2799 mesh->printTol = 1.0e-10; 2800 2801 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 2802 PetscFunctionReturn(0); 2803 } 2804 2805 /*@ 2806 DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram. 2807 2808 Collective 2809 2810 Input Parameter: 2811 . comm - The communicator for the DMPlex object 2812 2813 Output Parameter: 2814 . mesh - The DMPlex object 2815 2816 Level: beginner 2817 2818 @*/ 2819 PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh) 2820 { 2821 PetscErrorCode ierr; 2822 2823 PetscFunctionBegin; 2824 PetscValidPointer(mesh,2); 2825 ierr = DMCreate(comm, mesh);CHKERRQ(ierr); 2826 ierr = DMSetType(*mesh, DMPLEX);CHKERRQ(ierr); 2827 PetscFunctionReturn(0); 2828 } 2829 2830 /*@C 2831 DMPlexBuildFromCellListParallel - Build distributed DMPLEX topology from a list of vertices for each cell (common mesh generator output) 2832 2833 Input Parameters: 2834 + dm - The DM 2835 . numCells - The number of cells owned by this process 2836 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE 2837 . NVertices - The global number of vertices, or PETSC_DECIDE 2838 . numCorners - The number of vertices for each cell 2839 - cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 2840 2841 Output Parameter: 2842 . vertexSF - (Optional) SF describing complete vertex ownership 2843 2844 Notes: 2845 Two triangles sharing a face 2846 $ 2847 $ 2 2848 $ / | \ 2849 $ / | \ 2850 $ / | \ 2851 $ 0 0 | 1 3 2852 $ \ | / 2853 $ \ | / 2854 $ \ | / 2855 $ 1 2856 would have input 2857 $ numCells = 2, numVertices = 4 2858 $ cells = [0 1 2 1 3 2] 2859 $ 2860 which would result in the DMPlex 2861 $ 2862 $ 4 2863 $ / | \ 2864 $ / | \ 2865 $ / | \ 2866 $ 2 0 | 1 5 2867 $ \ | / 2868 $ \ | / 2869 $ \ | / 2870 $ 3 2871 2872 Vertices are implicitly numbered consecutively 0,...,NVertices. 2873 Each rank owns a chunk of numVertices consecutive vertices. 2874 If numVertices is PETSC_DECIDE, PETSc will distribute them as evenly as possible using PetscLayout. 2875 If both NVertices and numVertices are PETSC_DECIDE, NVertices is computed by PETSc as the maximum vertex index in cells + 1. 2876 If only NVertices is PETSC_DECIDE, it is computed as the sum of numVertices over all ranks. 2877 2878 The cell distribution is arbitrary non-overlapping, independent of the vertex distribution. 2879 2880 Not currently supported in Fortran. 2881 2882 Level: advanced 2883 2884 .seealso: DMPlexBuildFromCellList(), DMPlexCreateFromCellListParallelPetsc(), DMPlexBuildCoordinatesFromCellListParallel() 2885 @*/ 2886 PetscErrorCode DMPlexBuildFromCellListParallel(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt NVertices, PetscInt numCorners, const PetscInt cells[], PetscSF *vertexSF) 2887 { 2888 PetscSF sfPoint, sfVert; 2889 PetscLayout vLayout; 2890 PetscSFNode *vertexLocal, *vertexOwner, *remoteVertex; 2891 PetscInt numVerticesAdj, *verticesAdj, numVerticesGhost = 0, *localVertex, *cones, c, p, v, g, dim; 2892 PetscMPIInt rank, size; 2893 PetscErrorCode ierr; 2894 2895 PetscFunctionBegin; 2896 PetscValidLogicalCollectiveInt(dm,NVertices,4); 2897 ierr = PetscLogEventBegin(DMPLEX_BuildFromCellList,dm,0,0,0);CHKERRQ(ierr); 2898 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr); 2899 ierr = MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);CHKERRQ(ierr); 2900 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2901 /* Get/check global number of vertices */ 2902 { 2903 PetscInt NVerticesInCells, i; 2904 const PetscInt len = numCells * numCorners; 2905 2906 /* NVerticesInCells = max(cells) + 1 */ 2907 NVerticesInCells = PETSC_MIN_INT; 2908 for (i=0; i<len; i++) if (cells[i] > NVerticesInCells) NVerticesInCells = cells[i]; 2909 ++NVerticesInCells; 2910 ierr = MPI_Allreduce(MPI_IN_PLACE, &NVerticesInCells, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject) dm));CHKERRQ(ierr); 2911 2912 if (numVertices == PETSC_DECIDE && NVertices == PETSC_DECIDE) NVertices = NVerticesInCells; 2913 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); 2914 } 2915 /* Partition vertices */ 2916 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) dm), &vLayout);CHKERRQ(ierr); 2917 ierr = PetscLayoutSetLocalSize(vLayout, numVertices);CHKERRQ(ierr); 2918 ierr = PetscLayoutSetSize(vLayout, NVertices);CHKERRQ(ierr); 2919 ierr = PetscLayoutSetBlockSize(vLayout, 1);CHKERRQ(ierr); 2920 ierr = PetscLayoutSetUp(vLayout);CHKERRQ(ierr); 2921 /* Count locally unique vertices */ 2922 { 2923 PetscHSetI vhash; 2924 PetscInt off = 0; 2925 2926 ierr = PetscHSetICreate(&vhash);CHKERRQ(ierr); 2927 for (c = 0; c < numCells; ++c) { 2928 for (p = 0; p < numCorners; ++p) { 2929 ierr = PetscHSetIAdd(vhash, cells[c*numCorners+p]);CHKERRQ(ierr); 2930 } 2931 } 2932 ierr = PetscHSetIGetSize(vhash, &numVerticesAdj);CHKERRQ(ierr); 2933 ierr = PetscMalloc1(numVerticesAdj, &verticesAdj);CHKERRQ(ierr); 2934 ierr = PetscHSetIGetElems(vhash, &off, verticesAdj);CHKERRQ(ierr); 2935 ierr = PetscHSetIDestroy(&vhash);CHKERRQ(ierr); 2936 if (off != numVerticesAdj) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid number of local vertices %D should be %D", off, numVerticesAdj); 2937 } 2938 ierr = PetscSortInt(numVerticesAdj, verticesAdj);CHKERRQ(ierr); 2939 /* Create cones */ 2940 ierr = DMPlexSetChart(dm, 0, numCells+numVerticesAdj);CHKERRQ(ierr); 2941 for (c = 0; c < numCells; ++c) {ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr);} 2942 ierr = DMSetUp(dm);CHKERRQ(ierr); 2943 ierr = DMPlexGetCones(dm,&cones);CHKERRQ(ierr); 2944 for (c = 0; c < numCells; ++c) { 2945 for (p = 0; p < numCorners; ++p) { 2946 const PetscInt gv = cells[c*numCorners+p]; 2947 PetscInt lv; 2948 2949 /* Positions within verticesAdj form 0-based local vertex numbering; 2950 we need to shift it by numCells to get correct DAG points (cells go first) */ 2951 ierr = PetscFindInt(gv, numVerticesAdj, verticesAdj, &lv);CHKERRQ(ierr); 2952 if (lv < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not find global vertex %D in local connectivity", gv); 2953 cones[c*numCorners+p] = lv+numCells; 2954 } 2955 } 2956 /* Create SF for vertices */ 2957 ierr = PetscSFCreate(PetscObjectComm((PetscObject)dm), &sfVert);CHKERRQ(ierr); 2958 ierr = PetscObjectSetName((PetscObject) sfVert, "Vertex Ownership SF");CHKERRQ(ierr); 2959 ierr = PetscSFSetFromOptions(sfVert);CHKERRQ(ierr); 2960 ierr = PetscSFSetGraphLayout(sfVert, vLayout, numVerticesAdj, NULL, PETSC_OWN_POINTER, verticesAdj);CHKERRQ(ierr); 2961 ierr = PetscFree(verticesAdj);CHKERRQ(ierr); 2962 /* Build pointSF */ 2963 ierr = PetscMalloc2(numVerticesAdj, &vertexLocal, numVertices, &vertexOwner);CHKERRQ(ierr); 2964 for (v = 0; v < numVerticesAdj; ++v) {vertexLocal[v].index = v+numCells; vertexLocal[v].rank = rank;} 2965 for (v = 0; v < numVertices; ++v) {vertexOwner[v].index = -1; vertexOwner[v].rank = -1;} 2966 ierr = PetscSFReduceBegin(sfVert, MPIU_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);CHKERRQ(ierr); 2967 ierr = PetscSFReduceEnd(sfVert, MPIU_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);CHKERRQ(ierr); 2968 for (v = 0; v < numVertices; ++v) if (vertexOwner[v].rank < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Global vertex %D was unclaimed", v + vLayout->rstart); 2969 ierr = PetscSFBcastBegin(sfVert, MPIU_2INT, vertexOwner, vertexLocal);CHKERRQ(ierr); 2970 ierr = PetscSFBcastEnd(sfVert, MPIU_2INT, vertexOwner, vertexLocal);CHKERRQ(ierr); 2971 for (v = 0; v < numVerticesAdj; ++v) if (vertexLocal[v].rank != rank) ++numVerticesGhost; 2972 ierr = PetscMalloc1(numVerticesGhost, &localVertex);CHKERRQ(ierr); 2973 ierr = PetscMalloc1(numVerticesGhost, &remoteVertex);CHKERRQ(ierr); 2974 for (v = 0, g = 0; v < numVerticesAdj; ++v) { 2975 if (vertexLocal[v].rank != rank) { 2976 localVertex[g] = v+numCells; 2977 remoteVertex[g].index = vertexLocal[v].index; 2978 remoteVertex[g].rank = vertexLocal[v].rank; 2979 ++g; 2980 } 2981 } 2982 ierr = PetscFree2(vertexLocal, vertexOwner);CHKERRQ(ierr); 2983 if (g != numVerticesGhost) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of vertex ghosts %D should be %D", g, numVerticesGhost); 2984 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 2985 ierr = PetscObjectSetName((PetscObject) sfPoint, "point SF");CHKERRQ(ierr); 2986 ierr = PetscSFSetGraph(sfPoint, numCells+numVerticesAdj, numVerticesGhost, localVertex, PETSC_OWN_POINTER, remoteVertex, PETSC_OWN_POINTER);CHKERRQ(ierr); 2987 ierr = PetscLayoutDestroy(&vLayout);CHKERRQ(ierr); 2988 /* Fill in the rest of the topology structure */ 2989 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 2990 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 2991 ierr = PetscLogEventEnd(DMPLEX_BuildFromCellList,dm,0,0,0);CHKERRQ(ierr); 2992 if (vertexSF) *vertexSF = sfVert; 2993 else {ierr = PetscSFDestroy(&sfVert);CHKERRQ(ierr);} 2994 PetscFunctionReturn(0); 2995 } 2996 2997 /*@C 2998 DMPlexBuildCoordinatesFromCellListParallel - Build DM coordinates from a list of coordinates for each owned vertex (common mesh generator output) 2999 3000 Input Parameters: 3001 + dm - The DM 3002 . spaceDim - The spatial dimension used for coordinates 3003 . sfVert - SF describing complete vertex ownership 3004 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 3005 3006 Level: advanced 3007 3008 Notes: 3009 Not currently supported in Fortran. 3010 3011 .seealso: DMPlexBuildCoordinatesFromCellList(), DMPlexCreateFromCellListParallelPetsc(), DMPlexBuildFromCellListParallel() 3012 @*/ 3013 PetscErrorCode DMPlexBuildCoordinatesFromCellListParallel(DM dm, PetscInt spaceDim, PetscSF sfVert, const PetscReal vertexCoords[]) 3014 { 3015 PetscSection coordSection; 3016 Vec coordinates; 3017 PetscScalar *coords; 3018 PetscInt numVertices, numVerticesAdj, coordSize, v, vStart, vEnd; 3019 PetscErrorCode ierr; 3020 3021 PetscFunctionBegin; 3022 ierr = PetscLogEventBegin(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0);CHKERRQ(ierr); 3023 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 3024 if (vStart < 0 || vEnd < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM is not set up properly. DMPlexBuildFromCellList() should be called first."); 3025 ierr = DMSetCoordinateDim(dm, spaceDim);CHKERRQ(ierr); 3026 ierr = PetscSFGetGraph(sfVert, &numVertices, &numVerticesAdj, NULL, NULL);CHKERRQ(ierr); 3027 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); 3028 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 3029 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 3030 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 3031 ierr = PetscSectionSetChart(coordSection, vStart, vEnd);CHKERRQ(ierr); 3032 for (v = vStart; v < vEnd; ++v) { 3033 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 3034 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 3035 } 3036 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 3037 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 3038 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 3039 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 3040 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 3041 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 3042 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 3043 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 3044 { 3045 MPI_Datatype coordtype; 3046 3047 /* Need a temp buffer for coords if we have complex/single */ 3048 ierr = MPI_Type_contiguous(spaceDim, MPIU_SCALAR, &coordtype);CHKERRQ(ierr); 3049 ierr = MPI_Type_commit(&coordtype);CHKERRQ(ierr); 3050 #if defined(PETSC_USE_COMPLEX) 3051 { 3052 PetscScalar *svertexCoords; 3053 PetscInt i; 3054 ierr = PetscMalloc1(numVertices*spaceDim,&svertexCoords);CHKERRQ(ierr); 3055 for (i=0; i<numVertices*spaceDim; i++) svertexCoords[i] = vertexCoords[i]; 3056 ierr = PetscSFBcastBegin(sfVert, coordtype, svertexCoords, coords);CHKERRQ(ierr); 3057 ierr = PetscSFBcastEnd(sfVert, coordtype, svertexCoords, coords);CHKERRQ(ierr); 3058 ierr = PetscFree(svertexCoords);CHKERRQ(ierr); 3059 } 3060 #else 3061 ierr = PetscSFBcastBegin(sfVert, coordtype, vertexCoords, coords);CHKERRQ(ierr); 3062 ierr = PetscSFBcastEnd(sfVert, coordtype, vertexCoords, coords);CHKERRQ(ierr); 3063 #endif 3064 ierr = MPI_Type_free(&coordtype);CHKERRQ(ierr); 3065 } 3066 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 3067 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 3068 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 3069 ierr = PetscLogEventEnd(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0);CHKERRQ(ierr); 3070 PetscFunctionReturn(0); 3071 } 3072 3073 /*@ 3074 DMPlexCreateFromCellListParallelPetsc - Create distributed DMPLEX from a list of vertices for each cell (common mesh generator output) 3075 3076 Input Parameters: 3077 + comm - The communicator 3078 . dim - The topological dimension of the mesh 3079 . numCells - The number of cells owned by this process 3080 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE 3081 . NVertices - The global number of vertices, or PETSC_DECIDE 3082 . numCorners - The number of vertices for each cell 3083 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 3084 . cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 3085 . spaceDim - The spatial dimension used for coordinates 3086 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 3087 3088 Output Parameter: 3089 + dm - The DM 3090 - vertexSF - (Optional) SF describing complete vertex ownership 3091 3092 Notes: 3093 This function is just a convenient sequence of DMCreate(), DMSetType(), DMSetDimension(), 3094 DMPlexBuildFromCellListParallel(), DMPlexInterpolate(), DMPlexBuildCoordinatesFromCellListParallel() 3095 3096 See DMPlexBuildFromCellListParallel() for an example and details about the topology-related parameters. 3097 See DMPlexBuildCoordinatesFromCellListParallel() for details about the geometry-related parameters. 3098 3099 Level: intermediate 3100 3101 .seealso: DMPlexCreateFromCellListPetsc(), DMPlexBuildFromCellListParallel(), DMPlexBuildCoordinatesFromCellListParallel(), DMPlexCreateFromDAG(), DMPlexCreate() 3102 @*/ 3103 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) 3104 { 3105 PetscSF sfVert; 3106 PetscErrorCode ierr; 3107 3108 PetscFunctionBegin; 3109 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 3110 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 3111 PetscValidLogicalCollectiveInt(*dm, dim, 2); 3112 PetscValidLogicalCollectiveInt(*dm, spaceDim, 8); 3113 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 3114 ierr = DMPlexBuildFromCellListParallel(*dm, numCells, numVertices, NVertices, numCorners, cells, &sfVert);CHKERRQ(ierr); 3115 if (interpolate) { 3116 DM idm; 3117 3118 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 3119 ierr = DMDestroy(dm);CHKERRQ(ierr); 3120 *dm = idm; 3121 } 3122 ierr = DMPlexBuildCoordinatesFromCellListParallel(*dm, spaceDim, sfVert, vertexCoords);CHKERRQ(ierr); 3123 if (vertexSF) *vertexSF = sfVert; 3124 else {ierr = PetscSFDestroy(&sfVert);CHKERRQ(ierr);} 3125 PetscFunctionReturn(0); 3126 } 3127 3128 3129 /*@ 3130 DMPlexCreateFromCellListParallel - Deprecated, use DMPlexCreateFromCellListParallelPetsc() 3131 3132 Level: deprecated 3133 3134 .seealso: DMPlexCreateFromCellListParallelPetsc() 3135 @*/ 3136 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) 3137 { 3138 PetscErrorCode ierr; 3139 PetscInt i; 3140 PetscInt *pintCells; 3141 3142 PetscFunctionBegin; 3143 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)); 3144 if (sizeof(int) == sizeof(PetscInt)) { 3145 pintCells = (PetscInt *) cells; 3146 } else { 3147 ierr = PetscMalloc1(numCells*numCorners, &pintCells);CHKERRQ(ierr); 3148 for (i = 0; i < numCells*numCorners; i++) { 3149 pintCells[i] = (PetscInt) cells[i]; 3150 } 3151 } 3152 ierr = DMPlexCreateFromCellListParallelPetsc(comm, dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, pintCells, spaceDim, vertexCoords, vertexSF, dm);CHKERRQ(ierr); 3153 if (sizeof(int) != sizeof(PetscInt)) { 3154 ierr = PetscFree(pintCells);CHKERRQ(ierr); 3155 } 3156 PetscFunctionReturn(0); 3157 } 3158 3159 /*@C 3160 DMPlexBuildFromCellList - Build DMPLEX topology from a list of vertices for each cell (common mesh generator output) 3161 3162 Input Parameters: 3163 + dm - The DM 3164 . numCells - The number of cells owned by this process 3165 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE 3166 . numCorners - The number of vertices for each cell 3167 - cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 3168 3169 Level: advanced 3170 3171 Notes: 3172 Two triangles sharing a face 3173 $ 3174 $ 2 3175 $ / | \ 3176 $ / | \ 3177 $ / | \ 3178 $ 0 0 | 1 3 3179 $ \ | / 3180 $ \ | / 3181 $ \ | / 3182 $ 1 3183 would have input 3184 $ numCells = 2, numVertices = 4 3185 $ cells = [0 1 2 1 3 2] 3186 $ 3187 which would result in the DMPlex 3188 $ 3189 $ 4 3190 $ / | \ 3191 $ / | \ 3192 $ / | \ 3193 $ 2 0 | 1 5 3194 $ \ | / 3195 $ \ | / 3196 $ \ | / 3197 $ 3 3198 3199 If numVertices is PETSC_DECIDE, it is computed by PETSc as the maximum vertex index in cells + 1. 3200 3201 Not currently supported in Fortran. 3202 3203 .seealso: DMPlexBuildFromCellListParallel(), DMPlexBuildCoordinatesFromCellList(), DMPlexCreateFromCellListPetsc() 3204 @*/ 3205 PetscErrorCode DMPlexBuildFromCellList(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const PetscInt cells[]) 3206 { 3207 PetscInt *cones, c, p, dim; 3208 PetscErrorCode ierr; 3209 3210 PetscFunctionBegin; 3211 ierr = PetscLogEventBegin(DMPLEX_BuildFromCellList,dm,0,0,0);CHKERRQ(ierr); 3212 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3213 /* Get/check global number of vertices */ 3214 { 3215 PetscInt NVerticesInCells, i; 3216 const PetscInt len = numCells * numCorners; 3217 3218 /* NVerticesInCells = max(cells) + 1 */ 3219 NVerticesInCells = PETSC_MIN_INT; 3220 for (i=0; i<len; i++) if (cells[i] > NVerticesInCells) NVerticesInCells = cells[i]; 3221 ++NVerticesInCells; 3222 3223 if (numVertices == PETSC_DECIDE) numVertices = NVerticesInCells; 3224 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); 3225 } 3226 ierr = DMPlexSetChart(dm, 0, numCells+numVertices);CHKERRQ(ierr); 3227 for (c = 0; c < numCells; ++c) { 3228 ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr); 3229 } 3230 ierr = DMSetUp(dm);CHKERRQ(ierr); 3231 ierr = DMPlexGetCones(dm,&cones);CHKERRQ(ierr); 3232 for (c = 0; c < numCells; ++c) { 3233 for (p = 0; p < numCorners; ++p) { 3234 cones[c*numCorners+p] = cells[c*numCorners+p]+numCells; 3235 } 3236 } 3237 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 3238 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 3239 ierr = PetscLogEventEnd(DMPLEX_BuildFromCellList,dm,0,0,0);CHKERRQ(ierr); 3240 PetscFunctionReturn(0); 3241 } 3242 3243 /*@C 3244 DMPlexBuildCoordinatesFromCellList - Build DM coordinates from a list of coordinates for each owned vertex (common mesh generator output) 3245 3246 Input Parameters: 3247 + dm - The DM 3248 . spaceDim - The spatial dimension used for coordinates 3249 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 3250 3251 Level: advanced 3252 3253 Notes: 3254 Not currently supported in Fortran. 3255 3256 .seealso: DMPlexBuildCoordinatesFromCellListParallel(), DMPlexCreateFromCellListPetsc(), DMPlexBuildFromCellList() 3257 @*/ 3258 PetscErrorCode DMPlexBuildCoordinatesFromCellList(DM dm, PetscInt spaceDim, const PetscReal vertexCoords[]) 3259 { 3260 PetscSection coordSection; 3261 Vec coordinates; 3262 DM cdm; 3263 PetscScalar *coords; 3264 PetscInt v, vStart, vEnd, d; 3265 PetscErrorCode ierr; 3266 3267 PetscFunctionBegin; 3268 ierr = PetscLogEventBegin(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0);CHKERRQ(ierr); 3269 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 3270 if (vStart < 0 || vEnd < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM is not set up properly. DMPlexBuildFromCellList() should be called first."); 3271 ierr = DMSetCoordinateDim(dm, spaceDim);CHKERRQ(ierr); 3272 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 3273 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 3274 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 3275 ierr = PetscSectionSetChart(coordSection, vStart, vEnd);CHKERRQ(ierr); 3276 for (v = vStart; v < vEnd; ++v) { 3277 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 3278 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 3279 } 3280 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 3281 3282 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 3283 ierr = DMCreateLocalVector(cdm, &coordinates);CHKERRQ(ierr); 3284 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 3285 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 3286 ierr = VecGetArrayWrite(coordinates, &coords);CHKERRQ(ierr); 3287 for (v = 0; v < vEnd-vStart; ++v) { 3288 for (d = 0; d < spaceDim; ++d) { 3289 coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d]; 3290 } 3291 } 3292 ierr = VecRestoreArrayWrite(coordinates, &coords);CHKERRQ(ierr); 3293 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 3294 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 3295 ierr = PetscLogEventEnd(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0);CHKERRQ(ierr); 3296 PetscFunctionReturn(0); 3297 } 3298 3299 /*@ 3300 DMPlexCreateFromCellListPetsc - Create DMPLEX from a list of vertices for each cell (common mesh generator output) 3301 3302 Input Parameters: 3303 + comm - The communicator 3304 . dim - The topological dimension of the mesh 3305 . numCells - The number of cells 3306 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE 3307 . numCorners - The number of vertices for each cell 3308 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 3309 . cells - An array of numCells*numCorners numbers, the vertices for each cell 3310 . spaceDim - The spatial dimension used for coordinates 3311 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 3312 3313 Output Parameter: 3314 . dm - The DM 3315 3316 Notes: 3317 This function is just a convenient sequence of DMCreate(), DMSetType(), DMSetDimension(), DMPlexBuildFromCellList(), 3318 DMPlexInterpolate(), DMPlexBuildCoordinatesFromCellList() 3319 3320 See DMPlexBuildFromCellList() for an example and details about the topology-related parameters. 3321 See DMPlexBuildCoordinatesFromCellList() for details about the geometry-related parameters. 3322 3323 Level: intermediate 3324 3325 .seealso: DMPlexCreateFromCellListParallelPetsc(), DMPlexBuildFromCellList(), DMPlexBuildCoordinatesFromCellList(), DMPlexCreateFromDAG(), DMPlexCreate() 3326 @*/ 3327 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) 3328 { 3329 PetscErrorCode ierr; 3330 3331 PetscFunctionBegin; 3332 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."); 3333 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 3334 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 3335 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 3336 ierr = DMPlexBuildFromCellList(*dm, numCells, numVertices, numCorners, cells);CHKERRQ(ierr); 3337 if (interpolate) { 3338 DM idm; 3339 3340 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 3341 ierr = DMDestroy(dm);CHKERRQ(ierr); 3342 *dm = idm; 3343 } 3344 ierr = DMPlexBuildCoordinatesFromCellList(*dm, spaceDim, vertexCoords);CHKERRQ(ierr); 3345 PetscFunctionReturn(0); 3346 } 3347 3348 /*@ 3349 DMPlexCreateFromCellList - Deprecated, use DMPlexCreateFromCellListPetsc() 3350 3351 Level: deprecated 3352 3353 .seealso: DMPlexCreateFromCellListPetsc() 3354 @*/ 3355 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) 3356 { 3357 PetscErrorCode ierr; 3358 PetscInt i; 3359 PetscInt *pintCells; 3360 PetscReal *prealVC; 3361 3362 PetscFunctionBegin; 3363 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)); 3364 if (sizeof(int) == sizeof(PetscInt)) { 3365 pintCells = (PetscInt *) cells; 3366 } else { 3367 ierr = PetscMalloc1(numCells*numCorners, &pintCells);CHKERRQ(ierr); 3368 for (i = 0; i < numCells*numCorners; i++) { 3369 pintCells[i] = (PetscInt) cells[i]; 3370 } 3371 } 3372 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)); 3373 if (sizeof(double) == sizeof(PetscReal)) { 3374 prealVC = (PetscReal *) vertexCoords; 3375 } else { 3376 ierr = PetscMalloc1(numVertices*spaceDim, &prealVC);CHKERRQ(ierr); 3377 for (i = 0; i < numVertices*spaceDim; i++) { 3378 prealVC[i] = (PetscReal) vertexCoords[i]; 3379 } 3380 } 3381 ierr = DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, pintCells, spaceDim, prealVC, dm);CHKERRQ(ierr); 3382 if (sizeof(int) != sizeof(PetscInt)) { 3383 ierr = PetscFree(pintCells);CHKERRQ(ierr); 3384 } 3385 if (sizeof(double) != sizeof(PetscReal)) { 3386 ierr = PetscFree(prealVC);CHKERRQ(ierr); 3387 } 3388 PetscFunctionReturn(0); 3389 } 3390 3391 /*@ 3392 DMPlexCreateFromDAG - This takes as input the adjacency-list representation of the Directed Acyclic Graph (Hasse Diagram) encoding a mesh, and produces a DM 3393 3394 Input Parameters: 3395 + dm - The empty DM object, usually from DMCreate() and DMSetDimension() 3396 . depth - The depth of the DAG 3397 . numPoints - Array of size depth + 1 containing the number of points at each depth 3398 . coneSize - The cone size of each point 3399 . cones - The concatenation of the cone points for each point, the cone list must be oriented correctly for each point 3400 . coneOrientations - The orientation of each cone point 3401 - vertexCoords - An array of numPoints[0]*spacedim numbers representing the coordinates of each vertex, with spacedim the value set via DMSetCoordinateDim() 3402 3403 Output Parameter: 3404 . dm - The DM 3405 3406 Note: Two triangles sharing a face would have input 3407 $ depth = 1, numPoints = [4 2], coneSize = [3 3 0 0 0 0] 3408 $ cones = [2 3 4 3 5 4], coneOrientations = [0 0 0 0 0 0] 3409 $ vertexCoords = [-1.0 0.0 0.0 -1.0 0.0 1.0 1.0 0.0] 3410 $ 3411 which would result in the DMPlex 3412 $ 3413 $ 4 3414 $ / | \ 3415 $ / | \ 3416 $ / | \ 3417 $ 2 0 | 1 5 3418 $ \ | / 3419 $ \ | / 3420 $ \ | / 3421 $ 3 3422 $ 3423 $ Notice that all points are numbered consecutively, unlike DMPlexCreateFromCellListPetsc() 3424 3425 Level: advanced 3426 3427 .seealso: DMPlexCreateFromCellListPetsc(), DMPlexCreate() 3428 @*/ 3429 PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[]) 3430 { 3431 Vec coordinates; 3432 PetscSection coordSection; 3433 PetscScalar *coords; 3434 PetscInt coordSize, firstVertex = -1, pStart = 0, pEnd = 0, p, v, dim, dimEmbed, d, off; 3435 PetscErrorCode ierr; 3436 3437 PetscFunctionBegin; 3438 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3439 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 3440 if (dimEmbed < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Embedding dimension %D cannot be less than intrinsic dimension %d",dimEmbed,dim); 3441 for (d = 0; d <= depth; ++d) pEnd += numPoints[d]; 3442 ierr = DMPlexSetChart(dm, pStart, pEnd);CHKERRQ(ierr); 3443 for (p = pStart; p < pEnd; ++p) { 3444 ierr = DMPlexSetConeSize(dm, p, coneSize[p-pStart]);CHKERRQ(ierr); 3445 if (firstVertex < 0 && !coneSize[p - pStart]) { 3446 firstVertex = p - pStart; 3447 } 3448 } 3449 if (firstVertex < 0 && numPoints[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected %D vertices but could not find any", numPoints[0]); 3450 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 3451 for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) { 3452 ierr = DMPlexSetCone(dm, p, &cones[off]);CHKERRQ(ierr); 3453 ierr = DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);CHKERRQ(ierr); 3454 } 3455 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 3456 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 3457 /* Build coordinates */ 3458 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 3459 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 3460 ierr = PetscSectionSetFieldComponents(coordSection, 0, dimEmbed);CHKERRQ(ierr); 3461 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0]);CHKERRQ(ierr); 3462 for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) { 3463 ierr = PetscSectionSetDof(coordSection, v, dimEmbed);CHKERRQ(ierr); 3464 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed);CHKERRQ(ierr); 3465 } 3466 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 3467 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 3468 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 3469 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 3470 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 3471 ierr = VecSetBlockSize(coordinates, dimEmbed);CHKERRQ(ierr); 3472 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 3473 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 3474 for (v = 0; v < numPoints[0]; ++v) { 3475 PetscInt off; 3476 3477 ierr = PetscSectionGetOffset(coordSection, v+firstVertex, &off);CHKERRQ(ierr); 3478 for (d = 0; d < dimEmbed; ++d) { 3479 coords[off+d] = vertexCoords[v*dimEmbed+d]; 3480 } 3481 } 3482 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 3483 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 3484 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 3485 PetscFunctionReturn(0); 3486 } 3487 3488 /*@C 3489 DMPlexCreateCellVertexFromFile - Create a DMPlex mesh from a simple cell-vertex file. 3490 3491 + comm - The MPI communicator 3492 . filename - Name of the .dat file 3493 - interpolate - Create faces and edges in the mesh 3494 3495 Output Parameter: 3496 . dm - The DM object representing the mesh 3497 3498 Note: The format is the simplest possible: 3499 $ Ne 3500 $ v0 v1 ... vk 3501 $ Nv 3502 $ x y z marker 3503 3504 Level: beginner 3505 3506 .seealso: DMPlexCreateFromFile(), DMPlexCreateMedFromFile(), DMPlexCreateGmsh(), DMPlexCreate() 3507 @*/ 3508 PetscErrorCode DMPlexCreateCellVertexFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 3509 { 3510 DMLabel marker; 3511 PetscViewer viewer; 3512 Vec coordinates; 3513 PetscSection coordSection; 3514 PetscScalar *coords; 3515 char line[PETSC_MAX_PATH_LEN]; 3516 PetscInt dim = 3, cdim = 3, coordSize, v, c, d; 3517 PetscMPIInt rank; 3518 int snum, Nv, Nc; 3519 PetscErrorCode ierr; 3520 3521 PetscFunctionBegin; 3522 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 3523 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 3524 ierr = PetscViewerSetType(viewer, PETSCVIEWERASCII);CHKERRQ(ierr); 3525 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 3526 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 3527 if (!rank) { 3528 ierr = PetscViewerRead(viewer, line, 2, NULL, PETSC_STRING);CHKERRQ(ierr); 3529 snum = sscanf(line, "%d %d", &Nc, &Nv); 3530 if (snum != 2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 3531 } else { 3532 Nc = Nv = 0; 3533 } 3534 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 3535 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 3536 ierr = DMPlexSetChart(*dm, 0, Nc+Nv);CHKERRQ(ierr); 3537 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 3538 ierr = DMSetCoordinateDim(*dm, cdim);CHKERRQ(ierr); 3539 /* Read topology */ 3540 if (!rank) { 3541 PetscInt cone[8], corners = 8; 3542 int vbuf[8], v; 3543 3544 for (c = 0; c < Nc; ++c) {ierr = DMPlexSetConeSize(*dm, c, corners);CHKERRQ(ierr);} 3545 ierr = DMSetUp(*dm);CHKERRQ(ierr); 3546 for (c = 0; c < Nc; ++c) { 3547 ierr = PetscViewerRead(viewer, line, corners, NULL, PETSC_STRING);CHKERRQ(ierr); 3548 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]); 3549 if (snum != corners) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 3550 for (v = 0; v < corners; ++v) cone[v] = vbuf[v] + Nc; 3551 /* Hexahedra are inverted */ 3552 { 3553 PetscInt tmp = cone[1]; 3554 cone[1] = cone[3]; 3555 cone[3] = tmp; 3556 } 3557 ierr = DMPlexSetCone(*dm, c, cone);CHKERRQ(ierr); 3558 } 3559 } 3560 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 3561 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 3562 /* Read coordinates */ 3563 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 3564 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 3565 ierr = PetscSectionSetFieldComponents(coordSection, 0, cdim);CHKERRQ(ierr); 3566 ierr = PetscSectionSetChart(coordSection, Nc, Nc + Nv);CHKERRQ(ierr); 3567 for (v = Nc; v < Nc+Nv; ++v) { 3568 ierr = PetscSectionSetDof(coordSection, v, cdim);CHKERRQ(ierr); 3569 ierr = PetscSectionSetFieldDof(coordSection, v, 0, cdim);CHKERRQ(ierr); 3570 } 3571 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 3572 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 3573 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 3574 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 3575 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 3576 ierr = VecSetBlockSize(coordinates, cdim);CHKERRQ(ierr); 3577 ierr = VecSetType(coordinates, VECSTANDARD);CHKERRQ(ierr); 3578 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 3579 if (!rank) { 3580 double x[3]; 3581 int val; 3582 3583 ierr = DMCreateLabel(*dm, "marker");CHKERRQ(ierr); 3584 ierr = DMGetLabel(*dm, "marker", &marker);CHKERRQ(ierr); 3585 for (v = 0; v < Nv; ++v) { 3586 ierr = PetscViewerRead(viewer, line, 4, NULL, PETSC_STRING);CHKERRQ(ierr); 3587 snum = sscanf(line, "%lg %lg %lg %d", &x[0], &x[1], &x[2], &val); 3588 if (snum != 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 3589 for (d = 0; d < cdim; ++d) coords[v*cdim+d] = x[d]; 3590 if (val) {ierr = DMLabelSetValue(marker, v+Nc, val);CHKERRQ(ierr);} 3591 } 3592 } 3593 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 3594 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 3595 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 3596 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3597 if (interpolate) { 3598 DM idm; 3599 DMLabel bdlabel; 3600 3601 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 3602 ierr = DMDestroy(dm);CHKERRQ(ierr); 3603 *dm = idm; 3604 3605 ierr = DMGetLabel(*dm, "marker", &bdlabel);CHKERRQ(ierr); 3606 ierr = DMPlexMarkBoundaryFaces(*dm, PETSC_DETERMINE, bdlabel);CHKERRQ(ierr); 3607 ierr = DMPlexLabelComplete(*dm, bdlabel);CHKERRQ(ierr); 3608 } 3609 PetscFunctionReturn(0); 3610 } 3611 3612 /*@C 3613 DMPlexCreateFromFile - This takes a filename and produces a DM 3614 3615 Input Parameters: 3616 + comm - The communicator 3617 . filename - A file name 3618 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 3619 3620 Output Parameter: 3621 . dm - The DM 3622 3623 Options Database Keys: 3624 . -dm_plex_create_from_hdf5_xdmf - use the PETSC_VIEWER_HDF5_XDMF format for reading HDF5 3625 3626 Use -dm_plex_create_ prefix to pass options to the internal PetscViewer, e.g. 3627 $ -dm_plex_create_viewer_hdf5_collective 3628 3629 Level: beginner 3630 3631 .seealso: DMPlexCreateFromDAG(), DMPlexCreateFromCellListPetsc(), DMPlexCreate() 3632 @*/ 3633 PetscErrorCode DMPlexCreateFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 3634 { 3635 const char *extGmsh = ".msh"; 3636 const char *extGmsh2 = ".msh2"; 3637 const char *extGmsh4 = ".msh4"; 3638 const char *extCGNS = ".cgns"; 3639 const char *extExodus = ".exo"; 3640 const char *extGenesis = ".gen"; 3641 const char *extFluent = ".cas"; 3642 const char *extHDF5 = ".h5"; 3643 const char *extMed = ".med"; 3644 const char *extPLY = ".ply"; 3645 const char *extCV = ".dat"; 3646 size_t len; 3647 PetscBool isGmsh, isGmsh2, isGmsh4, isCGNS, isExodus, isGenesis, isFluent, isHDF5, isMed, isPLY, isCV; 3648 PetscMPIInt rank; 3649 PetscErrorCode ierr; 3650 3651 PetscFunctionBegin; 3652 PetscValidCharPointer(filename, 2); 3653 PetscValidPointer(dm, 4); 3654 ierr = DMInitializePackage();CHKERRQ(ierr); 3655 ierr = PetscLogEventBegin(DMPLEX_CreateFromFile,0,0,0,0);CHKERRQ(ierr); 3656 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 3657 ierr = PetscStrlen(filename, &len);CHKERRQ(ierr); 3658 if (!len) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Filename must be a valid path"); 3659 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGmsh, 4, &isGmsh);CHKERRQ(ierr); 3660 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extGmsh2, 5, &isGmsh2);CHKERRQ(ierr); 3661 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extGmsh4, 5, &isGmsh4);CHKERRQ(ierr); 3662 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extCGNS, 5, &isCGNS);CHKERRQ(ierr); 3663 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extExodus, 4, &isExodus);CHKERRQ(ierr); 3664 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGenesis, 4, &isGenesis);CHKERRQ(ierr); 3665 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extFluent, 4, &isFluent);CHKERRQ(ierr); 3666 ierr = PetscStrncmp(&filename[PetscMax(0,len-3)], extHDF5, 3, &isHDF5);CHKERRQ(ierr); 3667 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extMed, 4, &isMed);CHKERRQ(ierr); 3668 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extPLY, 4, &isPLY);CHKERRQ(ierr); 3669 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extCV, 4, &isCV);CHKERRQ(ierr); 3670 if (isGmsh || isGmsh2 || isGmsh4) { 3671 ierr = DMPlexCreateGmshFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3672 } else if (isCGNS) { 3673 ierr = DMPlexCreateCGNSFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3674 } else if (isExodus || isGenesis) { 3675 ierr = DMPlexCreateExodusFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3676 } else if (isFluent) { 3677 ierr = DMPlexCreateFluentFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3678 } else if (isHDF5) { 3679 PetscBool load_hdf5_xdmf = PETSC_FALSE; 3680 PetscViewer viewer; 3681 3682 /* PETSC_VIEWER_HDF5_XDMF is used if the filename ends with .xdmf.h5, or if -dm_plex_create_from_hdf5_xdmf option is present */ 3683 ierr = PetscStrncmp(&filename[PetscMax(0,len-8)], ".xdmf", 5, &load_hdf5_xdmf);CHKERRQ(ierr); 3684 ierr = PetscOptionsGetBool(NULL, NULL, "-dm_plex_create_from_hdf5_xdmf", &load_hdf5_xdmf, NULL);CHKERRQ(ierr); 3685 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 3686 ierr = PetscViewerSetType(viewer, PETSCVIEWERHDF5);CHKERRQ(ierr); 3687 ierr = PetscViewerSetOptionsPrefix(viewer, "dm_plex_create_");CHKERRQ(ierr); 3688 ierr = PetscViewerSetFromOptions(viewer);CHKERRQ(ierr); 3689 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 3690 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 3691 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 3692 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 3693 if (load_hdf5_xdmf) {ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_XDMF);CHKERRQ(ierr);} 3694 ierr = DMLoad(*dm, viewer);CHKERRQ(ierr); 3695 if (load_hdf5_xdmf) {ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);} 3696 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3697 3698 if (interpolate) { 3699 DM idm; 3700 3701 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 3702 ierr = DMDestroy(dm);CHKERRQ(ierr); 3703 *dm = idm; 3704 } 3705 } else if (isMed) { 3706 ierr = DMPlexCreateMedFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3707 } else if (isPLY) { 3708 ierr = DMPlexCreatePLYFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3709 } else if (isCV) { 3710 ierr = DMPlexCreateCellVertexFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 3711 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot load file %s: unrecognized extension", filename); 3712 ierr = PetscLogEventEnd(DMPLEX_CreateFromFile,0,0,0,0);CHKERRQ(ierr); 3713 PetscFunctionReturn(0); 3714 } 3715 3716 /*@ 3717 DMPlexCreateReferenceCellByType - Create a DMPLEX with the appropriate FEM reference cell 3718 3719 Collective 3720 3721 Input Parameters: 3722 + comm - The communicator 3723 - ct - The cell type of the reference cell 3724 3725 Output Parameter: 3726 . refdm - The reference cell 3727 3728 Level: intermediate 3729 3730 .seealso: DMPlexCreateReferenceCell(), DMPlexCreateBoxMesh() 3731 @*/ 3732 PetscErrorCode DMPlexCreateReferenceCellByType(MPI_Comm comm, DMPolytopeType ct, DM *refdm) 3733 { 3734 DM rdm; 3735 Vec coords; 3736 PetscErrorCode ierr; 3737 3738 PetscFunctionBegin; 3739 ierr = DMCreate(comm, &rdm);CHKERRQ(ierr); 3740 ierr = DMSetType(rdm, DMPLEX);CHKERRQ(ierr); 3741 switch (ct) { 3742 case DM_POLYTOPE_POINT: 3743 { 3744 PetscInt numPoints[1] = {1}; 3745 PetscInt coneSize[1] = {0}; 3746 PetscInt cones[1] = {0}; 3747 PetscInt coneOrientations[1] = {0}; 3748 PetscScalar vertexCoords[1] = {0.0}; 3749 3750 ierr = DMSetDimension(rdm, 0);CHKERRQ(ierr); 3751 ierr = DMPlexCreateFromDAG(rdm, 0, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3752 } 3753 break; 3754 case DM_POLYTOPE_SEGMENT: 3755 { 3756 PetscInt numPoints[2] = {2, 1}; 3757 PetscInt coneSize[3] = {2, 0, 0}; 3758 PetscInt cones[2] = {1, 2}; 3759 PetscInt coneOrientations[2] = {0, 0}; 3760 PetscScalar vertexCoords[2] = {-1.0, 1.0}; 3761 3762 ierr = DMSetDimension(rdm, 1);CHKERRQ(ierr); 3763 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3764 } 3765 break; 3766 case DM_POLYTOPE_TRIANGLE: 3767 { 3768 PetscInt numPoints[2] = {3, 1}; 3769 PetscInt coneSize[4] = {3, 0, 0, 0}; 3770 PetscInt cones[3] = {1, 2, 3}; 3771 PetscInt coneOrientations[3] = {0, 0, 0}; 3772 PetscScalar vertexCoords[6] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0}; 3773 3774 ierr = DMSetDimension(rdm, 2);CHKERRQ(ierr); 3775 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3776 } 3777 break; 3778 case DM_POLYTOPE_QUADRILATERAL: 3779 { 3780 PetscInt numPoints[2] = {4, 1}; 3781 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 3782 PetscInt cones[4] = {1, 2, 3, 4}; 3783 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 3784 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0}; 3785 3786 ierr = DMSetDimension(rdm, 2);CHKERRQ(ierr); 3787 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3788 } 3789 break; 3790 case DM_POLYTOPE_SEG_PRISM_TENSOR: 3791 { 3792 PetscInt numPoints[2] = {4, 1}; 3793 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 3794 PetscInt cones[4] = {1, 2, 3, 4}; 3795 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 3796 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0}; 3797 3798 ierr = DMSetDimension(rdm, 2);CHKERRQ(ierr); 3799 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3800 } 3801 break; 3802 case DM_POLYTOPE_TETRAHEDRON: 3803 { 3804 PetscInt numPoints[2] = {4, 1}; 3805 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 3806 PetscInt cones[4] = {1, 3, 2, 4}; 3807 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 3808 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}; 3809 3810 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 3811 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3812 } 3813 break; 3814 case DM_POLYTOPE_HEXAHEDRON: 3815 { 3816 PetscInt numPoints[2] = {8, 1}; 3817 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0}; 3818 PetscInt cones[8] = {1, 4, 3, 2, 5, 6, 7, 8}; 3819 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 3820 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, 3821 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0}; 3822 3823 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 3824 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3825 } 3826 break; 3827 case DM_POLYTOPE_TRI_PRISM: 3828 { 3829 PetscInt numPoints[2] = {6, 1}; 3830 PetscInt coneSize[7] = {6, 0, 0, 0, 0, 0, 0}; 3831 PetscInt cones[6] = {1, 3, 2, 4, 5, 6}; 3832 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 3833 PetscScalar vertexCoords[18] = {-1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 3834 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0}; 3835 3836 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 3837 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3838 } 3839 break; 3840 case DM_POLYTOPE_TRI_PRISM_TENSOR: 3841 { 3842 PetscInt numPoints[2] = {6, 1}; 3843 PetscInt coneSize[7] = {6, 0, 0, 0, 0, 0, 0}; 3844 PetscInt cones[6] = {1, 2, 3, 4, 5, 6}; 3845 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 3846 PetscScalar vertexCoords[18] = {-1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 3847 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0}; 3848 3849 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 3850 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3851 } 3852 break; 3853 case DM_POLYTOPE_QUAD_PRISM_TENSOR: 3854 { 3855 PetscInt numPoints[2] = {8, 1}; 3856 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0}; 3857 PetscInt cones[8] = {1, 2, 3, 4, 5, 6, 7, 8}; 3858 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 3859 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, 3860 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0}; 3861 3862 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 3863 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 3864 } 3865 break; 3866 default: SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Cannot create reference cell for cell type %s", DMPolytopeTypes[ct]); 3867 } 3868 { 3869 PetscInt Nv, v; 3870 3871 /* Must create the celltype label here so that we do not automatically try to compute the types */ 3872 ierr = DMCreateLabel(rdm, "celltype");CHKERRQ(ierr); 3873 ierr = DMPlexSetCellType(rdm, 0, ct);CHKERRQ(ierr); 3874 ierr = DMPlexGetChart(rdm, NULL, &Nv);CHKERRQ(ierr); 3875 for (v = 1; v < Nv; ++v) {ierr = DMPlexSetCellType(rdm, v, DM_POLYTOPE_POINT);CHKERRQ(ierr);} 3876 } 3877 ierr = DMPlexInterpolate(rdm, refdm);CHKERRQ(ierr); 3878 if (rdm->coordinateDM) { 3879 DM ncdm; 3880 PetscSection cs; 3881 PetscInt pEnd = -1; 3882 3883 ierr = DMGetLocalSection(rdm->coordinateDM, &cs);CHKERRQ(ierr); 3884 if (cs) {ierr = PetscSectionGetChart(cs, NULL, &pEnd);CHKERRQ(ierr);} 3885 if (pEnd >= 0) { 3886 ierr = DMClone(*refdm, &ncdm);CHKERRQ(ierr); 3887 ierr = DMCopyDisc(rdm->coordinateDM, ncdm);CHKERRQ(ierr); 3888 ierr = DMSetLocalSection(ncdm, cs);CHKERRQ(ierr); 3889 ierr = DMSetCoordinateDM(*refdm, ncdm);CHKERRQ(ierr); 3890 ierr = DMDestroy(&ncdm);CHKERRQ(ierr); 3891 } 3892 } 3893 ierr = DMGetCoordinatesLocal(rdm, &coords);CHKERRQ(ierr); 3894 if (coords) { 3895 ierr = DMSetCoordinatesLocal(*refdm, coords);CHKERRQ(ierr); 3896 } else { 3897 ierr = DMGetCoordinates(rdm, &coords);CHKERRQ(ierr); 3898 if (coords) {ierr = DMSetCoordinates(*refdm, coords);CHKERRQ(ierr);} 3899 } 3900 ierr = DMDestroy(&rdm);CHKERRQ(ierr); 3901 PetscFunctionReturn(0); 3902 } 3903 3904 /*@ 3905 DMPlexCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell 3906 3907 Collective 3908 3909 Input Parameters: 3910 + comm - The communicator 3911 . dim - The spatial dimension 3912 - simplex - Flag for simplex, otherwise use a tensor-product cell 3913 3914 Output Parameter: 3915 . refdm - The reference cell 3916 3917 Level: intermediate 3918 3919 .seealso: DMPlexCreateReferenceCellByType(), DMPlexCreateBoxMesh() 3920 @*/ 3921 PetscErrorCode DMPlexCreateReferenceCell(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *refdm) 3922 { 3923 PetscErrorCode ierr; 3924 3925 PetscFunctionBeginUser; 3926 switch (dim) { 3927 case 0: ierr = DMPlexCreateReferenceCellByType(comm, DM_POLYTOPE_POINT, refdm);CHKERRQ(ierr);break; 3928 case 1: ierr = DMPlexCreateReferenceCellByType(comm, DM_POLYTOPE_SEGMENT, refdm);CHKERRQ(ierr);break; 3929 case 2: 3930 if (simplex) {ierr = DMPlexCreateReferenceCellByType(comm, DM_POLYTOPE_TRIANGLE, refdm);CHKERRQ(ierr);} 3931 else {ierr = DMPlexCreateReferenceCellByType(comm, DM_POLYTOPE_QUADRILATERAL, refdm);CHKERRQ(ierr);} 3932 break; 3933 case 3: 3934 if (simplex) {ierr = DMPlexCreateReferenceCellByType(comm, DM_POLYTOPE_TETRAHEDRON, refdm);CHKERRQ(ierr);} 3935 else {ierr = DMPlexCreateReferenceCellByType(comm, DM_POLYTOPE_HEXAHEDRON, refdm);CHKERRQ(ierr);} 3936 break; 3937 default: 3938 SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Cannot create reference cell for dimension %D", dim); 3939 } 3940 PetscFunctionReturn(0); 3941 } 3942