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