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 EXTERN_C_BEGIN 551 #undef __FUNCT__ 552 #define __FUNCT__ "DMCreate_Plex" 553 PetscErrorCode DMCreate_Plex(DM dm) 554 { 555 DM_Plex *mesh; 556 PetscInt unit, d; 557 PetscErrorCode ierr; 558 559 PetscFunctionBegin; 560 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 561 ierr = PetscNewLog(dm, DM_Plex, &mesh);CHKERRQ(ierr); 562 dm->data = mesh; 563 564 mesh->refct = 1; 565 mesh->dim = 0; 566 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);CHKERRQ(ierr); 567 mesh->maxConeSize = 0; 568 mesh->cones = NULL; 569 mesh->coneOrientations = NULL; 570 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);CHKERRQ(ierr); 571 mesh->maxSupportSize = 0; 572 mesh->supports = NULL; 573 mesh->refinementUniform = PETSC_TRUE; 574 mesh->refinementLimit = -1.0; 575 576 mesh->facesTmp = NULL; 577 578 mesh->subpointMap = NULL; 579 580 for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0; 581 582 mesh->labels = NULL; 583 mesh->globalVertexNumbers = NULL; 584 mesh->globalCellNumbers = NULL; 585 for (d = 0; d < 8; ++d) mesh->hybridPointMax[d] = PETSC_DETERMINE; 586 mesh->vtkCellHeight = 0; 587 mesh->preallocCenterDim = -1; 588 589 mesh->integrateResidualFEM = NULL; 590 mesh->integrateJacobianActionFEM = NULL; 591 mesh->integrateJacobianFEM = NULL; 592 593 mesh->printSetValues = PETSC_FALSE; 594 mesh->printFEM = 0; 595 596 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 597 PetscFunctionReturn(0); 598 } 599 EXTERN_C_END 600 601 #undef __FUNCT__ 602 #define __FUNCT__ "DMPlexCreate" 603 /*@ 604 DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram. 605 606 Collective on MPI_Comm 607 608 Input Parameter: 609 . comm - The communicator for the DMPlex object 610 611 Output Parameter: 612 . mesh - The DMPlex object 613 614 Level: beginner 615 616 .keywords: DMPlex, create 617 @*/ 618 PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh) 619 { 620 PetscErrorCode ierr; 621 622 PetscFunctionBegin; 623 PetscValidPointer(mesh,2); 624 ierr = DMCreate(comm, mesh);CHKERRQ(ierr); 625 ierr = DMSetType(*mesh, DMPLEX);CHKERRQ(ierr); 626 PetscFunctionReturn(0); 627 } 628 629 #undef __FUNCT__ 630 #define __FUNCT__ "DMPlexClone" 631 /*@ 632 DMPlexClone - Creates a DMPlex object with the same mesh as the original. 633 634 Collective on MPI_Comm 635 636 Input Parameter: 637 . dm - The original DMPlex object 638 639 Output Parameter: 640 . newdm - The new DMPlex object 641 642 Level: beginner 643 644 .keywords: DMPlex, create 645 @*/ 646 PetscErrorCode DMPlexClone(DM dm, DM *newdm) 647 { 648 DM_Plex *mesh; 649 void *ctx; 650 PetscErrorCode ierr; 651 652 PetscFunctionBegin; 653 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 654 PetscValidPointer(newdm,2); 655 ierr = DMCreate(PetscObjectComm((PetscObject)dm), newdm);CHKERRQ(ierr); 656 ierr = PetscSFDestroy(&(*newdm)->sf);CHKERRQ(ierr); 657 ierr = PetscObjectReference((PetscObject) dm->sf);CHKERRQ(ierr); 658 (*newdm)->sf = dm->sf; 659 mesh = (DM_Plex*) dm->data; 660 mesh->refct++; 661 (*newdm)->data = mesh; 662 ierr = PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);CHKERRQ(ierr); 663 ierr = DMInitialize_Plex(*newdm);CHKERRQ(ierr); 664 ierr = DMGetApplicationContext(dm, &ctx);CHKERRQ(ierr); 665 ierr = DMSetApplicationContext(*newdm, ctx);CHKERRQ(ierr); 666 PetscFunctionReturn(0); 667 } 668