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