1 #define PETSCDM_DLL 2 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 3 #include <petscdmda.h> 4 #include <petscsf.h> 5 6 #undef __FUNCT__ 7 #define __FUNCT__ "DMPlexCreateDoublet" 8 /*@ 9 DMPlexCreateDoublet - Creates a mesh of two cells of the specified type, optionally with later refinement. 10 11 Collective on MPI_Comm 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 .keywords: DM, create 27 .seealso: DMSetType(), DMCreate() 28 @*/ 29 PetscErrorCode DMPlexCreateDoublet(MPI_Comm comm, PetscInt dim, PetscBool simplex, PetscBool interpolate, PetscBool refinementUniform, PetscReal refinementLimit, DM *newdm) 30 { 31 DM dm; 32 PetscInt p; 33 PetscMPIInt rank; 34 PetscErrorCode ierr; 35 36 PetscFunctionBegin; 37 ierr = DMCreate(comm, &dm);CHKERRQ(ierr); 38 ierr = DMSetType(dm, DMPLEX);CHKERRQ(ierr); 39 ierr = DMSetDimension(dm, dim);CHKERRQ(ierr); 40 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 41 switch (dim) { 42 case 2: 43 if (simplex) {ierr = PetscObjectSetName((PetscObject) dm, "triangular");CHKERRQ(ierr);} 44 else {ierr = PetscObjectSetName((PetscObject) dm, "quadrilateral");CHKERRQ(ierr);} 45 break; 46 case 3: 47 if (simplex) {ierr = PetscObjectSetName((PetscObject) dm, "tetrahedral");CHKERRQ(ierr);} 48 else {ierr = PetscObjectSetName((PetscObject) dm, "hexahedral");CHKERRQ(ierr);} 49 break; 50 default: 51 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %d", dim); 52 } 53 if (rank) { 54 PetscInt numPoints[2] = {0, 0}; 55 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 56 } else { 57 switch (dim) { 58 case 2: 59 if (simplex) { 60 PetscInt numPoints[2] = {4, 2}; 61 PetscInt coneSize[6] = {3, 3, 0, 0, 0, 0}; 62 PetscInt cones[6] = {2, 3, 4, 5, 4, 3}; 63 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 64 PetscScalar vertexCoords[8] = {-0.5, 0.5, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5}; 65 PetscInt markerPoints[8] = {2, 1, 3, 1, 4, 1, 5, 1}; 66 67 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 68 for (p = 0; p < 4; ++p) {ierr = DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 69 } else { 70 PetscInt numPoints[2] = {6, 2}; 71 PetscInt coneSize[8] = {4, 4, 0, 0, 0, 0, 0, 0}; 72 PetscInt cones[8] = {2, 3, 4, 5, 3, 6, 7, 4}; 73 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 74 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}; 75 76 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 77 } 78 break; 79 case 3: 80 if (simplex) { 81 PetscInt numPoints[2] = {5, 2}; 82 PetscInt coneSize[7] = {4, 4, 0, 0, 0, 0, 0}; 83 PetscInt cones[8] = {4, 3, 5, 2, 5, 3, 4, 6}; 84 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 85 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}; 86 PetscInt markerPoints[10] = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1}; 87 88 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 89 for (p = 0; p < 5; ++p) {ierr = DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);} 90 } else { 91 PetscInt numPoints[2] = {12, 2}; 92 PetscInt coneSize[14] = {8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 93 PetscInt cones[16] = {2, 3, 4, 5, 6, 7, 8, 9, 5, 4, 10, 11, 7, 12, 13, 8}; 94 PetscInt coneOrientations[16] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 95 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, 96 -1.0, -0.5, 0.5, 0.0, -0.5, 0.5, 0.0, 0.5, 0.5, -1.0, 0.5, 0.5, 97 1.0, 0.5, -0.5, 1.0, -0.5, -0.5, 1.0, -0.5, 0.5, 1.0, 0.5, 0.5}; 98 99 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 100 } 101 break; 102 default: 103 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %d", dim); 104 } 105 } 106 *newdm = dm; 107 if (refinementLimit > 0.0) { 108 DM rdm; 109 const char *name; 110 111 ierr = DMPlexSetRefinementUniform(*newdm, PETSC_FALSE);CHKERRQ(ierr); 112 ierr = DMPlexSetRefinementLimit(*newdm, refinementLimit);CHKERRQ(ierr); 113 ierr = DMRefine(*newdm, comm, &rdm);CHKERRQ(ierr); 114 ierr = PetscObjectGetName((PetscObject) *newdm, &name);CHKERRQ(ierr); 115 ierr = PetscObjectSetName((PetscObject) rdm, name);CHKERRQ(ierr); 116 ierr = DMDestroy(newdm);CHKERRQ(ierr); 117 *newdm = rdm; 118 } 119 if (interpolate) { 120 DM idm = NULL; 121 const char *name; 122 123 ierr = DMPlexInterpolate(*newdm, &idm);CHKERRQ(ierr); 124 ierr = PetscObjectGetName((PetscObject) *newdm, &name);CHKERRQ(ierr); 125 ierr = PetscObjectSetName((PetscObject) idm, name);CHKERRQ(ierr); 126 ierr = DMPlexCopyCoordinates(*newdm, idm);CHKERRQ(ierr); 127 ierr = DMCopyLabels(*newdm, idm);CHKERRQ(ierr); 128 ierr = DMDestroy(newdm);CHKERRQ(ierr); 129 *newdm = idm; 130 } 131 { 132 DM refinedMesh = NULL; 133 DM distributedMesh = NULL; 134 135 /* Distribute mesh over processes */ 136 ierr = DMPlexDistribute(*newdm, 0, NULL, &distributedMesh);CHKERRQ(ierr); 137 if (distributedMesh) { 138 ierr = DMDestroy(newdm);CHKERRQ(ierr); 139 *newdm = distributedMesh; 140 } 141 if (refinementUniform) { 142 ierr = DMPlexSetRefinementUniform(*newdm, refinementUniform);CHKERRQ(ierr); 143 ierr = DMRefine(*newdm, comm, &refinedMesh);CHKERRQ(ierr); 144 if (refinedMesh) { 145 ierr = DMDestroy(newdm);CHKERRQ(ierr); 146 *newdm = refinedMesh; 147 } 148 } 149 } 150 PetscFunctionReturn(0); 151 } 152 153 #undef __FUNCT__ 154 #define __FUNCT__ "DMPlexCreateSquareBoundary" 155 /*@ 156 DMPlexCreateSquareBoundary - Creates a 1D mesh the is the boundary of a square lattice. 157 158 Collective on MPI_Comm 159 160 Input Parameters: 161 + comm - The communicator for the DM object 162 . lower - The lower left corner coordinates 163 . upper - The upper right corner coordinates 164 - edges - The number of cells in each direction 165 166 Output Parameter: 167 . dm - The DM object 168 169 Note: Here is the numbering returned for 2 cells in each direction: 170 $ 18--5-17--4--16 171 $ | | | 172 $ 6 10 3 173 $ | | | 174 $ 19-11-20--9--15 175 $ | | | 176 $ 7 8 2 177 $ | | | 178 $ 12--0-13--1--14 179 180 Level: beginner 181 182 .keywords: DM, create 183 .seealso: DMPlexCreateBoxMesh(), DMPlexCreateCubeBoundary(), DMSetType(), DMCreate() 184 @*/ 185 PetscErrorCode DMPlexCreateSquareBoundary(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[]) 186 { 187 const PetscInt numVertices = (edges[0]+1)*(edges[1]+1); 188 const PetscInt numEdges = edges[0]*(edges[1]+1) + (edges[0]+1)*edges[1]; 189 PetscInt markerTop = 1; 190 PetscInt markerBottom = 1; 191 PetscInt markerRight = 1; 192 PetscInt markerLeft = 1; 193 PetscBool markerSeparate = PETSC_FALSE; 194 Vec coordinates; 195 PetscSection coordSection; 196 PetscScalar *coords; 197 PetscInt coordSize; 198 PetscMPIInt rank; 199 PetscInt v, vx, vy; 200 PetscErrorCode ierr; 201 202 PetscFunctionBegin; 203 ierr = PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);CHKERRQ(ierr); 204 if (markerSeparate) { 205 markerTop = 3; 206 markerBottom = 1; 207 markerRight = 2; 208 markerLeft = 4; 209 } 210 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 211 if (!rank) { 212 PetscInt e, ex, ey; 213 214 ierr = DMPlexSetChart(dm, 0, numEdges+numVertices);CHKERRQ(ierr); 215 for (e = 0; e < numEdges; ++e) { 216 ierr = DMPlexSetConeSize(dm, e, 2);CHKERRQ(ierr); 217 } 218 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 219 for (vx = 0; vx <= edges[0]; vx++) { 220 for (ey = 0; ey < edges[1]; ey++) { 221 PetscInt edge = vx*edges[1] + ey + edges[0]*(edges[1]+1); 222 PetscInt vertex = ey*(edges[0]+1) + vx + numEdges; 223 PetscInt cone[2]; 224 225 cone[0] = vertex; cone[1] = vertex+edges[0]+1; 226 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 227 if (vx == edges[0]) { 228 ierr = DMSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 229 ierr = DMSetLabelValue(dm, "marker", cone[0], markerRight);CHKERRQ(ierr); 230 if (ey == edges[1]-1) { 231 ierr = DMSetLabelValue(dm, "marker", cone[1], markerRight);CHKERRQ(ierr); 232 } 233 } else if (vx == 0) { 234 ierr = DMSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 235 ierr = DMSetLabelValue(dm, "marker", cone[0], markerLeft);CHKERRQ(ierr); 236 if (ey == edges[1]-1) { 237 ierr = DMSetLabelValue(dm, "marker", cone[1], markerLeft);CHKERRQ(ierr); 238 } 239 } 240 } 241 } 242 for (vy = 0; vy <= edges[1]; vy++) { 243 for (ex = 0; ex < edges[0]; ex++) { 244 PetscInt edge = vy*edges[0] + ex; 245 PetscInt vertex = vy*(edges[0]+1) + ex + numEdges; 246 PetscInt cone[2]; 247 248 cone[0] = vertex; cone[1] = vertex+1; 249 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 250 if (vy == edges[1]) { 251 ierr = DMSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 252 ierr = DMSetLabelValue(dm, "marker", cone[0], markerTop);CHKERRQ(ierr); 253 if (ex == edges[0]-1) { 254 ierr = DMSetLabelValue(dm, "marker", cone[1], markerTop);CHKERRQ(ierr); 255 } 256 } else if (vy == 0) { 257 ierr = DMSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 258 ierr = DMSetLabelValue(dm, "marker", cone[0], markerBottom);CHKERRQ(ierr); 259 if (ex == edges[0]-1) { 260 ierr = DMSetLabelValue(dm, "marker", cone[1], markerBottom);CHKERRQ(ierr); 261 } 262 } 263 } 264 } 265 } 266 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 267 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 268 /* Build coordinates */ 269 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 270 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 271 ierr = PetscSectionSetChart(coordSection, numEdges, numEdges + numVertices);CHKERRQ(ierr); 272 ierr = PetscSectionSetFieldComponents(coordSection, 0, 2);CHKERRQ(ierr); 273 for (v = numEdges; v < numEdges+numVertices; ++v) { 274 ierr = PetscSectionSetDof(coordSection, v, 2);CHKERRQ(ierr); 275 ierr = PetscSectionSetFieldDof(coordSection, v, 0, 2);CHKERRQ(ierr); 276 } 277 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 278 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 279 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 280 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 281 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 282 ierr = VecSetBlockSize(coordinates, 2);CHKERRQ(ierr); 283 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 284 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 285 for (vy = 0; vy <= edges[1]; ++vy) { 286 for (vx = 0; vx <= edges[0]; ++vx) { 287 coords[(vy*(edges[0]+1)+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/edges[0])*vx; 288 coords[(vy*(edges[0]+1)+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/edges[1])*vy; 289 } 290 } 291 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 292 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 293 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 294 PetscFunctionReturn(0); 295 } 296 297 #undef __FUNCT__ 298 #define __FUNCT__ "DMPlexCreateCubeBoundary" 299 /*@ 300 DMPlexCreateCubeBoundary - Creates a 2D mesh the is the boundary of a cubic lattice. 301 302 Collective on MPI_Comm 303 304 Input Parameters: 305 + comm - The communicator for the DM object 306 . lower - The lower left front corner coordinates 307 . upper - The upper right back corner coordinates 308 - edges - The number of cells in each direction 309 310 Output Parameter: 311 . dm - The DM object 312 313 Level: beginner 314 315 .keywords: DM, create 316 .seealso: DMPlexCreateBoxMesh(), DMPlexCreateSquareBoundary(), DMSetType(), DMCreate() 317 @*/ 318 PetscErrorCode DMPlexCreateCubeBoundary(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt faces[]) 319 { 320 PetscInt vertices[3], numVertices; 321 PetscInt numFaces = 2*faces[0]*faces[1] + 2*faces[1]*faces[2] + 2*faces[0]*faces[2]; 322 Vec coordinates; 323 PetscSection coordSection; 324 PetscScalar *coords; 325 PetscInt coordSize; 326 PetscMPIInt rank; 327 PetscInt v, vx, vy, vz; 328 PetscInt voffset, iface=0, cone[4]; 329 PetscErrorCode ierr; 330 331 PetscFunctionBegin; 332 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"); 333 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 334 vertices[0] = faces[0]+1; vertices[1] = faces[1]+1; vertices[2] = faces[2]+1; 335 numVertices = vertices[0]*vertices[1]*vertices[2]; 336 if (!rank) { 337 PetscInt f; 338 339 ierr = DMPlexSetChart(dm, 0, numFaces+numVertices);CHKERRQ(ierr); 340 for (f = 0; f < numFaces; ++f) { 341 ierr = DMPlexSetConeSize(dm, f, 4);CHKERRQ(ierr); 342 } 343 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 344 for (v = 0; v < numFaces+numVertices; ++v) { 345 ierr = DMSetLabelValue(dm, "marker", v, 1);CHKERRQ(ierr); 346 } 347 348 /* Side 0 (Top) */ 349 for (vy = 0; vy < faces[1]; vy++) { 350 for (vx = 0; vx < faces[0]; vx++) { 351 voffset = numFaces + vertices[0]*vertices[1]*(vertices[2]-1) + vy*vertices[0] + vx; 352 cone[0] = voffset; cone[1] = voffset+1; cone[2] = voffset+vertices[0]+1; cone[3] = voffset+vertices[0]; 353 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 354 iface++; 355 } 356 } 357 358 /* Side 1 (Bottom) */ 359 for (vy = 0; vy < faces[1]; vy++) { 360 for (vx = 0; vx < faces[0]; vx++) { 361 voffset = numFaces + vy*(faces[0]+1) + vx; 362 cone[0] = voffset+1; cone[1] = voffset; cone[2] = voffset+vertices[0]; cone[3] = voffset+vertices[0]+1; 363 ierr = DMPlexSetCone(dm, iface, cone);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 iface++; 375 } 376 } 377 378 /* Side 3 (Back) */ 379 for (vz = 0; vz < faces[2]; vz++) { 380 for (vx = 0; vx < faces[0]; vx++) { 381 voffset = numFaces + vz*vertices[0]*vertices[1] + vertices[0]*(vertices[1]-1) + vx; 382 cone[0] = voffset+vertices[0]*vertices[1]; cone[1] = voffset+vertices[0]*vertices[1]+1; 383 cone[2] = voffset+1; cone[3] = voffset; 384 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 385 iface++; 386 } 387 } 388 389 /* Side 4 (Left) */ 390 for (vz = 0; vz < faces[2]; vz++) { 391 for (vy = 0; vy < faces[1]; vy++) { 392 voffset = numFaces + vz*vertices[0]*vertices[1] + vy*vertices[0]; 393 cone[0] = voffset; cone[1] = voffset+vertices[0]*vertices[1]; 394 cone[2] = voffset+vertices[0]*vertices[1]+vertices[0]; cone[3] = voffset+vertices[0]; 395 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 396 iface++; 397 } 398 } 399 400 /* Side 5 (Right) */ 401 for (vz = 0; vz < faces[2]; vz++) { 402 for (vy = 0; vy < faces[1]; vy++) { 403 voffset = numFaces + vz*vertices[0]*vertices[1] + vy*vertices[0] + vx; 404 cone[0] = voffset+vertices[0]*vertices[1]; cone[1] = voffset; 405 cone[2] = voffset+vertices[0]; cone[3] = voffset+vertices[0]*vertices[1]+vertices[0]; 406 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 407 iface++; 408 } 409 } 410 } 411 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 412 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 413 /* Build coordinates */ 414 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 415 ierr = PetscSectionSetChart(coordSection, numFaces, numFaces + numVertices);CHKERRQ(ierr); 416 for (v = numFaces; v < numFaces+numVertices; ++v) { 417 ierr = PetscSectionSetDof(coordSection, v, 3);CHKERRQ(ierr); 418 } 419 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 420 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 421 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 422 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 423 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 424 ierr = VecSetBlockSize(coordinates, 3);CHKERRQ(ierr); 425 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 426 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 427 for (vz = 0; vz <= faces[2]; ++vz) { 428 for (vy = 0; vy <= faces[1]; ++vy) { 429 for (vx = 0; vx <= faces[0]; ++vx) { 430 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+0] = lower[0] + ((upper[0] - lower[0])/faces[0])*vx; 431 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+1] = lower[1] + ((upper[1] - lower[1])/faces[1])*vy; 432 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+2] = lower[2] + ((upper[2] - lower[2])/faces[2])*vz; 433 } 434 } 435 } 436 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 437 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 438 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 439 PetscFunctionReturn(0); 440 } 441 442 #undef __FUNCT__ 443 #define __FUNCT__ "DMPlexCreateCubeMesh_Internal" 444 static PetscErrorCode DMPlexCreateCubeMesh_Internal(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[], DMBoundaryType bdX, DMBoundaryType bdY, DMBoundaryType bdZ) 445 { 446 PetscInt markerTop = 1, faceMarkerTop = 1; 447 PetscInt markerBottom = 1, faceMarkerBottom = 1; 448 PetscInt markerFront = 1, faceMarkerFront = 1; 449 PetscInt markerBack = 1, faceMarkerBack = 1; 450 PetscInt markerRight = 1, faceMarkerRight = 1; 451 PetscInt markerLeft = 1, faceMarkerLeft = 1; 452 PetscInt dim; 453 PetscBool markerSeparate = PETSC_FALSE; 454 PetscMPIInt rank; 455 PetscErrorCode ierr; 456 457 PetscFunctionBegin; 458 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 459 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 460 ierr = DMCreateLabel(dm,"marker");CHKERRQ(ierr); 461 ierr = DMCreateLabel(dm,"Face Sets");CHKERRQ(ierr); 462 switch (dim) { 463 case 2: 464 faceMarkerTop = 3; 465 faceMarkerBottom = 1; 466 faceMarkerRight = 2; 467 faceMarkerLeft = 4; 468 break; 469 case 3: 470 faceMarkerBottom = 1; 471 faceMarkerTop = 2; 472 faceMarkerFront = 3; 473 faceMarkerBack = 4; 474 faceMarkerRight = 5; 475 faceMarkerLeft = 6; 476 break; 477 default: 478 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Dimension %d not supported",dim); 479 break; 480 } 481 ierr = PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);CHKERRQ(ierr); 482 if (markerSeparate) { 483 markerBottom = faceMarkerBottom; 484 markerTop = faceMarkerTop; 485 markerFront = faceMarkerFront; 486 markerBack = faceMarkerBack; 487 markerRight = faceMarkerRight; 488 markerLeft = faceMarkerLeft; 489 } 490 { 491 const PetscInt numXEdges = !rank ? edges[0] : 0; 492 const PetscInt numYEdges = !rank ? edges[1] : 0; 493 const PetscInt numZEdges = !rank ? edges[2] : 0; 494 const PetscInt numXVertices = !rank ? (bdX == DM_BOUNDARY_PERIODIC || bdX == DM_BOUNDARY_TWIST ? edges[0] : edges[0]+1) : 0; 495 const PetscInt numYVertices = !rank ? (bdY == DM_BOUNDARY_PERIODIC || bdY == DM_BOUNDARY_TWIST ? edges[1] : edges[1]+1) : 0; 496 const PetscInt numZVertices = !rank ? (bdZ == DM_BOUNDARY_PERIODIC || bdY == DM_BOUNDARY_TWIST ? edges[2] : edges[2]+1) : 0; 497 const PetscInt numCells = numXEdges*numYEdges*numZEdges; 498 const PetscInt numXFaces = numYEdges*numZEdges; 499 const PetscInt numYFaces = numXEdges*numZEdges; 500 const PetscInt numZFaces = numXEdges*numYEdges; 501 const PetscInt numTotXFaces = numXVertices*numXFaces; 502 const PetscInt numTotYFaces = numYVertices*numYFaces; 503 const PetscInt numTotZFaces = numZVertices*numZFaces; 504 const PetscInt numFaces = numTotXFaces + numTotYFaces + numTotZFaces; 505 const PetscInt numTotXEdges = numXEdges*numYVertices*numZVertices; 506 const PetscInt numTotYEdges = numYEdges*numXVertices*numZVertices; 507 const PetscInt numTotZEdges = numZEdges*numXVertices*numYVertices; 508 const PetscInt numVertices = numXVertices*numYVertices*numZVertices; 509 const PetscInt numEdges = numTotXEdges + numTotYEdges + numTotZEdges; 510 const PetscInt firstVertex = (dim == 2) ? numFaces : numCells; 511 const PetscInt firstXFace = (dim == 2) ? 0 : numCells + numVertices; 512 const PetscInt firstYFace = firstXFace + numTotXFaces; 513 const PetscInt firstZFace = firstYFace + numTotYFaces; 514 const PetscInt firstXEdge = numCells + numFaces + numVertices; 515 const PetscInt firstYEdge = firstXEdge + numTotXEdges; 516 const PetscInt firstZEdge = firstYEdge + numTotYEdges; 517 Vec coordinates; 518 PetscSection coordSection; 519 PetscScalar *coords; 520 PetscInt coordSize; 521 PetscInt v, vx, vy, vz; 522 PetscInt c, f, fx, fy, fz, e, ex, ey, ez; 523 524 ierr = DMPlexSetChart(dm, 0, numCells+numFaces+numEdges+numVertices);CHKERRQ(ierr); 525 for (c = 0; c < numCells; c++) { 526 ierr = DMPlexSetConeSize(dm, c, 6);CHKERRQ(ierr); 527 } 528 for (f = firstXFace; f < firstXFace+numFaces; ++f) { 529 ierr = DMPlexSetConeSize(dm, f, 4);CHKERRQ(ierr); 530 } 531 for (e = firstXEdge; e < firstXEdge+numEdges; ++e) { 532 ierr = DMPlexSetConeSize(dm, e, 2);CHKERRQ(ierr); 533 } 534 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 535 /* Build cells */ 536 for (fz = 0; fz < numZEdges; ++fz) { 537 for (fy = 0; fy < numYEdges; ++fy) { 538 for (fx = 0; fx < numXEdges; ++fx) { 539 PetscInt cell = (fz*numYEdges + fy)*numXEdges + fx; 540 PetscInt faceB = firstZFace + (fy*numXEdges+fx)*numZVertices + fz; 541 PetscInt faceT = firstZFace + (fy*numXEdges+fx)*numZVertices + ((fz+1)%numZVertices); 542 PetscInt faceF = firstYFace + (fz*numXEdges+fx)*numYVertices + fy; 543 PetscInt faceK = firstYFace + (fz*numXEdges+fx)*numYVertices + ((fy+1)%numYVertices); 544 PetscInt faceL = firstXFace + (fz*numYEdges+fy)*numXVertices + fx; 545 PetscInt faceR = firstXFace + (fz*numYEdges+fy)*numXVertices + ((fx+1)%numXVertices); 546 /* B, T, F, K, R, L */ 547 PetscInt ornt[6] = {-4, 0, 0, -1, 0, -4}; /* ??? */ 548 PetscInt cone[6]; 549 550 /* no boundary twisting in 3D */ 551 cone[0] = faceB; cone[1] = faceT; cone[2] = faceF; cone[3] = faceK; cone[4] = faceR; cone[5] = faceL; 552 ierr = DMPlexSetCone(dm, cell, cone);CHKERRQ(ierr); 553 ierr = DMPlexSetConeOrientation(dm, cell, ornt);CHKERRQ(ierr); 554 } 555 } 556 } 557 /* Build x faces */ 558 for (fz = 0; fz < numZEdges; ++fz) { 559 for (fy = 0; fy < numYEdges; ++fy) { 560 for (fx = 0; fx < numXVertices; ++fx) { 561 PetscInt face = firstXFace + (fz*numYEdges+fy)*numXVertices + fx; 562 PetscInt edgeL = firstZEdge + ( fy* numXVertices+fx)*numZEdges + fz; 563 PetscInt edgeR = firstZEdge + (((fy+1)%numYVertices)*numXVertices+fx)*numZEdges + fz; 564 PetscInt edgeB = firstYEdge + ( fz* numXVertices+fx)*numYEdges + fy; 565 PetscInt edgeT = firstYEdge + (((fz+1)%numZVertices)*numXVertices+fx)*numYEdges + fy; 566 PetscInt ornt[4] = {0, 0, -2, -2}; 567 PetscInt cone[4]; 568 569 if (dim == 3) { 570 /* markers */ 571 if (bdX != DM_BOUNDARY_PERIODIC) { 572 if (fx == numXVertices-1) { 573 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerRight);CHKERRQ(ierr); 574 ierr = DMSetLabelValue(dm, "marker", face, markerRight);CHKERRQ(ierr); 575 } 576 else if (fx == 0) { 577 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerLeft);CHKERRQ(ierr); 578 ierr = DMSetLabelValue(dm, "marker", face, markerLeft);CHKERRQ(ierr); 579 } 580 } 581 } 582 cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL; 583 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr); 584 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr); 585 } 586 } 587 } 588 /* Build y faces */ 589 for (fz = 0; fz < numZEdges; ++fz) { 590 for (fx = 0; fx < numYEdges; ++fx) { 591 for (fy = 0; fy < numYVertices; ++fy) { 592 PetscInt face = firstYFace + (fz*numXEdges+fx)*numYVertices + fy; 593 PetscInt edgeL = firstZEdge + (fy*numXVertices+ fx )*numZEdges + fz; 594 PetscInt edgeR = firstZEdge + (fy*numXVertices+((fx+1)%numXVertices))*numZEdges + fz; 595 PetscInt edgeB = firstXEdge + ( fz *numYVertices+fy)*numXEdges + fx; 596 PetscInt edgeT = firstXEdge + (((fz+1)%numZVertices)*numYVertices+fy)*numXEdges + fx; 597 PetscInt ornt[4] = {0, 0, -2, -2}; 598 PetscInt cone[4]; 599 600 if (dim == 3) { 601 /* markers */ 602 if (bdY != DM_BOUNDARY_PERIODIC) { 603 if (fy == numYVertices-1) { 604 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerBack);CHKERRQ(ierr); 605 ierr = DMSetLabelValue(dm, "marker", face, markerBack);CHKERRQ(ierr); 606 } 607 else if (fy == 0) { 608 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerFront);CHKERRQ(ierr); 609 ierr = DMSetLabelValue(dm, "marker", face, markerFront);CHKERRQ(ierr); 610 } 611 } 612 } 613 cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL; 614 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr); 615 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr); 616 } 617 } 618 } 619 /* Build z faces */ 620 for (fy = 0; fy < numYEdges; ++fy) { 621 for (fx = 0; fx < numXEdges; ++fx) { 622 for (fz = 0; fz < numZVertices; fz++) { 623 PetscInt face = firstZFace + (fy*numXEdges+fx)*numZVertices + fz; 624 PetscInt edgeL = firstYEdge + (fz*numXVertices+ fx )*numYEdges + fy; 625 PetscInt edgeR = firstYEdge + (fz*numXVertices+((fx+1)%numXVertices))*numYEdges + fy; 626 PetscInt edgeB = firstXEdge + (fz*numYVertices+ fy )*numXEdges + fx; 627 PetscInt edgeT = firstXEdge + (fz*numYVertices+((fy+1)%numYVertices))*numXEdges + fx; 628 PetscInt ornt[4] = {0, 0, -2, -2}; 629 PetscInt cone[4]; 630 631 if (dim == 2) { 632 if (bdX == DM_BOUNDARY_TWIST && fx == numXEdges-1) {edgeR += numYEdges-1-2*fy; ornt[1] = -2;} 633 if (bdY == DM_BOUNDARY_TWIST && fy == numYEdges-1) {edgeT += numXEdges-1-2*fx; ornt[2] = 0;} 634 } 635 else { 636 /* markers */ 637 if (bdZ != DM_BOUNDARY_PERIODIC) { 638 if (fz == numZVertices-1) { 639 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerTop);CHKERRQ(ierr); 640 ierr = DMSetLabelValue(dm, "marker", face, markerTop);CHKERRQ(ierr); 641 } 642 else if (fz == 0) { 643 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerBottom);CHKERRQ(ierr); 644 ierr = DMSetLabelValue(dm, "marker", face, markerBottom);CHKERRQ(ierr); 645 } 646 } 647 } 648 cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL; 649 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr); 650 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr); 651 } 652 } 653 } 654 /* Build Z edges*/ 655 for (vy = 0; vy < numYVertices; vy++) { 656 for (vx = 0; vx < numXVertices; vx++) { 657 for (ez = 0; ez < numZEdges; ez++) { 658 const PetscInt edge = firstZEdge + (vy*numXVertices+vx)*numZEdges + ez; 659 const PetscInt vertexB = firstVertex + ( ez *numYVertices+vy)*numXVertices + vx; 660 const PetscInt vertexT = firstVertex + (((ez+1)%numZVertices)*numYVertices+vy)*numXVertices + vx; 661 PetscInt cone[2]; 662 663 if (dim == 3) { 664 if (bdX != DM_BOUNDARY_PERIODIC) { 665 if (vx == numXVertices-1) { 666 ierr = DMSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 667 } 668 else if (vx == 0) { 669 ierr = DMSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 670 } 671 } 672 if (bdY != DM_BOUNDARY_PERIODIC) { 673 if (vy == numYVertices-1) { 674 ierr = DMSetLabelValue(dm, "marker", edge, markerBack);CHKERRQ(ierr); 675 } 676 else if (vy == 0) { 677 ierr = DMSetLabelValue(dm, "marker", edge, markerFront);CHKERRQ(ierr); 678 } 679 } 680 } 681 cone[0] = vertexB; cone[1] = vertexT; 682 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 683 } 684 } 685 } 686 /* Build Y edges*/ 687 for (vz = 0; vz < numZVertices; vz++) { 688 for (vx = 0; vx < numXVertices; vx++) { 689 for (ey = 0; ey < numYEdges; ey++) { 690 const PetscInt nextv = (dim == 2 && bdY == DM_BOUNDARY_TWIST && ey == numYEdges-1) ? (numXVertices-vx-1) : (vz*numYVertices+((ey+1)%numYVertices))*numXVertices + vx; 691 const PetscInt edge = firstYEdge + (vz*numXVertices+vx)*numYEdges + ey; 692 const PetscInt vertexF = firstVertex + (vz*numYVertices+ey)*numXVertices + vx; 693 const PetscInt vertexK = firstVertex + nextv; 694 PetscInt cone[2]; 695 696 cone[0] = vertexF; cone[1] = vertexK; 697 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 698 if (dim == 2) { 699 if ((bdX != DM_BOUNDARY_PERIODIC) && (bdX != DM_BOUNDARY_TWIST)) { 700 if (vx == numXVertices-1) { 701 ierr = DMSetLabelValue(dm, "Face Sets", edge, faceMarkerRight);CHKERRQ(ierr); 702 ierr = DMSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 703 ierr = DMSetLabelValue(dm, "marker", cone[0], markerRight);CHKERRQ(ierr); 704 if (ey == numYEdges-1) { 705 ierr = DMSetLabelValue(dm, "marker", cone[1], markerRight);CHKERRQ(ierr); 706 } 707 } 708 else if (vx == 0) { 709 ierr = DMSetLabelValue(dm, "Face Sets", edge, faceMarkerLeft);CHKERRQ(ierr); 710 ierr = DMSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 711 ierr = DMSetLabelValue(dm, "marker", cone[0], markerLeft);CHKERRQ(ierr); 712 if (ey == numYEdges-1) { 713 ierr = DMSetLabelValue(dm, "marker", cone[1], markerLeft);CHKERRQ(ierr); 714 } 715 } 716 } 717 } 718 else { 719 if (bdX != DM_BOUNDARY_PERIODIC) { 720 if (vx == numXVertices-1) { 721 ierr = DMSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 722 } 723 else if (vx == 0) { 724 ierr = DMSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 725 } 726 } 727 if (bdZ != DM_BOUNDARY_PERIODIC) { 728 if (vz == numZVertices-1) { 729 ierr = DMSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 730 } 731 else if (vz == 0) { 732 ierr = DMSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 733 } 734 } 735 } 736 } 737 } 738 } 739 /* Build X edges*/ 740 for (vz = 0; vz < numZVertices; vz++) { 741 for (vy = 0; vy < numYVertices; vy++) { 742 for (ex = 0; ex < numXEdges; ex++) { 743 const PetscInt nextv = (dim == 2 && bdX == DM_BOUNDARY_TWIST && ex == numXEdges-1) ? (numYVertices-vy-1)*numXVertices : (vz*numYVertices+vy)*numXVertices + (ex+1)%numXVertices; 744 const PetscInt edge = firstXEdge + (vz*numYVertices+vy)*numXEdges + ex; 745 const PetscInt vertexL = firstVertex + (vz*numYVertices+vy)*numXVertices + ex; 746 const PetscInt vertexR = firstVertex + nextv; 747 PetscInt cone[2]; 748 749 cone[0] = vertexL; cone[1] = vertexR; 750 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 751 if (dim == 2) { 752 if ((bdY != DM_BOUNDARY_PERIODIC) && (bdY != DM_BOUNDARY_TWIST)) { 753 if (vy == numYVertices-1) { 754 ierr = DMSetLabelValue(dm, "Face Sets", edge, faceMarkerTop);CHKERRQ(ierr); 755 ierr = DMSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 756 ierr = DMSetLabelValue(dm, "marker", cone[0], markerTop);CHKERRQ(ierr); 757 if (ex == numXEdges-1) { 758 ierr = DMSetLabelValue(dm, "marker", cone[1], markerTop);CHKERRQ(ierr); 759 } 760 } 761 else if (vy == 0) { 762 ierr = DMSetLabelValue(dm, "Face Sets", edge, faceMarkerBottom);CHKERRQ(ierr); 763 ierr = DMSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 764 ierr = DMSetLabelValue(dm, "marker", cone[0], markerBottom);CHKERRQ(ierr); 765 if (ex == numXEdges-1) { 766 ierr = DMSetLabelValue(dm, "marker", cone[1], markerBottom);CHKERRQ(ierr); 767 } 768 } 769 } 770 } 771 else { 772 if (bdY != DM_BOUNDARY_PERIODIC) { 773 if (vy == numYVertices-1) { 774 ierr = DMSetLabelValue(dm, "marker", edge, markerBack);CHKERRQ(ierr); 775 } 776 else if (vy == 0) { 777 ierr = DMSetLabelValue(dm, "marker", edge, markerFront);CHKERRQ(ierr); 778 } 779 } 780 if (bdZ != DM_BOUNDARY_PERIODIC) { 781 if (vz == numZVertices-1) { 782 ierr = DMSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 783 } 784 else if (vz == 0) { 785 ierr = DMSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 786 } 787 } 788 } 789 } 790 } 791 } 792 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 793 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 794 /* Build coordinates */ 795 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 796 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 797 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 798 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVertices);CHKERRQ(ierr); 799 for (v = firstVertex; v < firstVertex+numVertices; ++v) { 800 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 801 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 802 } 803 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 804 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 805 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 806 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 807 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 808 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 809 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 810 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 811 for (vz = 0; vz < numZVertices; ++vz) { 812 for (vy = 0; vy < numYVertices; ++vy) { 813 for (vx = 0; vx < numXVertices; ++vx) { 814 coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+0] = lower[0] + ((upper[0] - lower[0])/numXEdges)*vx; 815 coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+1] = lower[1] + ((upper[1] - lower[1])/numYEdges)*vy; 816 if (dim == 3) { 817 coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+2] = lower[2] + ((upper[2] - lower[2])/numZEdges)*vz; 818 } 819 } 820 } 821 } 822 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 823 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 824 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 825 } 826 PetscFunctionReturn(0); 827 } 828 829 #undef __FUNCT__ 830 #define __FUNCT__ "DMPlexCreateSquareMesh" 831 /*@ 832 DMPlexCreateSquareMesh - Creates a 2D mesh for a square lattice. 833 834 Collective on MPI_Comm 835 836 Input Parameters: 837 + comm - The communicator for the DM object 838 . lower - The lower left corner coordinates 839 . upper - The upper right corner coordinates 840 . edges - The number of cells in each direction 841 . bdX - The boundary type for the X direction 842 - bdY - The boundary type for the Y direction 843 844 Output Parameter: 845 . dm - The DM object 846 847 Note: Here is the numbering returned for 2 cells in each direction: 848 $ 22--8-23--9--24 849 $ | | | 850 $ 13 2 14 3 15 851 $ | | | 852 $ 19--6-20--7--21 853 $ | | | 854 $ 10 0 11 1 12 855 $ | | | 856 $ 16--4-17--5--18 857 858 Level: beginner 859 860 .keywords: DM, create 861 .seealso: DMPlexCreateBoxMesh(), DMPlexCreateSquareBoundary(), DMPlexCreateCubeBoundary(), DMSetType(), DMCreate() 862 @*/ 863 PetscErrorCode DMPlexCreateSquareMesh(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[], DMBoundaryType bdX, DMBoundaryType bdY) 864 { 865 PetscReal lower3[3], upper3[3]; 866 PetscInt edges3[3]; 867 PetscErrorCode ierr; 868 869 PetscFunctionBegin; 870 lower3[0] = lower[0]; lower3[1] = lower[1]; lower3[2] = 0.; 871 upper3[0] = upper[0]; upper3[1] = upper[1]; upper3[2] = 0.; 872 edges3[0] = edges[0]; edges3[1] = edges[1]; edges3[2] = 0; 873 ierr = DMPlexCreateCubeMesh_Internal(dm, lower3, upper3, edges3, bdX, bdY, DM_BOUNDARY_NONE);CHKERRQ(ierr); 874 PetscFunctionReturn(0); 875 } 876 877 #undef __FUNCT__ 878 #define __FUNCT__ "DMPlexCreateBoxMesh" 879 /*@ 880 DMPlexCreateBoxMesh - Creates a mesh on the tensor product of unit intervals (box) using simplices. 881 882 Collective on MPI_Comm 883 884 Input Parameters: 885 + comm - The communicator for the DM object 886 . dim - The spatial dimension 887 . numFaces - Number of faces per dimension 888 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 889 890 Output Parameter: 891 . dm - The DM object 892 893 Level: beginner 894 895 .keywords: DM, create 896 .seealso: DMPlexCreateHexBoxMesh(), DMSetType(), DMCreate() 897 @*/ 898 PetscErrorCode DMPlexCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscInt numFaces, PetscBool interpolate, DM *dm) 899 { 900 DM boundary; 901 PetscErrorCode ierr; 902 903 PetscFunctionBegin; 904 PetscValidPointer(dm, 4); 905 ierr = DMCreate(comm, &boundary);CHKERRQ(ierr); 906 PetscValidLogicalCollectiveInt(boundary,dim,2); 907 ierr = DMSetType(boundary, DMPLEX);CHKERRQ(ierr); 908 ierr = DMSetDimension(boundary, dim-1);CHKERRQ(ierr); 909 ierr = DMSetCoordinateDim(boundary, dim);CHKERRQ(ierr); 910 switch (dim) { 911 case 2: 912 { 913 PetscReal lower[2] = {0.0, 0.0}; 914 PetscReal upper[2] = {1.0, 1.0}; 915 PetscInt edges[2]; 916 917 edges[0] = numFaces; edges[1] = numFaces; 918 ierr = DMPlexCreateSquareBoundary(boundary, lower, upper, edges);CHKERRQ(ierr); 919 break; 920 } 921 case 3: 922 { 923 PetscReal lower[3] = {0.0, 0.0, 0.0}; 924 PetscReal upper[3] = {1.0, 1.0, 1.0}; 925 PetscInt faces[3]; 926 927 faces[0] = numFaces; faces[1] = numFaces; faces[2] = numFaces; 928 ierr = DMPlexCreateCubeBoundary(boundary, lower, upper, faces);CHKERRQ(ierr); 929 break; 930 } 931 default: 932 SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim); 933 } 934 ierr = DMPlexGenerate(boundary, NULL, interpolate, dm);CHKERRQ(ierr); 935 ierr = DMDestroy(&boundary);CHKERRQ(ierr); 936 PetscFunctionReturn(0); 937 } 938 939 #undef __FUNCT__ 940 #define __FUNCT__ "DMPlexCreateHexBoxMesh" 941 /*@ 942 DMPlexCreateHexBoxMesh - Creates a mesh on the tensor product of unit intervals (box) using hexahedra. 943 944 Collective on MPI_Comm 945 946 Input Parameters: 947 + comm - The communicator for the DM object 948 . dim - The spatial dimension 949 . periodicX - The boundary type for the X direction 950 . periodicY - The boundary type for the Y direction 951 . periodicZ - The boundary type for the Z direction 952 - cells - The number of cells in each direction 953 954 Output Parameter: 955 . dm - The DM object 956 957 Level: beginner 958 959 .keywords: DM, create 960 .seealso: DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 961 @*/ 962 PetscErrorCode DMPlexCreateHexBoxMesh(MPI_Comm comm, PetscInt dim, const PetscInt cells[], DMBoundaryType periodicX, DMBoundaryType periodicY, DMBoundaryType periodicZ, DM *dm) 963 { 964 PetscErrorCode ierr; 965 966 PetscFunctionBegin; 967 PetscValidPointer(dm, 7); 968 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 969 PetscValidLogicalCollectiveInt(*dm,dim,2); 970 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 971 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 972 switch (dim) { 973 case 2: 974 { 975 PetscReal lower[2] = {0.0, 0.0}; 976 PetscReal upper[2] = {1.0, 1.0}; 977 978 ierr = DMPlexCreateSquareMesh(*dm, lower, upper, cells, periodicX, periodicY);CHKERRQ(ierr); 979 break; 980 } 981 case 3: 982 { 983 PetscReal lower[3] = {0.0, 0.0, 0.0}; 984 PetscReal upper[3] = {1.0, 1.0, 1.0}; 985 986 ierr = DMPlexCreateCubeMesh_Internal(*dm, lower, upper, cells, periodicX, periodicY, periodicZ);CHKERRQ(ierr); 987 break; 988 } 989 default: 990 SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim); 991 } 992 PetscFunctionReturn(0); 993 } 994 995 #undef __FUNCT__ 996 #define __FUNCT__ "DMPlexCreateQuadSphereMesh" 997 /*@ 998 DMPlexCreateQuadSphereMesh - Creates a mesh on the sphere using quadrilaterals. 999 1000 Collective on MPI_Comm 1001 1002 Input Parameters: 1003 . comm - The communicator for the DM object 1004 1005 Output Parameter: 1006 . dm - The DM object 1007 1008 Level: beginner 1009 1010 .keywords: DM, create 1011 .seealso: DMPlexCreateBoxMesh(), DMPlexCreateHexBoxMesh(), DMSetType(), DMCreate() 1012 @*/ 1013 PetscErrorCode DMPlexCreateQuadSphereMesh(MPI_Comm comm, DM *dm) 1014 { 1015 /* 1016 12-21--13 1017 | | 1018 25 4 24 1019 | | 1020 12-25--9-16--8-24--13 1021 | | | | 1022 23 5 17 0 15 3 22 1023 | | | | 1024 10-20--6-14--7-19--11 1025 | | 1026 20 1 19 1027 | | 1028 10-18--11 1029 | | 1030 23 2 22 1031 | | 1032 12-21--13 1033 */ 1034 const PetscInt dim = 2, embedDim = 3; 1035 const PetscReal d = 1.0/PetscSqrtReal(3.0); 1036 PetscMPIInt rank; 1037 PetscErrorCode ierr; 1038 1039 PetscFunctionBegin; 1040 PetscValidPointer(dm, 3); 1041 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1042 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1043 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1044 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) *dm), &rank);CHKERRQ(ierr); 1045 { 1046 PetscSection coordSection; 1047 Vec coordinates; 1048 PetscScalar *coords; 1049 PetscInt coordSize; 1050 const PetscInt numCells = !rank ? 6 : 0; 1051 const PetscInt numEdges = !rank ? 12 : 0; 1052 const PetscInt numVerts = !rank ? 8 : 0; 1053 const PetscInt firstVertex = numCells; 1054 const PetscInt firstEdge = numCells + numVerts; 1055 PetscInt cone[4], ornt[4]; 1056 PetscInt c, e, v; 1057 1058 /* Build Topology */ 1059 ierr = DMPlexSetChart(*dm, 0, numCells+numEdges+numVerts);CHKERRQ(ierr); 1060 for (c = 0; c < numCells; c++) { 1061 ierr = DMPlexSetConeSize(*dm, c, 4);CHKERRQ(ierr); 1062 } 1063 for (e = firstEdge; e < firstEdge+numEdges; ++e) { 1064 ierr = DMPlexSetConeSize(*dm, e, 2);CHKERRQ(ierr); 1065 } 1066 ierr = DMSetUp(*dm);CHKERRQ(ierr); /* Allocate space for cones */ 1067 /* Cell 0 */ 1068 cone[0] = 14; cone[1] = 15; cone[2] = 16; cone[3] = 17; 1069 ierr = DMPlexSetCone(*dm, 0, cone);CHKERRQ(ierr); 1070 ornt[0] = 0; ornt[1] = 0; ornt[2] = 0; ornt[3] = 0; 1071 ierr = DMPlexSetConeOrientation(*dm, 0, ornt);CHKERRQ(ierr); 1072 /* Cell 1 */ 1073 cone[0] = 18; cone[1] = 19; cone[2] = 14; cone[3] = 20; 1074 ierr = DMPlexSetCone(*dm, 1, cone);CHKERRQ(ierr); 1075 ornt[0] = 0; ornt[1] = 0; ornt[2] = -2; ornt[3] = 0; 1076 ierr = DMPlexSetConeOrientation(*dm, 1, ornt);CHKERRQ(ierr); 1077 /* Cell 2 */ 1078 cone[0] = 21; cone[1] = 22; cone[2] = 18; cone[3] = 23; 1079 ierr = DMPlexSetCone(*dm, 2, cone);CHKERRQ(ierr); 1080 ornt[0] = 0; ornt[1] = 0; ornt[2] = -2; ornt[3] = 0; 1081 ierr = DMPlexSetConeOrientation(*dm, 2, ornt);CHKERRQ(ierr); 1082 /* Cell 3 */ 1083 cone[0] = 19; cone[1] = 22; cone[2] = 24; cone[3] = 15; 1084 ierr = DMPlexSetCone(*dm, 3, cone);CHKERRQ(ierr); 1085 ornt[0] = -2; ornt[1] = -2; ornt[2] = 0; ornt[3] = -2; 1086 ierr = DMPlexSetConeOrientation(*dm, 3, ornt);CHKERRQ(ierr); 1087 /* Cell 4 */ 1088 cone[0] = 16; cone[1] = 24; cone[2] = 21; cone[3] = 25; 1089 ierr = DMPlexSetCone(*dm, 4, cone);CHKERRQ(ierr); 1090 ornt[0] = -2; ornt[1] = -2; ornt[2] = -2; ornt[3] = 0; 1091 ierr = DMPlexSetConeOrientation(*dm, 4, ornt);CHKERRQ(ierr); 1092 /* Cell 5 */ 1093 cone[0] = 20; cone[1] = 17; cone[2] = 25; cone[3] = 23; 1094 ierr = DMPlexSetCone(*dm, 5, cone);CHKERRQ(ierr); 1095 ornt[0] = -2; ornt[1] = -2; ornt[2] = -2; ornt[3] = -2; 1096 ierr = DMPlexSetConeOrientation(*dm, 5, ornt);CHKERRQ(ierr); 1097 /* Edges */ 1098 cone[0] = 6; cone[1] = 7; 1099 ierr = DMPlexSetCone(*dm, 14, cone);CHKERRQ(ierr); 1100 cone[0] = 7; cone[1] = 8; 1101 ierr = DMPlexSetCone(*dm, 15, cone);CHKERRQ(ierr); 1102 cone[0] = 8; cone[1] = 9; 1103 ierr = DMPlexSetCone(*dm, 16, cone);CHKERRQ(ierr); 1104 cone[0] = 9; cone[1] = 6; 1105 ierr = DMPlexSetCone(*dm, 17, cone);CHKERRQ(ierr); 1106 cone[0] = 10; cone[1] = 11; 1107 ierr = DMPlexSetCone(*dm, 18, cone);CHKERRQ(ierr); 1108 cone[0] = 11; cone[1] = 7; 1109 ierr = DMPlexSetCone(*dm, 19, cone);CHKERRQ(ierr); 1110 cone[0] = 6; cone[1] = 10; 1111 ierr = DMPlexSetCone(*dm, 20, cone);CHKERRQ(ierr); 1112 cone[0] = 12; cone[1] = 13; 1113 ierr = DMPlexSetCone(*dm, 21, cone);CHKERRQ(ierr); 1114 cone[0] = 13; cone[1] = 11; 1115 ierr = DMPlexSetCone(*dm, 22, cone);CHKERRQ(ierr); 1116 cone[0] = 10; cone[1] = 12; 1117 ierr = DMPlexSetCone(*dm, 23, cone);CHKERRQ(ierr); 1118 cone[0] = 13; cone[1] = 8; 1119 ierr = DMPlexSetCone(*dm, 24, cone);CHKERRQ(ierr); 1120 cone[0] = 12; cone[1] = 9; 1121 ierr = DMPlexSetCone(*dm, 25, cone);CHKERRQ(ierr); 1122 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1123 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1124 /* Build coordinates */ 1125 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 1126 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1127 ierr = PetscSectionSetFieldComponents(coordSection, 0, embedDim);CHKERRQ(ierr); 1128 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVerts);CHKERRQ(ierr); 1129 for (v = firstVertex; v < firstVertex+numVerts; ++v) { 1130 ierr = PetscSectionSetDof(coordSection, v, embedDim);CHKERRQ(ierr); 1131 ierr = PetscSectionSetFieldDof(coordSection, v, 0, embedDim);CHKERRQ(ierr); 1132 } 1133 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1134 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1135 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 1136 ierr = VecSetBlockSize(coordinates, embedDim);CHKERRQ(ierr); 1137 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1138 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1139 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1140 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1141 coords[0*embedDim+0] = -d; coords[0*embedDim+1] = d; coords[0*embedDim+2] = -d; 1142 coords[1*embedDim+0] = d; coords[1*embedDim+1] = d; coords[1*embedDim+2] = -d; 1143 coords[2*embedDim+0] = d; coords[2*embedDim+1] = -d; coords[2*embedDim+2] = -d; 1144 coords[3*embedDim+0] = -d; coords[3*embedDim+1] = -d; coords[3*embedDim+2] = -d; 1145 coords[4*embedDim+0] = -d; coords[4*embedDim+1] = d; coords[4*embedDim+2] = d; 1146 coords[5*embedDim+0] = d; coords[5*embedDim+1] = d; coords[5*embedDim+2] = d; 1147 coords[6*embedDim+0] = -d; coords[6*embedDim+1] = -d; coords[6*embedDim+2] = d; 1148 coords[7*embedDim+0] = d; coords[7*embedDim+1] = -d; coords[7*embedDim+2] = d; 1149 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1150 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 1151 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1152 } 1153 PetscFunctionReturn(0); 1154 } 1155 1156 /* External function declarations here */ 1157 extern PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling); 1158 extern PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, Mat *mat); 1159 extern PetscErrorCode DMCreateDefaultSection_Plex(DM dm); 1160 extern PetscErrorCode DMCreateDefaultConstraints_Plex(DM dm); 1161 extern PetscErrorCode DMCreateMatrix_Plex(DM dm, Mat *J); 1162 extern PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm); 1163 extern PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *dmRefined); 1164 extern PetscErrorCode DMCoarsen_Plex(DM dm, MPI_Comm comm, DM *dmCoarsened); 1165 extern PetscErrorCode DMAdaptLabel_Plex(DM dm, DMLabel adaptLabel, DM *dmRefined); 1166 extern PetscErrorCode DMRefineHierarchy_Plex(DM dm, PetscInt nlevels, DM dmRefined[]); 1167 extern PetscErrorCode DMCoarsenHierarchy_Plex(DM dm, PetscInt nlevels, DM dmCoarsened[]); 1168 extern PetscErrorCode DMClone_Plex(DM dm, DM *newdm); 1169 extern PetscErrorCode DMSetUp_Plex(DM dm); 1170 extern PetscErrorCode DMDestroy_Plex(DM dm); 1171 extern PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer); 1172 extern PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer); 1173 extern PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm); 1174 extern PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, DMPointLocationType ltype, PetscSF cellSF); 1175 extern PetscErrorCode DMProjectFunctionLocal_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,InsertMode,Vec); 1176 extern PetscErrorCode DMProjectFunctionLabelLocal_Plex(DM,PetscReal,DMLabel,PetscInt,const PetscInt[],PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,InsertMode,Vec); 1177 extern PetscErrorCode DMProjectFieldLocal_Plex(DM,Vec,void (**)(PetscInt,PetscInt,PetscInt,const PetscInt[],const PetscInt[],const PetscScalar[],const PetscScalar[],const PetscScalar[],const PetscInt[],const PetscInt[],const PetscScalar[],const PetscScalar[],const PetscScalar[],PetscReal,const PetscReal[],PetscScalar[]),InsertMode,Vec); 1178 extern PetscErrorCode DMComputeL2Diff_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,Vec,PetscReal *); 1179 extern PetscErrorCode DMComputeL2GradientDiff_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[], const PetscReal[],PetscInt,PetscScalar *,void *),void **,Vec,const PetscReal [],PetscReal *); 1180 extern PetscErrorCode DMComputeL2FieldDiff_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,Vec,PetscReal *); 1181 1182 #undef __FUNCT__ 1183 #define __FUNCT__ "DMPlexReplace_Static" 1184 /* Replace dm with the contents of dmNew 1185 - Share the DM_Plex structure 1186 - Share the coordinates 1187 - Share the SF 1188 */ 1189 static PetscErrorCode DMPlexReplace_Static(DM dm, DM dmNew) 1190 { 1191 PetscSF sf; 1192 DM coordDM, coarseDM; 1193 Vec coords; 1194 const PetscReal *maxCell, *L; 1195 const DMBoundaryType *bd; 1196 PetscErrorCode ierr; 1197 1198 PetscFunctionBegin; 1199 ierr = DMGetPointSF(dmNew, &sf);CHKERRQ(ierr); 1200 ierr = DMSetPointSF(dm, sf);CHKERRQ(ierr); 1201 ierr = DMGetCoordinateDM(dmNew, &coordDM);CHKERRQ(ierr); 1202 ierr = DMGetCoordinatesLocal(dmNew, &coords);CHKERRQ(ierr); 1203 ierr = DMSetCoordinateDM(dm, coordDM);CHKERRQ(ierr); 1204 ierr = DMSetCoordinatesLocal(dm, coords);CHKERRQ(ierr); 1205 ierr = DMGetPeriodicity(dm, &maxCell, &L, &bd);CHKERRQ(ierr); 1206 if (L) {ierr = DMSetPeriodicity(dmNew, maxCell, L, bd);CHKERRQ(ierr);} 1207 ierr = DMDestroy_Plex(dm);CHKERRQ(ierr); 1208 dm->data = dmNew->data; 1209 ((DM_Plex *) dmNew->data)->refct++; 1210 dmNew->labels->refct++; 1211 if (!--(dm->labels->refct)) { 1212 DMLabelLink next = dm->labels->next; 1213 1214 /* destroy the labels */ 1215 while (next) { 1216 DMLabelLink tmp = next->next; 1217 1218 ierr = DMLabelDestroy(&next->label);CHKERRQ(ierr); 1219 ierr = PetscFree(next);CHKERRQ(ierr); 1220 next = tmp; 1221 } 1222 ierr = PetscFree(dm->labels);CHKERRQ(ierr); 1223 } 1224 dm->labels = dmNew->labels; 1225 dm->depthLabel = dmNew->depthLabel; 1226 ierr = DMGetCoarseDM(dmNew,&coarseDM);CHKERRQ(ierr); 1227 ierr = DMSetCoarseDM(dm,coarseDM);CHKERRQ(ierr); 1228 PetscFunctionReturn(0); 1229 } 1230 1231 #undef __FUNCT__ 1232 #define __FUNCT__ "DMPlexSwap_Static" 1233 /* Swap dm with the contents of dmNew 1234 - Swap the DM_Plex structure 1235 - Swap the coordinates 1236 - Swap the point PetscSF 1237 */ 1238 static PetscErrorCode DMPlexSwap_Static(DM dmA, DM dmB) 1239 { 1240 DM coordDMA, coordDMB; 1241 Vec coordsA, coordsB; 1242 PetscSF sfA, sfB; 1243 void *tmp; 1244 DMLabelLinkList listTmp; 1245 DMLabel depthTmp; 1246 PetscErrorCode ierr; 1247 1248 PetscFunctionBegin; 1249 ierr = DMGetPointSF(dmA, &sfA);CHKERRQ(ierr); 1250 ierr = DMGetPointSF(dmB, &sfB);CHKERRQ(ierr); 1251 ierr = PetscObjectReference((PetscObject) sfA);CHKERRQ(ierr); 1252 ierr = DMSetPointSF(dmA, sfB);CHKERRQ(ierr); 1253 ierr = DMSetPointSF(dmB, sfA);CHKERRQ(ierr); 1254 ierr = PetscObjectDereference((PetscObject) sfA);CHKERRQ(ierr); 1255 1256 ierr = DMGetCoordinateDM(dmA, &coordDMA);CHKERRQ(ierr); 1257 ierr = DMGetCoordinateDM(dmB, &coordDMB);CHKERRQ(ierr); 1258 ierr = PetscObjectReference((PetscObject) coordDMA);CHKERRQ(ierr); 1259 ierr = DMSetCoordinateDM(dmA, coordDMB);CHKERRQ(ierr); 1260 ierr = DMSetCoordinateDM(dmB, coordDMA);CHKERRQ(ierr); 1261 ierr = PetscObjectDereference((PetscObject) coordDMA);CHKERRQ(ierr); 1262 1263 ierr = DMGetCoordinatesLocal(dmA, &coordsA);CHKERRQ(ierr); 1264 ierr = DMGetCoordinatesLocal(dmB, &coordsB);CHKERRQ(ierr); 1265 ierr = PetscObjectReference((PetscObject) coordsA);CHKERRQ(ierr); 1266 ierr = DMSetCoordinatesLocal(dmA, coordsB);CHKERRQ(ierr); 1267 ierr = DMSetCoordinatesLocal(dmB, coordsA);CHKERRQ(ierr); 1268 ierr = PetscObjectDereference((PetscObject) coordsA);CHKERRQ(ierr); 1269 1270 tmp = dmA->data; 1271 dmA->data = dmB->data; 1272 dmB->data = tmp; 1273 listTmp = dmA->labels; 1274 dmA->labels = dmB->labels; 1275 dmB->labels = listTmp; 1276 depthTmp = dmA->depthLabel; 1277 dmA->depthLabel = dmB->depthLabel; 1278 dmB->depthLabel = depthTmp; 1279 PetscFunctionReturn(0); 1280 } 1281 1282 #undef __FUNCT__ 1283 #define __FUNCT__ "DMSetFromOptions_NonRefinement_Plex" 1284 PetscErrorCode DMSetFromOptions_NonRefinement_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 1285 { 1286 DM_Plex *mesh = (DM_Plex*) dm->data; 1287 PetscErrorCode ierr; 1288 1289 PetscFunctionBegin; 1290 /* Handle viewing */ 1291 ierr = PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMView", PETSC_FALSE, &mesh->printSetValues, NULL);CHKERRQ(ierr); 1292 ierr = PetscOptionsInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMView", 0, &mesh->printFEM, NULL);CHKERRQ(ierr); 1293 ierr = PetscOptionsReal("-dm_plex_print_tol", "Tolerance for FEM output", "DMView", mesh->printTol, &mesh->printTol, NULL);CHKERRQ(ierr); 1294 /* Point Location */ 1295 ierr = PetscOptionsBool("-dm_plex_hash_location", "Use grid hashing for point location", "DMView", PETSC_FALSE, &mesh->useHashLocation, NULL);CHKERRQ(ierr); 1296 /* Generation and remeshing */ 1297 ierr = PetscOptionsBool("-dm_plex_remesh_bd", "Allow changes to the boundary on remeshing", "DMView", PETSC_FALSE, &mesh->remeshBd, NULL);CHKERRQ(ierr); 1298 /* Projection behavior */ 1299 ierr = PetscOptionsInt("-dm_plex_max_projection_height", "Maxmimum mesh point height used to project locally", "DMPlexSetMaxProjectionHeight", 0, &mesh->maxProjectionHeight, NULL);CHKERRQ(ierr); 1300 ierr = PetscOptionsBool("-dm_plex_regular_refinement", "Use special nested projection algorithm for regular refinement", "DMPlexSetRegularRefinement", mesh->regularRefinement, &mesh->regularRefinement, NULL);CHKERRQ(ierr); 1301 PetscFunctionReturn(0); 1302 } 1303 1304 #undef __FUNCT__ 1305 #define __FUNCT__ "DMSetFromOptions_Plex" 1306 PetscErrorCode DMSetFromOptions_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 1307 { 1308 PetscInt refine = 0, coarsen = 0, r; 1309 PetscBool isHierarchy; 1310 PetscErrorCode ierr; 1311 1312 PetscFunctionBegin; 1313 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1314 ierr = PetscOptionsHead(PetscOptionsObject,"DMPlex Options");CHKERRQ(ierr); 1315 /* Handle DMPlex refinement */ 1316 ierr = PetscOptionsInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL);CHKERRQ(ierr); 1317 ierr = PetscOptionsInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy);CHKERRQ(ierr); 1318 if (refine) {ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE);CHKERRQ(ierr);} 1319 if (refine && isHierarchy) { 1320 DM *dms, coarseDM; 1321 1322 ierr = DMGetCoarseDM(dm, &coarseDM);CHKERRQ(ierr); 1323 ierr = PetscObjectReference((PetscObject)coarseDM);CHKERRQ(ierr); 1324 ierr = PetscMalloc1(refine,&dms);CHKERRQ(ierr); 1325 ierr = DMRefineHierarchy(dm, refine, dms);CHKERRQ(ierr); 1326 /* Total hack since we do not pass in a pointer */ 1327 ierr = DMPlexSwap_Static(dm, dms[refine-1]);CHKERRQ(ierr); 1328 if (refine == 1) { 1329 ierr = DMSetCoarseDM(dm, dms[0]);CHKERRQ(ierr); 1330 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr); 1331 } else { 1332 ierr = DMSetCoarseDM(dm, dms[refine-2]);CHKERRQ(ierr); 1333 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr); 1334 ierr = DMSetCoarseDM(dms[0], dms[refine-1]);CHKERRQ(ierr); 1335 ierr = DMPlexSetRegularRefinement(dms[0], PETSC_TRUE);CHKERRQ(ierr); 1336 } 1337 ierr = DMSetCoarseDM(dms[refine-1], coarseDM);CHKERRQ(ierr); 1338 ierr = PetscObjectDereference((PetscObject)coarseDM);CHKERRQ(ierr); 1339 /* Free DMs */ 1340 for (r = 0; r < refine; ++r) { 1341 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr); 1342 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 1343 } 1344 ierr = PetscFree(dms);CHKERRQ(ierr); 1345 } else { 1346 for (r = 0; r < refine; ++r) { 1347 DM refinedMesh; 1348 1349 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 1350 ierr = DMRefine(dm, PetscObjectComm((PetscObject) dm), &refinedMesh);CHKERRQ(ierr); 1351 /* Total hack since we do not pass in a pointer */ 1352 ierr = DMPlexReplace_Static(dm, refinedMesh);CHKERRQ(ierr); 1353 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 1354 ierr = DMDestroy(&refinedMesh);CHKERRQ(ierr); 1355 } 1356 } 1357 /* Handle DMPlex coarsening */ 1358 ierr = PetscOptionsInt("-dm_coarsen", "Coarsen the mesh", "DMCreate", coarsen, &coarsen, NULL);CHKERRQ(ierr); 1359 ierr = PetscOptionsInt("-dm_coarsen_hierarchy", "The number of coarsenings", "DMCreate", coarsen, &coarsen, &isHierarchy);CHKERRQ(ierr); 1360 if (coarsen && isHierarchy) { 1361 DM *dms; 1362 1363 ierr = PetscMalloc1(coarsen, &dms);CHKERRQ(ierr); 1364 ierr = DMCoarsenHierarchy(dm, coarsen, dms);CHKERRQ(ierr); 1365 /* Free DMs */ 1366 for (r = 0; r < coarsen; ++r) { 1367 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr); 1368 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 1369 } 1370 ierr = PetscFree(dms);CHKERRQ(ierr); 1371 } else { 1372 for (r = 0; r < coarsen; ++r) { 1373 DM coarseMesh; 1374 1375 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 1376 ierr = DMCoarsen(dm, PetscObjectComm((PetscObject) dm), &coarseMesh);CHKERRQ(ierr); 1377 /* Total hack since we do not pass in a pointer */ 1378 ierr = DMPlexReplace_Static(dm, coarseMesh);CHKERRQ(ierr); 1379 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 1380 ierr = DMDestroy(&coarseMesh);CHKERRQ(ierr); 1381 } 1382 } 1383 /* Handle */ 1384 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 1385 ierr = PetscOptionsTail();CHKERRQ(ierr); 1386 PetscFunctionReturn(0); 1387 } 1388 1389 #undef __FUNCT__ 1390 #define __FUNCT__ "DMCreateGlobalVector_Plex" 1391 static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec) 1392 { 1393 PetscErrorCode ierr; 1394 1395 PetscFunctionBegin; 1396 ierr = DMCreateGlobalVector_Section_Private(dm,vec);CHKERRQ(ierr); 1397 /* ierr = VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM);CHKERRQ(ierr); */ 1398 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex);CHKERRQ(ierr); 1399 ierr = VecSetOperation(*vec, VECOP_VIEWNATIVE, (void (*)(void)) VecView_Plex_Native);CHKERRQ(ierr); 1400 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex);CHKERRQ(ierr); 1401 ierr = VecSetOperation(*vec, VECOP_LOADNATIVE, (void (*)(void)) VecLoad_Plex_Native);CHKERRQ(ierr); 1402 PetscFunctionReturn(0); 1403 } 1404 1405 #undef __FUNCT__ 1406 #define __FUNCT__ "DMCreateLocalVector_Plex" 1407 static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec) 1408 { 1409 PetscErrorCode ierr; 1410 1411 PetscFunctionBegin; 1412 ierr = DMCreateLocalVector_Section_Private(dm,vec);CHKERRQ(ierr); 1413 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex_Local);CHKERRQ(ierr); 1414 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex_Local);CHKERRQ(ierr); 1415 PetscFunctionReturn(0); 1416 } 1417 1418 #undef __FUNCT__ 1419 #define __FUNCT__ "DMGetDimPoints_Plex" 1420 static PetscErrorCode DMGetDimPoints_Plex(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd) 1421 { 1422 PetscInt depth, d; 1423 PetscErrorCode ierr; 1424 1425 PetscFunctionBegin; 1426 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1427 if (depth == 1) { 1428 ierr = DMGetDimension(dm, &d);CHKERRQ(ierr); 1429 if (dim == 0) {ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr);} 1430 else if (dim == d) {ierr = DMPlexGetDepthStratum(dm, 1, pStart, pEnd);CHKERRQ(ierr);} 1431 else {*pStart = 0; *pEnd = 0;} 1432 } else { 1433 ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr); 1434 } 1435 PetscFunctionReturn(0); 1436 } 1437 1438 #undef __FUNCT__ 1439 #define __FUNCT__ "DMGetNeighors_Plex" 1440 static PetscErrorCode DMGetNeighors_Plex(DM dm, PetscInt *nranks, const PetscMPIInt *ranks[]) 1441 { 1442 PetscSF sf; 1443 PetscErrorCode ierr; 1444 1445 PetscFunctionBegin; 1446 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 1447 ierr = PetscSFGetRanks(sf, nranks, ranks, NULL, NULL, NULL);CHKERRQ(ierr); 1448 PetscFunctionReturn(0); 1449 } 1450 1451 #undef __FUNCT__ 1452 #define __FUNCT__ "DMInitialize_Plex" 1453 PetscErrorCode DMInitialize_Plex(DM dm) 1454 { 1455 PetscErrorCode ierr; 1456 1457 PetscFunctionBegin; 1458 dm->ops->view = DMView_Plex; 1459 dm->ops->load = DMLoad_Plex; 1460 dm->ops->setfromoptions = DMSetFromOptions_Plex; 1461 dm->ops->clone = DMClone_Plex; 1462 dm->ops->setup = DMSetUp_Plex; 1463 dm->ops->createdefaultsection = DMCreateDefaultSection_Plex; 1464 dm->ops->createdefaultconstraints = DMCreateDefaultConstraints_Plex; 1465 dm->ops->createglobalvector = DMCreateGlobalVector_Plex; 1466 dm->ops->createlocalvector = DMCreateLocalVector_Plex; 1467 dm->ops->getlocaltoglobalmapping = NULL; 1468 dm->ops->createfieldis = NULL; 1469 dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex; 1470 dm->ops->getcoloring = NULL; 1471 dm->ops->creatematrix = DMCreateMatrix_Plex; 1472 dm->ops->createinterpolation = DMCreateInterpolation_Plex; 1473 dm->ops->getaggregates = NULL; 1474 dm->ops->getinjection = DMCreateInjection_Plex; 1475 dm->ops->refine = DMRefine_Plex; 1476 dm->ops->coarsen = DMCoarsen_Plex; 1477 dm->ops->refinehierarchy = DMRefineHierarchy_Plex; 1478 dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Plex; 1479 dm->ops->globaltolocalbegin = NULL; 1480 dm->ops->globaltolocalend = NULL; 1481 dm->ops->localtoglobalbegin = NULL; 1482 dm->ops->localtoglobalend = NULL; 1483 dm->ops->destroy = DMDestroy_Plex; 1484 dm->ops->createsubdm = DMCreateSubDM_Plex; 1485 dm->ops->getdimpoints = DMGetDimPoints_Plex; 1486 dm->ops->locatepoints = DMLocatePoints_Plex; 1487 dm->ops->projectfunctionlocal = DMProjectFunctionLocal_Plex; 1488 dm->ops->projectfunctionlabellocal = DMProjectFunctionLabelLocal_Plex; 1489 dm->ops->projectfieldlocal = DMProjectFieldLocal_Plex; 1490 dm->ops->computel2diff = DMComputeL2Diff_Plex; 1491 dm->ops->computel2gradientdiff = DMComputeL2GradientDiff_Plex; 1492 dm->ops->computel2fielddiff = DMComputeL2FieldDiff_Plex; 1493 dm->ops->getneighbors = DMGetNeighors_Plex; 1494 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMAdaptLabel_C",DMAdaptLabel_Plex);CHKERRQ(ierr); 1495 PetscFunctionReturn(0); 1496 } 1497 1498 #undef __FUNCT__ 1499 #define __FUNCT__ "DMClone_Plex" 1500 PetscErrorCode DMClone_Plex(DM dm, DM *newdm) 1501 { 1502 DM_Plex *mesh = (DM_Plex *) dm->data; 1503 PetscErrorCode ierr; 1504 1505 PetscFunctionBegin; 1506 mesh->refct++; 1507 (*newdm)->data = mesh; 1508 ierr = PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);CHKERRQ(ierr); 1509 ierr = DMInitialize_Plex(*newdm);CHKERRQ(ierr); 1510 PetscFunctionReturn(0); 1511 } 1512 1513 /*MC 1514 DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram. 1515 In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is 1516 specified by a PetscSection object. Ownership in the global representation is determined by 1517 ownership of the underlying DMPlex points. This is specified by another PetscSection object. 1518 1519 Level: intermediate 1520 1521 .seealso: DMType, DMPlexCreate(), DMCreate(), DMSetType() 1522 M*/ 1523 1524 #undef __FUNCT__ 1525 #define __FUNCT__ "DMCreate_Plex" 1526 PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm) 1527 { 1528 DM_Plex *mesh; 1529 PetscInt unit, d; 1530 PetscErrorCode ierr; 1531 1532 PetscFunctionBegin; 1533 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1534 ierr = PetscNewLog(dm,&mesh);CHKERRQ(ierr); 1535 dm->dim = 0; 1536 dm->data = mesh; 1537 1538 mesh->refct = 1; 1539 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);CHKERRQ(ierr); 1540 mesh->maxConeSize = 0; 1541 mesh->cones = NULL; 1542 mesh->coneOrientations = NULL; 1543 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);CHKERRQ(ierr); 1544 mesh->maxSupportSize = 0; 1545 mesh->supports = NULL; 1546 mesh->refinementUniform = PETSC_TRUE; 1547 mesh->refinementLimit = -1.0; 1548 1549 mesh->facesTmp = NULL; 1550 1551 mesh->tetgenOpts = NULL; 1552 mesh->triangleOpts = NULL; 1553 ierr = PetscPartitionerCreate(PetscObjectComm((PetscObject)dm), &mesh->partitioner);CHKERRQ(ierr); 1554 ierr = PetscPartitionerSetTypeFromOptions_Internal(mesh->partitioner);CHKERRQ(ierr); 1555 mesh->remeshBd = PETSC_FALSE; 1556 1557 mesh->subpointMap = NULL; 1558 1559 for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0; 1560 1561 mesh->regularRefinement = PETSC_FALSE; 1562 mesh->depthState = -1; 1563 mesh->globalVertexNumbers = NULL; 1564 mesh->globalCellNumbers = NULL; 1565 mesh->anchorSection = NULL; 1566 mesh->anchorIS = NULL; 1567 mesh->createanchors = NULL; 1568 mesh->computeanchormatrix = NULL; 1569 mesh->parentSection = NULL; 1570 mesh->parents = NULL; 1571 mesh->childIDs = NULL; 1572 mesh->childSection = NULL; 1573 mesh->children = NULL; 1574 mesh->referenceTree = NULL; 1575 mesh->getchildsymmetry = NULL; 1576 for (d = 0; d < 8; ++d) mesh->hybridPointMax[d] = PETSC_DETERMINE; 1577 mesh->vtkCellHeight = 0; 1578 mesh->useCone = PETSC_FALSE; 1579 mesh->useClosure = PETSC_TRUE; 1580 mesh->useAnchors = PETSC_FALSE; 1581 1582 mesh->maxProjectionHeight = 0; 1583 1584 mesh->printSetValues = PETSC_FALSE; 1585 mesh->printFEM = 0; 1586 mesh->printTol = 1.0e-10; 1587 1588 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 1589 PetscFunctionReturn(0); 1590 } 1591 1592 #undef __FUNCT__ 1593 #define __FUNCT__ "DMPlexCreate" 1594 /*@ 1595 DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram. 1596 1597 Collective on MPI_Comm 1598 1599 Input Parameter: 1600 . comm - The communicator for the DMPlex object 1601 1602 Output Parameter: 1603 . mesh - The DMPlex object 1604 1605 Level: beginner 1606 1607 .keywords: DMPlex, create 1608 @*/ 1609 PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh) 1610 { 1611 PetscErrorCode ierr; 1612 1613 PetscFunctionBegin; 1614 PetscValidPointer(mesh,2); 1615 ierr = DMCreate(comm, mesh);CHKERRQ(ierr); 1616 ierr = DMSetType(*mesh, DMPLEX);CHKERRQ(ierr); 1617 PetscFunctionReturn(0); 1618 } 1619 1620 #undef __FUNCT__ 1621 #define __FUNCT__ "DMPlexBuildFromCellList_Parallel_Private" 1622 /* 1623 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 1624 */ 1625 static PetscErrorCode DMPlexBuildFromCellList_Parallel_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[], PetscSF *sfVert) 1626 { 1627 PetscSF sfPoint; 1628 PetscLayout vLayout; 1629 PetscHashI vhash; 1630 PetscSFNode *remoteVerticesAdj, *vertexLocal, *vertexOwner, *remoteVertex; 1631 const PetscInt *vrange; 1632 PetscInt numVerticesAdj, off, *verticesAdj, numVerticesGhost = 0, *localVertex, *cone, c, p, v, g; 1633 PETSC_UNUSED PetscHashIIter ret, iter; 1634 PetscMPIInt rank, numProcs; 1635 PetscErrorCode ierr; 1636 1637 PetscFunctionBegin; 1638 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr); 1639 ierr = MPI_Comm_size(PetscObjectComm((PetscObject) dm), &numProcs);CHKERRQ(ierr); 1640 /* Partition vertices */ 1641 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) dm), &vLayout);CHKERRQ(ierr); 1642 ierr = PetscLayoutSetLocalSize(vLayout, numVertices);CHKERRQ(ierr); 1643 ierr = PetscLayoutSetBlockSize(vLayout, 1);CHKERRQ(ierr); 1644 ierr = PetscLayoutSetUp(vLayout);CHKERRQ(ierr); 1645 ierr = PetscLayoutGetRanges(vLayout, &vrange);CHKERRQ(ierr); 1646 /* Count vertices and map them to procs */ 1647 PetscHashICreate(vhash); 1648 for (c = 0; c < numCells; ++c) { 1649 for (p = 0; p < numCorners; ++p) { 1650 PetscHashIPut(vhash, cells[c*numCorners+p], ret, iter); 1651 } 1652 } 1653 PetscHashISize(vhash, numVerticesAdj); 1654 ierr = PetscMalloc1(numVerticesAdj, &verticesAdj);CHKERRQ(ierr); 1655 off = 0; ierr = PetscHashIGetKeys(vhash, &off, verticesAdj);CHKERRQ(ierr); 1656 if (off != numVerticesAdj) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid number of local vertices %D should be %D", off, numVerticesAdj); 1657 ierr = PetscMalloc1(numVerticesAdj, &remoteVerticesAdj);CHKERRQ(ierr); 1658 for (v = 0; v < numVerticesAdj; ++v) { 1659 const PetscInt gv = verticesAdj[v]; 1660 PetscInt vrank; 1661 1662 ierr = PetscFindInt(gv, numProcs+1, vrange, &vrank);CHKERRQ(ierr); 1663 vrank = vrank < 0 ? -(vrank+2) : vrank; 1664 remoteVerticesAdj[v].index = gv - vrange[vrank]; 1665 remoteVerticesAdj[v].rank = vrank; 1666 } 1667 /* Create cones */ 1668 ierr = DMPlexSetChart(dm, 0, numCells+numVerticesAdj);CHKERRQ(ierr); 1669 for (c = 0; c < numCells; ++c) {ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr);} 1670 ierr = DMSetUp(dm);CHKERRQ(ierr); 1671 ierr = DMGetWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 1672 for (c = 0; c < numCells; ++c) { 1673 for (p = 0; p < numCorners; ++p) { 1674 const PetscInt gv = cells[c*numCorners+p]; 1675 PetscInt lv; 1676 1677 ierr = PetscFindInt(gv, numVerticesAdj, verticesAdj, &lv);CHKERRQ(ierr); 1678 if (lv < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not find global vertex %D in local connectivity", gv); 1679 cone[p] = lv+numCells; 1680 } 1681 ierr = DMPlexSetCone(dm, c, cone);CHKERRQ(ierr); 1682 } 1683 ierr = DMRestoreWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 1684 /* Create SF for vertices */ 1685 ierr = PetscSFCreate(PetscObjectComm((PetscObject)dm), sfVert);CHKERRQ(ierr); 1686 ierr = PetscObjectSetName((PetscObject) *sfVert, "Vertex Ownership SF");CHKERRQ(ierr); 1687 ierr = PetscSFSetFromOptions(*sfVert);CHKERRQ(ierr); 1688 ierr = PetscSFSetGraph(*sfVert, numVertices, numVerticesAdj, NULL, PETSC_OWN_POINTER, remoteVerticesAdj, PETSC_OWN_POINTER);CHKERRQ(ierr); 1689 ierr = PetscFree(verticesAdj);CHKERRQ(ierr); 1690 /* Build pointSF */ 1691 ierr = PetscMalloc2(numVerticesAdj, &vertexLocal, numVertices, &vertexOwner);CHKERRQ(ierr); 1692 for (v = 0; v < numVerticesAdj; ++v) {vertexLocal[v].index = v+numCells; vertexLocal[v].rank = rank;} 1693 for (v = 0; v < numVertices; ++v) {vertexOwner[v].index = -1; vertexOwner[v].rank = -1;} 1694 ierr = PetscSFReduceBegin(*sfVert, MPI_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);CHKERRQ(ierr); 1695 ierr = PetscSFReduceEnd(*sfVert, MPI_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);CHKERRQ(ierr); 1696 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); 1697 ierr = PetscSFBcastBegin(*sfVert, MPI_2INT, vertexOwner, vertexLocal);CHKERRQ(ierr); 1698 ierr = PetscSFBcastEnd(*sfVert, MPI_2INT, vertexOwner, vertexLocal);CHKERRQ(ierr); 1699 for (v = 0; v < numVerticesAdj; ++v) if (vertexLocal[v].rank != rank) ++numVerticesGhost; 1700 ierr = PetscMalloc1(numVerticesGhost, &localVertex);CHKERRQ(ierr); 1701 ierr = PetscMalloc1(numVerticesGhost, &remoteVertex);CHKERRQ(ierr); 1702 for (v = 0, g = 0; v < numVerticesAdj; ++v) { 1703 if (vertexLocal[v].rank != rank) { 1704 localVertex[g] = v+numCells; 1705 remoteVertex[g].index = vertexLocal[v].index; 1706 remoteVertex[g].rank = vertexLocal[v].rank; 1707 ++g; 1708 } 1709 } 1710 ierr = PetscFree2(vertexLocal, vertexOwner);CHKERRQ(ierr); 1711 if (g != numVerticesGhost) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of vertex ghosts %D should be %D", g, numVerticesGhost); 1712 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr); 1713 ierr = PetscObjectSetName((PetscObject) sfPoint, "point SF");CHKERRQ(ierr); 1714 ierr = PetscSFSetGraph(sfPoint, numCells+numVerticesAdj, numVerticesGhost, localVertex, PETSC_OWN_POINTER, remoteVertex, PETSC_OWN_POINTER);CHKERRQ(ierr); 1715 ierr = PetscLayoutDestroy(&vLayout);CHKERRQ(ierr); 1716 PetscHashIDestroy(vhash); 1717 /* Fill in the rest of the topology structure */ 1718 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 1719 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 1720 PetscFunctionReturn(0); 1721 } 1722 1723 #undef __FUNCT__ 1724 #define __FUNCT__ "DMPlexBuildCoordinates_Parallel_Private" 1725 /* 1726 This takes as input the coordinates for each owned vertex 1727 */ 1728 static PetscErrorCode DMPlexBuildCoordinates_Parallel_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscSF sfVert, const double vertexCoords[]) 1729 { 1730 PetscSection coordSection; 1731 Vec coordinates; 1732 PetscScalar *coords; 1733 PetscInt numVertices, numVerticesAdj, coordSize, v; 1734 PetscErrorCode ierr; 1735 1736 PetscFunctionBegin; 1737 ierr = PetscSFGetGraph(sfVert, &numVertices, &numVerticesAdj, NULL, NULL);CHKERRQ(ierr); 1738 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1739 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1740 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 1741 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVerticesAdj);CHKERRQ(ierr); 1742 for (v = numCells; v < numCells+numVerticesAdj; ++v) { 1743 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 1744 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 1745 } 1746 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1747 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1748 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 1749 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 1750 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1751 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1752 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1753 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1754 { 1755 MPI_Datatype coordtype; 1756 1757 /* Need a temp buffer for coords if we have complex/single */ 1758 ierr = MPI_Type_contiguous(spaceDim, MPI_DOUBLE, &coordtype);CHKERRQ(ierr); 1759 ierr = MPI_Type_commit(&coordtype);CHKERRQ(ierr); 1760 ierr = PetscSFBcastBegin(sfVert, coordtype, vertexCoords, coords);CHKERRQ(ierr); 1761 ierr = PetscSFBcastEnd(sfVert, coordtype, vertexCoords, coords);CHKERRQ(ierr); 1762 ierr = MPI_Type_free(&coordtype);CHKERRQ(ierr); 1763 } 1764 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1765 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 1766 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1767 PetscFunctionReturn(0); 1768 } 1769 1770 #undef __FUNCT__ 1771 #define __FUNCT__ "DMPlexCreateFromCellListParallel" 1772 /*@C 1773 DMPlexCreateFromCellListParallel - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM 1774 1775 Input Parameters: 1776 + comm - The communicator 1777 . dim - The topological dimension of the mesh 1778 . numCells - The number of cells owned by this process 1779 . numVertices - The number of vertices owned by this process 1780 . numCorners - The number of vertices for each cell 1781 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 1782 . cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 1783 . spaceDim - The spatial dimension used for coordinates 1784 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 1785 1786 Output Parameter: 1787 . dm - The DM 1788 1789 Note: Two triangles sharing a face 1790 $ 1791 $ 2 1792 $ / | \ 1793 $ / | \ 1794 $ / | \ 1795 $ 0 0 | 1 3 1796 $ \ | / 1797 $ \ | / 1798 $ \ | / 1799 $ 1 1800 would have input 1801 $ numCells = 2, numVertices = 4 1802 $ cells = [0 1 2 1 3 2] 1803 $ 1804 which would result in the DMPlex 1805 $ 1806 $ 4 1807 $ / | \ 1808 $ / | \ 1809 $ / | \ 1810 $ 2 0 | 1 5 1811 $ \ | / 1812 $ \ | / 1813 $ \ | / 1814 $ 3 1815 1816 Level: beginner 1817 1818 .seealso: DMPlexCreateFromCellList(), DMPlexCreateFromDAG(), DMPlexCreate() 1819 @*/ 1820 PetscErrorCode DMPlexCreateFromCellListParallel(MPI_Comm comm, PetscInt dim, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, PetscBool interpolate, const int cells[], PetscInt spaceDim, const double vertexCoords[], DM *dm) 1821 { 1822 PetscSF sfVert; 1823 PetscErrorCode ierr; 1824 1825 PetscFunctionBegin; 1826 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1827 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1828 PetscValidLogicalCollectiveInt(*dm, dim, 2); 1829 PetscValidLogicalCollectiveInt(*dm, spaceDim, 8); 1830 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1831 ierr = DMPlexBuildFromCellList_Parallel_Private(*dm, numCells, numVertices, numCorners, cells, &sfVert);CHKERRQ(ierr); 1832 if (interpolate) { 1833 DM idm = NULL; 1834 1835 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1836 ierr = DMDestroy(dm);CHKERRQ(ierr); 1837 *dm = idm; 1838 } 1839 ierr = DMPlexBuildCoordinates_Parallel_Private(*dm, spaceDim, numCells, sfVert, vertexCoords);CHKERRQ(ierr); 1840 ierr = PetscSFDestroy(&sfVert);CHKERRQ(ierr); 1841 PetscFunctionReturn(0); 1842 } 1843 1844 #undef __FUNCT__ 1845 #define __FUNCT__ "DMPlexBuildFromCellList_Private" 1846 /* 1847 This takes as input the common mesh generator output, a list of the vertices for each cell 1848 */ 1849 static PetscErrorCode DMPlexBuildFromCellList_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[]) 1850 { 1851 PetscInt *cone, c, p; 1852 PetscErrorCode ierr; 1853 1854 PetscFunctionBegin; 1855 ierr = DMPlexSetChart(dm, 0, numCells+numVertices);CHKERRQ(ierr); 1856 for (c = 0; c < numCells; ++c) { 1857 ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr); 1858 } 1859 ierr = DMSetUp(dm);CHKERRQ(ierr); 1860 ierr = DMGetWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 1861 for (c = 0; c < numCells; ++c) { 1862 for (p = 0; p < numCorners; ++p) { 1863 cone[p] = cells[c*numCorners+p]+numCells; 1864 } 1865 ierr = DMPlexSetCone(dm, c, cone);CHKERRQ(ierr); 1866 } 1867 ierr = DMRestoreWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 1868 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 1869 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 1870 PetscFunctionReturn(0); 1871 } 1872 1873 #undef __FUNCT__ 1874 #define __FUNCT__ "DMPlexBuildCoordinates_Private" 1875 /* 1876 This takes as input the coordinates for each vertex 1877 */ 1878 static PetscErrorCode DMPlexBuildCoordinates_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, const double vertexCoords[]) 1879 { 1880 PetscSection coordSection; 1881 Vec coordinates; 1882 DM cdm; 1883 PetscScalar *coords; 1884 PetscInt v, d; 1885 PetscErrorCode ierr; 1886 1887 PetscFunctionBegin; 1888 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1889 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1890 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 1891 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVertices);CHKERRQ(ierr); 1892 for (v = numCells; v < numCells+numVertices; ++v) { 1893 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 1894 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 1895 } 1896 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1897 1898 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 1899 ierr = DMCreateLocalVector(cdm, &coordinates);CHKERRQ(ierr); 1900 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 1901 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1902 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1903 for (v = 0; v < numVertices; ++v) { 1904 for (d = 0; d < spaceDim; ++d) { 1905 coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d]; 1906 } 1907 } 1908 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1909 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 1910 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1911 PetscFunctionReturn(0); 1912 } 1913 1914 #undef __FUNCT__ 1915 #define __FUNCT__ "DMPlexCreateFromCellList" 1916 /*@C 1917 DMPlexCreateFromCellList - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM 1918 1919 Input Parameters: 1920 + comm - The communicator 1921 . dim - The topological dimension of the mesh 1922 . numCells - The number of cells 1923 . numVertices - The number of vertices 1924 . numCorners - The number of vertices for each cell 1925 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 1926 . cells - An array of numCells*numCorners numbers, the vertices for each cell 1927 . spaceDim - The spatial dimension used for coordinates 1928 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 1929 1930 Output Parameter: 1931 . dm - The DM 1932 1933 Note: Two triangles sharing a face 1934 $ 1935 $ 2 1936 $ / | \ 1937 $ / | \ 1938 $ / | \ 1939 $ 0 0 | 1 3 1940 $ \ | / 1941 $ \ | / 1942 $ \ | / 1943 $ 1 1944 would have input 1945 $ numCells = 2, numVertices = 4 1946 $ cells = [0 1 2 1 3 2] 1947 $ 1948 which would result in the DMPlex 1949 $ 1950 $ 4 1951 $ / | \ 1952 $ / | \ 1953 $ / | \ 1954 $ 2 0 | 1 5 1955 $ \ | / 1956 $ \ | / 1957 $ \ | / 1958 $ 3 1959 1960 Level: beginner 1961 1962 .seealso: DMPlexCreateFromDAG(), DMPlexCreate() 1963 @*/ 1964 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) 1965 { 1966 PetscErrorCode ierr; 1967 1968 PetscFunctionBegin; 1969 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1970 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1971 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1972 ierr = DMPlexBuildFromCellList_Private(*dm, numCells, numVertices, numCorners, cells);CHKERRQ(ierr); 1973 if (interpolate) { 1974 DM idm = NULL; 1975 1976 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1977 ierr = DMDestroy(dm);CHKERRQ(ierr); 1978 *dm = idm; 1979 } 1980 ierr = DMPlexBuildCoordinates_Private(*dm, spaceDim, numCells, numVertices, vertexCoords);CHKERRQ(ierr); 1981 PetscFunctionReturn(0); 1982 } 1983 1984 #undef __FUNCT__ 1985 #define __FUNCT__ "DMPlexCreateFromDAG" 1986 /*@ 1987 DMPlexCreateFromDAG - This takes as input the adjacency-list representation of the Directed Acyclic Graph (Hasse Diagram) encoding a mesh, and produces a DM 1988 1989 Input Parameters: 1990 + dm - The empty DM object, usually from DMCreate() and DMSetDimension() 1991 . depth - The depth of the DAG 1992 . numPoints - The number of points at each depth 1993 . coneSize - The cone size of each point 1994 . cones - The concatenation of the cone points for each point, the cone list must be oriented correctly for each point 1995 . coneOrientations - The orientation of each cone point 1996 - vertexCoords - An array of numVertices*dim numbers, the coordinates of each vertex 1997 1998 Output Parameter: 1999 . dm - The DM 2000 2001 Note: Two triangles sharing a face would have input 2002 $ depth = 1, numPoints = [4 2], coneSize = [3 3 0 0 0 0] 2003 $ cones = [2 3 4 3 5 4], coneOrientations = [0 0 0 0 0 0] 2004 $ vertexCoords = [-1.0 0.0 0.0 -1.0 0.0 1.0 1.0 0.0] 2005 $ 2006 which would result in the DMPlex 2007 $ 2008 $ 4 2009 $ / | \ 2010 $ / | \ 2011 $ / | \ 2012 $ 2 0 | 1 5 2013 $ \ | / 2014 $ \ | / 2015 $ \ | / 2016 $ 3 2017 $ 2018 $ Notice that all points are numbered consecutively, unlikely DMPlexCreateFromCellList() 2019 2020 Level: advanced 2021 2022 .seealso: DMPlexCreateFromCellList(), DMPlexCreate() 2023 @*/ 2024 PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[]) 2025 { 2026 Vec coordinates; 2027 PetscSection coordSection; 2028 PetscScalar *coords; 2029 PetscInt coordSize, firstVertex = -1, pStart = 0, pEnd = 0, p, v, dim, dimEmbed, d, off; 2030 PetscErrorCode ierr; 2031 2032 PetscFunctionBegin; 2033 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2034 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 2035 if (dimEmbed < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Embedding dimension %d cannot be less than intrinsic dimension %d",dimEmbed,dim); 2036 for (d = 0; d <= depth; ++d) pEnd += numPoints[d]; 2037 ierr = DMPlexSetChart(dm, pStart, pEnd);CHKERRQ(ierr); 2038 for (p = pStart; p < pEnd; ++p) { 2039 ierr = DMPlexSetConeSize(dm, p, coneSize[p-pStart]);CHKERRQ(ierr); 2040 if (firstVertex < 0 && !coneSize[p - pStart]) { 2041 firstVertex = p - pStart; 2042 } 2043 } 2044 if (firstVertex < 0 && numPoints[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected %d vertices but could not find any", numPoints[0]); 2045 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 2046 for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) { 2047 ierr = DMPlexSetCone(dm, p, &cones[off]);CHKERRQ(ierr); 2048 ierr = DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);CHKERRQ(ierr); 2049 } 2050 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 2051 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 2052 /* Build coordinates */ 2053 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2054 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 2055 ierr = PetscSectionSetFieldComponents(coordSection, 0, dimEmbed);CHKERRQ(ierr); 2056 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0]);CHKERRQ(ierr); 2057 for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) { 2058 ierr = PetscSectionSetDof(coordSection, v, dimEmbed);CHKERRQ(ierr); 2059 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed);CHKERRQ(ierr); 2060 } 2061 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 2062 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 2063 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 2064 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 2065 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 2066 ierr = VecSetBlockSize(coordinates, dimEmbed);CHKERRQ(ierr); 2067 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 2068 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2069 for (v = 0; v < numPoints[0]; ++v) { 2070 PetscInt off; 2071 2072 ierr = PetscSectionGetOffset(coordSection, v+firstVertex, &off);CHKERRQ(ierr); 2073 for (d = 0; d < dimEmbed; ++d) { 2074 coords[off+d] = vertexCoords[v*dimEmbed+d]; 2075 } 2076 } 2077 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2078 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 2079 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 2080 PetscFunctionReturn(0); 2081 } 2082 2083 #undef __FUNCT__ 2084 #define __FUNCT__ "DMPlexCreateFromFile" 2085 /*@C 2086 DMPlexCreateFromFile - This takes a filename and produces a DM 2087 2088 Input Parameters: 2089 + comm - The communicator 2090 . filename - A file name 2091 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 2092 2093 Output Parameter: 2094 . dm - The DM 2095 2096 Level: beginner 2097 2098 .seealso: DMPlexCreateFromDAG(), DMPlexCreateFromCellList(), DMPlexCreate() 2099 @*/ 2100 PetscErrorCode DMPlexCreateFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 2101 { 2102 const char *extGmsh = ".msh"; 2103 const char *extCGNS = ".cgns"; 2104 const char *extExodus = ".exo"; 2105 const char *extFluent = ".cas"; 2106 const char *extHDF5 = ".h5"; 2107 size_t len; 2108 PetscBool isGmsh, isCGNS, isExodus, isFluent, isHDF5; 2109 PetscMPIInt rank; 2110 PetscErrorCode ierr; 2111 2112 PetscFunctionBegin; 2113 PetscValidPointer(filename, 2); 2114 PetscValidPointer(dm, 4); 2115 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 2116 ierr = PetscStrlen(filename, &len);CHKERRQ(ierr); 2117 if (!len) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Filename must be a valid path"); 2118 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGmsh, 4, &isGmsh);CHKERRQ(ierr); 2119 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extCGNS, 5, &isCGNS);CHKERRQ(ierr); 2120 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extExodus, 4, &isExodus);CHKERRQ(ierr); 2121 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extFluent, 4, &isFluent);CHKERRQ(ierr); 2122 ierr = PetscStrncmp(&filename[PetscMax(0,len-3)], extHDF5, 3, &isHDF5);CHKERRQ(ierr); 2123 if (isGmsh) { 2124 ierr = DMPlexCreateGmshFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2125 } else if (isCGNS) { 2126 ierr = DMPlexCreateCGNSFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2127 } else if (isExodus) { 2128 ierr = DMPlexCreateExodusFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2129 } else if (isFluent) { 2130 ierr = DMPlexCreateFluentFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 2131 } else if (isHDF5) { 2132 PetscViewer viewer; 2133 2134 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 2135 ierr = PetscViewerSetType(viewer, PETSCVIEWERHDF5);CHKERRQ(ierr); 2136 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 2137 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 2138 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 2139 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 2140 ierr = DMLoad(*dm, viewer);CHKERRQ(ierr); 2141 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 2142 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot load file %s: unrecognized extension", filename); 2143 PetscFunctionReturn(0); 2144 } 2145 2146 #undef __FUNCT__ 2147 #define __FUNCT__ "DMPlexCreateReferenceCell" 2148 /*@ 2149 DMPlexCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell 2150 2151 Collective on comm 2152 2153 Input Parameters: 2154 + comm - The communicator 2155 . dim - The spatial dimension 2156 - simplex - Flag for simplex, otherwise use a tensor-product cell 2157 2158 Output Parameter: 2159 . refdm - The reference cell 2160 2161 Level: intermediate 2162 2163 .keywords: reference cell 2164 .seealso: 2165 @*/ 2166 PetscErrorCode DMPlexCreateReferenceCell(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *refdm) 2167 { 2168 DM rdm; 2169 Vec coords; 2170 PetscErrorCode ierr; 2171 2172 PetscFunctionBeginUser; 2173 ierr = DMCreate(comm, &rdm);CHKERRQ(ierr); 2174 ierr = DMSetType(rdm, DMPLEX);CHKERRQ(ierr); 2175 ierr = DMSetDimension(rdm, dim);CHKERRQ(ierr); 2176 switch (dim) { 2177 case 0: 2178 { 2179 PetscInt numPoints[1] = {1}; 2180 PetscInt coneSize[1] = {0}; 2181 PetscInt cones[1] = {0}; 2182 PetscInt coneOrientations[1] = {0}; 2183 PetscScalar vertexCoords[1] = {0.0}; 2184 2185 ierr = DMPlexCreateFromDAG(rdm, 0, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2186 } 2187 break; 2188 case 1: 2189 { 2190 PetscInt numPoints[2] = {2, 1}; 2191 PetscInt coneSize[3] = {2, 0, 0}; 2192 PetscInt cones[2] = {1, 2}; 2193 PetscInt coneOrientations[2] = {0, 0}; 2194 PetscScalar vertexCoords[2] = {-1.0, 1.0}; 2195 2196 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2197 } 2198 break; 2199 case 2: 2200 if (simplex) { 2201 PetscInt numPoints[2] = {3, 1}; 2202 PetscInt coneSize[4] = {3, 0, 0, 0}; 2203 PetscInt cones[3] = {1, 2, 3}; 2204 PetscInt coneOrientations[3] = {0, 0, 0}; 2205 PetscScalar vertexCoords[6] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0}; 2206 2207 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2208 } else { 2209 PetscInt numPoints[2] = {4, 1}; 2210 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 2211 PetscInt cones[4] = {1, 2, 3, 4}; 2212 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 2213 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0}; 2214 2215 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2216 } 2217 break; 2218 case 3: 2219 if (simplex) { 2220 PetscInt numPoints[2] = {4, 1}; 2221 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 2222 PetscInt cones[4] = {1, 3, 2, 4}; 2223 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 2224 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}; 2225 2226 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2227 } else { 2228 PetscInt numPoints[2] = {8, 1}; 2229 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0}; 2230 PetscInt cones[8] = {1, 4, 3, 2, 5, 6, 7, 8}; 2231 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 2232 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, 2233 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0}; 2234 2235 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2236 } 2237 break; 2238 default: 2239 SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Cannot create reference cell for dimension %d", dim); 2240 } 2241 *refdm = NULL; 2242 ierr = DMPlexInterpolate(rdm, refdm);CHKERRQ(ierr); 2243 if (rdm->coordinateDM) { 2244 DM ncdm; 2245 PetscSection cs; 2246 PetscInt pEnd = -1; 2247 2248 ierr = DMGetDefaultSection(rdm->coordinateDM, &cs);CHKERRQ(ierr); 2249 if (cs) {ierr = PetscSectionGetChart(cs, NULL, &pEnd);CHKERRQ(ierr);} 2250 if (pEnd >= 0) { 2251 ierr = DMClone(rdm->coordinateDM, &ncdm);CHKERRQ(ierr); 2252 ierr = DMSetDefaultSection(ncdm, cs);CHKERRQ(ierr); 2253 ierr = DMSetCoordinateDM(*refdm, ncdm);CHKERRQ(ierr); 2254 ierr = DMDestroy(&ncdm);CHKERRQ(ierr); 2255 } 2256 } 2257 ierr = DMGetCoordinatesLocal(rdm, &coords);CHKERRQ(ierr); 2258 if (coords) { 2259 ierr = DMSetCoordinatesLocal(*refdm, coords);CHKERRQ(ierr); 2260 } else { 2261 ierr = DMGetCoordinates(rdm, &coords);CHKERRQ(ierr); 2262 if (coords) {ierr = DMSetCoordinates(*refdm, coords);CHKERRQ(ierr);} 2263 } 2264 ierr = DMDestroy(&rdm);CHKERRQ(ierr); 2265 PetscFunctionReturn(0); 2266 } 2267