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__ "DMSetFromOptions_Plex" 8 PetscErrorCode DMSetFromOptions_Plex(DM dm) 9 { 10 DM_Plex *mesh = (DM_Plex*) dm->data; 11 PetscErrorCode ierr; 12 13 PetscFunctionBegin; 14 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 15 ierr = PetscOptionsHead("DMPlex Options");CHKERRQ(ierr); 16 /* Handle DMPlex refinement */ 17 /* Handle associated vectors */ 18 /* Handle viewing */ 19 ierr = PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMView", PETSC_FALSE, &mesh->printSetValues, NULL);CHKERRQ(ierr); 20 ierr = PetscOptionsInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMView", 0, &mesh->printFEM, NULL);CHKERRQ(ierr); 21 ierr = PetscOptionsReal("-dm_plex_print_tol", "Tolerance for FEM output", "DMView", PETSC_FALSE, &mesh->printTol, NULL);CHKERRQ(ierr); 22 ierr = PetscOptionsTail();CHKERRQ(ierr); 23 PetscFunctionReturn(0); 24 } 25 26 #undef __FUNCT__ 27 #define __FUNCT__ "DMPlexCreateSquareBoundary" 28 /* 29 Simple square boundary: 30 31 18--5-17--4--16 32 | | | 33 6 10 3 34 | | | 35 19-11-20--9--15 36 | | | 37 7 8 2 38 | | | 39 12--0-13--1--14 40 */ 41 PetscErrorCode DMPlexCreateSquareBoundary(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[]) 42 { 43 PetscInt numVertices = (edges[0]+1)*(edges[1]+1); 44 PetscInt numEdges = edges[0]*(edges[1]+1) + (edges[0]+1)*edges[1]; 45 PetscInt markerTop = 1; 46 PetscInt markerBottom = 1; 47 PetscInt markerRight = 1; 48 PetscInt markerLeft = 1; 49 PetscBool markerSeparate = PETSC_FALSE; 50 Vec coordinates; 51 PetscSection coordSection; 52 PetscScalar *coords; 53 PetscInt coordSize; 54 PetscMPIInt rank; 55 PetscInt v, vx, vy; 56 PetscErrorCode ierr; 57 58 PetscFunctionBegin; 59 ierr = PetscOptionsGetBool(((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);CHKERRQ(ierr); 60 if (markerSeparate) { 61 markerTop = 1; 62 markerBottom = 0; 63 markerRight = 0; 64 markerLeft = 0; 65 } 66 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 67 if (!rank) { 68 PetscInt e, ex, ey; 69 70 ierr = DMPlexSetChart(dm, 0, numEdges+numVertices);CHKERRQ(ierr); 71 for (e = 0; e < numEdges; ++e) { 72 ierr = DMPlexSetConeSize(dm, e, 2);CHKERRQ(ierr); 73 } 74 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 75 for (vx = 0; vx <= edges[0]; vx++) { 76 for (ey = 0; ey < edges[1]; ey++) { 77 PetscInt edge = vx*edges[1] + ey + edges[0]*(edges[1]+1); 78 PetscInt vertex = ey*(edges[0]+1) + vx + numEdges; 79 PetscInt cone[2]; 80 81 cone[0] = vertex; cone[1] = vertex+edges[0]+1; 82 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 83 if (vx == edges[0]) { 84 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 85 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerRight);CHKERRQ(ierr); 86 if (ey == edges[1]-1) { 87 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerRight);CHKERRQ(ierr); 88 } 89 } else if (vx == 0) { 90 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 91 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerLeft);CHKERRQ(ierr); 92 if (ey == edges[1]-1) { 93 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerLeft);CHKERRQ(ierr); 94 } 95 } 96 } 97 } 98 for (vy = 0; vy <= edges[1]; vy++) { 99 for (ex = 0; ex < edges[0]; ex++) { 100 PetscInt edge = vy*edges[0] + ex; 101 PetscInt vertex = vy*(edges[0]+1) + ex + numEdges; 102 PetscInt cone[2]; 103 104 cone[0] = vertex; cone[1] = vertex+1; 105 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 106 if (vy == edges[1]) { 107 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 108 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerTop);CHKERRQ(ierr); 109 if (ex == edges[0]-1) { 110 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerTop);CHKERRQ(ierr); 111 } 112 } else if (vy == 0) { 113 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 114 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerBottom);CHKERRQ(ierr); 115 if (ex == edges[0]-1) { 116 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerBottom);CHKERRQ(ierr); 117 } 118 } 119 } 120 } 121 } 122 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 123 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 124 /* Build coordinates */ 125 ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 126 ierr = PetscSectionSetChart(coordSection, numEdges, numEdges + numVertices);CHKERRQ(ierr); 127 for (v = numEdges; v < numEdges+numVertices; ++v) { 128 ierr = PetscSectionSetDof(coordSection, v, 2);CHKERRQ(ierr); 129 } 130 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 131 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 132 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 133 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 134 ierr = VecSetFromOptions(coordinates);CHKERRQ(ierr); 135 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 136 for (vy = 0; vy <= edges[1]; ++vy) { 137 for (vx = 0; vx <= edges[0]; ++vx) { 138 coords[(vy*(edges[0]+1)+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/edges[0])*vx; 139 coords[(vy*(edges[0]+1)+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/edges[1])*vy; 140 } 141 } 142 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 143 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 144 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 145 PetscFunctionReturn(0); 146 } 147 148 #undef __FUNCT__ 149 #define __FUNCT__ "DMPlexCreateCubeBoundary" 150 /* 151 Simple cubic boundary: 152 153 2-------3 154 /| /| 155 6-------7 | 156 | | | | 157 | 0-----|-1 158 |/ |/ 159 4-------5 160 */ 161 PetscErrorCode DMPlexCreateCubeBoundary(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt faces[]) 162 { 163 PetscInt numVertices = (faces[0]+1)*(faces[1]+1)*(faces[2]+1); 164 PetscInt numFaces = 6; 165 Vec coordinates; 166 PetscSection coordSection; 167 PetscScalar *coords; 168 PetscInt coordSize; 169 PetscMPIInt rank; 170 PetscInt v, vx, vy, vz; 171 PetscErrorCode ierr; 172 173 PetscFunctionBegin; 174 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"); 175 if ((faces[0] > 1) || (faces[1] > 1) || (faces[2] > 1)) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Currently can't handle more than 1 face per side"); 176 ierr = PetscMalloc(numVertices*2 * sizeof(PetscReal), &coords);CHKERRQ(ierr); 177 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 178 if (!rank) { 179 PetscInt f; 180 181 ierr = DMPlexSetChart(dm, 0, numFaces+numVertices);CHKERRQ(ierr); 182 for (f = 0; f < numFaces; ++f) { 183 ierr = DMPlexSetConeSize(dm, f, 4);CHKERRQ(ierr); 184 } 185 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 186 for (v = 0; v < numFaces+numVertices; ++v) { 187 ierr = DMPlexSetLabelValue(dm, "marker", v, 1);CHKERRQ(ierr); 188 } 189 { /* Side 0 (Front) */ 190 PetscInt cone[4]; 191 cone[0] = numFaces+4; cone[1] = numFaces+5; cone[2] = numFaces+7; cone[3] = numFaces+6; 192 ierr = DMPlexSetCone(dm, 0, cone);CHKERRQ(ierr); 193 } 194 { /* Side 1 (Back) */ 195 PetscInt cone[4]; 196 cone[0] = numFaces+1; cone[1] = numFaces+0; cone[2] = numFaces+2; cone[3] = numFaces+3; 197 ierr = DMPlexSetCone(dm, 1, cone);CHKERRQ(ierr); 198 } 199 { /* Side 2 (Bottom) */ 200 PetscInt cone[4]; 201 cone[0] = numFaces+0; cone[1] = numFaces+1; cone[2] = numFaces+5; cone[3] = numFaces+4; 202 ierr = DMPlexSetCone(dm, 2, cone);CHKERRQ(ierr); 203 } 204 { /* Side 3 (Top) */ 205 PetscInt cone[4]; 206 cone[0] = numFaces+6; cone[1] = numFaces+7; cone[2] = numFaces+3; cone[3] = numFaces+2; 207 ierr = DMPlexSetCone(dm, 3, cone);CHKERRQ(ierr); 208 } 209 { /* Side 4 (Left) */ 210 PetscInt cone[4]; 211 cone[0] = numFaces+0; cone[1] = numFaces+4; cone[2] = numFaces+6; cone[3] = numFaces+2; 212 ierr = DMPlexSetCone(dm, 4, cone);CHKERRQ(ierr); 213 } 214 { /* Side 5 (Right) */ 215 PetscInt cone[4]; 216 cone[0] = numFaces+5; cone[1] = numFaces+1; cone[2] = numFaces+3; cone[3] = numFaces+7; 217 ierr = DMPlexSetCone(dm, 5, cone);CHKERRQ(ierr); 218 } 219 } 220 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 221 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 222 /* Build coordinates */ 223 ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 224 ierr = PetscSectionSetChart(coordSection, numFaces, numFaces + numVertices);CHKERRQ(ierr); 225 for (v = numFaces; v < numFaces+numVertices; ++v) { 226 ierr = PetscSectionSetDof(coordSection, v, 3);CHKERRQ(ierr); 227 } 228 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 229 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 230 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 231 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 232 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 233 ierr = VecSetFromOptions(coordinates);CHKERRQ(ierr); 234 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 235 for (vz = 0; vz <= faces[2]; ++vz) { 236 for (vy = 0; vy <= faces[1]; ++vy) { 237 for (vx = 0; vx <= faces[0]; ++vx) { 238 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+0] = lower[0] + ((upper[0] - lower[0])/faces[0])*vx; 239 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+1] = lower[1] + ((upper[1] - lower[1])/faces[1])*vy; 240 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+2] = lower[2] + ((upper[2] - lower[2])/faces[2])*vz; 241 } 242 } 243 } 244 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 245 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 246 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 247 PetscFunctionReturn(0); 248 } 249 250 #undef __FUNCT__ 251 #define __FUNCT__ "DMPlexCreateSquareMesh" 252 /* 253 Simple square mesh: 254 255 22--8-23--9--24 256 | | | 257 13 2 14 3 15 258 | | | 259 19--6-20--7--21 260 | | | 261 10 0 11 1 12 262 | | | 263 16--4-17--5--18 264 */ 265 PetscErrorCode DMPlexCreateSquareMesh(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[]) 266 { 267 PetscInt markerTop = 1; 268 PetscInt markerBottom = 1; 269 PetscInt markerRight = 1; 270 PetscInt markerLeft = 1; 271 PetscBool markerSeparate = PETSC_FALSE; 272 PetscMPIInt rank; 273 PetscErrorCode ierr; 274 275 PetscFunctionBegin; 276 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 277 ierr = PetscOptionsGetBool(((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);CHKERRQ(ierr); 278 if (markerSeparate) { 279 markerTop = 3; 280 markerBottom = 1; 281 markerRight = 2; 282 markerLeft = 4; 283 } 284 { 285 const PetscInt numXEdges = !rank ? edges[0] : 0; 286 const PetscInt numYEdges = !rank ? edges[1] : 0; 287 const PetscInt numXVertices = !rank ? edges[0]+1 : 0; 288 const PetscInt numYVertices = !rank ? edges[1]+1 : 0; 289 const PetscInt numTotXEdges = numXEdges*numYVertices; 290 const PetscInt numTotYEdges = numYEdges*numXVertices; 291 const PetscInt numVertices = numXVertices*numYVertices; 292 const PetscInt numEdges = numTotXEdges + numTotYEdges; 293 const PetscInt numFaces = numXEdges*numYEdges; 294 const PetscInt firstVertex = numFaces; 295 const PetscInt firstXEdge = numFaces + numVertices; 296 const PetscInt firstYEdge = numFaces + numVertices + numTotXEdges; 297 Vec coordinates; 298 PetscSection coordSection; 299 PetscScalar *coords; 300 PetscInt coordSize; 301 PetscInt v, vx, vy; 302 PetscInt f, fx, fy, e, ex, ey; 303 304 ierr = DMPlexSetChart(dm, 0, numFaces+numEdges+numVertices);CHKERRQ(ierr); 305 for (f = 0; f < numFaces; ++f) { 306 ierr = DMPlexSetConeSize(dm, f, 4);CHKERRQ(ierr); 307 } 308 for (e = firstXEdge; e < firstXEdge+numEdges; ++e) { 309 ierr = DMPlexSetConeSize(dm, e, 2);CHKERRQ(ierr); 310 } 311 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 312 /* Build faces */ 313 for (fy = 0; fy < numYEdges; fy++) { 314 for (fx = 0; fx < numXEdges; fx++) { 315 const PetscInt face = fy*numXEdges + fx; 316 const PetscInt edgeL = firstYEdge + fx*numYEdges + fy; 317 const PetscInt edgeB = firstXEdge + fy*numXEdges + fx; 318 const PetscInt ornt[4] = {0, 0, -2, -2}; 319 PetscInt cone[4]; 320 321 cone[0] = edgeB; cone[1] = edgeL+numYEdges; cone[2] = edgeB+numXEdges; cone[3] = edgeL; 322 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr); 323 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr); 324 } 325 } 326 /* Build Y edges*/ 327 for (vx = 0; vx < numXVertices; vx++) { 328 for (ey = 0; ey < numYEdges; ey++) { 329 const PetscInt edge = firstYEdge + vx*numYEdges + ey; 330 const PetscInt vertex = firstVertex + ey*numXVertices + vx; 331 PetscInt cone[2]; 332 333 cone[0] = vertex; cone[1] = vertex+numXVertices; 334 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 335 if (vx == numXVertices-1) { 336 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 337 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerRight);CHKERRQ(ierr); 338 if (ey == numYEdges-1) { 339 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerRight);CHKERRQ(ierr); 340 } 341 } else if (vx == 0) { 342 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 343 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerLeft);CHKERRQ(ierr); 344 if (ey == numYEdges-1) { 345 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerLeft);CHKERRQ(ierr); 346 } 347 } 348 } 349 } 350 /* Build X edges*/ 351 for (vy = 0; vy < numYVertices; vy++) { 352 for (ex = 0; ex < numXEdges; ex++) { 353 const PetscInt edge = firstXEdge + vy*numXEdges + ex; 354 const PetscInt vertex = firstVertex + vy*numXVertices + ex; 355 PetscInt cone[2]; 356 357 cone[0] = vertex; cone[1] = vertex+1; 358 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 359 if (vy == numYVertices-1) { 360 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 361 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerTop);CHKERRQ(ierr); 362 if (ex == numXEdges-1) { 363 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerTop);CHKERRQ(ierr); 364 } 365 } else if (vy == 0) { 366 ierr = DMPlexSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 367 ierr = DMPlexSetLabelValue(dm, "marker", cone[0], markerBottom);CHKERRQ(ierr); 368 if (ex == numXEdges-1) { 369 ierr = DMPlexSetLabelValue(dm, "marker", cone[1], markerBottom);CHKERRQ(ierr); 370 } 371 } 372 } 373 } 374 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 375 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 376 /* Build coordinates */ 377 ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 378 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVertices);CHKERRQ(ierr); 379 for (v = firstVertex; v < firstVertex+numVertices; ++v) { 380 ierr = PetscSectionSetDof(coordSection, v, 2);CHKERRQ(ierr); 381 } 382 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 383 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 384 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 385 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 386 ierr = VecSetFromOptions(coordinates);CHKERRQ(ierr); 387 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 388 for (vy = 0; vy < numYVertices; ++vy) { 389 for (vx = 0; vx < numXVertices; ++vx) { 390 coords[(vy*numXVertices+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/numXEdges)*vx; 391 coords[(vy*numXVertices+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/numYEdges)*vy; 392 } 393 } 394 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 395 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 396 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 397 } 398 PetscFunctionReturn(0); 399 } 400 401 #undef __FUNCT__ 402 #define __FUNCT__ "DMPlexCreateBoxMesh" 403 /*@ 404 DMPlexCreateBoxMesh - Creates a mesh on the tensor product of unit intervals (box). 405 406 Collective on MPI_Comm 407 408 Input Parameters: 409 + comm - The communicator for the DM object 410 . dim - The spatial dimension 411 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 412 413 Output Parameter: 414 . dm - The DM object 415 416 Level: beginner 417 418 .keywords: DM, create 419 .seealso: DMSetType(), DMCreate() 420 @*/ 421 PetscErrorCode DMPlexCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscBool interpolate, DM *dm) 422 { 423 DM boundary; 424 PetscErrorCode ierr; 425 426 PetscFunctionBegin; 427 PetscValidPointer(dm, 4); 428 ierr = DMCreate(comm, &boundary);CHKERRQ(ierr); 429 PetscValidLogicalCollectiveInt(boundary,dim,2); 430 ierr = DMSetType(boundary, DMPLEX);CHKERRQ(ierr); 431 ierr = DMPlexSetDimension(boundary, dim-1);CHKERRQ(ierr); 432 switch (dim) { 433 case 2: 434 { 435 PetscReal lower[2] = {0.0, 0.0}; 436 PetscReal upper[2] = {1.0, 1.0}; 437 PetscInt edges[2] = {2, 2}; 438 439 ierr = DMPlexCreateSquareBoundary(boundary, lower, upper, edges);CHKERRQ(ierr); 440 break; 441 } 442 case 3: 443 { 444 PetscReal lower[3] = {0.0, 0.0, 0.0}; 445 PetscReal upper[3] = {1.0, 1.0, 1.0}; 446 PetscInt faces[3] = {1, 1, 1}; 447 448 ierr = DMPlexCreateCubeBoundary(boundary, lower, upper, faces);CHKERRQ(ierr); 449 break; 450 } 451 default: 452 SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim); 453 } 454 ierr = DMPlexGenerate(boundary, NULL, interpolate, dm);CHKERRQ(ierr); 455 ierr = DMDestroy(&boundary);CHKERRQ(ierr); 456 PetscFunctionReturn(0); 457 } 458 459 #undef __FUNCT__ 460 #define __FUNCT__ "DMPlexCreateHexBoxMesh" 461 PetscErrorCode DMPlexCreateHexBoxMesh(MPI_Comm comm, PetscInt dim, const PetscInt cells[], DM *dm) 462 { 463 PetscErrorCode ierr; 464 465 PetscFunctionBegin; 466 PetscValidPointer(dm, 4); 467 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 468 PetscValidLogicalCollectiveInt(*dm,dim,2); 469 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 470 ierr = DMPlexSetDimension(*dm, dim);CHKERRQ(ierr); 471 switch (dim) { 472 case 2: 473 { 474 PetscReal lower[2] = {0.0, 0.0}; 475 PetscReal upper[2] = {1.0, 1.0}; 476 477 ierr = DMPlexCreateSquareMesh(*dm, lower, upper, cells);CHKERRQ(ierr); 478 break; 479 } 480 #if 0 481 case 3: 482 { 483 PetscReal lower[3] = {0.0, 0.0, 0.0}; 484 PetscReal upper[3] = {1.0, 1.0, 1.0}; 485 486 ierr = DMPlexCreateCubeMesh(boundary, lower, upper, cells);CHKERRQ(ierr); 487 break; 488 } 489 #endif 490 default: 491 SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim); 492 } 493 PetscFunctionReturn(0); 494 } 495 496 /* External function declarations here */ 497 extern PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling); 498 extern PetscErrorCode DMCreateMatrix_Plex(DM dm, MatType mtype, Mat *J); 499 extern PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm); 500 extern PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *dmRefined); 501 extern PetscErrorCode DMClone_Plex(DM dm, DM *newdm); 502 extern PetscErrorCode DMSetUp_Plex(DM dm); 503 extern PetscErrorCode DMDestroy_Plex(DM dm); 504 extern PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer); 505 extern PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm); 506 extern PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, IS *cellIS); 507 508 #undef __FUNCT__ 509 #define __FUNCT__ "DMCreateGlobalVector_Plex" 510 static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec) 511 { 512 PetscErrorCode ierr; 513 514 PetscFunctionBegin; 515 ierr = DMCreateGlobalVector_Section_Private(dm,vec);CHKERRQ(ierr); 516 /* ierr = VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM);CHKERRQ(ierr); */ 517 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void))VecView_Plex);CHKERRQ(ierr); 518 PetscFunctionReturn(0); 519 } 520 521 #undef __FUNCT__ 522 #define __FUNCT__ "DMCreateLocalVector_Plex" 523 static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec) 524 { 525 PetscErrorCode ierr; 526 527 PetscFunctionBegin; 528 ierr = DMCreateLocalVector_Section_Private(dm,vec);CHKERRQ(ierr); 529 ierr = VecSetOperation(*vec, VECOP_VIEW, (void(*)(void)) VecView_Plex_Local);CHKERRQ(ierr); 530 PetscFunctionReturn(0); 531 } 532 533 #undef __FUNCT__ 534 #define __FUNCT__ "DMInitialize_Plex" 535 PetscErrorCode DMInitialize_Plex(DM dm) 536 { 537 PetscErrorCode ierr; 538 539 PetscFunctionBegin; 540 ierr = PetscStrallocpy(VECSTANDARD, (char**)&dm->vectype);CHKERRQ(ierr); 541 542 dm->ops->view = DMView_Plex; 543 dm->ops->setfromoptions = DMSetFromOptions_Plex; 544 dm->ops->clone = DMClone_Plex; 545 dm->ops->setup = DMSetUp_Plex; 546 dm->ops->createglobalvector = DMCreateGlobalVector_Plex; 547 dm->ops->createlocalvector = DMCreateLocalVector_Plex; 548 dm->ops->getlocaltoglobalmapping = NULL; 549 dm->ops->getlocaltoglobalmappingblock = NULL; 550 dm->ops->createfieldis = NULL; 551 dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex; 552 dm->ops->getcoloring = 0; 553 dm->ops->creatematrix = DMCreateMatrix_Plex; 554 dm->ops->createinterpolation = 0; 555 dm->ops->getaggregates = 0; 556 dm->ops->getinjection = 0; 557 dm->ops->refine = DMRefine_Plex; 558 dm->ops->coarsen = 0; 559 dm->ops->refinehierarchy = 0; 560 dm->ops->coarsenhierarchy = 0; 561 dm->ops->globaltolocalbegin = NULL; 562 dm->ops->globaltolocalend = NULL; 563 dm->ops->localtoglobalbegin = NULL; 564 dm->ops->localtoglobalend = NULL; 565 dm->ops->destroy = DMDestroy_Plex; 566 dm->ops->createsubdm = DMCreateSubDM_Plex; 567 dm->ops->locatepoints = DMLocatePoints_Plex; 568 PetscFunctionReturn(0); 569 } 570 571 #undef __FUNCT__ 572 #define __FUNCT__ "DMClone_Plex" 573 PetscErrorCode DMClone_Plex(DM dm, DM *newdm) 574 { 575 DM_Plex *mesh = (DM_Plex *) dm->data; 576 PetscErrorCode ierr; 577 578 PetscFunctionBegin; 579 mesh->refct++; 580 (*newdm)->data = mesh; 581 ierr = PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);CHKERRQ(ierr); 582 ierr = DMInitialize_Plex(*newdm);CHKERRQ(ierr); 583 PetscFunctionReturn(0); 584 } 585 586 /*MC 587 DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram. 588 In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is 589 specified by a PetscSection object. Ownership in the global representation is determined by 590 ownership of the underlying DMPlex points. This is specified by another PetscSection object. 591 592 Level: intermediate 593 594 .seealso: DMType, DMPlexCreate(), DMCreate(), DMSetType() 595 M*/ 596 597 #undef __FUNCT__ 598 #define __FUNCT__ "DMCreate_Plex" 599 PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm) 600 { 601 DM_Plex *mesh; 602 PetscInt unit, d; 603 PetscErrorCode ierr; 604 605 PetscFunctionBegin; 606 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 607 ierr = PetscNewLog(dm, DM_Plex, &mesh);CHKERRQ(ierr); 608 dm->data = mesh; 609 610 mesh->refct = 1; 611 mesh->dim = 0; 612 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);CHKERRQ(ierr); 613 mesh->maxConeSize = 0; 614 mesh->cones = NULL; 615 mesh->coneOrientations = NULL; 616 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);CHKERRQ(ierr); 617 mesh->maxSupportSize = 0; 618 mesh->supports = NULL; 619 mesh->refinementUniform = PETSC_TRUE; 620 mesh->refinementLimit = -1.0; 621 622 mesh->facesTmp = NULL; 623 624 mesh->subpointMap = NULL; 625 626 for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0; 627 628 mesh->labels = NULL; 629 mesh->depthLabel = NULL; 630 mesh->globalVertexNumbers = NULL; 631 mesh->globalCellNumbers = NULL; 632 for (d = 0; d < 8; ++d) mesh->hybridPointMax[d] = PETSC_DETERMINE; 633 mesh->vtkCellHeight = 0; 634 mesh->preallocCenterDim = -1; 635 636 mesh->printSetValues = PETSC_FALSE; 637 mesh->printFEM = 0; 638 mesh->printTol = 1.0e-10; 639 640 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 641 PetscFunctionReturn(0); 642 } 643 644 #undef __FUNCT__ 645 #define __FUNCT__ "DMPlexCreate" 646 /*@ 647 DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram. 648 649 Collective on MPI_Comm 650 651 Input Parameter: 652 . comm - The communicator for the DMPlex object 653 654 Output Parameter: 655 . mesh - The DMPlex object 656 657 Level: beginner 658 659 .keywords: DMPlex, create 660 @*/ 661 PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh) 662 { 663 PetscErrorCode ierr; 664 665 PetscFunctionBegin; 666 PetscValidPointer(mesh,2); 667 ierr = DMCreate(comm, mesh);CHKERRQ(ierr); 668 ierr = DMSetType(*mesh, DMPLEX);CHKERRQ(ierr); 669 PetscFunctionReturn(0); 670 } 671 672 #undef __FUNCT__ 673 #define __FUNCT__ "DMPlexBuildFromCellList_Private" 674 /* 675 This takes as input the common mesh generator output, a list of the vertices for each cell 676 */ 677 PetscErrorCode DMPlexBuildFromCellList_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[]) 678 { 679 PetscInt *cone, c, p; 680 PetscErrorCode ierr; 681 682 PetscFunctionBegin; 683 ierr = DMPlexSetChart(dm, 0, numCells+numVertices);CHKERRQ(ierr); 684 for (c = 0; c < numCells; ++c) { 685 ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr); 686 } 687 ierr = DMSetUp(dm);CHKERRQ(ierr); 688 ierr = DMGetWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 689 for (c = 0; c < numCells; ++c) { 690 for (p = 0; p < numCorners; ++p) { 691 cone[p] = cells[c*numCorners+p]+numCells; 692 } 693 ierr = DMPlexSetCone(dm, c, cone);CHKERRQ(ierr); 694 } 695 ierr = DMRestoreWorkArray(dm, numCorners, PETSC_INT, &cone);CHKERRQ(ierr); 696 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 697 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 698 PetscFunctionReturn(0); 699 } 700 701 #undef __FUNCT__ 702 #define __FUNCT__ "DMPlexBuildCoordinates_Private" 703 /* 704 This takes as input the coordinates for each vertex 705 */ 706 PetscErrorCode DMPlexBuildCoordinates_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, const double vertexCoords[]) 707 { 708 PetscSection coordSection; 709 Vec coordinates; 710 PetscScalar *coords; 711 PetscInt coordSize, v, d; 712 PetscErrorCode ierr; 713 714 PetscFunctionBegin; 715 ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 716 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 717 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 718 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVertices);CHKERRQ(ierr); 719 for (v = numCells; v < numCells+numVertices; ++v) { 720 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 721 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 722 } 723 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 724 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 725 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 726 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 727 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 728 ierr = VecSetFromOptions(coordinates);CHKERRQ(ierr); 729 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 730 for (v = 0; v < numVertices; ++v) { 731 for (d = 0; d < spaceDim; ++d) { 732 coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d]; 733 } 734 } 735 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 736 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 737 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 738 PetscFunctionReturn(0); 739 } 740 741 #undef __FUNCT__ 742 #define __FUNCT__ "DMPlexCreateFromCellList" 743 /*@C 744 DMPlexCreateFromCellList - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM 745 746 Input Parameters: 747 + comm - The communicator 748 . dim - The topological dimension of the mesh 749 . numCells - The number of cells 750 . numVertices - The number of vertices 751 . numCorners - The number of vertices for each cell 752 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 753 . cells - An array of numCells*numCorners numbers, the vertices for each cell 754 . spaceDim - The spatial dimension used for coordinates 755 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 756 757 Output Parameter: 758 . dm - The DM 759 760 Note: Two triangles sharing a face 761 $ 762 $ 2 763 $ / | \ 764 $ / | \ 765 $ / | \ 766 $ 0 0 | 1 3 767 $ \ | / 768 $ \ | / 769 $ \ | / 770 $ 1 771 would have input 772 $ numCells = 2, numVertices = 4 773 $ cells = [0 1 2 1 3 2] 774 $ 775 which would result in the DMPlex 776 $ 777 $ 4 778 $ / | \ 779 $ / | \ 780 $ / | \ 781 $ 2 0 | 1 5 782 $ \ | / 783 $ \ | / 784 $ \ | / 785 $ 3 786 787 Level: beginner 788 789 .seealso: DMPlexCreate() 790 @*/ 791 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) 792 { 793 PetscErrorCode ierr; 794 795 PetscFunctionBegin; 796 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 797 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 798 ierr = DMPlexSetDimension(*dm, dim);CHKERRQ(ierr); 799 ierr = DMPlexBuildFromCellList_Private(*dm, numCells, numVertices, numCorners, cells);CHKERRQ(ierr); 800 if (interpolate) { 801 DM idm; 802 803 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 804 ierr = DMDestroy(dm);CHKERRQ(ierr); 805 *dm = idm; 806 } 807 ierr = DMPlexBuildCoordinates_Private(*dm, spaceDim, numCells, numVertices, vertexCoords);CHKERRQ(ierr); 808 PetscFunctionReturn(0); 809 } 810 811 #undef __FUNCT__ 812 #define __FUNCT__ "DMPlexCreateFromDAG" 813 /* 814 This takes as input the raw Hasse Diagram data 815 */ 816 PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[]) 817 { 818 Vec coordinates; 819 PetscSection coordSection; 820 PetscScalar *coords; 821 PetscInt coordSize, firstVertex = numPoints[depth], pStart = 0, pEnd = 0, p, v, dim, d, off; 822 PetscErrorCode ierr; 823 824 PetscFunctionBegin; 825 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 826 for (d = 0; d <= depth; ++d) pEnd += numPoints[d]; 827 ierr = DMPlexSetChart(dm, pStart, pEnd);CHKERRQ(ierr); 828 for (p = pStart; p < pEnd; ++p) { 829 ierr = DMPlexSetConeSize(dm, p, coneSize[p-pStart]);CHKERRQ(ierr); 830 } 831 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 832 for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) { 833 ierr = DMPlexSetCone(dm, p, &cones[off]);CHKERRQ(ierr); 834 ierr = DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);CHKERRQ(ierr); 835 } 836 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 837 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 838 /* Build coordinates */ 839 ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 840 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 841 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 842 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0]);CHKERRQ(ierr); 843 for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) { 844 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 845 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 846 } 847 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 848 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 849 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 850 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 851 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 852 ierr = VecSetFromOptions(coordinates);CHKERRQ(ierr); 853 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 854 for (v = 0; v < numPoints[0]; ++v) { 855 PetscInt off; 856 857 ierr = PetscSectionGetOffset(coordSection, v+firstVertex, &off);CHKERRQ(ierr); 858 for (d = 0; d < dim; ++d) { 859 coords[off+d] = vertexCoords[v*dim+d]; 860 } 861 } 862 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 863 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 864 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 865 PetscFunctionReturn(0); 866 } 867