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