1 #define PETSCDM_DLL 2 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 3 #include <petsc/private/hashseti.h> /*I "petscdmplex.h" I*/ 4 #include <petscsf.h> 5 6 PetscLogEvent DMPLEX_CreateFromFile, DMPLEX_BuildFromCellList, DMPLEX_BuildCoordinatesFromCellList; 7 8 /* External function declarations here */ 9 static PetscErrorCode DMInitialize_Plex(DM dm); 10 11 /* This copies internal things in the Plex structure that we generally want when making a new, related Plex */ 12 PetscErrorCode DMPlexCopy_Internal(DM dmin, PetscBool copyPeriodicity, DM dmout) 13 { 14 const DMBoundaryType *bd; 15 const PetscReal *maxCell, *L; 16 PetscBool isper, dist; 17 PetscErrorCode ierr; 18 19 PetscFunctionBegin; 20 if (copyPeriodicity) { 21 ierr = DMGetPeriodicity(dmin, &isper, &maxCell, &L, &bd);CHKERRQ(ierr); 22 ierr = DMSetPeriodicity(dmout, isper, maxCell, L, bd);CHKERRQ(ierr); 23 } 24 ierr = DMPlexDistributeGetDefault(dmin, &dist);CHKERRQ(ierr); 25 ierr = DMPlexDistributeSetDefault(dmout, dist);CHKERRQ(ierr); 26 ((DM_Plex *) dmout->data)->useHashLocation = ((DM_Plex *) dmin->data)->useHashLocation; 27 PetscFunctionReturn(0); 28 } 29 30 /* Replace dm with the contents of ndm, and then destroy ndm 31 - Share the DM_Plex structure 32 - Share the coordinates 33 - Share the SF 34 */ 35 static PetscErrorCode DMPlexReplace_Static(DM dm, DM *ndm) 36 { 37 PetscSF sf; 38 DM dmNew = *ndm, coordDM, coarseDM; 39 Vec coords; 40 PetscBool isper; 41 const PetscReal *maxCell, *L; 42 const DMBoundaryType *bd; 43 PetscInt dim, cdim; 44 PetscErrorCode ierr; 45 46 PetscFunctionBegin; 47 if (dm == dmNew) { 48 ierr = DMDestroy(ndm);CHKERRQ(ierr); 49 PetscFunctionReturn(0); 50 } 51 dm->setupcalled = dmNew->setupcalled; 52 ierr = DMGetDimension(dmNew, &dim);CHKERRQ(ierr); 53 ierr = DMSetDimension(dm, dim);CHKERRQ(ierr); 54 ierr = DMGetCoordinateDim(dmNew, &cdim);CHKERRQ(ierr); 55 ierr = DMSetCoordinateDim(dm, cdim);CHKERRQ(ierr); 56 ierr = DMGetPointSF(dmNew, &sf);CHKERRQ(ierr); 57 ierr = DMSetPointSF(dm, sf);CHKERRQ(ierr); 58 ierr = DMGetCoordinateDM(dmNew, &coordDM);CHKERRQ(ierr); 59 ierr = DMGetCoordinatesLocal(dmNew, &coords);CHKERRQ(ierr); 60 ierr = DMSetCoordinateDM(dm, coordDM);CHKERRQ(ierr); 61 ierr = DMSetCoordinatesLocal(dm, coords);CHKERRQ(ierr); 62 /* Do not want to create the coordinate field if it does not already exist, so do not call DMGetCoordinateField() */ 63 ierr = DMFieldDestroy(&dm->coordinateField);CHKERRQ(ierr); 64 dm->coordinateField = dmNew->coordinateField; 65 ((DM_Plex *) dmNew->data)->coordFunc = ((DM_Plex *) dm->data)->coordFunc; 66 ierr = DMGetPeriodicity(dmNew, &isper, &maxCell, &L, &bd);CHKERRQ(ierr); 67 ierr = DMSetPeriodicity(dm, isper, maxCell, L, bd);CHKERRQ(ierr); 68 ierr = DMDestroy_Plex(dm);CHKERRQ(ierr); 69 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 70 dm->data = dmNew->data; 71 ((DM_Plex *) dmNew->data)->refct++; 72 ierr = DMDestroyLabelLinkList_Internal(dm);CHKERRQ(ierr); 73 ierr = DMCopyLabels(dmNew, dm, PETSC_OWN_POINTER, PETSC_TRUE, DM_COPY_LABELS_FAIL);CHKERRQ(ierr); 74 ierr = DMGetCoarseDM(dmNew,&coarseDM);CHKERRQ(ierr); 75 ierr = DMSetCoarseDM(dm,coarseDM);CHKERRQ(ierr); 76 ierr = DMDestroy(ndm);CHKERRQ(ierr); 77 PetscFunctionReturn(0); 78 } 79 80 /* Swap dm with the contents of dmNew 81 - Swap the DM_Plex structure 82 - Swap the coordinates 83 - Swap the point PetscSF 84 */ 85 static PetscErrorCode DMPlexSwap_Static(DM dmA, DM dmB) 86 { 87 DM coordDMA, coordDMB; 88 Vec coordsA, coordsB; 89 PetscSF sfA, sfB; 90 DMField fieldTmp; 91 void *tmp; 92 DMLabelLink listTmp; 93 DMLabel depthTmp; 94 PetscInt tmpI; 95 PetscErrorCode ierr; 96 97 PetscFunctionBegin; 98 if (dmA == dmB) PetscFunctionReturn(0); 99 ierr = DMGetPointSF(dmA, &sfA);CHKERRQ(ierr); 100 ierr = DMGetPointSF(dmB, &sfB);CHKERRQ(ierr); 101 ierr = PetscObjectReference((PetscObject) sfA);CHKERRQ(ierr); 102 ierr = DMSetPointSF(dmA, sfB);CHKERRQ(ierr); 103 ierr = DMSetPointSF(dmB, sfA);CHKERRQ(ierr); 104 ierr = PetscObjectDereference((PetscObject) sfA);CHKERRQ(ierr); 105 106 ierr = DMGetCoordinateDM(dmA, &coordDMA);CHKERRQ(ierr); 107 ierr = DMGetCoordinateDM(dmB, &coordDMB);CHKERRQ(ierr); 108 ierr = PetscObjectReference((PetscObject) coordDMA);CHKERRQ(ierr); 109 ierr = DMSetCoordinateDM(dmA, coordDMB);CHKERRQ(ierr); 110 ierr = DMSetCoordinateDM(dmB, coordDMA);CHKERRQ(ierr); 111 ierr = PetscObjectDereference((PetscObject) coordDMA);CHKERRQ(ierr); 112 113 ierr = DMGetCoordinatesLocal(dmA, &coordsA);CHKERRQ(ierr); 114 ierr = DMGetCoordinatesLocal(dmB, &coordsB);CHKERRQ(ierr); 115 ierr = PetscObjectReference((PetscObject) coordsA);CHKERRQ(ierr); 116 ierr = DMSetCoordinatesLocal(dmA, coordsB);CHKERRQ(ierr); 117 ierr = DMSetCoordinatesLocal(dmB, coordsA);CHKERRQ(ierr); 118 ierr = PetscObjectDereference((PetscObject) coordsA);CHKERRQ(ierr); 119 120 fieldTmp = dmA->coordinateField; 121 dmA->coordinateField = dmB->coordinateField; 122 dmB->coordinateField = fieldTmp; 123 tmp = dmA->data; 124 dmA->data = dmB->data; 125 dmB->data = tmp; 126 listTmp = dmA->labels; 127 dmA->labels = dmB->labels; 128 dmB->labels = listTmp; 129 depthTmp = dmA->depthLabel; 130 dmA->depthLabel = dmB->depthLabel; 131 dmB->depthLabel = depthTmp; 132 depthTmp = dmA->celltypeLabel; 133 dmA->celltypeLabel = dmB->celltypeLabel; 134 dmB->celltypeLabel = depthTmp; 135 tmpI = dmA->levelup; 136 dmA->levelup = dmB->levelup; 137 dmB->levelup = tmpI; 138 PetscFunctionReturn(0); 139 } 140 141 static PetscErrorCode DMPlexInterpolateInPlace_Internal(DM dm) 142 { 143 DM idm; 144 PetscErrorCode ierr; 145 146 PetscFunctionBegin; 147 ierr = DMPlexInterpolate(dm, &idm);CHKERRQ(ierr); 148 ierr = DMPlexCopyCoordinates(dm, idm);CHKERRQ(ierr); 149 ierr = DMPlexReplace_Static(dm, &idm);CHKERRQ(ierr); 150 PetscFunctionReturn(0); 151 } 152 153 /*@C 154 DMPlexCreateCoordinateSpace - Creates a finite element space for the coordinates 155 156 Collective 157 158 Input Parameters: 159 + DM - The DM 160 . degree - The degree of the finite element 161 - coordFunc - An optional function to map new points from refinement to the surface 162 163 Level: advanced 164 165 .seealso: PetscFECreateLagrange(), DMGetCoordinateDM() 166 @*/ 167 PetscErrorCode DMPlexCreateCoordinateSpace(DM dm, PetscInt degree, PetscPointFunc coordFunc) 168 { 169 DM_Plex *mesh = (DM_Plex *) dm->data; 170 DM cdm; 171 PetscDS cds; 172 PetscFE fe; 173 PetscClassId id; 174 PetscErrorCode ierr; 175 176 PetscFunctionBegin; 177 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 178 ierr = DMGetDS(cdm, &cds);CHKERRQ(ierr); 179 ierr = PetscDSGetDiscretization(cds, 0, (PetscObject *) &fe);CHKERRQ(ierr); 180 ierr = PetscObjectGetClassId((PetscObject) fe, &id);CHKERRQ(ierr); 181 if (id != PETSCFE_CLASSID) { 182 PetscBool simplex; 183 PetscInt dim, dE, qorder; 184 185 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 186 ierr = DMGetCoordinateDim(dm, &dE);CHKERRQ(ierr); 187 ierr = DMPlexIsSimplex(dm, &simplex);CHKERRQ(ierr); 188 qorder = degree; 189 ierr = PetscObjectOptionsBegin((PetscObject) cdm);CHKERRQ(ierr); 190 ierr = PetscOptionsBoundedInt("-coord_dm_default_quadrature_order", "Quadrature order is one less than quadrature points per edge", "DMPlexCreateCoordinateSpace", qorder, &qorder, NULL, 0);CHKERRQ(ierr); 191 ierr = PetscOptionsEnd();CHKERRQ(ierr); 192 ierr = PetscFECreateLagrange(PETSC_COMM_SELF, dim, dE, simplex, degree, qorder, &fe);CHKERRQ(ierr); 193 ierr = DMSetField(cdm, 0, NULL, (PetscObject) fe);CHKERRQ(ierr); 194 ierr = DMCreateDS(cdm);CHKERRQ(ierr); 195 ierr = DMProjectCoordinates(dm, fe);CHKERRQ(ierr); 196 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 197 } 198 mesh->coordFunc = coordFunc; 199 PetscFunctionReturn(0); 200 } 201 202 /*@ 203 DMPlexCreateDoublet - Creates a mesh of two cells of the specified type, optionally with later refinement. 204 205 Collective 206 207 Input Parameters: 208 + comm - The communicator for the DM object 209 . dim - The spatial dimension 210 . simplex - Flag for simplicial cells, otherwise they are tensor product cells 211 . interpolate - Flag to create intermediate mesh pieces (edges, faces) 212 - refinementLimit - A nonzero number indicates the largest admissible volume for a refined cell 213 214 Output Parameter: 215 . dm - The DM object 216 217 Level: beginner 218 219 .seealso: DMSetType(), DMCreate() 220 @*/ 221 PetscErrorCode DMPlexCreateDoublet(MPI_Comm comm, PetscInt dim, PetscBool simplex, PetscBool interpolate, PetscReal refinementLimit, DM *newdm) 222 { 223 DM dm; 224 PetscMPIInt rank; 225 PetscErrorCode ierr; 226 227 PetscFunctionBegin; 228 ierr = DMCreate(comm, &dm);CHKERRQ(ierr); 229 ierr = DMSetType(dm, DMPLEX);CHKERRQ(ierr); 230 ierr = DMSetDimension(dm, dim);CHKERRQ(ierr); 231 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 232 switch (dim) { 233 case 2: 234 if (simplex) {ierr = PetscObjectSetName((PetscObject) dm, "triangular");CHKERRQ(ierr);} 235 else {ierr = PetscObjectSetName((PetscObject) dm, "quadrilateral");CHKERRQ(ierr);} 236 break; 237 case 3: 238 if (simplex) {ierr = PetscObjectSetName((PetscObject) dm, "tetrahedral");CHKERRQ(ierr);} 239 else {ierr = PetscObjectSetName((PetscObject) dm, "hexahedral");CHKERRQ(ierr);} 240 break; 241 default: 242 SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %D", dim); 243 } 244 if (rank) { 245 PetscInt numPoints[2] = {0, 0}; 246 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 247 } else { 248 switch (dim) { 249 case 2: 250 if (simplex) { 251 PetscInt numPoints[2] = {4, 2}; 252 PetscInt coneSize[6] = {3, 3, 0, 0, 0, 0}; 253 PetscInt cones[6] = {2, 3, 4, 5, 4, 3}; 254 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 255 PetscScalar vertexCoords[8] = {-0.5, 0.5, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5}; 256 257 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 258 } else { 259 PetscInt numPoints[2] = {6, 2}; 260 PetscInt coneSize[8] = {4, 4, 0, 0, 0, 0, 0, 0}; 261 PetscInt cones[8] = {2, 3, 4, 5, 3, 6, 7, 4}; 262 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 263 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}; 264 265 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 266 } 267 break; 268 case 3: 269 if (simplex) { 270 PetscInt numPoints[2] = {5, 2}; 271 PetscInt coneSize[7] = {4, 4, 0, 0, 0, 0, 0}; 272 PetscInt cones[8] = {4, 3, 5, 2, 5, 3, 4, 6}; 273 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 274 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}; 275 276 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 277 } else { 278 PetscInt numPoints[2] = {12, 2}; 279 PetscInt coneSize[14] = {8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 280 PetscInt cones[16] = {2, 3, 4, 5, 6, 7, 8, 9, 5, 4, 10, 11, 7, 12, 13, 8}; 281 PetscInt coneOrientations[16] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 282 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, 283 -1.0, -0.5, 0.5, 0.0, -0.5, 0.5, 0.0, 0.5, 0.5, -1.0, 0.5, 0.5, 284 1.0, 0.5, -0.5, 1.0, -0.5, -0.5, 1.0, -0.5, 0.5, 1.0, 0.5, 0.5}; 285 286 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 287 } 288 break; 289 default: 290 SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %D", dim); 291 } 292 } 293 *newdm = dm; 294 if (refinementLimit > 0.0) { 295 DM rdm; 296 const char *name; 297 298 ierr = DMPlexSetRefinementUniform(*newdm, PETSC_FALSE);CHKERRQ(ierr); 299 ierr = DMPlexSetRefinementLimit(*newdm, refinementLimit);CHKERRQ(ierr); 300 ierr = DMRefine(*newdm, comm, &rdm);CHKERRQ(ierr); 301 ierr = PetscObjectGetName((PetscObject) *newdm, &name);CHKERRQ(ierr); 302 ierr = PetscObjectSetName((PetscObject) rdm, name);CHKERRQ(ierr); 303 ierr = DMDestroy(newdm);CHKERRQ(ierr); 304 *newdm = rdm; 305 } 306 if (interpolate) { 307 DM idm; 308 309 ierr = DMPlexInterpolate(*newdm, &idm);CHKERRQ(ierr); 310 ierr = DMDestroy(newdm);CHKERRQ(ierr); 311 *newdm = idm; 312 } 313 PetscFunctionReturn(0); 314 } 315 316 static PetscErrorCode DMPlexCreateBoxSurfaceMesh_Tensor_1D_Internal(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[]) 317 { 318 const PetscInt numVertices = 2; 319 PetscInt markerRight = 1; 320 PetscInt markerLeft = 1; 321 PetscBool markerSeparate = PETSC_FALSE; 322 Vec coordinates; 323 PetscSection coordSection; 324 PetscScalar *coords; 325 PetscInt coordSize; 326 PetscMPIInt rank; 327 PetscInt cdim = 1, v; 328 PetscErrorCode ierr; 329 330 PetscFunctionBegin; 331 ierr = PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);CHKERRQ(ierr); 332 if (markerSeparate) { 333 markerRight = 2; 334 markerLeft = 1; 335 } 336 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRMPI(ierr); 337 if (!rank) { 338 ierr = DMPlexSetChart(dm, 0, numVertices);CHKERRQ(ierr); 339 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 340 ierr = DMSetLabelValue(dm, "marker", 0, markerLeft);CHKERRQ(ierr); 341 ierr = DMSetLabelValue(dm, "marker", 1, markerRight);CHKERRQ(ierr); 342 } 343 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 344 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 345 /* Build coordinates */ 346 ierr = DMSetCoordinateDim(dm, cdim);CHKERRQ(ierr); 347 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 348 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 349 ierr = PetscSectionSetChart(coordSection, 0, numVertices);CHKERRQ(ierr); 350 ierr = PetscSectionSetFieldComponents(coordSection, 0, cdim);CHKERRQ(ierr); 351 for (v = 0; v < numVertices; ++v) { 352 ierr = PetscSectionSetDof(coordSection, v, cdim);CHKERRQ(ierr); 353 ierr = PetscSectionSetFieldDof(coordSection, v, 0, cdim);CHKERRQ(ierr); 354 } 355 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 356 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 357 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 358 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 359 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 360 ierr = VecSetBlockSize(coordinates, cdim);CHKERRQ(ierr); 361 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 362 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 363 coords[0] = lower[0]; 364 coords[1] = upper[0]; 365 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 366 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 367 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 368 PetscFunctionReturn(0); 369 } 370 371 static PetscErrorCode DMPlexCreateBoxSurfaceMesh_Tensor_2D_Internal(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[]) 372 { 373 const PetscInt numVertices = (edges[0]+1)*(edges[1]+1); 374 const PetscInt numEdges = edges[0]*(edges[1]+1) + (edges[0]+1)*edges[1]; 375 PetscInt markerTop = 1; 376 PetscInt markerBottom = 1; 377 PetscInt markerRight = 1; 378 PetscInt markerLeft = 1; 379 PetscBool markerSeparate = PETSC_FALSE; 380 Vec coordinates; 381 PetscSection coordSection; 382 PetscScalar *coords; 383 PetscInt coordSize; 384 PetscMPIInt rank; 385 PetscInt v, vx, vy; 386 PetscErrorCode ierr; 387 388 PetscFunctionBegin; 389 ierr = PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);CHKERRQ(ierr); 390 if (markerSeparate) { 391 markerTop = 3; 392 markerBottom = 1; 393 markerRight = 2; 394 markerLeft = 4; 395 } 396 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRMPI(ierr); 397 if (rank == 0) { 398 PetscInt e, ex, ey; 399 400 ierr = DMPlexSetChart(dm, 0, numEdges+numVertices);CHKERRQ(ierr); 401 for (e = 0; e < numEdges; ++e) { 402 ierr = DMPlexSetConeSize(dm, e, 2);CHKERRQ(ierr); 403 } 404 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 405 for (vx = 0; vx <= edges[0]; vx++) { 406 for (ey = 0; ey < edges[1]; ey++) { 407 PetscInt edge = vx*edges[1] + ey + edges[0]*(edges[1]+1); 408 PetscInt vertex = ey*(edges[0]+1) + vx + numEdges; 409 PetscInt cone[2]; 410 411 cone[0] = vertex; cone[1] = vertex+edges[0]+1; 412 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 413 if (vx == edges[0]) { 414 ierr = DMSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 415 ierr = DMSetLabelValue(dm, "marker", cone[0], markerRight);CHKERRQ(ierr); 416 if (ey == edges[1]-1) { 417 ierr = DMSetLabelValue(dm, "marker", cone[1], markerRight);CHKERRQ(ierr); 418 ierr = DMSetLabelValue(dm, "Face Sets", cone[1], markerRight);CHKERRQ(ierr); 419 } 420 } else if (vx == 0) { 421 ierr = DMSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 422 ierr = DMSetLabelValue(dm, "marker", cone[0], markerLeft);CHKERRQ(ierr); 423 if (ey == edges[1]-1) { 424 ierr = DMSetLabelValue(dm, "marker", cone[1], markerLeft);CHKERRQ(ierr); 425 ierr = DMSetLabelValue(dm, "Face Sets", cone[1], markerLeft);CHKERRQ(ierr); 426 } 427 } 428 } 429 } 430 for (vy = 0; vy <= edges[1]; vy++) { 431 for (ex = 0; ex < edges[0]; ex++) { 432 PetscInt edge = vy*edges[0] + ex; 433 PetscInt vertex = vy*(edges[0]+1) + ex + numEdges; 434 PetscInt cone[2]; 435 436 cone[0] = vertex; cone[1] = vertex+1; 437 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 438 if (vy == edges[1]) { 439 ierr = DMSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 440 ierr = DMSetLabelValue(dm, "marker", cone[0], markerTop);CHKERRQ(ierr); 441 if (ex == edges[0]-1) { 442 ierr = DMSetLabelValue(dm, "marker", cone[1], markerTop);CHKERRQ(ierr); 443 ierr = DMSetLabelValue(dm, "Face Sets", cone[1], markerTop);CHKERRQ(ierr); 444 } 445 } else if (vy == 0) { 446 ierr = DMSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 447 ierr = DMSetLabelValue(dm, "marker", cone[0], markerBottom);CHKERRQ(ierr); 448 if (ex == edges[0]-1) { 449 ierr = DMSetLabelValue(dm, "marker", cone[1], markerBottom);CHKERRQ(ierr); 450 ierr = DMSetLabelValue(dm, "Face Sets", cone[1], markerBottom);CHKERRQ(ierr); 451 } 452 } 453 } 454 } 455 } 456 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 457 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 458 /* Build coordinates */ 459 ierr = DMSetCoordinateDim(dm, 2);CHKERRQ(ierr); 460 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 461 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 462 ierr = PetscSectionSetChart(coordSection, numEdges, numEdges + numVertices);CHKERRQ(ierr); 463 ierr = PetscSectionSetFieldComponents(coordSection, 0, 2);CHKERRQ(ierr); 464 for (v = numEdges; v < numEdges+numVertices; ++v) { 465 ierr = PetscSectionSetDof(coordSection, v, 2);CHKERRQ(ierr); 466 ierr = PetscSectionSetFieldDof(coordSection, v, 0, 2);CHKERRQ(ierr); 467 } 468 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 469 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 470 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 471 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 472 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 473 ierr = VecSetBlockSize(coordinates, 2);CHKERRQ(ierr); 474 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 475 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 476 for (vy = 0; vy <= edges[1]; ++vy) { 477 for (vx = 0; vx <= edges[0]; ++vx) { 478 coords[(vy*(edges[0]+1)+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/edges[0])*vx; 479 coords[(vy*(edges[0]+1)+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/edges[1])*vy; 480 } 481 } 482 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 483 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 484 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 485 PetscFunctionReturn(0); 486 } 487 488 static PetscErrorCode DMPlexCreateBoxSurfaceMesh_Tensor_3D_Internal(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt faces[]) 489 { 490 PetscInt vertices[3], numVertices; 491 PetscInt numFaces = 2*faces[0]*faces[1] + 2*faces[1]*faces[2] + 2*faces[0]*faces[2]; 492 Vec coordinates; 493 PetscSection coordSection; 494 PetscScalar *coords; 495 PetscInt coordSize; 496 PetscMPIInt rank; 497 PetscInt v, vx, vy, vz; 498 PetscInt voffset, iface=0, cone[4]; 499 PetscErrorCode ierr; 500 501 PetscFunctionBegin; 502 PetscCheckFalse((faces[0] < 1) || (faces[1] < 1) || (faces[2] < 1),PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Must have at least 1 face per side"); 503 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRMPI(ierr); 504 vertices[0] = faces[0]+1; vertices[1] = faces[1]+1; vertices[2] = faces[2]+1; 505 numVertices = vertices[0]*vertices[1]*vertices[2]; 506 if (rank == 0) { 507 PetscInt f; 508 509 ierr = DMPlexSetChart(dm, 0, numFaces+numVertices);CHKERRQ(ierr); 510 for (f = 0; f < numFaces; ++f) { 511 ierr = DMPlexSetConeSize(dm, f, 4);CHKERRQ(ierr); 512 } 513 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 514 515 /* Side 0 (Top) */ 516 for (vy = 0; vy < faces[1]; vy++) { 517 for (vx = 0; vx < faces[0]; vx++) { 518 voffset = numFaces + vertices[0]*vertices[1]*(vertices[2]-1) + vy*vertices[0] + vx; 519 cone[0] = voffset; cone[1] = voffset+1; cone[2] = voffset+vertices[0]+1; cone[3] = voffset+vertices[0]; 520 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 521 ierr = DMSetLabelValue(dm, "marker", iface, 1);CHKERRQ(ierr); 522 ierr = DMSetLabelValue(dm, "marker", voffset+0, 1);CHKERRQ(ierr); 523 ierr = DMSetLabelValue(dm, "marker", voffset+1, 1);CHKERRQ(ierr); 524 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]+0, 1);CHKERRQ(ierr); 525 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]+1, 1);CHKERRQ(ierr); 526 iface++; 527 } 528 } 529 530 /* Side 1 (Bottom) */ 531 for (vy = 0; vy < faces[1]; vy++) { 532 for (vx = 0; vx < faces[0]; vx++) { 533 voffset = numFaces + vy*(faces[0]+1) + vx; 534 cone[0] = voffset+1; cone[1] = voffset; cone[2] = voffset+vertices[0]; cone[3] = voffset+vertices[0]+1; 535 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 536 ierr = DMSetLabelValue(dm, "marker", iface, 1);CHKERRQ(ierr); 537 ierr = DMSetLabelValue(dm, "marker", voffset+0, 1);CHKERRQ(ierr); 538 ierr = DMSetLabelValue(dm, "marker", voffset+1, 1);CHKERRQ(ierr); 539 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]+0, 1);CHKERRQ(ierr); 540 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]+1, 1);CHKERRQ(ierr); 541 iface++; 542 } 543 } 544 545 /* Side 2 (Front) */ 546 for (vz = 0; vz < faces[2]; vz++) { 547 for (vx = 0; vx < faces[0]; vx++) { 548 voffset = numFaces + vz*vertices[0]*vertices[1] + vx; 549 cone[0] = voffset; cone[1] = voffset+1; cone[2] = voffset+vertices[0]*vertices[1]+1; cone[3] = voffset+vertices[0]*vertices[1]; 550 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 551 ierr = DMSetLabelValue(dm, "marker", iface, 1);CHKERRQ(ierr); 552 ierr = DMSetLabelValue(dm, "marker", voffset+0, 1);CHKERRQ(ierr); 553 ierr = DMSetLabelValue(dm, "marker", voffset+1, 1);CHKERRQ(ierr); 554 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+0, 1);CHKERRQ(ierr); 555 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+1, 1);CHKERRQ(ierr); 556 iface++; 557 } 558 } 559 560 /* Side 3 (Back) */ 561 for (vz = 0; vz < faces[2]; vz++) { 562 for (vx = 0; vx < faces[0]; vx++) { 563 voffset = numFaces + vz*vertices[0]*vertices[1] + vertices[0]*(vertices[1]-1) + vx; 564 cone[0] = voffset+vertices[0]*vertices[1]; cone[1] = voffset+vertices[0]*vertices[1]+1; 565 cone[2] = voffset+1; cone[3] = voffset; 566 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 567 ierr = DMSetLabelValue(dm, "marker", iface, 1);CHKERRQ(ierr); 568 ierr = DMSetLabelValue(dm, "marker", voffset+0, 1);CHKERRQ(ierr); 569 ierr = DMSetLabelValue(dm, "marker", voffset+1, 1);CHKERRQ(ierr); 570 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+0, 1);CHKERRQ(ierr); 571 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+1, 1);CHKERRQ(ierr); 572 iface++; 573 } 574 } 575 576 /* Side 4 (Left) */ 577 for (vz = 0; vz < faces[2]; vz++) { 578 for (vy = 0; vy < faces[1]; vy++) { 579 voffset = numFaces + vz*vertices[0]*vertices[1] + vy*vertices[0]; 580 cone[0] = voffset; cone[1] = voffset+vertices[0]*vertices[1]; 581 cone[2] = voffset+vertices[0]*vertices[1]+vertices[0]; cone[3] = voffset+vertices[0]; 582 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 583 ierr = DMSetLabelValue(dm, "marker", iface, 1);CHKERRQ(ierr); 584 ierr = DMSetLabelValue(dm, "marker", voffset+0, 1);CHKERRQ(ierr); 585 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]+0, 1);CHKERRQ(ierr); 586 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[1]+0, 1);CHKERRQ(ierr); 587 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+vertices[0], 1);CHKERRQ(ierr); 588 iface++; 589 } 590 } 591 592 /* Side 5 (Right) */ 593 for (vz = 0; vz < faces[2]; vz++) { 594 for (vy = 0; vy < faces[1]; vy++) { 595 voffset = numFaces + vz*vertices[0]*vertices[1] + vy*vertices[0] + faces[0]; 596 cone[0] = voffset+vertices[0]*vertices[1]; cone[1] = voffset; 597 cone[2] = voffset+vertices[0]; cone[3] = voffset+vertices[0]*vertices[1]+vertices[0]; 598 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr); 599 ierr = DMSetLabelValue(dm, "marker", iface, 1);CHKERRQ(ierr); 600 ierr = DMSetLabelValue(dm, "marker", voffset+0, 1);CHKERRQ(ierr); 601 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]+0, 1);CHKERRQ(ierr); 602 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+0, 1);CHKERRQ(ierr); 603 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+vertices[0], 1);CHKERRQ(ierr); 604 iface++; 605 } 606 } 607 } 608 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 609 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 610 /* Build coordinates */ 611 ierr = DMSetCoordinateDim(dm, 3);CHKERRQ(ierr); 612 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 613 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 614 ierr = PetscSectionSetChart(coordSection, numFaces, numFaces + numVertices);CHKERRQ(ierr); 615 ierr = PetscSectionSetFieldComponents(coordSection, 0, 3);CHKERRQ(ierr); 616 for (v = numFaces; v < numFaces+numVertices; ++v) { 617 ierr = PetscSectionSetDof(coordSection, v, 3);CHKERRQ(ierr); 618 ierr = PetscSectionSetFieldDof(coordSection, v, 0, 3);CHKERRQ(ierr); 619 } 620 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 621 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 622 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 623 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 624 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 625 ierr = VecSetBlockSize(coordinates, 3);CHKERRQ(ierr); 626 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 627 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 628 for (vz = 0; vz <= faces[2]; ++vz) { 629 for (vy = 0; vy <= faces[1]; ++vy) { 630 for (vx = 0; vx <= faces[0]; ++vx) { 631 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+0] = lower[0] + ((upper[0] - lower[0])/faces[0])*vx; 632 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+1] = lower[1] + ((upper[1] - lower[1])/faces[1])*vy; 633 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+2] = lower[2] + ((upper[2] - lower[2])/faces[2])*vz; 634 } 635 } 636 } 637 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 638 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 639 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 640 PetscFunctionReturn(0); 641 } 642 643 static PetscErrorCode DMPlexCreateBoxSurfaceMesh_Internal(DM dm, PetscInt dim, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], PetscBool interpolate) 644 { 645 PetscErrorCode ierr; 646 647 PetscFunctionBegin; 648 PetscValidLogicalCollectiveInt(dm, dim, 2); 649 ierr = DMSetDimension(dm, dim-1);CHKERRQ(ierr); 650 ierr = DMSetCoordinateDim(dm, dim);CHKERRQ(ierr); 651 switch (dim) { 652 case 1: ierr = DMPlexCreateBoxSurfaceMesh_Tensor_1D_Internal(dm, lower, upper, faces);CHKERRQ(ierr);break; 653 case 2: ierr = DMPlexCreateBoxSurfaceMesh_Tensor_2D_Internal(dm, lower, upper, faces);CHKERRQ(ierr);break; 654 case 3: ierr = DMPlexCreateBoxSurfaceMesh_Tensor_3D_Internal(dm, lower, upper, faces);CHKERRQ(ierr);break; 655 default: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Dimension not supported: %D", dim); 656 } 657 if (interpolate) {ierr = DMPlexInterpolateInPlace_Internal(dm);CHKERRQ(ierr);} 658 PetscFunctionReturn(0); 659 } 660 661 /*@C 662 DMPlexCreateBoxSurfaceMesh - Creates a mesh on the surface of the tensor product of unit intervals (box) using tensor cells (hexahedra). 663 664 Collective 665 666 Input Parameters: 667 + comm - The communicator for the DM object 668 . dim - The spatial dimension of the box, so the resulting mesh is has dimension dim-1 669 . faces - Number of faces per dimension, or NULL for (1,) in 1D and (2, 2) in 2D and (1, 1, 1) in 3D 670 . lower - The lower left corner, or NULL for (0, 0, 0) 671 . upper - The upper right corner, or NULL for (1, 1, 1) 672 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 673 674 Output Parameter: 675 . dm - The DM object 676 677 Level: beginner 678 679 .seealso: DMSetFromOptions(), DMPlexCreateBoxMesh(), DMPlexCreateFromFile(), DMSetType(), DMCreate() 680 @*/ 681 PetscErrorCode DMPlexCreateBoxSurfaceMesh(MPI_Comm comm, PetscInt dim, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], PetscBool interpolate, DM *dm) 682 { 683 PetscInt fac[3] = {1, 1, 1}; 684 PetscReal low[3] = {0, 0, 0}; 685 PetscReal upp[3] = {1, 1, 1}; 686 PetscErrorCode ierr; 687 688 PetscFunctionBegin; 689 ierr = DMCreate(comm,dm);CHKERRQ(ierr); 690 ierr = DMSetType(*dm,DMPLEX);CHKERRQ(ierr); 691 ierr = DMPlexCreateBoxSurfaceMesh_Internal(*dm, dim, faces ? faces : fac, lower ? lower : low, upper ? upper : upp, interpolate);CHKERRQ(ierr); 692 PetscFunctionReturn(0); 693 } 694 695 static PetscErrorCode DMPlexCreateLineMesh_Internal(DM dm,PetscInt segments,PetscReal lower,PetscReal upper,DMBoundaryType bd) 696 { 697 PetscInt i,fStart,fEnd,numCells = 0,numVerts = 0; 698 PetscInt numPoints[2],*coneSize,*cones,*coneOrientations; 699 PetscScalar *vertexCoords; 700 PetscReal L,maxCell; 701 PetscBool markerSeparate = PETSC_FALSE; 702 PetscInt markerLeft = 1, faceMarkerLeft = 1; 703 PetscInt markerRight = 1, faceMarkerRight = 2; 704 PetscBool wrap = (bd == DM_BOUNDARY_PERIODIC || bd == DM_BOUNDARY_TWIST) ? PETSC_TRUE : PETSC_FALSE; 705 PetscMPIInt rank; 706 PetscErrorCode ierr; 707 708 PetscFunctionBegin; 709 PetscValidPointer(dm,1); 710 711 ierr = DMSetDimension(dm,1);CHKERRQ(ierr); 712 ierr = DMCreateLabel(dm,"marker");CHKERRQ(ierr); 713 ierr = DMCreateLabel(dm,"Face Sets");CHKERRQ(ierr); 714 715 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm),&rank);CHKERRMPI(ierr); 716 if (rank == 0) numCells = segments; 717 if (rank == 0) numVerts = segments + (wrap ? 0 : 1); 718 719 numPoints[0] = numVerts ; numPoints[1] = numCells; 720 ierr = PetscMalloc4(numCells+numVerts,&coneSize,numCells*2,&cones,numCells+numVerts,&coneOrientations,numVerts,&vertexCoords);CHKERRQ(ierr); 721 ierr = PetscArrayzero(coneOrientations,numCells+numVerts);CHKERRQ(ierr); 722 for (i = 0; i < numCells; ++i) { coneSize[i] = 2; } 723 for (i = 0; i < numVerts; ++i) { coneSize[numCells+i] = 0; } 724 for (i = 0; i < numCells; ++i) { cones[2*i] = numCells + i%numVerts; cones[2*i+1] = numCells + (i+1)%numVerts; } 725 for (i = 0; i < numVerts; ++i) { vertexCoords[i] = lower + (upper-lower)*((PetscReal)i/(PetscReal)numCells); } 726 ierr = DMPlexCreateFromDAG(dm,1,numPoints,coneSize,cones,coneOrientations,vertexCoords);CHKERRQ(ierr); 727 ierr = PetscFree4(coneSize,cones,coneOrientations,vertexCoords);CHKERRQ(ierr); 728 729 ierr = PetscOptionsGetBool(((PetscObject)dm)->options,((PetscObject)dm)->prefix,"-dm_plex_separate_marker",&markerSeparate,NULL);CHKERRQ(ierr); 730 if (markerSeparate) { markerLeft = faceMarkerLeft; markerRight = faceMarkerRight;} 731 if (!wrap && rank == 0) { 732 ierr = DMPlexGetHeightStratum(dm,1,&fStart,&fEnd);CHKERRQ(ierr); 733 ierr = DMSetLabelValue(dm,"marker",fStart,markerLeft);CHKERRQ(ierr); 734 ierr = DMSetLabelValue(dm,"marker",fEnd-1,markerRight);CHKERRQ(ierr); 735 ierr = DMSetLabelValue(dm,"Face Sets",fStart,faceMarkerLeft);CHKERRQ(ierr); 736 ierr = DMSetLabelValue(dm,"Face Sets",fEnd-1,faceMarkerRight);CHKERRQ(ierr); 737 } 738 if (wrap) { 739 L = upper - lower; 740 maxCell = (PetscReal)1.1*(L/(PetscReal)PetscMax(1,segments)); 741 ierr = DMSetPeriodicity(dm,PETSC_TRUE,&maxCell,&L,&bd);CHKERRQ(ierr); 742 } 743 ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE);CHKERRQ(ierr); 744 PetscFunctionReturn(0); 745 } 746 747 static PetscErrorCode DMPlexCreateBoxMesh_Simplex_Internal(DM dm, PetscInt dim, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool interpolate) 748 { 749 DM boundary, vol; 750 PetscInt i; 751 PetscErrorCode ierr; 752 753 PetscFunctionBegin; 754 PetscValidPointer(dm, 1); 755 for (i = 0; i < dim; ++i) PetscCheckFalse(periodicity[i] != DM_BOUNDARY_NONE,PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Periodicity is not supported for simplex meshes"); 756 ierr = DMCreate(PetscObjectComm((PetscObject) dm), &boundary);CHKERRQ(ierr); 757 ierr = DMSetType(boundary, DMPLEX);CHKERRQ(ierr); 758 ierr = DMPlexCreateBoxSurfaceMesh_Internal(boundary, dim, faces, lower, upper, PETSC_FALSE);CHKERRQ(ierr); 759 ierr = DMPlexGenerate(boundary, NULL, interpolate, &vol);CHKERRQ(ierr); 760 ierr = DMPlexCopy_Internal(dm, PETSC_TRUE, vol);CHKERRQ(ierr); 761 ierr = DMPlexReplace_Static(dm, &vol);CHKERRQ(ierr); 762 ierr = DMDestroy(&boundary);CHKERRQ(ierr); 763 PetscFunctionReturn(0); 764 } 765 766 static PetscErrorCode DMPlexCreateCubeMesh_Internal(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[], DMBoundaryType bdX, DMBoundaryType bdY, DMBoundaryType bdZ) 767 { 768 DMLabel cutLabel = NULL; 769 PetscInt markerTop = 1, faceMarkerTop = 1; 770 PetscInt markerBottom = 1, faceMarkerBottom = 1; 771 PetscInt markerFront = 1, faceMarkerFront = 1; 772 PetscInt markerBack = 1, faceMarkerBack = 1; 773 PetscInt markerRight = 1, faceMarkerRight = 1; 774 PetscInt markerLeft = 1, faceMarkerLeft = 1; 775 PetscInt dim; 776 PetscBool markerSeparate = PETSC_FALSE, cutMarker = PETSC_FALSE; 777 PetscMPIInt rank; 778 PetscErrorCode ierr; 779 780 PetscFunctionBegin; 781 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 782 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRMPI(ierr); 783 ierr = DMCreateLabel(dm,"marker");CHKERRQ(ierr); 784 ierr = DMCreateLabel(dm,"Face Sets");CHKERRQ(ierr); 785 ierr = PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_periodic_cut", &cutMarker, NULL);CHKERRQ(ierr); 786 if (bdX == DM_BOUNDARY_PERIODIC || bdX == DM_BOUNDARY_TWIST || 787 bdY == DM_BOUNDARY_PERIODIC || bdY == DM_BOUNDARY_TWIST || 788 bdZ == DM_BOUNDARY_PERIODIC || bdZ == DM_BOUNDARY_TWIST) { 789 790 if (cutMarker) {ierr = DMCreateLabel(dm, "periodic_cut");CHKERRQ(ierr); ierr = DMGetLabel(dm, "periodic_cut", &cutLabel);CHKERRQ(ierr);} 791 } 792 switch (dim) { 793 case 2: 794 faceMarkerTop = 3; 795 faceMarkerBottom = 1; 796 faceMarkerRight = 2; 797 faceMarkerLeft = 4; 798 break; 799 case 3: 800 faceMarkerBottom = 1; 801 faceMarkerTop = 2; 802 faceMarkerFront = 3; 803 faceMarkerBack = 4; 804 faceMarkerRight = 5; 805 faceMarkerLeft = 6; 806 break; 807 default: 808 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Dimension %D not supported",dim); 809 } 810 ierr = PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);CHKERRQ(ierr); 811 if (markerSeparate) { 812 markerBottom = faceMarkerBottom; 813 markerTop = faceMarkerTop; 814 markerFront = faceMarkerFront; 815 markerBack = faceMarkerBack; 816 markerRight = faceMarkerRight; 817 markerLeft = faceMarkerLeft; 818 } 819 { 820 const PetscInt numXEdges = rank == 0 ? edges[0] : 0; 821 const PetscInt numYEdges = rank == 0 ? edges[1] : 0; 822 const PetscInt numZEdges = rank == 0 ? edges[2] : 0; 823 const PetscInt numXVertices = rank == 0 ? (bdX == DM_BOUNDARY_PERIODIC || bdX == DM_BOUNDARY_TWIST ? edges[0] : edges[0]+1) : 0; 824 const PetscInt numYVertices = rank == 0 ? (bdY == DM_BOUNDARY_PERIODIC || bdY == DM_BOUNDARY_TWIST ? edges[1] : edges[1]+1) : 0; 825 const PetscInt numZVertices = rank == 0 ? (bdZ == DM_BOUNDARY_PERIODIC || bdZ == DM_BOUNDARY_TWIST ? edges[2] : edges[2]+1) : 0; 826 const PetscInt numCells = numXEdges*numYEdges*numZEdges; 827 const PetscInt numXFaces = numYEdges*numZEdges; 828 const PetscInt numYFaces = numXEdges*numZEdges; 829 const PetscInt numZFaces = numXEdges*numYEdges; 830 const PetscInt numTotXFaces = numXVertices*numXFaces; 831 const PetscInt numTotYFaces = numYVertices*numYFaces; 832 const PetscInt numTotZFaces = numZVertices*numZFaces; 833 const PetscInt numFaces = numTotXFaces + numTotYFaces + numTotZFaces; 834 const PetscInt numTotXEdges = numXEdges*numYVertices*numZVertices; 835 const PetscInt numTotYEdges = numYEdges*numXVertices*numZVertices; 836 const PetscInt numTotZEdges = numZEdges*numXVertices*numYVertices; 837 const PetscInt numVertices = numXVertices*numYVertices*numZVertices; 838 const PetscInt numEdges = numTotXEdges + numTotYEdges + numTotZEdges; 839 const PetscInt firstVertex = (dim == 2) ? numFaces : numCells; 840 const PetscInt firstXFace = (dim == 2) ? 0 : numCells + numVertices; 841 const PetscInt firstYFace = firstXFace + numTotXFaces; 842 const PetscInt firstZFace = firstYFace + numTotYFaces; 843 const PetscInt firstXEdge = numCells + numFaces + numVertices; 844 const PetscInt firstYEdge = firstXEdge + numTotXEdges; 845 const PetscInt firstZEdge = firstYEdge + numTotYEdges; 846 Vec coordinates; 847 PetscSection coordSection; 848 PetscScalar *coords; 849 PetscInt coordSize; 850 PetscInt v, vx, vy, vz; 851 PetscInt c, f, fx, fy, fz, e, ex, ey, ez; 852 853 ierr = DMPlexSetChart(dm, 0, numCells+numFaces+numEdges+numVertices);CHKERRQ(ierr); 854 for (c = 0; c < numCells; c++) { 855 ierr = DMPlexSetConeSize(dm, c, 6);CHKERRQ(ierr); 856 } 857 for (f = firstXFace; f < firstXFace+numFaces; ++f) { 858 ierr = DMPlexSetConeSize(dm, f, 4);CHKERRQ(ierr); 859 } 860 for (e = firstXEdge; e < firstXEdge+numEdges; ++e) { 861 ierr = DMPlexSetConeSize(dm, e, 2);CHKERRQ(ierr); 862 } 863 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 864 /* Build cells */ 865 for (fz = 0; fz < numZEdges; ++fz) { 866 for (fy = 0; fy < numYEdges; ++fy) { 867 for (fx = 0; fx < numXEdges; ++fx) { 868 PetscInt cell = (fz*numYEdges + fy)*numXEdges + fx; 869 PetscInt faceB = firstZFace + (fy*numXEdges+fx)*numZVertices + fz; 870 PetscInt faceT = firstZFace + (fy*numXEdges+fx)*numZVertices + ((fz+1)%numZVertices); 871 PetscInt faceF = firstYFace + (fz*numXEdges+fx)*numYVertices + fy; 872 PetscInt faceK = firstYFace + (fz*numXEdges+fx)*numYVertices + ((fy+1)%numYVertices); 873 PetscInt faceL = firstXFace + (fz*numYEdges+fy)*numXVertices + fx; 874 PetscInt faceR = firstXFace + (fz*numYEdges+fy)*numXVertices + ((fx+1)%numXVertices); 875 /* B, T, F, K, R, L */ 876 PetscInt ornt[6] = {-2, 0, 0, -3, 0, -2}; /* ??? */ 877 PetscInt cone[6]; 878 879 /* no boundary twisting in 3D */ 880 cone[0] = faceB; cone[1] = faceT; cone[2] = faceF; cone[3] = faceK; cone[4] = faceR; cone[5] = faceL; 881 ierr = DMPlexSetCone(dm, cell, cone);CHKERRQ(ierr); 882 ierr = DMPlexSetConeOrientation(dm, cell, ornt);CHKERRQ(ierr); 883 if (bdX != DM_BOUNDARY_NONE && fx == numXEdges-1 && cutLabel) {ierr = DMLabelSetValue(cutLabel, cell, 2);CHKERRQ(ierr);} 884 if (bdY != DM_BOUNDARY_NONE && fy == numYEdges-1 && cutLabel) {ierr = DMLabelSetValue(cutLabel, cell, 2);CHKERRQ(ierr);} 885 if (bdZ != DM_BOUNDARY_NONE && fz == numZEdges-1 && cutLabel) {ierr = DMLabelSetValue(cutLabel, cell, 2);CHKERRQ(ierr);} 886 } 887 } 888 } 889 /* Build x faces */ 890 for (fz = 0; fz < numZEdges; ++fz) { 891 for (fy = 0; fy < numYEdges; ++fy) { 892 for (fx = 0; fx < numXVertices; ++fx) { 893 PetscInt face = firstXFace + (fz*numYEdges+fy) *numXVertices+fx; 894 PetscInt edgeL = firstZEdge + (fy *numXVertices+fx)*numZEdges + fz; 895 PetscInt edgeR = firstZEdge + (((fy+1)%numYVertices)*numXVertices+fx)*numZEdges + fz; 896 PetscInt edgeB = firstYEdge + (fz *numXVertices+fx)*numYEdges + fy; 897 PetscInt edgeT = firstYEdge + (((fz+1)%numZVertices)*numXVertices+fx)*numYEdges + fy; 898 PetscInt ornt[4] = {0, 0, -1, -1}; 899 PetscInt cone[4]; 900 901 if (dim == 3) { 902 /* markers */ 903 if (bdX != DM_BOUNDARY_PERIODIC) { 904 if (fx == numXVertices-1) { 905 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerRight);CHKERRQ(ierr); 906 ierr = DMSetLabelValue(dm, "marker", face, markerRight);CHKERRQ(ierr); 907 } 908 else if (fx == 0) { 909 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerLeft);CHKERRQ(ierr); 910 ierr = DMSetLabelValue(dm, "marker", face, markerLeft);CHKERRQ(ierr); 911 } 912 } 913 } 914 cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL; 915 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr); 916 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr); 917 } 918 } 919 } 920 /* Build y faces */ 921 for (fz = 0; fz < numZEdges; ++fz) { 922 for (fx = 0; fx < numXEdges; ++fx) { 923 for (fy = 0; fy < numYVertices; ++fy) { 924 PetscInt face = firstYFace + (fz*numXEdges+fx)*numYVertices + fy; 925 PetscInt edgeL = firstZEdge + (fy*numXVertices+ fx)*numZEdges + fz; 926 PetscInt edgeR = firstZEdge + (fy*numXVertices+((fx+1)%numXVertices))*numZEdges + fz; 927 PetscInt edgeB = firstXEdge + (fz *numYVertices+fy)*numXEdges + fx; 928 PetscInt edgeT = firstXEdge + (((fz+1)%numZVertices)*numYVertices+fy)*numXEdges + fx; 929 PetscInt ornt[4] = {0, 0, -1, -1}; 930 PetscInt cone[4]; 931 932 if (dim == 3) { 933 /* markers */ 934 if (bdY != DM_BOUNDARY_PERIODIC) { 935 if (fy == numYVertices-1) { 936 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerBack);CHKERRQ(ierr); 937 ierr = DMSetLabelValue(dm, "marker", face, markerBack);CHKERRQ(ierr); 938 } 939 else if (fy == 0) { 940 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerFront);CHKERRQ(ierr); 941 ierr = DMSetLabelValue(dm, "marker", face, markerFront);CHKERRQ(ierr); 942 } 943 } 944 } 945 cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL; 946 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr); 947 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr); 948 } 949 } 950 } 951 /* Build z faces */ 952 for (fy = 0; fy < numYEdges; ++fy) { 953 for (fx = 0; fx < numXEdges; ++fx) { 954 for (fz = 0; fz < numZVertices; fz++) { 955 PetscInt face = firstZFace + (fy*numXEdges+fx)*numZVertices + fz; 956 PetscInt edgeL = firstYEdge + (fz*numXVertices+ fx)*numYEdges + fy; 957 PetscInt edgeR = firstYEdge + (fz*numXVertices+((fx+1)%numXVertices))*numYEdges + fy; 958 PetscInt edgeB = firstXEdge + (fz*numYVertices+ fy)*numXEdges + fx; 959 PetscInt edgeT = firstXEdge + (fz*numYVertices+((fy+1)%numYVertices))*numXEdges + fx; 960 PetscInt ornt[4] = {0, 0, -1, -1}; 961 PetscInt cone[4]; 962 963 if (dim == 2) { 964 if (bdX == DM_BOUNDARY_TWIST && fx == numXEdges-1) {edgeR += numYEdges-1-2*fy; ornt[1] = -1;} 965 if (bdY == DM_BOUNDARY_TWIST && fy == numYEdges-1) {edgeT += numXEdges-1-2*fx; ornt[2] = 0;} 966 if (bdX != DM_BOUNDARY_NONE && fx == numXEdges-1 && cutLabel) {ierr = DMLabelSetValue(cutLabel, face, 2);CHKERRQ(ierr);} 967 if (bdY != DM_BOUNDARY_NONE && fy == numYEdges-1 && cutLabel) {ierr = DMLabelSetValue(cutLabel, face, 2);CHKERRQ(ierr);} 968 } else { 969 /* markers */ 970 if (bdZ != DM_BOUNDARY_PERIODIC) { 971 if (fz == numZVertices-1) { 972 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerTop);CHKERRQ(ierr); 973 ierr = DMSetLabelValue(dm, "marker", face, markerTop);CHKERRQ(ierr); 974 } 975 else if (fz == 0) { 976 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerBottom);CHKERRQ(ierr); 977 ierr = DMSetLabelValue(dm, "marker", face, markerBottom);CHKERRQ(ierr); 978 } 979 } 980 } 981 cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL; 982 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr); 983 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr); 984 } 985 } 986 } 987 /* Build Z edges*/ 988 for (vy = 0; vy < numYVertices; vy++) { 989 for (vx = 0; vx < numXVertices; vx++) { 990 for (ez = 0; ez < numZEdges; ez++) { 991 const PetscInt edge = firstZEdge + (vy*numXVertices+vx)*numZEdges + ez; 992 const PetscInt vertexB = firstVertex + (ez *numYVertices+vy)*numXVertices + vx; 993 const PetscInt vertexT = firstVertex + (((ez+1)%numZVertices)*numYVertices+vy)*numXVertices + vx; 994 PetscInt cone[2]; 995 996 if (dim == 3) { 997 if (bdX != DM_BOUNDARY_PERIODIC) { 998 if (vx == numXVertices-1) { 999 ierr = DMSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 1000 } 1001 else if (vx == 0) { 1002 ierr = DMSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 1003 } 1004 } 1005 if (bdY != DM_BOUNDARY_PERIODIC) { 1006 if (vy == numYVertices-1) { 1007 ierr = DMSetLabelValue(dm, "marker", edge, markerBack);CHKERRQ(ierr); 1008 } 1009 else if (vy == 0) { 1010 ierr = DMSetLabelValue(dm, "marker", edge, markerFront);CHKERRQ(ierr); 1011 } 1012 } 1013 } 1014 cone[0] = vertexB; cone[1] = vertexT; 1015 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 1016 } 1017 } 1018 } 1019 /* Build Y edges*/ 1020 for (vz = 0; vz < numZVertices; vz++) { 1021 for (vx = 0; vx < numXVertices; vx++) { 1022 for (ey = 0; ey < numYEdges; ey++) { 1023 const PetscInt nextv = (dim == 2 && bdY == DM_BOUNDARY_TWIST && ey == numYEdges-1) ? (numXVertices-vx-1) : (vz*numYVertices+((ey+1)%numYVertices))*numXVertices + vx; 1024 const PetscInt edge = firstYEdge + (vz*numXVertices+vx)*numYEdges + ey; 1025 const PetscInt vertexF = firstVertex + (vz*numYVertices+ey)*numXVertices + vx; 1026 const PetscInt vertexK = firstVertex + nextv; 1027 PetscInt cone[2]; 1028 1029 cone[0] = vertexF; cone[1] = vertexK; 1030 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 1031 if (dim == 2) { 1032 if ((bdX != DM_BOUNDARY_PERIODIC) && (bdX != DM_BOUNDARY_TWIST)) { 1033 if (vx == numXVertices-1) { 1034 ierr = DMSetLabelValue(dm, "Face Sets", edge, faceMarkerRight);CHKERRQ(ierr); 1035 ierr = DMSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 1036 ierr = DMSetLabelValue(dm, "marker", cone[0], markerRight);CHKERRQ(ierr); 1037 if (ey == numYEdges-1) { 1038 ierr = DMSetLabelValue(dm, "marker", cone[1], markerRight);CHKERRQ(ierr); 1039 } 1040 } else if (vx == 0) { 1041 ierr = DMSetLabelValue(dm, "Face Sets", edge, faceMarkerLeft);CHKERRQ(ierr); 1042 ierr = DMSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 1043 ierr = DMSetLabelValue(dm, "marker", cone[0], markerLeft);CHKERRQ(ierr); 1044 if (ey == numYEdges-1) { 1045 ierr = DMSetLabelValue(dm, "marker", cone[1], markerLeft);CHKERRQ(ierr); 1046 } 1047 } 1048 } else { 1049 if (vx == 0 && cutLabel) { 1050 ierr = DMLabelSetValue(cutLabel, edge, 1);CHKERRQ(ierr); 1051 ierr = DMLabelSetValue(cutLabel, cone[0], 1);CHKERRQ(ierr); 1052 if (ey == numYEdges-1) { 1053 ierr = DMLabelSetValue(cutLabel, cone[1], 1);CHKERRQ(ierr); 1054 } 1055 } 1056 } 1057 } else { 1058 if (bdX != DM_BOUNDARY_PERIODIC) { 1059 if (vx == numXVertices-1) { 1060 ierr = DMSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr); 1061 } else if (vx == 0) { 1062 ierr = DMSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr); 1063 } 1064 } 1065 if (bdZ != DM_BOUNDARY_PERIODIC) { 1066 if (vz == numZVertices-1) { 1067 ierr = DMSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 1068 } else if (vz == 0) { 1069 ierr = DMSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 1070 } 1071 } 1072 } 1073 } 1074 } 1075 } 1076 /* Build X edges*/ 1077 for (vz = 0; vz < numZVertices; vz++) { 1078 for (vy = 0; vy < numYVertices; vy++) { 1079 for (ex = 0; ex < numXEdges; ex++) { 1080 const PetscInt nextv = (dim == 2 && bdX == DM_BOUNDARY_TWIST && ex == numXEdges-1) ? (numYVertices-vy-1)*numXVertices : (vz*numYVertices+vy)*numXVertices + (ex+1)%numXVertices; 1081 const PetscInt edge = firstXEdge + (vz*numYVertices+vy)*numXEdges + ex; 1082 const PetscInt vertexL = firstVertex + (vz*numYVertices+vy)*numXVertices + ex; 1083 const PetscInt vertexR = firstVertex + nextv; 1084 PetscInt cone[2]; 1085 1086 cone[0] = vertexL; cone[1] = vertexR; 1087 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr); 1088 if (dim == 2) { 1089 if ((bdY != DM_BOUNDARY_PERIODIC) && (bdY != DM_BOUNDARY_TWIST)) { 1090 if (vy == numYVertices-1) { 1091 ierr = DMSetLabelValue(dm, "Face Sets", edge, faceMarkerTop);CHKERRQ(ierr); 1092 ierr = DMSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 1093 ierr = DMSetLabelValue(dm, "marker", cone[0], markerTop);CHKERRQ(ierr); 1094 if (ex == numXEdges-1) { 1095 ierr = DMSetLabelValue(dm, "marker", cone[1], markerTop);CHKERRQ(ierr); 1096 } 1097 } else if (vy == 0) { 1098 ierr = DMSetLabelValue(dm, "Face Sets", edge, faceMarkerBottom);CHKERRQ(ierr); 1099 ierr = DMSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 1100 ierr = DMSetLabelValue(dm, "marker", cone[0], markerBottom);CHKERRQ(ierr); 1101 if (ex == numXEdges-1) { 1102 ierr = DMSetLabelValue(dm, "marker", cone[1], markerBottom);CHKERRQ(ierr); 1103 } 1104 } 1105 } else { 1106 if (vy == 0 && cutLabel) { 1107 ierr = DMLabelSetValue(cutLabel, edge, 1);CHKERRQ(ierr); 1108 ierr = DMLabelSetValue(cutLabel, cone[0], 1);CHKERRQ(ierr); 1109 if (ex == numXEdges-1) { 1110 ierr = DMLabelSetValue(cutLabel, cone[1], 1);CHKERRQ(ierr); 1111 } 1112 } 1113 } 1114 } else { 1115 if (bdY != DM_BOUNDARY_PERIODIC) { 1116 if (vy == numYVertices-1) { 1117 ierr = DMSetLabelValue(dm, "marker", edge, markerBack);CHKERRQ(ierr); 1118 } 1119 else if (vy == 0) { 1120 ierr = DMSetLabelValue(dm, "marker", edge, markerFront);CHKERRQ(ierr); 1121 } 1122 } 1123 if (bdZ != DM_BOUNDARY_PERIODIC) { 1124 if (vz == numZVertices-1) { 1125 ierr = DMSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr); 1126 } 1127 else if (vz == 0) { 1128 ierr = DMSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr); 1129 } 1130 } 1131 } 1132 } 1133 } 1134 } 1135 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 1136 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 1137 /* Build coordinates */ 1138 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1139 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1140 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 1141 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVertices);CHKERRQ(ierr); 1142 for (v = firstVertex; v < firstVertex+numVertices; ++v) { 1143 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 1144 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 1145 } 1146 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1147 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1148 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 1149 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1150 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1151 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 1152 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1153 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1154 for (vz = 0; vz < numZVertices; ++vz) { 1155 for (vy = 0; vy < numYVertices; ++vy) { 1156 for (vx = 0; vx < numXVertices; ++vx) { 1157 coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+0] = lower[0] + ((upper[0] - lower[0])/numXEdges)*vx; 1158 coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+1] = lower[1] + ((upper[1] - lower[1])/numYEdges)*vy; 1159 if (dim == 3) { 1160 coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+2] = lower[2] + ((upper[2] - lower[2])/numZEdges)*vz; 1161 } 1162 } 1163 } 1164 } 1165 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1166 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 1167 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1168 } 1169 PetscFunctionReturn(0); 1170 } 1171 1172 static PetscErrorCode DMPlexCreateBoxMesh_Tensor_Internal(DM dm, PetscInt dim, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[]) 1173 { 1174 DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE}; 1175 PetscInt fac[3] = {0, 0, 0}, d; 1176 PetscErrorCode ierr; 1177 1178 PetscFunctionBegin; 1179 PetscValidPointer(dm, 1); 1180 PetscValidLogicalCollectiveInt(dm, dim, 2); 1181 ierr = DMSetDimension(dm, dim);CHKERRQ(ierr); 1182 for (d = 0; d < dim; ++d) {fac[d] = faces[d]; bdt[d] = periodicity[d];} 1183 ierr = DMPlexCreateCubeMesh_Internal(dm, lower, upper, fac, bdt[0], bdt[1], bdt[2]);CHKERRQ(ierr); 1184 if (periodicity[0] == DM_BOUNDARY_PERIODIC || periodicity[0] == DM_BOUNDARY_TWIST || 1185 periodicity[1] == DM_BOUNDARY_PERIODIC || periodicity[1] == DM_BOUNDARY_TWIST || 1186 (dim > 2 && (periodicity[2] == DM_BOUNDARY_PERIODIC || periodicity[2] == DM_BOUNDARY_TWIST))) { 1187 PetscReal L[3]; 1188 PetscReal maxCell[3]; 1189 1190 for (d = 0; d < dim; ++d) { 1191 L[d] = upper[d] - lower[d]; 1192 maxCell[d] = 1.1 * (L[d] / PetscMax(1, faces[d])); 1193 } 1194 ierr = DMSetPeriodicity(dm, PETSC_TRUE, maxCell, L, periodicity);CHKERRQ(ierr); 1195 } 1196 ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE);CHKERRQ(ierr); 1197 PetscFunctionReturn(0); 1198 } 1199 1200 static PetscErrorCode DMPlexCreateBoxMesh_Internal(DM dm, PetscInt dim, PetscBool simplex, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool interpolate) 1201 { 1202 PetscErrorCode ierr; 1203 1204 PetscFunctionBegin; 1205 if (dim == 1) {ierr = DMPlexCreateLineMesh_Internal(dm, faces[0], lower[0], upper[0], periodicity[0]);CHKERRQ(ierr);} 1206 else if (simplex) {ierr = DMPlexCreateBoxMesh_Simplex_Internal(dm, dim, faces, lower, upper, periodicity, interpolate);CHKERRQ(ierr);} 1207 else {ierr = DMPlexCreateBoxMesh_Tensor_Internal(dm, dim, faces, lower, upper, periodicity);CHKERRQ(ierr);} 1208 if (!interpolate && dim > 1 && !simplex) { 1209 DM udm; 1210 1211 ierr = DMPlexUninterpolate(dm, &udm);CHKERRQ(ierr); 1212 ierr = DMPlexCopyCoordinates(dm, udm);CHKERRQ(ierr); 1213 ierr = DMPlexReplace_Static(dm, &udm);CHKERRQ(ierr); 1214 } 1215 PetscFunctionReturn(0); 1216 } 1217 1218 /*@C 1219 DMPlexCreateBoxMesh - Creates a mesh on the tensor product of unit intervals (box) using simplices or tensor cells (hexahedra). 1220 1221 Collective 1222 1223 Input Parameters: 1224 + comm - The communicator for the DM object 1225 . dim - The spatial dimension 1226 . simplex - PETSC_TRUE for simplices, PETSC_FALSE for tensor cells 1227 . faces - Number of faces per dimension, or NULL for (1,) in 1D and (2, 2) in 2D and (1, 1, 1) in 3D 1228 . lower - The lower left corner, or NULL for (0, 0, 0) 1229 . upper - The upper right corner, or NULL for (1, 1, 1) 1230 . periodicity - The boundary type for the X,Y,Z direction, or NULL for DM_BOUNDARY_NONE 1231 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 1232 1233 Output Parameter: 1234 . dm - The DM object 1235 1236 Note: If you want to customize this mesh using options, you just need to 1237 $ DMCreate(comm, &dm); 1238 $ DMSetType(dm, DMPLEX); 1239 $ DMSetFromOptions(dm); 1240 and use the options on the DMSetFromOptions() page. 1241 1242 Here is the numbering returned for 2 faces in each direction for tensor cells: 1243 $ 10---17---11---18----12 1244 $ | | | 1245 $ | | | 1246 $ 20 2 22 3 24 1247 $ | | | 1248 $ | | | 1249 $ 7---15----8---16----9 1250 $ | | | 1251 $ | | | 1252 $ 19 0 21 1 23 1253 $ | | | 1254 $ | | | 1255 $ 4---13----5---14----6 1256 1257 and for simplicial cells 1258 1259 $ 14----8---15----9----16 1260 $ |\ 5 |\ 7 | 1261 $ | \ | \ | 1262 $ 13 2 14 3 15 1263 $ | 4 \ | 6 \ | 1264 $ | \ | \ | 1265 $ 11----6---12----7----13 1266 $ |\ |\ | 1267 $ | \ 1 | \ 3 | 1268 $ 10 0 11 1 12 1269 $ | 0 \ | 2 \ | 1270 $ | \ | \ | 1271 $ 8----4----9----5----10 1272 1273 Level: beginner 1274 1275 .seealso: DMSetFromOptions(), DMPlexCreateFromFile(), DMPlexCreateHexCylinderMesh(), DMSetType(), DMCreate() 1276 @*/ 1277 PetscErrorCode DMPlexCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscBool simplex, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool interpolate, DM *dm) 1278 { 1279 PetscInt fac[3] = {1, 1, 1}; 1280 PetscReal low[3] = {0, 0, 0}; 1281 PetscReal upp[3] = {1, 1, 1}; 1282 DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE}; 1283 PetscErrorCode ierr; 1284 1285 PetscFunctionBegin; 1286 ierr = DMCreate(comm,dm);CHKERRQ(ierr); 1287 ierr = DMSetType(*dm,DMPLEX);CHKERRQ(ierr); 1288 ierr = DMPlexCreateBoxMesh_Internal(*dm, dim, simplex, faces ? faces : fac, lower ? lower : low, upper ? upper : upp, periodicity ? periodicity : bdt, interpolate);CHKERRQ(ierr); 1289 PetscFunctionReturn(0); 1290 } 1291 1292 static PetscErrorCode DMPlexCreateWedgeBoxMesh_Internal(DM dm, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[]) 1293 { 1294 DM bdm, vol; 1295 PetscInt i; 1296 PetscErrorCode ierr; 1297 1298 PetscFunctionBegin; 1299 for (i = 0; i < 3; ++i) PetscCheckFalse(periodicity[i] != DM_BOUNDARY_NONE,PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Periodicity not yet supported"); 1300 ierr = DMCreate(PetscObjectComm((PetscObject) dm), &bdm);CHKERRQ(ierr); 1301 ierr = DMSetType(bdm, DMPLEX);CHKERRQ(ierr); 1302 ierr = DMSetDimension(bdm, 2);CHKERRQ(ierr); 1303 ierr = DMPlexCreateBoxMesh_Simplex_Internal(bdm, 2, faces, lower, upper, periodicity, PETSC_TRUE);CHKERRQ(ierr); 1304 ierr = DMPlexExtrude(bdm, faces[2], upper[2] - lower[2], PETSC_TRUE, PETSC_FALSE, NULL, NULL, &vol);CHKERRQ(ierr); 1305 ierr = DMDestroy(&bdm);CHKERRQ(ierr); 1306 ierr = DMPlexReplace_Static(dm, &vol);CHKERRQ(ierr); 1307 if (lower[2] != 0.0) { 1308 Vec v; 1309 PetscScalar *x; 1310 PetscInt cDim, n; 1311 1312 ierr = DMGetCoordinatesLocal(dm, &v);CHKERRQ(ierr); 1313 ierr = VecGetBlockSize(v, &cDim);CHKERRQ(ierr); 1314 ierr = VecGetLocalSize(v, &n);CHKERRQ(ierr); 1315 ierr = VecGetArray(v, &x);CHKERRQ(ierr); 1316 x += cDim; 1317 for (i = 0; i < n; i += cDim) x[i] += lower[2]; 1318 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1319 ierr = DMSetCoordinatesLocal(dm, v);CHKERRQ(ierr); 1320 } 1321 PetscFunctionReturn(0); 1322 } 1323 1324 /*@ 1325 DMPlexCreateWedgeBoxMesh - Creates a 3-D mesh tesselating the (x,y) plane and extruding in the third direction using wedge cells. 1326 1327 Collective 1328 1329 Input Parameters: 1330 + comm - The communicator for the DM object 1331 . faces - Number of faces per dimension, or NULL for (1, 1, 1) 1332 . lower - The lower left corner, or NULL for (0, 0, 0) 1333 . upper - The upper right corner, or NULL for (1, 1, 1) 1334 . periodicity - The boundary type for the X,Y,Z direction, or NULL for DM_BOUNDARY_NONE 1335 . orderHeight - If PETSC_TRUE, orders the extruded cells in the height first. Otherwise, orders the cell on the layers first 1336 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 1337 1338 Output Parameter: 1339 . dm - The DM object 1340 1341 Level: beginner 1342 1343 .seealso: DMPlexCreateHexCylinderMesh(), DMPlexCreateWedgeCylinderMesh(), DMExtrude(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 1344 @*/ 1345 PetscErrorCode DMPlexCreateWedgeBoxMesh(MPI_Comm comm, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool orderHeight, PetscBool interpolate, DM *dm) 1346 { 1347 PetscInt fac[3] = {1, 1, 1}; 1348 PetscReal low[3] = {0, 0, 0}; 1349 PetscReal upp[3] = {1, 1, 1}; 1350 DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE}; 1351 PetscErrorCode ierr; 1352 1353 PetscFunctionBegin; 1354 ierr = DMCreate(comm,dm);CHKERRQ(ierr); 1355 ierr = DMSetType(*dm,DMPLEX);CHKERRQ(ierr); 1356 ierr = DMPlexCreateWedgeBoxMesh_Internal(*dm, faces ? faces : fac, lower ? lower : low, upper ? upper : upp, periodicity ? periodicity : bdt);CHKERRQ(ierr); 1357 if (!interpolate) { 1358 DM udm; 1359 1360 ierr = DMPlexUninterpolate(*dm, &udm);CHKERRQ(ierr); 1361 ierr = DMPlexReplace_Static(*dm, &udm);CHKERRQ(ierr); 1362 } 1363 PetscFunctionReturn(0); 1364 } 1365 1366 /*@C 1367 DMPlexSetOptionsPrefix - Sets the prefix used for searching for all DM options in the database. 1368 1369 Logically Collective on dm 1370 1371 Input Parameters: 1372 + dm - the DM context 1373 - prefix - the prefix to prepend to all option names 1374 1375 Notes: 1376 A hyphen (-) must NOT be given at the beginning of the prefix name. 1377 The first character of all runtime options is AUTOMATICALLY the hyphen. 1378 1379 Level: advanced 1380 1381 .seealso: SNESSetFromOptions() 1382 @*/ 1383 PetscErrorCode DMPlexSetOptionsPrefix(DM dm, const char prefix[]) 1384 { 1385 DM_Plex *mesh = (DM_Plex *) dm->data; 1386 PetscErrorCode ierr; 1387 1388 PetscFunctionBegin; 1389 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1390 ierr = PetscObjectSetOptionsPrefix((PetscObject) dm, prefix);CHKERRQ(ierr); 1391 ierr = PetscObjectSetOptionsPrefix((PetscObject) mesh->partitioner, prefix);CHKERRQ(ierr); 1392 PetscFunctionReturn(0); 1393 } 1394 1395 /* Remap geometry to cylinder 1396 TODO: This only works for a single refinement, then it is broken 1397 1398 Interior square: Linear interpolation is correct 1399 The other cells all have vertices on rays from the origin. We want to uniformly expand the spacing 1400 such that the last vertex is on the unit circle. So the closest and farthest vertices are at distance 1401 1402 phi = arctan(y/x) 1403 d_close = sqrt(1/8 + 1/4 sin^2(phi)) 1404 d_far = sqrt(1/2 + sin^2(phi)) 1405 1406 so we remap them using 1407 1408 x_new = x_close + (x - x_close) (1 - d_close) / (d_far - d_close) 1409 y_new = y_close + (y - y_close) (1 - d_close) / (d_far - d_close) 1410 1411 If pi/4 < phi < 3pi/4 or -3pi/4 < phi < -pi/4, then we switch x and y. 1412 */ 1413 static void snapToCylinder(PetscInt dim, PetscInt Nf, PetscInt NfAux, 1414 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 1415 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 1416 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 1417 { 1418 const PetscReal dis = 1.0/PetscSqrtReal(2.0); 1419 const PetscReal ds2 = 0.5*dis; 1420 1421 if ((PetscAbsScalar(u[0]) <= ds2) && (PetscAbsScalar(u[1]) <= ds2)) { 1422 f0[0] = u[0]; 1423 f0[1] = u[1]; 1424 } else { 1425 PetscReal phi, sinp, cosp, dc, df, x, y, xc, yc; 1426 1427 x = PetscRealPart(u[0]); 1428 y = PetscRealPart(u[1]); 1429 phi = PetscAtan2Real(y, x); 1430 sinp = PetscSinReal(phi); 1431 cosp = PetscCosReal(phi); 1432 if ((PetscAbsReal(phi) > PETSC_PI/4.0) && (PetscAbsReal(phi) < 3.0*PETSC_PI/4.0)) { 1433 dc = PetscAbsReal(ds2/sinp); 1434 df = PetscAbsReal(dis/sinp); 1435 xc = ds2*x/PetscAbsReal(y); 1436 yc = ds2*PetscSignReal(y); 1437 } else { 1438 dc = PetscAbsReal(ds2/cosp); 1439 df = PetscAbsReal(dis/cosp); 1440 xc = ds2*PetscSignReal(x); 1441 yc = ds2*y/PetscAbsReal(x); 1442 } 1443 f0[0] = xc + (u[0] - xc)*(1.0 - dc)/(df - dc); 1444 f0[1] = yc + (u[1] - yc)*(1.0 - dc)/(df - dc); 1445 } 1446 f0[2] = u[2]; 1447 } 1448 1449 static PetscErrorCode DMPlexCreateHexCylinderMesh_Internal(DM dm, DMBoundaryType periodicZ) 1450 { 1451 const PetscInt dim = 3; 1452 PetscInt numCells, numVertices; 1453 PetscMPIInt rank; 1454 PetscErrorCode ierr; 1455 1456 PetscFunctionBegin; 1457 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRMPI(ierr); 1458 ierr = DMSetDimension(dm, dim);CHKERRQ(ierr); 1459 /* Create topology */ 1460 { 1461 PetscInt cone[8], c; 1462 1463 numCells = rank == 0 ? 5 : 0; 1464 numVertices = rank == 0 ? 16 : 0; 1465 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1466 numCells *= 3; 1467 numVertices = rank == 0 ? 24 : 0; 1468 } 1469 ierr = DMPlexSetChart(dm, 0, numCells+numVertices);CHKERRQ(ierr); 1470 for (c = 0; c < numCells; c++) {ierr = DMPlexSetConeSize(dm, c, 8);CHKERRQ(ierr);} 1471 ierr = DMSetUp(dm);CHKERRQ(ierr); 1472 if (rank == 0) { 1473 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1474 cone[0] = 15; cone[1] = 18; cone[2] = 17; cone[3] = 16; 1475 cone[4] = 31; cone[5] = 32; cone[6] = 33; cone[7] = 34; 1476 ierr = DMPlexSetCone(dm, 0, cone);CHKERRQ(ierr); 1477 cone[0] = 16; cone[1] = 17; cone[2] = 24; cone[3] = 23; 1478 cone[4] = 32; cone[5] = 36; cone[6] = 37; cone[7] = 33; /* 22 25 26 21 */ 1479 ierr = DMPlexSetCone(dm, 1, cone);CHKERRQ(ierr); 1480 cone[0] = 18; cone[1] = 27; cone[2] = 24; cone[3] = 17; 1481 cone[4] = 34; cone[5] = 33; cone[6] = 37; cone[7] = 38; 1482 ierr = DMPlexSetCone(dm, 2, cone);CHKERRQ(ierr); 1483 cone[0] = 29; cone[1] = 27; cone[2] = 18; cone[3] = 15; 1484 cone[4] = 35; cone[5] = 31; cone[6] = 34; cone[7] = 38; 1485 ierr = DMPlexSetCone(dm, 3, cone);CHKERRQ(ierr); 1486 cone[0] = 29; cone[1] = 15; cone[2] = 16; cone[3] = 23; 1487 cone[4] = 35; cone[5] = 36; cone[6] = 32; cone[7] = 31; 1488 ierr = DMPlexSetCone(dm, 4, cone);CHKERRQ(ierr); 1489 1490 cone[0] = 31; cone[1] = 34; cone[2] = 33; cone[3] = 32; 1491 cone[4] = 19; cone[5] = 22; cone[6] = 21; cone[7] = 20; 1492 ierr = DMPlexSetCone(dm, 5, cone);CHKERRQ(ierr); 1493 cone[0] = 32; cone[1] = 33; cone[2] = 37; cone[3] = 36; 1494 cone[4] = 22; cone[5] = 25; cone[6] = 26; cone[7] = 21; 1495 ierr = DMPlexSetCone(dm, 6, cone);CHKERRQ(ierr); 1496 cone[0] = 34; cone[1] = 38; cone[2] = 37; cone[3] = 33; 1497 cone[4] = 20; cone[5] = 21; cone[6] = 26; cone[7] = 28; 1498 ierr = DMPlexSetCone(dm, 7, cone);CHKERRQ(ierr); 1499 cone[0] = 35; cone[1] = 38; cone[2] = 34; cone[3] = 31; 1500 cone[4] = 30; cone[5] = 19; cone[6] = 20; cone[7] = 28; 1501 ierr = DMPlexSetCone(dm, 8, cone);CHKERRQ(ierr); 1502 cone[0] = 35; cone[1] = 31; cone[2] = 32; cone[3] = 36; 1503 cone[4] = 30; cone[5] = 25; cone[6] = 22; cone[7] = 19; 1504 ierr = DMPlexSetCone(dm, 9, cone);CHKERRQ(ierr); 1505 1506 cone[0] = 19; cone[1] = 20; cone[2] = 21; cone[3] = 22; 1507 cone[4] = 15; cone[5] = 16; cone[6] = 17; cone[7] = 18; 1508 ierr = DMPlexSetCone(dm, 10, cone);CHKERRQ(ierr); 1509 cone[0] = 22; cone[1] = 21; cone[2] = 26; cone[3] = 25; 1510 cone[4] = 16; cone[5] = 23; cone[6] = 24; cone[7] = 17; 1511 ierr = DMPlexSetCone(dm, 11, cone);CHKERRQ(ierr); 1512 cone[0] = 20; cone[1] = 28; cone[2] = 26; cone[3] = 21; 1513 cone[4] = 18; cone[5] = 17; cone[6] = 24; cone[7] = 27; 1514 ierr = DMPlexSetCone(dm, 12, cone);CHKERRQ(ierr); 1515 cone[0] = 30; cone[1] = 28; cone[2] = 20; cone[3] = 19; 1516 cone[4] = 29; cone[5] = 15; cone[6] = 18; cone[7] = 27; 1517 ierr = DMPlexSetCone(dm, 13, cone);CHKERRQ(ierr); 1518 cone[0] = 30; cone[1] = 19; cone[2] = 22; cone[3] = 25; 1519 cone[4] = 29; cone[5] = 23; cone[6] = 16; cone[7] = 15; 1520 ierr = DMPlexSetCone(dm, 14, cone);CHKERRQ(ierr); 1521 } else { 1522 cone[0] = 5; cone[1] = 8; cone[2] = 7; cone[3] = 6; 1523 cone[4] = 9; cone[5] = 12; cone[6] = 11; cone[7] = 10; 1524 ierr = DMPlexSetCone(dm, 0, cone);CHKERRQ(ierr); 1525 cone[0] = 6; cone[1] = 7; cone[2] = 14; cone[3] = 13; 1526 cone[4] = 12; cone[5] = 15; cone[6] = 16; cone[7] = 11; 1527 ierr = DMPlexSetCone(dm, 1, cone);CHKERRQ(ierr); 1528 cone[0] = 8; cone[1] = 17; cone[2] = 14; cone[3] = 7; 1529 cone[4] = 10; cone[5] = 11; cone[6] = 16; cone[7] = 18; 1530 ierr = DMPlexSetCone(dm, 2, cone);CHKERRQ(ierr); 1531 cone[0] = 19; cone[1] = 17; cone[2] = 8; cone[3] = 5; 1532 cone[4] = 20; cone[5] = 9; cone[6] = 10; cone[7] = 18; 1533 ierr = DMPlexSetCone(dm, 3, cone);CHKERRQ(ierr); 1534 cone[0] = 19; cone[1] = 5; cone[2] = 6; cone[3] = 13; 1535 cone[4] = 20; cone[5] = 15; cone[6] = 12; cone[7] = 9; 1536 ierr = DMPlexSetCone(dm, 4, cone);CHKERRQ(ierr); 1537 } 1538 } 1539 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 1540 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 1541 } 1542 /* Create cube geometry */ 1543 { 1544 Vec coordinates; 1545 PetscSection coordSection; 1546 PetscScalar *coords; 1547 PetscInt coordSize, v; 1548 const PetscReal dis = 1.0/PetscSqrtReal(2.0); 1549 const PetscReal ds2 = dis/2.0; 1550 1551 /* Build coordinates */ 1552 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1553 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1554 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 1555 ierr = PetscSectionSetChart(coordSection, numCells, numCells+numVertices);CHKERRQ(ierr); 1556 for (v = numCells; v < numCells+numVertices; ++v) { 1557 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 1558 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 1559 } 1560 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1561 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1562 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 1563 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1564 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1565 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 1566 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1567 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1568 if (rank == 0) { 1569 coords[0*dim+0] = -ds2; coords[0*dim+1] = -ds2; coords[0*dim+2] = 0.0; 1570 coords[1*dim+0] = ds2; coords[1*dim+1] = -ds2; coords[1*dim+2] = 0.0; 1571 coords[2*dim+0] = ds2; coords[2*dim+1] = ds2; coords[2*dim+2] = 0.0; 1572 coords[3*dim+0] = -ds2; coords[3*dim+1] = ds2; coords[3*dim+2] = 0.0; 1573 coords[4*dim+0] = -ds2; coords[4*dim+1] = -ds2; coords[4*dim+2] = 1.0; 1574 coords[5*dim+0] = -ds2; coords[5*dim+1] = ds2; coords[5*dim+2] = 1.0; 1575 coords[6*dim+0] = ds2; coords[6*dim+1] = ds2; coords[6*dim+2] = 1.0; 1576 coords[7*dim+0] = ds2; coords[7*dim+1] = -ds2; coords[7*dim+2] = 1.0; 1577 coords[ 8*dim+0] = dis; coords[ 8*dim+1] = -dis; coords[ 8*dim+2] = 0.0; 1578 coords[ 9*dim+0] = dis; coords[ 9*dim+1] = dis; coords[ 9*dim+2] = 0.0; 1579 coords[10*dim+0] = dis; coords[10*dim+1] = -dis; coords[10*dim+2] = 1.0; 1580 coords[11*dim+0] = dis; coords[11*dim+1] = dis; coords[11*dim+2] = 1.0; 1581 coords[12*dim+0] = -dis; coords[12*dim+1] = dis; coords[12*dim+2] = 0.0; 1582 coords[13*dim+0] = -dis; coords[13*dim+1] = dis; coords[13*dim+2] = 1.0; 1583 coords[14*dim+0] = -dis; coords[14*dim+1] = -dis; coords[14*dim+2] = 0.0; 1584 coords[15*dim+0] = -dis; coords[15*dim+1] = -dis; coords[15*dim+2] = 1.0; 1585 if (periodicZ == DM_BOUNDARY_PERIODIC) { 1586 /* 15 31 19 */ coords[16*dim+0] = -ds2; coords[16*dim+1] = -ds2; coords[16*dim+2] = 0.5; 1587 /* 16 32 22 */ coords[17*dim+0] = ds2; coords[17*dim+1] = -ds2; coords[17*dim+2] = 0.5; 1588 /* 17 33 21 */ coords[18*dim+0] = ds2; coords[18*dim+1] = ds2; coords[18*dim+2] = 0.5; 1589 /* 18 34 20 */ coords[19*dim+0] = -ds2; coords[19*dim+1] = ds2; coords[19*dim+2] = 0.5; 1590 /* 29 35 30 */ coords[20*dim+0] = -dis; coords[20*dim+1] = -dis; coords[20*dim+2] = 0.5; 1591 /* 23 36 25 */ coords[21*dim+0] = dis; coords[21*dim+1] = -dis; coords[21*dim+2] = 0.5; 1592 /* 24 37 26 */ coords[22*dim+0] = dis; coords[22*dim+1] = dis; coords[22*dim+2] = 0.5; 1593 /* 27 38 28 */ coords[23*dim+0] = -dis; coords[23*dim+1] = dis; coords[23*dim+2] = 0.5; 1594 } 1595 } 1596 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1597 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 1598 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1599 } 1600 /* Create periodicity */ 1601 if (periodicZ == DM_BOUNDARY_PERIODIC || periodicZ == DM_BOUNDARY_TWIST) { 1602 PetscReal L[3]; 1603 PetscReal maxCell[3]; 1604 DMBoundaryType bdType[3]; 1605 PetscReal lower[3] = {0.0, 0.0, 0.0}; 1606 PetscReal upper[3] = {1.0, 1.0, 1.5}; 1607 PetscInt i, numZCells = 3; 1608 1609 bdType[0] = DM_BOUNDARY_NONE; 1610 bdType[1] = DM_BOUNDARY_NONE; 1611 bdType[2] = periodicZ; 1612 for (i = 0; i < dim; i++) { 1613 L[i] = upper[i] - lower[i]; 1614 maxCell[i] = 1.1 * (L[i] / numZCells); 1615 } 1616 ierr = DMSetPeriodicity(dm, PETSC_TRUE, maxCell, L, bdType);CHKERRQ(ierr); 1617 } 1618 { 1619 DM cdm; 1620 PetscDS cds; 1621 PetscScalar c[2] = {1.0, 1.0}; 1622 1623 ierr = DMPlexCreateCoordinateSpace(dm, 1, snapToCylinder);CHKERRQ(ierr); 1624 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 1625 ierr = DMGetDS(cdm, &cds);CHKERRQ(ierr); 1626 ierr = PetscDSSetConstants(cds, 2, c);CHKERRQ(ierr); 1627 } 1628 /* Wait for coordinate creation before doing in-place modification */ 1629 ierr = DMPlexInterpolateInPlace_Internal(dm);CHKERRQ(ierr); 1630 PetscFunctionReturn(0); 1631 } 1632 1633 /*@ 1634 DMPlexCreateHexCylinderMesh - Creates a mesh on the tensor product of the unit interval with the circle (cylinder) using hexahedra. 1635 1636 Collective 1637 1638 Input Parameters: 1639 + comm - The communicator for the DM object 1640 - periodicZ - The boundary type for the Z direction 1641 1642 Output Parameter: 1643 . dm - The DM object 1644 1645 Note: 1646 Here is the output numbering looking from the bottom of the cylinder: 1647 $ 17-----14 1648 $ | | 1649 $ | 2 | 1650 $ | | 1651 $ 17-----8-----7-----14 1652 $ | | | | 1653 $ | 3 | 0 | 1 | 1654 $ | | | | 1655 $ 19-----5-----6-----13 1656 $ | | 1657 $ | 4 | 1658 $ | | 1659 $ 19-----13 1660 $ 1661 $ and up through the top 1662 $ 1663 $ 18-----16 1664 $ | | 1665 $ | 2 | 1666 $ | | 1667 $ 18----10----11-----16 1668 $ | | | | 1669 $ | 3 | 0 | 1 | 1670 $ | | | | 1671 $ 20-----9----12-----15 1672 $ | | 1673 $ | 4 | 1674 $ | | 1675 $ 20-----15 1676 1677 Level: beginner 1678 1679 .seealso: DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 1680 @*/ 1681 PetscErrorCode DMPlexCreateHexCylinderMesh(MPI_Comm comm, DMBoundaryType periodicZ, DM *dm) 1682 { 1683 PetscErrorCode ierr; 1684 1685 PetscFunctionBegin; 1686 PetscValidPointer(dm, 3); 1687 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1688 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1689 ierr = DMPlexCreateHexCylinderMesh_Internal(*dm, periodicZ);CHKERRQ(ierr); 1690 PetscFunctionReturn(0); 1691 } 1692 1693 static PetscErrorCode DMPlexCreateWedgeCylinderMesh_Internal(DM dm, PetscInt n, PetscBool interpolate) 1694 { 1695 const PetscInt dim = 3; 1696 PetscInt numCells, numVertices, v; 1697 PetscMPIInt rank; 1698 PetscErrorCode ierr; 1699 1700 PetscFunctionBegin; 1701 PetscCheckFalse(n < 0,PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Number of wedges %D cannot be negative", n); 1702 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRMPI(ierr); 1703 ierr = DMSetDimension(dm, dim);CHKERRQ(ierr); 1704 /* Must create the celltype label here so that we do not automatically try to compute the types */ 1705 ierr = DMCreateLabel(dm, "celltype");CHKERRQ(ierr); 1706 /* Create topology */ 1707 { 1708 PetscInt cone[6], c; 1709 1710 numCells = rank == 0 ? n : 0; 1711 numVertices = rank == 0 ? 2*(n+1) : 0; 1712 ierr = DMPlexSetChart(dm, 0, numCells+numVertices);CHKERRQ(ierr); 1713 for (c = 0; c < numCells; c++) {ierr = DMPlexSetConeSize(dm, c, 6);CHKERRQ(ierr);} 1714 ierr = DMSetUp(dm);CHKERRQ(ierr); 1715 for (c = 0; c < numCells; c++) { 1716 cone[0] = c+n*1; cone[1] = (c+1)%n+n*1; cone[2] = 0+3*n; 1717 cone[3] = c+n*2; cone[4] = (c+1)%n+n*2; cone[5] = 1+3*n; 1718 ierr = DMPlexSetCone(dm, c, cone);CHKERRQ(ierr); 1719 ierr = DMPlexSetCellType(dm, c, DM_POLYTOPE_TRI_PRISM_TENSOR);CHKERRQ(ierr); 1720 } 1721 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 1722 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 1723 } 1724 for (v = numCells; v < numCells+numVertices; ++v) { 1725 ierr = DMPlexSetCellType(dm, v, DM_POLYTOPE_POINT);CHKERRQ(ierr); 1726 } 1727 /* Create cylinder geometry */ 1728 { 1729 Vec coordinates; 1730 PetscSection coordSection; 1731 PetscScalar *coords; 1732 PetscInt coordSize, c; 1733 1734 /* Build coordinates */ 1735 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1736 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1737 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 1738 ierr = PetscSectionSetChart(coordSection, numCells, numCells+numVertices);CHKERRQ(ierr); 1739 for (v = numCells; v < numCells+numVertices; ++v) { 1740 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 1741 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 1742 } 1743 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1744 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1745 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 1746 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1747 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1748 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 1749 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 1750 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1751 for (c = 0; c < numCells; c++) { 1752 coords[(c+0*n)*dim+0] = PetscCosReal(2.0*c*PETSC_PI/n); coords[(c+0*n)*dim+1] = PetscSinReal(2.0*c*PETSC_PI/n); coords[(c+0*n)*dim+2] = 1.0; 1753 coords[(c+1*n)*dim+0] = PetscCosReal(2.0*c*PETSC_PI/n); coords[(c+1*n)*dim+1] = PetscSinReal(2.0*c*PETSC_PI/n); coords[(c+1*n)*dim+2] = 0.0; 1754 } 1755 if (rank == 0) { 1756 coords[(2*n+0)*dim+0] = 0.0; coords[(2*n+0)*dim+1] = 0.0; coords[(2*n+0)*dim+2] = 1.0; 1757 coords[(2*n+1)*dim+0] = 0.0; coords[(2*n+1)*dim+1] = 0.0; coords[(2*n+1)*dim+2] = 0.0; 1758 } 1759 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1760 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 1761 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1762 } 1763 /* Interpolate */ 1764 if (interpolate) {ierr = DMPlexInterpolateInPlace_Internal(dm);CHKERRQ(ierr);} 1765 PetscFunctionReturn(0); 1766 } 1767 1768 /*@ 1769 DMPlexCreateWedgeCylinderMesh - Creates a mesh on the tensor product of the unit interval with the circle (cylinder) using wedges. 1770 1771 Collective 1772 1773 Input Parameters: 1774 + comm - The communicator for the DM object 1775 . n - The number of wedges around the origin 1776 - interpolate - Create edges and faces 1777 1778 Output Parameter: 1779 . dm - The DM object 1780 1781 Level: beginner 1782 1783 .seealso: DMPlexCreateHexCylinderMesh(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 1784 @*/ 1785 PetscErrorCode DMPlexCreateWedgeCylinderMesh(MPI_Comm comm, PetscInt n, PetscBool interpolate, DM *dm) 1786 { 1787 PetscErrorCode ierr; 1788 1789 PetscFunctionBegin; 1790 PetscValidPointer(dm, 4); 1791 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1792 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1793 ierr = DMPlexCreateWedgeCylinderMesh_Internal(*dm, n, interpolate);CHKERRQ(ierr); 1794 PetscFunctionReturn(0); 1795 } 1796 1797 static inline PetscReal DiffNormReal(PetscInt dim, const PetscReal x[], const PetscReal y[]) 1798 { 1799 PetscReal prod = 0.0; 1800 PetscInt i; 1801 for (i = 0; i < dim; ++i) prod += PetscSqr(x[i] - y[i]); 1802 return PetscSqrtReal(prod); 1803 } 1804 static inline PetscReal DotReal(PetscInt dim, const PetscReal x[], const PetscReal y[]) 1805 { 1806 PetscReal prod = 0.0; 1807 PetscInt i; 1808 for (i = 0; i < dim; ++i) prod += x[i]*y[i]; 1809 return prod; 1810 } 1811 1812 /* The first constant is the sphere radius */ 1813 static void snapToSphere(PetscInt dim, PetscInt Nf, PetscInt NfAux, 1814 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 1815 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 1816 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 1817 { 1818 PetscReal r = PetscRealPart(constants[0]); 1819 PetscReal norm2 = 0.0, fac; 1820 PetscInt n = uOff[1] - uOff[0], d; 1821 1822 for (d = 0; d < n; ++d) norm2 += PetscSqr(PetscRealPart(u[d])); 1823 fac = r/PetscSqrtReal(norm2); 1824 for (d = 0; d < n; ++d) f0[d] = u[d]*fac; 1825 } 1826 1827 static PetscErrorCode DMPlexCreateSphereMesh_Internal(DM dm, PetscInt dim, PetscBool simplex, PetscReal R) 1828 { 1829 const PetscInt embedDim = dim+1; 1830 PetscSection coordSection; 1831 Vec coordinates; 1832 PetscScalar *coords; 1833 PetscReal *coordsIn; 1834 PetscInt numCells, numEdges, numVerts, firstVertex, v, firstEdge, coordSize, d, c, e; 1835 PetscMPIInt rank; 1836 PetscErrorCode ierr; 1837 1838 PetscFunctionBegin; 1839 PetscValidLogicalCollectiveBool(dm, simplex, 3); 1840 ierr = DMSetDimension(dm, dim);CHKERRQ(ierr); 1841 ierr = DMSetCoordinateDim(dm, dim+1);CHKERRQ(ierr); 1842 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRMPI(ierr); 1843 switch (dim) { 1844 case 2: 1845 if (simplex) { 1846 const PetscReal radius = PetscSqrtReal(1 + PETSC_PHI*PETSC_PHI)/(1.0 + PETSC_PHI); 1847 const PetscReal edgeLen = 2.0/(1.0 + PETSC_PHI) * (R/radius); 1848 const PetscInt degree = 5; 1849 PetscReal vertex[3] = {0.0, 1.0/(1.0 + PETSC_PHI), PETSC_PHI/(1.0 + PETSC_PHI)}; 1850 PetscInt s[3] = {1, 1, 1}; 1851 PetscInt cone[3]; 1852 PetscInt *graph, p, i, j, k; 1853 1854 vertex[0] *= R/radius; vertex[1] *= R/radius; vertex[2] *= R/radius; 1855 numCells = rank == 0 ? 20 : 0; 1856 numVerts = rank == 0 ? 12 : 0; 1857 firstVertex = numCells; 1858 /* Use icosahedron, which for a R-sphere has coordinates which are all cyclic permutations of 1859 1860 (0, \pm 1/\phi+1, \pm \phi/\phi+1) 1861 1862 where \phi^2 - \phi - 1 = 0, meaning \phi is the golden ratio \frac{1 + \sqrt{5}}{2}. The edge 1863 length is then given by 2/(1+\phi) = 2 * 0.38197 = 0.76393. 1864 */ 1865 /* Construct vertices */ 1866 ierr = PetscCalloc1(numVerts * embedDim, &coordsIn);CHKERRQ(ierr); 1867 if (rank == 0) { 1868 for (p = 0, i = 0; p < embedDim; ++p) { 1869 for (s[1] = -1; s[1] < 2; s[1] += 2) { 1870 for (s[2] = -1; s[2] < 2; s[2] += 2) { 1871 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[(d+p)%embedDim]*vertex[(d+p)%embedDim]; 1872 ++i; 1873 } 1874 } 1875 } 1876 } 1877 /* Construct graph */ 1878 ierr = PetscCalloc1(numVerts * numVerts, &graph);CHKERRQ(ierr); 1879 for (i = 0; i < numVerts; ++i) { 1880 for (j = 0, k = 0; j < numVerts; ++j) { 1881 if (PetscAbsReal(DiffNormReal(embedDim, &coordsIn[i*embedDim], &coordsIn[j*embedDim]) - edgeLen) < PETSC_SMALL) {graph[i*numVerts+j] = 1; ++k;} 1882 } 1883 PetscCheckFalse(k != degree,PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "Invalid icosahedron, vertex %D degree %D != %D", i, k, degree); 1884 } 1885 /* Build Topology */ 1886 ierr = DMPlexSetChart(dm, 0, numCells+numVerts);CHKERRQ(ierr); 1887 for (c = 0; c < numCells; c++) { 1888 ierr = DMPlexSetConeSize(dm, c, embedDim);CHKERRQ(ierr); 1889 } 1890 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 1891 /* Cells */ 1892 for (i = 0, c = 0; i < numVerts; ++i) { 1893 for (j = 0; j < i; ++j) { 1894 for (k = 0; k < j; ++k) { 1895 if (graph[i*numVerts+j] && graph[j*numVerts+k] && graph[k*numVerts+i]) { 1896 cone[0] = firstVertex+i; cone[1] = firstVertex+j; cone[2] = firstVertex+k; 1897 /* Check orientation */ 1898 { 1899 const PetscInt epsilon[3][3][3] = {{{0, 0, 0}, {0, 0, 1}, {0, -1, 0}}, {{0, 0, -1}, {0, 0, 0}, {1, 0, 0}}, {{0, 1, 0}, {-1, 0, 0}, {0, 0, 0}}}; 1900 PetscReal normal[3]; 1901 PetscInt e, f; 1902 1903 for (d = 0; d < embedDim; ++d) { 1904 normal[d] = 0.0; 1905 for (e = 0; e < embedDim; ++e) { 1906 for (f = 0; f < embedDim; ++f) { 1907 normal[d] += epsilon[d][e][f]*(coordsIn[j*embedDim+e] - coordsIn[i*embedDim+e])*(coordsIn[k*embedDim+f] - coordsIn[i*embedDim+f]); 1908 } 1909 } 1910 } 1911 if (DotReal(embedDim, normal, &coordsIn[i*embedDim]) < 0) {PetscInt tmp = cone[1]; cone[1] = cone[2]; cone[2] = tmp;} 1912 } 1913 ierr = DMPlexSetCone(dm, c++, cone);CHKERRQ(ierr); 1914 } 1915 } 1916 } 1917 } 1918 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 1919 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 1920 ierr = PetscFree(graph);CHKERRQ(ierr); 1921 } else { 1922 /* 1923 12-21--13 1924 | | 1925 25 4 24 1926 | | 1927 12-25--9-16--8-24--13 1928 | | | | 1929 23 5 17 0 15 3 22 1930 | | | | 1931 10-20--6-14--7-19--11 1932 | | 1933 20 1 19 1934 | | 1935 10-18--11 1936 | | 1937 23 2 22 1938 | | 1939 12-21--13 1940 */ 1941 PetscInt cone[4], ornt[4]; 1942 1943 numCells = rank == 0 ? 6 : 0; 1944 numEdges = rank == 0 ? 12 : 0; 1945 numVerts = rank == 0 ? 8 : 0; 1946 firstVertex = numCells; 1947 firstEdge = numCells + numVerts; 1948 /* Build Topology */ 1949 ierr = DMPlexSetChart(dm, 0, numCells+numEdges+numVerts);CHKERRQ(ierr); 1950 for (c = 0; c < numCells; c++) { 1951 ierr = DMPlexSetConeSize(dm, c, 4);CHKERRQ(ierr); 1952 } 1953 for (e = firstEdge; e < firstEdge+numEdges; ++e) { 1954 ierr = DMPlexSetConeSize(dm, e, 2);CHKERRQ(ierr); 1955 } 1956 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 1957 if (rank == 0) { 1958 /* Cell 0 */ 1959 cone[0] = 14; cone[1] = 15; cone[2] = 16; cone[3] = 17; 1960 ierr = DMPlexSetCone(dm, 0, cone);CHKERRQ(ierr); 1961 ornt[0] = 0; ornt[1] = 0; ornt[2] = 0; ornt[3] = 0; 1962 ierr = DMPlexSetConeOrientation(dm, 0, ornt);CHKERRQ(ierr); 1963 /* Cell 1 */ 1964 cone[0] = 18; cone[1] = 19; cone[2] = 14; cone[3] = 20; 1965 ierr = DMPlexSetCone(dm, 1, cone);CHKERRQ(ierr); 1966 ornt[0] = 0; ornt[1] = 0; ornt[2] = -1; ornt[3] = 0; 1967 ierr = DMPlexSetConeOrientation(dm, 1, ornt);CHKERRQ(ierr); 1968 /* Cell 2 */ 1969 cone[0] = 21; cone[1] = 22; cone[2] = 18; cone[3] = 23; 1970 ierr = DMPlexSetCone(dm, 2, cone);CHKERRQ(ierr); 1971 ornt[0] = 0; ornt[1] = 0; ornt[2] = -1; ornt[3] = 0; 1972 ierr = DMPlexSetConeOrientation(dm, 2, ornt);CHKERRQ(ierr); 1973 /* Cell 3 */ 1974 cone[0] = 19; cone[1] = 22; cone[2] = 24; cone[3] = 15; 1975 ierr = DMPlexSetCone(dm, 3, cone);CHKERRQ(ierr); 1976 ornt[0] = -1; ornt[1] = -1; ornt[2] = 0; ornt[3] = -1; 1977 ierr = DMPlexSetConeOrientation(dm, 3, ornt);CHKERRQ(ierr); 1978 /* Cell 4 */ 1979 cone[0] = 16; cone[1] = 24; cone[2] = 21; cone[3] = 25; 1980 ierr = DMPlexSetCone(dm, 4, cone);CHKERRQ(ierr); 1981 ornt[0] = -1; ornt[1] = -1; ornt[2] = -1; ornt[3] = 0; 1982 ierr = DMPlexSetConeOrientation(dm, 4, ornt);CHKERRQ(ierr); 1983 /* Cell 5 */ 1984 cone[0] = 20; cone[1] = 17; cone[2] = 25; cone[3] = 23; 1985 ierr = DMPlexSetCone(dm, 5, cone);CHKERRQ(ierr); 1986 ornt[0] = -1; ornt[1] = -1; ornt[2] = -1; ornt[3] = -1; 1987 ierr = DMPlexSetConeOrientation(dm, 5, ornt);CHKERRQ(ierr); 1988 /* Edges */ 1989 cone[0] = 6; cone[1] = 7; 1990 ierr = DMPlexSetCone(dm, 14, cone);CHKERRQ(ierr); 1991 cone[0] = 7; cone[1] = 8; 1992 ierr = DMPlexSetCone(dm, 15, cone);CHKERRQ(ierr); 1993 cone[0] = 8; cone[1] = 9; 1994 ierr = DMPlexSetCone(dm, 16, cone);CHKERRQ(ierr); 1995 cone[0] = 9; cone[1] = 6; 1996 ierr = DMPlexSetCone(dm, 17, cone);CHKERRQ(ierr); 1997 cone[0] = 10; cone[1] = 11; 1998 ierr = DMPlexSetCone(dm, 18, cone);CHKERRQ(ierr); 1999 cone[0] = 11; cone[1] = 7; 2000 ierr = DMPlexSetCone(dm, 19, cone);CHKERRQ(ierr); 2001 cone[0] = 6; cone[1] = 10; 2002 ierr = DMPlexSetCone(dm, 20, cone);CHKERRQ(ierr); 2003 cone[0] = 12; cone[1] = 13; 2004 ierr = DMPlexSetCone(dm, 21, cone);CHKERRQ(ierr); 2005 cone[0] = 13; cone[1] = 11; 2006 ierr = DMPlexSetCone(dm, 22, cone);CHKERRQ(ierr); 2007 cone[0] = 10; cone[1] = 12; 2008 ierr = DMPlexSetCone(dm, 23, cone);CHKERRQ(ierr); 2009 cone[0] = 13; cone[1] = 8; 2010 ierr = DMPlexSetCone(dm, 24, cone);CHKERRQ(ierr); 2011 cone[0] = 12; cone[1] = 9; 2012 ierr = DMPlexSetCone(dm, 25, cone);CHKERRQ(ierr); 2013 } 2014 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 2015 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 2016 /* Build coordinates */ 2017 ierr = PetscCalloc1(numVerts * embedDim, &coordsIn);CHKERRQ(ierr); 2018 if (rank == 0) { 2019 coordsIn[0*embedDim+0] = -R; coordsIn[0*embedDim+1] = R; coordsIn[0*embedDim+2] = -R; 2020 coordsIn[1*embedDim+0] = R; coordsIn[1*embedDim+1] = R; coordsIn[1*embedDim+2] = -R; 2021 coordsIn[2*embedDim+0] = R; coordsIn[2*embedDim+1] = -R; coordsIn[2*embedDim+2] = -R; 2022 coordsIn[3*embedDim+0] = -R; coordsIn[3*embedDim+1] = -R; coordsIn[3*embedDim+2] = -R; 2023 coordsIn[4*embedDim+0] = -R; coordsIn[4*embedDim+1] = R; coordsIn[4*embedDim+2] = R; 2024 coordsIn[5*embedDim+0] = R; coordsIn[5*embedDim+1] = R; coordsIn[5*embedDim+2] = R; 2025 coordsIn[6*embedDim+0] = -R; coordsIn[6*embedDim+1] = -R; coordsIn[6*embedDim+2] = R; 2026 coordsIn[7*embedDim+0] = R; coordsIn[7*embedDim+1] = -R; coordsIn[7*embedDim+2] = R; 2027 } 2028 } 2029 break; 2030 case 3: 2031 if (simplex) { 2032 const PetscReal edgeLen = 1.0/PETSC_PHI; 2033 PetscReal vertexA[4] = {0.5, 0.5, 0.5, 0.5}; 2034 PetscReal vertexB[4] = {1.0, 0.0, 0.0, 0.0}; 2035 PetscReal vertexC[4] = {0.5, 0.5*PETSC_PHI, 0.5/PETSC_PHI, 0.0}; 2036 const PetscInt degree = 12; 2037 PetscInt s[4] = {1, 1, 1}; 2038 PetscInt evenPerm[12][4] = {{0, 1, 2, 3}, {0, 2, 3, 1}, {0, 3, 1, 2}, {1, 0, 3, 2}, {1, 2, 0, 3}, {1, 3, 2, 0}, 2039 {2, 0, 1, 3}, {2, 1, 3, 0}, {2, 3, 0, 1}, {3, 0, 2, 1}, {3, 1, 0, 2}, {3, 2, 1, 0}}; 2040 PetscInt cone[4]; 2041 PetscInt *graph, p, i, j, k, l; 2042 2043 vertexA[0] *= R; vertexA[1] *= R; vertexA[2] *= R; vertexA[3] *= R; 2044 vertexB[0] *= R; vertexB[1] *= R; vertexB[2] *= R; vertexB[3] *= R; 2045 vertexC[0] *= R; vertexC[1] *= R; vertexC[2] *= R; vertexC[3] *= R; 2046 numCells = rank == 0 ? 600 : 0; 2047 numVerts = rank == 0 ? 120 : 0; 2048 firstVertex = numCells; 2049 /* Use the 600-cell, which for a unit sphere has coordinates which are 2050 2051 1/2 (\pm 1, \pm 1, \pm 1, \pm 1) 16 2052 (\pm 1, 0, 0, 0) all cyclic permutations 8 2053 1/2 (\pm 1, \pm phi, \pm 1/phi, 0) all even permutations 96 2054 2055 where \phi^2 - \phi - 1 = 0, meaning \phi is the golden ratio \frac{1 + \sqrt{5}}{2}. The edge 2056 length is then given by 1/\phi = 0.61803. 2057 2058 http://buzzard.pugetsound.edu/sage-practice/ch03s03.html 2059 http://mathworld.wolfram.com/600-Cell.html 2060 */ 2061 /* Construct vertices */ 2062 ierr = PetscCalloc1(numVerts * embedDim, &coordsIn);CHKERRQ(ierr); 2063 i = 0; 2064 if (rank == 0) { 2065 for (s[0] = -1; s[0] < 2; s[0] += 2) { 2066 for (s[1] = -1; s[1] < 2; s[1] += 2) { 2067 for (s[2] = -1; s[2] < 2; s[2] += 2) { 2068 for (s[3] = -1; s[3] < 2; s[3] += 2) { 2069 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[d]*vertexA[d]; 2070 ++i; 2071 } 2072 } 2073 } 2074 } 2075 for (p = 0; p < embedDim; ++p) { 2076 s[1] = s[2] = s[3] = 1; 2077 for (s[0] = -1; s[0] < 2; s[0] += 2) { 2078 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[(d+p)%embedDim]*vertexB[(d+p)%embedDim]; 2079 ++i; 2080 } 2081 } 2082 for (p = 0; p < 12; ++p) { 2083 s[3] = 1; 2084 for (s[0] = -1; s[0] < 2; s[0] += 2) { 2085 for (s[1] = -1; s[1] < 2; s[1] += 2) { 2086 for (s[2] = -1; s[2] < 2; s[2] += 2) { 2087 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[evenPerm[p][d]]*vertexC[evenPerm[p][d]]; 2088 ++i; 2089 } 2090 } 2091 } 2092 } 2093 } 2094 PetscCheckFalse(i != numVerts,PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "Invalid 600-cell, vertices %D != %D", i, numVerts); 2095 /* Construct graph */ 2096 ierr = PetscCalloc1(numVerts * numVerts, &graph);CHKERRQ(ierr); 2097 for (i = 0; i < numVerts; ++i) { 2098 for (j = 0, k = 0; j < numVerts; ++j) { 2099 if (PetscAbsReal(DiffNormReal(embedDim, &coordsIn[i*embedDim], &coordsIn[j*embedDim]) - edgeLen) < PETSC_SMALL) {graph[i*numVerts+j] = 1; ++k;} 2100 } 2101 PetscCheckFalse(k != degree,PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "Invalid 600-cell, vertex %D degree %D != %D", i, k, degree); 2102 } 2103 /* Build Topology */ 2104 ierr = DMPlexSetChart(dm, 0, numCells+numVerts);CHKERRQ(ierr); 2105 for (c = 0; c < numCells; c++) { 2106 ierr = DMPlexSetConeSize(dm, c, embedDim);CHKERRQ(ierr); 2107 } 2108 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 2109 /* Cells */ 2110 if (rank == 0) { 2111 for (i = 0, c = 0; i < numVerts; ++i) { 2112 for (j = 0; j < i; ++j) { 2113 for (k = 0; k < j; ++k) { 2114 for (l = 0; l < k; ++l) { 2115 if (graph[i*numVerts+j] && graph[j*numVerts+k] && graph[k*numVerts+i] && 2116 graph[l*numVerts+i] && graph[l*numVerts+j] && graph[l*numVerts+k]) { 2117 cone[0] = firstVertex+i; cone[1] = firstVertex+j; cone[2] = firstVertex+k; cone[3] = firstVertex+l; 2118 /* Check orientation: https://ef.gy/linear-algebra:normal-vectors-in-higher-dimensional-spaces */ 2119 { 2120 const PetscInt epsilon[4][4][4][4] = {{{{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 2121 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 1}, { 0, 0, -1, 0}}, 2122 {{0, 0, 0, 0}, { 0, 0, 0, -1}, { 0, 0, 0, 0}, { 0, 1, 0, 0}}, 2123 {{0, 0, 0, 0}, { 0, 0, 1, 0}, { 0, -1, 0, 0}, { 0, 0, 0, 0}}}, 2124 2125 {{{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, -1}, { 0, 0, 1, 0}}, 2126 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 2127 {{0, 0, 0, 1}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, {-1, 0, 0, 0}}, 2128 {{0, 0, -1, 0}, { 0, 0, 0, 0}, { 1, 0, 0, 0}, { 0, 0, 0, 0}}}, 2129 2130 {{{0, 0, 0, 0}, { 0, 0, 0, 1}, { 0, 0, 0, 0}, { 0, -1, 0, 0}}, 2131 {{0, 0, 0, -1}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 1, 0, 0, 0}}, 2132 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 2133 {{0, 1, 0, 0}, {-1, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}}, 2134 2135 {{{0, 0, 0, 0}, { 0, 0, -1, 0}, { 0, 1, 0, 0}, { 0, 0, 0, 0}}, 2136 {{0, 0, 1, 0}, { 0, 0, 0, 0}, {-1, 0, 0, 0}, { 0, 0, 0, 0}}, 2137 {{0, -1, 0, 0}, { 1, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}, 2138 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}}}; 2139 PetscReal normal[4]; 2140 PetscInt e, f, g; 2141 2142 for (d = 0; d < embedDim; ++d) { 2143 normal[d] = 0.0; 2144 for (e = 0; e < embedDim; ++e) { 2145 for (f = 0; f < embedDim; ++f) { 2146 for (g = 0; g < embedDim; ++g) { 2147 normal[d] += epsilon[d][e][f][g]*(coordsIn[j*embedDim+e] - coordsIn[i*embedDim+e])*(coordsIn[k*embedDim+f] - coordsIn[i*embedDim+f])*(coordsIn[l*embedDim+f] - coordsIn[i*embedDim+f]); 2148 } 2149 } 2150 } 2151 } 2152 if (DotReal(embedDim, normal, &coordsIn[i*embedDim]) < 0) {PetscInt tmp = cone[1]; cone[1] = cone[2]; cone[2] = tmp;} 2153 } 2154 ierr = DMPlexSetCone(dm, c++, cone);CHKERRQ(ierr); 2155 } 2156 } 2157 } 2158 } 2159 } 2160 } 2161 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 2162 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 2163 ierr = PetscFree(graph);CHKERRQ(ierr); 2164 break; 2165 } 2166 default: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Unsupported dimension for sphere: %D", dim); 2167 } 2168 /* Create coordinates */ 2169 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2170 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 2171 ierr = PetscSectionSetFieldComponents(coordSection, 0, embedDim);CHKERRQ(ierr); 2172 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVerts);CHKERRQ(ierr); 2173 for (v = firstVertex; v < firstVertex+numVerts; ++v) { 2174 ierr = PetscSectionSetDof(coordSection, v, embedDim);CHKERRQ(ierr); 2175 ierr = PetscSectionSetFieldDof(coordSection, v, 0, embedDim);CHKERRQ(ierr); 2176 } 2177 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 2178 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 2179 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 2180 ierr = VecSetBlockSize(coordinates, embedDim);CHKERRQ(ierr); 2181 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 2182 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 2183 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 2184 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2185 for (v = 0; v < numVerts; ++v) for (d = 0; d < embedDim; ++d) {coords[v*embedDim+d] = coordsIn[v*embedDim+d];} 2186 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2187 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 2188 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 2189 ierr = PetscFree(coordsIn);CHKERRQ(ierr); 2190 { 2191 DM cdm; 2192 PetscDS cds; 2193 PetscScalar c = R; 2194 2195 ierr = DMPlexCreateCoordinateSpace(dm, 1, snapToSphere);CHKERRQ(ierr); 2196 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 2197 ierr = DMGetDS(cdm, &cds);CHKERRQ(ierr); 2198 ierr = PetscDSSetConstants(cds, 1, &c);CHKERRQ(ierr); 2199 } 2200 /* Wait for coordinate creation before doing in-place modification */ 2201 if (simplex) {ierr = DMPlexInterpolateInPlace_Internal(dm);CHKERRQ(ierr);} 2202 PetscFunctionReturn(0); 2203 } 2204 2205 /*@ 2206 DMPlexCreateSphereMesh - Creates a mesh on the d-dimensional sphere, S^d. 2207 2208 Collective 2209 2210 Input Parameters: 2211 + comm - The communicator for the DM object 2212 . dim - The dimension 2213 . simplex - Use simplices, or tensor product cells 2214 - R - The radius 2215 2216 Output Parameter: 2217 . dm - The DM object 2218 2219 Level: beginner 2220 2221 .seealso: DMPlexCreateBallMesh(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 2222 @*/ 2223 PetscErrorCode DMPlexCreateSphereMesh(MPI_Comm comm, PetscInt dim, PetscBool simplex, PetscReal R, DM *dm) 2224 { 2225 PetscErrorCode ierr; 2226 2227 PetscFunctionBegin; 2228 PetscValidPointer(dm, 5); 2229 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 2230 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 2231 ierr = DMPlexCreateSphereMesh_Internal(*dm, dim, simplex, R);CHKERRQ(ierr); 2232 PetscFunctionReturn(0); 2233 } 2234 2235 static PetscErrorCode DMPlexCreateBallMesh_Internal(DM dm, PetscInt dim, PetscReal R) 2236 { 2237 DM sdm, vol; 2238 DMLabel bdlabel; 2239 PetscErrorCode ierr; 2240 2241 PetscFunctionBegin; 2242 ierr = DMCreate(PetscObjectComm((PetscObject) dm), &sdm);CHKERRQ(ierr); 2243 ierr = DMSetType(sdm, DMPLEX);CHKERRQ(ierr); 2244 ierr = PetscObjectSetOptionsPrefix((PetscObject) sdm, "bd_");CHKERRQ(ierr); 2245 ierr = DMPlexCreateSphereMesh_Internal(sdm, dim-1, PETSC_TRUE, R);CHKERRQ(ierr); 2246 ierr = DMSetFromOptions(sdm);CHKERRQ(ierr); 2247 ierr = DMViewFromOptions(sdm, NULL, "-dm_view");CHKERRQ(ierr); 2248 ierr = DMPlexGenerate(sdm, NULL, PETSC_TRUE, &vol);CHKERRQ(ierr); 2249 ierr = DMDestroy(&sdm);CHKERRQ(ierr); 2250 ierr = DMPlexReplace_Static(dm, &vol);CHKERRQ(ierr); 2251 ierr = DMCreateLabel(dm, "marker");CHKERRQ(ierr); 2252 ierr = DMGetLabel(dm, "marker", &bdlabel);CHKERRQ(ierr); 2253 ierr = DMPlexMarkBoundaryFaces(dm, PETSC_DETERMINE, bdlabel);CHKERRQ(ierr); 2254 ierr = DMPlexLabelComplete(dm, bdlabel);CHKERRQ(ierr); 2255 PetscFunctionReturn(0); 2256 } 2257 2258 /*@ 2259 DMPlexCreateBallMesh - Creates a simplex mesh on the d-dimensional ball, B^d. 2260 2261 Collective 2262 2263 Input Parameters: 2264 + comm - The communicator for the DM object 2265 . dim - The dimension 2266 - R - The radius 2267 2268 Output Parameter: 2269 . dm - The DM object 2270 2271 Options Database Keys: 2272 - bd_dm_refine - This will refine the surface mesh preserving the sphere geometry 2273 2274 Level: beginner 2275 2276 .seealso: DMPlexCreateSphereMesh(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate() 2277 @*/ 2278 PetscErrorCode DMPlexCreateBallMesh(MPI_Comm comm, PetscInt dim, PetscReal R, DM *dm) 2279 { 2280 PetscErrorCode ierr; 2281 2282 PetscFunctionBegin; 2283 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 2284 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 2285 ierr = DMPlexCreateBallMesh_Internal(*dm, dim, R);CHKERRQ(ierr); 2286 PetscFunctionReturn(0); 2287 } 2288 2289 static PetscErrorCode DMPlexCreateReferenceCell_Internal(DM rdm, DMPolytopeType ct) 2290 { 2291 PetscErrorCode ierr; 2292 2293 PetscFunctionBegin; 2294 switch (ct) { 2295 case DM_POLYTOPE_POINT: 2296 { 2297 PetscInt numPoints[1] = {1}; 2298 PetscInt coneSize[1] = {0}; 2299 PetscInt cones[1] = {0}; 2300 PetscInt coneOrientations[1] = {0}; 2301 PetscScalar vertexCoords[1] = {0.0}; 2302 2303 ierr = DMSetDimension(rdm, 0);CHKERRQ(ierr); 2304 ierr = DMPlexCreateFromDAG(rdm, 0, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2305 } 2306 break; 2307 case DM_POLYTOPE_SEGMENT: 2308 { 2309 PetscInt numPoints[2] = {2, 1}; 2310 PetscInt coneSize[3] = {2, 0, 0}; 2311 PetscInt cones[2] = {1, 2}; 2312 PetscInt coneOrientations[2] = {0, 0}; 2313 PetscScalar vertexCoords[2] = {-1.0, 1.0}; 2314 2315 ierr = DMSetDimension(rdm, 1);CHKERRQ(ierr); 2316 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2317 } 2318 break; 2319 case DM_POLYTOPE_POINT_PRISM_TENSOR: 2320 { 2321 PetscInt numPoints[2] = {2, 1}; 2322 PetscInt coneSize[3] = {2, 0, 0}; 2323 PetscInt cones[2] = {1, 2}; 2324 PetscInt coneOrientations[2] = {0, 0}; 2325 PetscScalar vertexCoords[2] = {-1.0, 1.0}; 2326 2327 ierr = DMSetDimension(rdm, 1);CHKERRQ(ierr); 2328 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2329 } 2330 break; 2331 case DM_POLYTOPE_TRIANGLE: 2332 { 2333 PetscInt numPoints[2] = {3, 1}; 2334 PetscInt coneSize[4] = {3, 0, 0, 0}; 2335 PetscInt cones[3] = {1, 2, 3}; 2336 PetscInt coneOrientations[3] = {0, 0, 0}; 2337 PetscScalar vertexCoords[6] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0}; 2338 2339 ierr = DMSetDimension(rdm, 2);CHKERRQ(ierr); 2340 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2341 } 2342 break; 2343 case DM_POLYTOPE_QUADRILATERAL: 2344 { 2345 PetscInt numPoints[2] = {4, 1}; 2346 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 2347 PetscInt cones[4] = {1, 2, 3, 4}; 2348 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 2349 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0}; 2350 2351 ierr = DMSetDimension(rdm, 2);CHKERRQ(ierr); 2352 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2353 } 2354 break; 2355 case DM_POLYTOPE_SEG_PRISM_TENSOR: 2356 { 2357 PetscInt numPoints[2] = {4, 1}; 2358 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 2359 PetscInt cones[4] = {1, 2, 3, 4}; 2360 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 2361 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0}; 2362 2363 ierr = DMSetDimension(rdm, 2);CHKERRQ(ierr); 2364 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2365 } 2366 break; 2367 case DM_POLYTOPE_TETRAHEDRON: 2368 { 2369 PetscInt numPoints[2] = {4, 1}; 2370 PetscInt coneSize[5] = {4, 0, 0, 0, 0}; 2371 PetscInt cones[4] = {1, 2, 3, 4}; 2372 PetscInt coneOrientations[4] = {0, 0, 0, 0}; 2373 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}; 2374 2375 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 2376 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2377 } 2378 break; 2379 case DM_POLYTOPE_HEXAHEDRON: 2380 { 2381 PetscInt numPoints[2] = {8, 1}; 2382 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0}; 2383 PetscInt cones[8] = {1, 2, 3, 4, 5, 6, 7, 8}; 2384 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 2385 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, 2386 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0}; 2387 2388 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 2389 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2390 } 2391 break; 2392 case DM_POLYTOPE_TRI_PRISM: 2393 { 2394 PetscInt numPoints[2] = {6, 1}; 2395 PetscInt coneSize[7] = {6, 0, 0, 0, 0, 0, 0}; 2396 PetscInt cones[6] = {1, 2, 3, 4, 5, 6}; 2397 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 2398 PetscScalar vertexCoords[18] = {-1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 2399 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0}; 2400 2401 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 2402 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2403 } 2404 break; 2405 case DM_POLYTOPE_TRI_PRISM_TENSOR: 2406 { 2407 PetscInt numPoints[2] = {6, 1}; 2408 PetscInt coneSize[7] = {6, 0, 0, 0, 0, 0, 0}; 2409 PetscInt cones[6] = {1, 2, 3, 4, 5, 6}; 2410 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 2411 PetscScalar vertexCoords[18] = {-1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 2412 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0}; 2413 2414 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 2415 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2416 } 2417 break; 2418 case DM_POLYTOPE_QUAD_PRISM_TENSOR: 2419 { 2420 PetscInt numPoints[2] = {8, 1}; 2421 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0}; 2422 PetscInt cones[8] = {1, 2, 3, 4, 5, 6, 7, 8}; 2423 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 2424 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, 2425 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0}; 2426 2427 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 2428 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2429 } 2430 break; 2431 case DM_POLYTOPE_PYRAMID: 2432 { 2433 PetscInt numPoints[2] = {5, 1}; 2434 PetscInt coneSize[6] = {5, 0, 0, 0, 0, 0}; 2435 PetscInt cones[5] = {1, 2, 3, 4, 5}; 2436 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 2437 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, 2438 0.0, 0.0, 1.0}; 2439 2440 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr); 2441 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 2442 } 2443 break; 2444 default: SETERRQ(PetscObjectComm((PetscObject) rdm), PETSC_ERR_ARG_WRONG, "Cannot create reference cell for cell type %s", DMPolytopeTypes[ct]); 2445 } 2446 { 2447 PetscInt Nv, v; 2448 2449 /* Must create the celltype label here so that we do not automatically try to compute the types */ 2450 ierr = DMCreateLabel(rdm, "celltype");CHKERRQ(ierr); 2451 ierr = DMPlexSetCellType(rdm, 0, ct);CHKERRQ(ierr); 2452 ierr = DMPlexGetChart(rdm, NULL, &Nv);CHKERRQ(ierr); 2453 for (v = 1; v < Nv; ++v) {ierr = DMPlexSetCellType(rdm, v, DM_POLYTOPE_POINT);CHKERRQ(ierr);} 2454 } 2455 ierr = DMPlexInterpolateInPlace_Internal(rdm);CHKERRQ(ierr); 2456 ierr = PetscObjectSetName((PetscObject) rdm, DMPolytopeTypes[ct]);CHKERRQ(ierr); 2457 PetscFunctionReturn(0); 2458 } 2459 2460 /*@ 2461 DMPlexCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell 2462 2463 Collective 2464 2465 Input Parameters: 2466 + comm - The communicator 2467 - ct - The cell type of the reference cell 2468 2469 Output Parameter: 2470 . refdm - The reference cell 2471 2472 Level: intermediate 2473 2474 .seealso: DMPlexCreateReferenceCell(), DMPlexCreateBoxMesh() 2475 @*/ 2476 PetscErrorCode DMPlexCreateReferenceCell(MPI_Comm comm, DMPolytopeType ct, DM *refdm) 2477 { 2478 PetscErrorCode ierr; 2479 2480 PetscFunctionBegin; 2481 ierr = DMCreate(comm, refdm);CHKERRQ(ierr); 2482 ierr = DMSetType(*refdm, DMPLEX);CHKERRQ(ierr); 2483 ierr = DMPlexCreateReferenceCell_Internal(*refdm, ct);CHKERRQ(ierr); 2484 PetscFunctionReturn(0); 2485 } 2486 2487 static PetscErrorCode DMPlexCreateBoundaryLabel_Private(DM dm, const char name[]) 2488 { 2489 DM plex; 2490 DMLabel label; 2491 PetscBool hasLabel; 2492 PetscErrorCode ierr; 2493 2494 PetscFunctionBeginUser; 2495 ierr = DMHasLabel(dm, name, &hasLabel);CHKERRQ(ierr); 2496 if (hasLabel) PetscFunctionReturn(0); 2497 ierr = DMCreateLabel(dm, name);CHKERRQ(ierr); 2498 ierr = DMGetLabel(dm, name, &label);CHKERRQ(ierr); 2499 ierr = DMConvert(dm, DMPLEX, &plex);CHKERRQ(ierr); 2500 ierr = DMPlexMarkBoundaryFaces(plex, 1, label);CHKERRQ(ierr); 2501 ierr = DMDestroy(&plex);CHKERRQ(ierr); 2502 PetscFunctionReturn(0); 2503 } 2504 2505 const char * const DMPlexShapes[] = {"box", "box_surface", "ball", "sphere", "cylinder", "unknown", "DMPlexShape", "DM_SHAPE_", NULL}; 2506 2507 static PetscErrorCode DMPlexCreateFromOptions_Internal(PetscOptionItems *PetscOptionsObject, PetscBool *useCoordSpace, DM dm) 2508 { 2509 DMPlexShape shape = DM_SHAPE_BOX; 2510 DMPolytopeType cell = DM_POLYTOPE_TRIANGLE; 2511 PetscInt dim = 2; 2512 PetscBool simplex = PETSC_TRUE, interpolate = PETSC_TRUE, adjCone = PETSC_FALSE, adjClosure = PETSC_TRUE, refDomain = PETSC_FALSE; 2513 PetscBool flg, flg2, fflg, bdfflg, nameflg; 2514 MPI_Comm comm; 2515 char filename[PETSC_MAX_PATH_LEN] = "<unspecified>"; 2516 char bdFilename[PETSC_MAX_PATH_LEN] = "<unspecified>"; 2517 char plexname[PETSC_MAX_PATH_LEN] = ""; 2518 PetscErrorCode ierr; 2519 2520 PetscFunctionBegin; 2521 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr); 2522 /* TODO Turn this into a registration interface */ 2523 ierr = PetscOptionsString("-dm_plex_filename", "File containing a mesh", "DMPlexCreateFromFile", filename, filename, sizeof(filename), &fflg);CHKERRQ(ierr); 2524 ierr = PetscOptionsString("-dm_plex_boundary_filename", "File containing a mesh boundary", "DMPlexCreateFromFile", bdFilename, bdFilename, sizeof(bdFilename), &bdfflg);CHKERRQ(ierr); 2525 ierr = PetscOptionsString("-dm_plex_name", "Name of the mesh in the file", "DMPlexCreateFromFile", plexname, plexname, sizeof(plexname), &nameflg);CHKERRQ(ierr); 2526 ierr = PetscOptionsEnum("-dm_plex_cell", "Cell shape", "", DMPolytopeTypes, (PetscEnum) cell, (PetscEnum *) &cell, NULL);CHKERRQ(ierr); 2527 ierr = PetscOptionsBool("-dm_plex_reference_cell_domain", "Use a reference cell domain", "", refDomain, &refDomain, NULL);CHKERRQ(ierr); 2528 ierr = PetscOptionsEnum("-dm_plex_shape", "Shape for built-in mesh", "", DMPlexShapes, (PetscEnum) shape, (PetscEnum *) &shape, &flg);CHKERRQ(ierr); 2529 ierr = PetscOptionsBoundedInt("-dm_plex_dim", "Topological dimension of the mesh", "DMGetDimension", dim, &dim, &flg, 0);CHKERRQ(ierr); 2530 PetscCheckFalse((dim < 0) || (dim > 3),comm, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %D should be in [1, 3]", dim); 2531 ierr = PetscOptionsBool("-dm_plex_simplex", "Mesh cell shape", "", simplex, &simplex, &flg);CHKERRQ(ierr); 2532 ierr = PetscOptionsBool("-dm_plex_interpolate", "Flag to create edges and faces automatically", "", interpolate, &interpolate, &flg);CHKERRQ(ierr); 2533 ierr = PetscOptionsBool("-dm_plex_adj_cone", "Set adjacency direction", "DMSetBasicAdjacency", adjCone, &adjCone, &flg);CHKERRQ(ierr); 2534 ierr = PetscOptionsBool("-dm_plex_adj_closure", "Set adjacency size", "DMSetBasicAdjacency", adjClosure, &adjClosure, &flg2);CHKERRQ(ierr); 2535 if (flg || flg2) {ierr = DMSetBasicAdjacency(dm, adjCone, adjClosure);CHKERRQ(ierr);} 2536 2537 switch (cell) { 2538 case DM_POLYTOPE_POINT: 2539 case DM_POLYTOPE_SEGMENT: 2540 case DM_POLYTOPE_POINT_PRISM_TENSOR: 2541 case DM_POLYTOPE_TRIANGLE: 2542 case DM_POLYTOPE_QUADRILATERAL: 2543 case DM_POLYTOPE_TETRAHEDRON: 2544 case DM_POLYTOPE_HEXAHEDRON: 2545 *useCoordSpace = PETSC_TRUE;break; 2546 default: *useCoordSpace = PETSC_FALSE;break; 2547 } 2548 2549 if (fflg) { 2550 DM dmnew; 2551 2552 ierr = DMPlexCreateFromFile(PetscObjectComm((PetscObject) dm), filename, plexname, interpolate, &dmnew);CHKERRQ(ierr); 2553 ierr = DMPlexCopy_Internal(dm, PETSC_FALSE, dmnew);CHKERRQ(ierr); 2554 ierr = DMPlexReplace_Static(dm, &dmnew);CHKERRQ(ierr); 2555 } else if (refDomain) { 2556 ierr = DMPlexCreateReferenceCell_Internal(dm, cell);CHKERRQ(ierr); 2557 } else if (bdfflg) { 2558 DM bdm, dmnew; 2559 2560 ierr = DMPlexCreateFromFile(PetscObjectComm((PetscObject) dm), bdFilename, plexname, interpolate, &bdm);CHKERRQ(ierr); 2561 ierr = PetscObjectSetOptionsPrefix((PetscObject) bdm, "bd_");CHKERRQ(ierr); 2562 ierr = DMSetFromOptions(bdm);CHKERRQ(ierr); 2563 ierr = DMPlexGenerate(bdm, NULL, interpolate, &dmnew);CHKERRQ(ierr); 2564 ierr = DMDestroy(&bdm);CHKERRQ(ierr); 2565 ierr = DMPlexCopy_Internal(dm, PETSC_FALSE, dmnew);CHKERRQ(ierr); 2566 ierr = DMPlexReplace_Static(dm, &dmnew);CHKERRQ(ierr); 2567 } else { 2568 ierr = PetscObjectSetName((PetscObject) dm, DMPlexShapes[shape]);CHKERRQ(ierr); 2569 switch (shape) { 2570 case DM_SHAPE_BOX: 2571 { 2572 PetscInt faces[3] = {0, 0, 0}; 2573 PetscReal lower[3] = {0, 0, 0}; 2574 PetscReal upper[3] = {1, 1, 1}; 2575 DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE}; 2576 PetscInt i, n; 2577 2578 n = dim; 2579 for (i = 0; i < dim; ++i) faces[i] = (dim == 1 ? 1 : 4-dim); 2580 ierr = PetscOptionsIntArray("-dm_plex_box_faces", "Number of faces along each dimension", "", faces, &n, &flg);CHKERRQ(ierr); 2581 n = 3; 2582 ierr = PetscOptionsRealArray("-dm_plex_box_lower", "Lower left corner of box", "", lower, &n, &flg);CHKERRQ(ierr); 2583 PetscCheckFalse(flg && (n != dim),comm, PETSC_ERR_ARG_SIZ, "Lower box point had %D values, should have been %D", n, dim); 2584 n = 3; 2585 ierr = PetscOptionsRealArray("-dm_plex_box_upper", "Upper right corner of box", "", upper, &n, &flg);CHKERRQ(ierr); 2586 PetscCheckFalse(flg && (n != dim),comm, PETSC_ERR_ARG_SIZ, "Upper box point had %D values, should have been %D", n, dim); 2587 n = 3; 2588 ierr = PetscOptionsEnumArray("-dm_plex_box_bd", "Boundary type for each dimension", "", DMBoundaryTypes, (PetscEnum *) bdt, &n, &flg);CHKERRQ(ierr); 2589 PetscCheckFalse(flg && (n != dim),comm, PETSC_ERR_ARG_SIZ, "Box boundary types had %D values, should have been %D", n, dim); 2590 switch (cell) { 2591 case DM_POLYTOPE_TRI_PRISM_TENSOR: 2592 ierr = DMPlexCreateWedgeBoxMesh_Internal(dm, faces, lower, upper, bdt);CHKERRQ(ierr); 2593 if (!interpolate) { 2594 DM udm; 2595 2596 ierr = DMPlexUninterpolate(dm, &udm);CHKERRQ(ierr); 2597 ierr = DMPlexReplace_Static(dm, &udm);CHKERRQ(ierr); 2598 } 2599 break; 2600 default: 2601 ierr = DMPlexCreateBoxMesh_Internal(dm, dim, simplex, faces, lower, upper, bdt, interpolate);CHKERRQ(ierr); 2602 break; 2603 } 2604 } 2605 break; 2606 case DM_SHAPE_BOX_SURFACE: 2607 { 2608 PetscInt faces[3] = {0, 0, 0}; 2609 PetscReal lower[3] = {0, 0, 0}; 2610 PetscReal upper[3] = {1, 1, 1}; 2611 PetscInt i, n; 2612 2613 n = dim+1; 2614 for (i = 0; i < dim+1; ++i) faces[i] = (dim+1 == 1 ? 1 : 4-(dim+1)); 2615 ierr = PetscOptionsIntArray("-dm_plex_box_faces", "Number of faces along each dimension", "", faces, &n, &flg);CHKERRQ(ierr); 2616 n = 3; 2617 ierr = PetscOptionsRealArray("-dm_plex_box_lower", "Lower left corner of box", "", lower, &n, &flg);CHKERRQ(ierr); 2618 PetscCheckFalse(flg && (n != dim+1),comm, PETSC_ERR_ARG_SIZ, "Lower box point had %D values, should have been %D", n, dim+1); 2619 n = 3; 2620 ierr = PetscOptionsRealArray("-dm_plex_box_upper", "Upper right corner of box", "", upper, &n, &flg);CHKERRQ(ierr); 2621 PetscCheckFalse(flg && (n != dim+1),comm, PETSC_ERR_ARG_SIZ, "Upper box point had %D values, should have been %D", n, dim+1); 2622 ierr = DMPlexCreateBoxSurfaceMesh_Internal(dm, dim+1, faces, lower, upper, interpolate);CHKERRQ(ierr); 2623 } 2624 break; 2625 case DM_SHAPE_SPHERE: 2626 { 2627 PetscReal R = 1.0; 2628 2629 ierr = PetscOptionsReal("-dm_plex_sphere_radius", "Radius of the sphere", "", R, &R, &flg);CHKERRQ(ierr); 2630 ierr = DMPlexCreateSphereMesh_Internal(dm, dim, simplex, R);CHKERRQ(ierr); 2631 } 2632 break; 2633 case DM_SHAPE_BALL: 2634 { 2635 PetscReal R = 1.0; 2636 2637 ierr = PetscOptionsReal("-dm_plex_ball_radius", "Radius of the ball", "", R, &R, &flg);CHKERRQ(ierr); 2638 ierr = DMPlexCreateBallMesh_Internal(dm, dim, R);CHKERRQ(ierr); 2639 } 2640 break; 2641 case DM_SHAPE_CYLINDER: 2642 { 2643 DMBoundaryType bdt = DM_BOUNDARY_NONE; 2644 PetscInt Nw = 6; 2645 2646 ierr = PetscOptionsEnum("-dm_plex_cylinder_bd", "Boundary type in the z direction", "", DMBoundaryTypes, (PetscEnum) bdt, (PetscEnum *) &bdt, NULL);CHKERRQ(ierr); 2647 ierr = PetscOptionsInt("-dm_plex_cylinder_num_wedges", "Number of wedges around the cylinder", "", Nw, &Nw, NULL);CHKERRQ(ierr); 2648 switch (cell) { 2649 case DM_POLYTOPE_TRI_PRISM_TENSOR: 2650 ierr = DMPlexCreateWedgeCylinderMesh_Internal(dm, Nw, interpolate);CHKERRQ(ierr); 2651 break; 2652 default: 2653 ierr = DMPlexCreateHexCylinderMesh_Internal(dm, bdt);CHKERRQ(ierr); 2654 break; 2655 } 2656 } 2657 break; 2658 default: SETERRQ(comm, PETSC_ERR_SUP, "Domain shape %s is unsupported", DMPlexShapes[shape]); 2659 } 2660 } 2661 ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE);CHKERRQ(ierr); 2662 if (!((PetscObject)dm)->name && nameflg) { 2663 ierr = PetscObjectSetName((PetscObject)dm, plexname);CHKERRQ(ierr); 2664 } 2665 PetscFunctionReturn(0); 2666 } 2667 2668 PetscErrorCode DMSetFromOptions_NonRefinement_Plex(PetscOptionItems *PetscOptionsObject, DM dm) 2669 { 2670 DM_Plex *mesh = (DM_Plex*) dm->data; 2671 PetscBool flg; 2672 char bdLabel[PETSC_MAX_PATH_LEN]; 2673 PetscErrorCode ierr; 2674 2675 PetscFunctionBegin; 2676 /* Handle viewing */ 2677 ierr = PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMPlexMatSetClosure", PETSC_FALSE, &mesh->printSetValues, NULL);CHKERRQ(ierr); 2678 ierr = PetscOptionsBoundedInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMPlexSNESComputeResidualFEM", 0, &mesh->printFEM, NULL,0);CHKERRQ(ierr); 2679 ierr = PetscOptionsReal("-dm_plex_print_tol", "Tolerance for FEM output", "DMPlexSNESComputeResidualFEM", mesh->printTol, &mesh->printTol, NULL);CHKERRQ(ierr); 2680 ierr = PetscOptionsBoundedInt("-dm_plex_print_l2", "Debug output level all L2 diff computations", "DMComputeL2Diff", 0, &mesh->printL2, NULL,0);CHKERRQ(ierr); 2681 ierr = DMMonitorSetFromOptions(dm, "-dm_plex_monitor_throughput", "Monitor the simulation throughput", "DMPlexMonitorThroughput", DMPlexMonitorThroughput, NULL, &flg);CHKERRQ(ierr); 2682 if (flg) {ierr = PetscLogDefaultBegin();CHKERRQ(ierr);} 2683 /* Labeling */ 2684 ierr = PetscOptionsString("-dm_plex_boundary_label", "Label to mark the mesh boundary", "", bdLabel, bdLabel, sizeof(bdLabel), &flg);CHKERRQ(ierr); 2685 if (flg) {ierr = DMPlexCreateBoundaryLabel_Private(dm, bdLabel);CHKERRQ(ierr);} 2686 /* Point Location */ 2687 ierr = PetscOptionsBool("-dm_plex_hash_location", "Use grid hashing for point location", "DMInterpolate", PETSC_FALSE, &mesh->useHashLocation, NULL);CHKERRQ(ierr); 2688 /* Partitioning and distribution */ 2689 ierr = PetscOptionsBool("-dm_plex_partition_balance", "Attempt to evenly divide points on partition boundary between processes", "DMPlexSetPartitionBalance", PETSC_FALSE, &mesh->partitionBalance, NULL);CHKERRQ(ierr); 2690 /* Generation and remeshing */ 2691 ierr = PetscOptionsBool("-dm_plex_remesh_bd", "Allow changes to the boundary on remeshing", "DMAdapt", PETSC_FALSE, &mesh->remeshBd, NULL);CHKERRQ(ierr); 2692 /* Projection behavior */ 2693 ierr = PetscOptionsBoundedInt("-dm_plex_max_projection_height", "Maxmimum mesh point height used to project locally", "DMPlexSetMaxProjectionHeight", 0, &mesh->maxProjectionHeight, NULL,0);CHKERRQ(ierr); 2694 ierr = PetscOptionsBool("-dm_plex_regular_refinement", "Use special nested projection algorithm for regular refinement", "DMPlexSetRegularRefinement", mesh->regularRefinement, &mesh->regularRefinement, NULL);CHKERRQ(ierr); 2695 /* Checking structure */ 2696 { 2697 PetscBool flg = PETSC_FALSE, flg2 = PETSC_FALSE, all = PETSC_FALSE; 2698 2699 ierr = PetscOptionsBool("-dm_plex_check_all", "Perform all checks", NULL, PETSC_FALSE, &all, &flg2);CHKERRQ(ierr); 2700 ierr = PetscOptionsBool("-dm_plex_check_symmetry", "Check that the adjacency information in the mesh is symmetric", "DMPlexCheckSymmetry", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr); 2701 if (all || (flg && flg2)) {ierr = DMPlexCheckSymmetry(dm);CHKERRQ(ierr);} 2702 ierr = PetscOptionsBool("-dm_plex_check_skeleton", "Check that each cell has the correct number of vertices (only for homogeneous simplex or tensor meshes)", "DMPlexCheckSkeleton", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr); 2703 if (all || (flg && flg2)) {ierr = DMPlexCheckSkeleton(dm, 0);CHKERRQ(ierr);} 2704 ierr = PetscOptionsBool("-dm_plex_check_faces", "Check that the faces of each cell give a vertex order this is consistent with what we expect from the cell type", "DMPlexCheckFaces", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr); 2705 if (all || (flg && flg2)) {ierr = DMPlexCheckFaces(dm, 0);CHKERRQ(ierr);} 2706 ierr = PetscOptionsBool("-dm_plex_check_geometry", "Check that cells have positive volume", "DMPlexCheckGeometry", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr); 2707 if (all || (flg && flg2)) {ierr = DMPlexCheckGeometry(dm);CHKERRQ(ierr);} 2708 ierr = PetscOptionsBool("-dm_plex_check_pointsf", "Check some necessary conditions for PointSF", "DMPlexCheckPointSF", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr); 2709 if (all || (flg && flg2)) {ierr = DMPlexCheckPointSF(dm);CHKERRQ(ierr);} 2710 ierr = PetscOptionsBool("-dm_plex_check_interface_cones", "Check points on inter-partition interfaces have conforming order of cone points", "DMPlexCheckInterfaceCones", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr); 2711 if (all || (flg && flg2)) {ierr = DMPlexCheckInterfaceCones(dm);CHKERRQ(ierr);} 2712 ierr = PetscOptionsBool("-dm_plex_check_cell_shape", "Check cell shape", "DMPlexCheckCellShape", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr); 2713 if (flg && flg2) {ierr = DMPlexCheckCellShape(dm, PETSC_TRUE, PETSC_DETERMINE);CHKERRQ(ierr);} 2714 } 2715 { 2716 PetscReal scale = 1.0; 2717 2718 ierr = PetscOptionsReal("-dm_plex_scale", "Scale factor for mesh coordinates", "DMPlexScale", scale, &scale, &flg);CHKERRQ(ierr); 2719 if (flg) { 2720 Vec coordinates, coordinatesLocal; 2721 2722 ierr = DMGetCoordinates(dm, &coordinates);CHKERRQ(ierr); 2723 ierr = DMGetCoordinatesLocal(dm, &coordinatesLocal);CHKERRQ(ierr); 2724 ierr = VecScale(coordinates, scale);CHKERRQ(ierr); 2725 ierr = VecScale(coordinatesLocal, scale);CHKERRQ(ierr); 2726 } 2727 } 2728 ierr = PetscPartitionerSetFromOptions(mesh->partitioner);CHKERRQ(ierr); 2729 PetscFunctionReturn(0); 2730 } 2731 2732 static PetscErrorCode DMSetFromOptions_Plex(PetscOptionItems *PetscOptionsObject,DM dm) 2733 { 2734 PetscFunctionList ordlist; 2735 char oname[256]; 2736 PetscReal volume = -1.0; 2737 PetscInt prerefine = 0, refine = 0, r, coarsen = 0, overlap = 0, extLayers = 0, dim; 2738 PetscBool uniformOrig, created = PETSC_FALSE, uniform = PETSC_TRUE, distribute, interpolate = PETSC_TRUE, coordSpace = PETSC_TRUE, remap = PETSC_TRUE, ghostCells = PETSC_FALSE, isHierarchy, ignoreModel = PETSC_FALSE, flg; 2739 PetscErrorCode ierr; 2740 2741 PetscFunctionBegin; 2742 PetscValidHeaderSpecific(dm, DM_CLASSID, 2); 2743 ierr = PetscOptionsHead(PetscOptionsObject,"DMPlex Options");CHKERRQ(ierr); 2744 /* Handle automatic creation */ 2745 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2746 if (dim < 0) {ierr = DMPlexCreateFromOptions_Internal(PetscOptionsObject, &coordSpace, dm);CHKERRQ(ierr);created = PETSC_TRUE;} 2747 /* Handle interpolation before distribution */ 2748 ierr = PetscOptionsBool("-dm_plex_interpolate_pre", "Flag to interpolate mesh before distribution", "", interpolate, &interpolate, &flg);CHKERRQ(ierr); 2749 if (flg) { 2750 DMPlexInterpolatedFlag interpolated; 2751 2752 ierr = DMPlexIsInterpolated(dm, &interpolated);CHKERRQ(ierr); 2753 if (interpolated == DMPLEX_INTERPOLATED_FULL && !interpolate) { 2754 DM udm; 2755 2756 ierr = DMPlexUninterpolate(dm, &udm);CHKERRQ(ierr); 2757 ierr = DMPlexReplace_Static(dm, &udm);CHKERRQ(ierr); 2758 } else if (interpolated != DMPLEX_INTERPOLATED_FULL && interpolate) { 2759 DM idm; 2760 2761 ierr = DMPlexInterpolate(dm, &idm);CHKERRQ(ierr); 2762 ierr = DMPlexReplace_Static(dm, &idm);CHKERRQ(ierr); 2763 } 2764 } 2765 /* Handle DMPlex refinement before distribution */ 2766 ierr = PetscOptionsBool("-dm_refine_ignore_model", "Flag to ignore the geometry model when refining", "DMCreate", ignoreModel, &ignoreModel, &flg);CHKERRQ(ierr); 2767 if (flg) {((DM_Plex *) dm->data)->ignoreModel = ignoreModel;} 2768 ierr = DMPlexGetRefinementUniform(dm, &uniformOrig);CHKERRQ(ierr); 2769 ierr = PetscOptionsBoundedInt("-dm_refine_pre", "The number of refinements before distribution", "DMCreate", prerefine, &prerefine, NULL,0);CHKERRQ(ierr); 2770 ierr = PetscOptionsBool("-dm_refine_remap_pre", "Flag to control coordinate remapping", "DMCreate", remap, &remap, NULL);CHKERRQ(ierr); 2771 ierr = PetscOptionsBool("-dm_refine_uniform_pre", "Flag for uniform refinement before distribution", "DMCreate", uniform, &uniform, &flg);CHKERRQ(ierr); 2772 if (flg) {ierr = DMPlexSetRefinementUniform(dm, uniform);CHKERRQ(ierr);} 2773 ierr = PetscOptionsReal("-dm_refine_volume_limit_pre", "The maximum cell volume after refinement before distribution", "DMCreate", volume, &volume, &flg);CHKERRQ(ierr); 2774 if (flg) { 2775 ierr = DMPlexSetRefinementUniform(dm, PETSC_FALSE);CHKERRQ(ierr); 2776 ierr = DMPlexSetRefinementLimit(dm, volume);CHKERRQ(ierr); 2777 prerefine = PetscMax(prerefine, 1); 2778 } 2779 for (r = 0; r < prerefine; ++r) { 2780 DM rdm; 2781 PetscPointFunc coordFunc = ((DM_Plex*) dm->data)->coordFunc; 2782 2783 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2784 ierr = DMRefine(dm, PetscObjectComm((PetscObject) dm), &rdm);CHKERRQ(ierr); 2785 ierr = DMPlexReplace_Static(dm, &rdm);CHKERRQ(ierr); 2786 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2787 if (coordFunc && remap) { 2788 ierr = DMPlexRemapGeometry(dm, 0.0, coordFunc);CHKERRQ(ierr); 2789 ((DM_Plex*) dm->data)->coordFunc = coordFunc; 2790 } 2791 } 2792 ierr = DMPlexSetRefinementUniform(dm, uniformOrig);CHKERRQ(ierr); 2793 /* Handle DMPlex extrusion before distribution */ 2794 ierr = PetscOptionsBoundedInt("-dm_extrude", "The number of layers to extrude", "", extLayers, &extLayers, NULL, 0);CHKERRQ(ierr); 2795 if (extLayers) { 2796 DM edm; 2797 2798 ierr = DMExtrude(dm, extLayers, &edm);CHKERRQ(ierr); 2799 ierr = DMPlexReplace_Static(dm, &edm);CHKERRQ(ierr); 2800 ((DM_Plex *) dm->data)->coordFunc = NULL; 2801 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2802 extLayers = 0; 2803 } 2804 /* Handle DMPlex reordering before distribution */ 2805 ierr = MatGetOrderingList(&ordlist);CHKERRQ(ierr); 2806 ierr = PetscOptionsFList("-dm_plex_reorder", "Set mesh reordering type", "DMPlexGetOrdering", ordlist, MATORDERINGNATURAL, oname, sizeof(oname), &flg);CHKERRQ(ierr); 2807 if (flg) { 2808 DM pdm; 2809 IS perm; 2810 2811 ierr = DMPlexGetOrdering(dm, oname, NULL, &perm);CHKERRQ(ierr); 2812 ierr = DMPlexPermute(dm, perm, &pdm);CHKERRQ(ierr); 2813 ierr = ISDestroy(&perm);CHKERRQ(ierr); 2814 ierr = DMPlexReplace_Static(dm, &pdm);CHKERRQ(ierr); 2815 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2816 } 2817 /* Handle DMPlex distribution */ 2818 ierr = DMPlexDistributeGetDefault(dm, &distribute);CHKERRQ(ierr); 2819 ierr = PetscOptionsBool("-dm_distribute", "Flag to redistribute a mesh among processes", "DMCreate", distribute, &distribute, NULL);CHKERRQ(ierr); 2820 ierr = PetscOptionsBoundedInt("-dm_distribute_overlap", "The size of the overlap halo", "DMCreate", overlap, &overlap, NULL, 0);CHKERRQ(ierr); 2821 if (distribute) { 2822 DM pdm = NULL; 2823 PetscPartitioner part; 2824 2825 ierr = DMPlexGetPartitioner(dm, &part);CHKERRQ(ierr); 2826 ierr = PetscPartitionerSetFromOptions(part);CHKERRQ(ierr); 2827 ierr = DMPlexDistribute(dm, overlap, NULL, &pdm);CHKERRQ(ierr); 2828 if (pdm) { 2829 ierr = DMPlexReplace_Static(dm, &pdm);CHKERRQ(ierr); 2830 } 2831 } 2832 /* Create coordinate space */ 2833 if (created) { 2834 DM_Plex *mesh = (DM_Plex *) dm->data; 2835 PetscInt degree = 1; 2836 PetscBool periodic, flg; 2837 2838 ierr = PetscOptionsBool("-dm_coord_space", "Use an FEM space for coordinates", "", coordSpace, &coordSpace, &flg);CHKERRQ(ierr); 2839 ierr = PetscOptionsInt("-dm_coord_petscspace_degree", "FEM degree for coordinate space", "", degree, °ree, NULL);CHKERRQ(ierr); 2840 if (coordSpace) {ierr = DMPlexCreateCoordinateSpace(dm, degree, mesh->coordFunc);CHKERRQ(ierr);} 2841 if (flg && !coordSpace) { 2842 DM cdm; 2843 PetscDS cds; 2844 PetscObject obj; 2845 PetscClassId id; 2846 2847 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 2848 ierr = DMGetDS(cdm, &cds);CHKERRQ(ierr); 2849 ierr = PetscDSGetDiscretization(cds, 0, &obj);CHKERRQ(ierr); 2850 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 2851 if (id == PETSCFE_CLASSID) { 2852 PetscContainer dummy; 2853 2854 ierr = PetscContainerCreate(PETSC_COMM_SELF, &dummy);CHKERRQ(ierr); 2855 ierr = PetscObjectSetName((PetscObject) dummy, "coordinates");CHKERRQ(ierr); 2856 ierr = DMSetField(cdm, 0, NULL, (PetscObject) dummy);CHKERRQ(ierr); 2857 ierr = PetscContainerDestroy(&dummy);CHKERRQ(ierr); 2858 ierr = DMClearDS(cdm);CHKERRQ(ierr); 2859 } 2860 mesh->coordFunc = NULL; 2861 } 2862 ierr = DMLocalizeCoordinates(dm);CHKERRQ(ierr); 2863 ierr = DMGetPeriodicity(dm, &periodic, NULL, NULL, NULL);CHKERRQ(ierr); 2864 if (periodic) {ierr = DMSetPeriodicity(dm, PETSC_TRUE, NULL, NULL, NULL);CHKERRQ(ierr);} 2865 } 2866 /* Handle DMPlex refinement */ 2867 remap = PETSC_TRUE; 2868 ierr = PetscOptionsBoundedInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL,0);CHKERRQ(ierr); 2869 ierr = PetscOptionsBool("-dm_refine_remap", "Flag to control coordinate remapping", "DMCreate", remap, &remap, NULL);CHKERRQ(ierr); 2870 ierr = PetscOptionsBoundedInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy,0);CHKERRQ(ierr); 2871 if (refine) {ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE);CHKERRQ(ierr);} 2872 if (refine && isHierarchy) { 2873 DM *dms, coarseDM; 2874 2875 ierr = DMGetCoarseDM(dm, &coarseDM);CHKERRQ(ierr); 2876 ierr = PetscObjectReference((PetscObject)coarseDM);CHKERRQ(ierr); 2877 ierr = PetscMalloc1(refine,&dms);CHKERRQ(ierr); 2878 ierr = DMRefineHierarchy(dm, refine, dms);CHKERRQ(ierr); 2879 /* Total hack since we do not pass in a pointer */ 2880 ierr = DMPlexSwap_Static(dm, dms[refine-1]);CHKERRQ(ierr); 2881 if (refine == 1) { 2882 ierr = DMSetCoarseDM(dm, dms[0]);CHKERRQ(ierr); 2883 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr); 2884 } else { 2885 ierr = DMSetCoarseDM(dm, dms[refine-2]);CHKERRQ(ierr); 2886 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr); 2887 ierr = DMSetCoarseDM(dms[0], dms[refine-1]);CHKERRQ(ierr); 2888 ierr = DMPlexSetRegularRefinement(dms[0], PETSC_TRUE);CHKERRQ(ierr); 2889 } 2890 ierr = DMSetCoarseDM(dms[refine-1], coarseDM);CHKERRQ(ierr); 2891 ierr = PetscObjectDereference((PetscObject)coarseDM);CHKERRQ(ierr); 2892 /* Free DMs */ 2893 for (r = 0; r < refine; ++r) { 2894 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr); 2895 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 2896 } 2897 ierr = PetscFree(dms);CHKERRQ(ierr); 2898 } else { 2899 for (r = 0; r < refine; ++r) { 2900 DM rdm; 2901 PetscPointFunc coordFunc = ((DM_Plex*) dm->data)->coordFunc; 2902 2903 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2904 ierr = DMRefine(dm, PetscObjectComm((PetscObject) dm), &rdm);CHKERRQ(ierr); 2905 /* Total hack since we do not pass in a pointer */ 2906 ierr = DMPlexReplace_Static(dm, &rdm);CHKERRQ(ierr); 2907 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2908 if (coordFunc && remap) { 2909 ierr = DMPlexRemapGeometry(dm, 0.0, coordFunc);CHKERRQ(ierr); 2910 ((DM_Plex*) dm->data)->coordFunc = coordFunc; 2911 } 2912 } 2913 } 2914 /* Handle DMPlex coarsening */ 2915 ierr = PetscOptionsBoundedInt("-dm_coarsen", "Coarsen the mesh", "DMCreate", coarsen, &coarsen, NULL,0);CHKERRQ(ierr); 2916 ierr = PetscOptionsBoundedInt("-dm_coarsen_hierarchy", "The number of coarsenings", "DMCreate", coarsen, &coarsen, &isHierarchy,0);CHKERRQ(ierr); 2917 if (coarsen && isHierarchy) { 2918 DM *dms; 2919 2920 ierr = PetscMalloc1(coarsen, &dms);CHKERRQ(ierr); 2921 ierr = DMCoarsenHierarchy(dm, coarsen, dms);CHKERRQ(ierr); 2922 /* Free DMs */ 2923 for (r = 0; r < coarsen; ++r) { 2924 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr); 2925 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr); 2926 } 2927 ierr = PetscFree(dms);CHKERRQ(ierr); 2928 } else { 2929 for (r = 0; r < coarsen; ++r) { 2930 DM cdm; 2931 PetscPointFunc coordFunc = ((DM_Plex*) dm->data)->coordFunc; 2932 2933 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2934 ierr = DMCoarsen(dm, PetscObjectComm((PetscObject) dm), &cdm);CHKERRQ(ierr); 2935 /* Total hack since we do not pass in a pointer */ 2936 ierr = DMPlexReplace_Static(dm, &cdm);CHKERRQ(ierr); 2937 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2938 if (coordFunc) { 2939 ierr = DMPlexRemapGeometry(dm, 0.0, coordFunc);CHKERRQ(ierr); 2940 ((DM_Plex*) dm->data)->coordFunc = coordFunc; 2941 } 2942 } 2943 } 2944 /* Handle ghost cells */ 2945 ierr = PetscOptionsBool("-dm_plex_create_fv_ghost_cells", "Flag to create finite volume ghost cells on the boundary", "DMCreate", ghostCells, &ghostCells, NULL);CHKERRQ(ierr); 2946 if (ghostCells) { 2947 DM gdm; 2948 char lname[PETSC_MAX_PATH_LEN]; 2949 2950 lname[0] = '\0'; 2951 ierr = PetscOptionsString("-dm_plex_fv_ghost_cells_label", "Label name for ghost cells boundary", "DMCreate", lname, lname, sizeof(lname), &flg);CHKERRQ(ierr); 2952 ierr = DMPlexConstructGhostCells(dm, flg ? lname : NULL, NULL, &gdm);CHKERRQ(ierr); 2953 ierr = DMPlexReplace_Static(dm, &gdm);CHKERRQ(ierr); 2954 } 2955 /* Handle */ 2956 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr); 2957 ierr = PetscOptionsTail();CHKERRQ(ierr); 2958 PetscFunctionReturn(0); 2959 } 2960 2961 static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec) 2962 { 2963 PetscErrorCode ierr; 2964 2965 PetscFunctionBegin; 2966 ierr = DMCreateGlobalVector_Section_Private(dm,vec);CHKERRQ(ierr); 2967 /* ierr = VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM);CHKERRQ(ierr); */ 2968 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex);CHKERRQ(ierr); 2969 ierr = VecSetOperation(*vec, VECOP_VIEWNATIVE, (void (*)(void)) VecView_Plex_Native);CHKERRQ(ierr); 2970 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex);CHKERRQ(ierr); 2971 ierr = VecSetOperation(*vec, VECOP_LOADNATIVE, (void (*)(void)) VecLoad_Plex_Native);CHKERRQ(ierr); 2972 PetscFunctionReturn(0); 2973 } 2974 2975 static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec) 2976 { 2977 PetscErrorCode ierr; 2978 2979 PetscFunctionBegin; 2980 ierr = DMCreateLocalVector_Section_Private(dm,vec);CHKERRQ(ierr); 2981 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex_Local);CHKERRQ(ierr); 2982 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex_Local);CHKERRQ(ierr); 2983 PetscFunctionReturn(0); 2984 } 2985 2986 static PetscErrorCode DMGetDimPoints_Plex(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd) 2987 { 2988 PetscInt depth, d; 2989 PetscErrorCode ierr; 2990 2991 PetscFunctionBegin; 2992 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 2993 if (depth == 1) { 2994 ierr = DMGetDimension(dm, &d);CHKERRQ(ierr); 2995 if (dim == 0) {ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr);} 2996 else if (dim == d) {ierr = DMPlexGetDepthStratum(dm, 1, pStart, pEnd);CHKERRQ(ierr);} 2997 else {*pStart = 0; *pEnd = 0;} 2998 } else { 2999 ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr); 3000 } 3001 PetscFunctionReturn(0); 3002 } 3003 3004 static PetscErrorCode DMGetNeighbors_Plex(DM dm, PetscInt *nranks, const PetscMPIInt *ranks[]) 3005 { 3006 PetscSF sf; 3007 PetscInt niranks, njranks, n; 3008 const PetscMPIInt *iranks, *jranks; 3009 DM_Plex *data = (DM_Plex*) dm->data; 3010 PetscErrorCode ierr; 3011 3012 PetscFunctionBegin; 3013 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 3014 if (!data->neighbors) { 3015 ierr = PetscSFSetUp(sf);CHKERRQ(ierr); 3016 ierr = PetscSFGetRootRanks(sf, &njranks, &jranks, NULL, NULL, NULL);CHKERRQ(ierr); 3017 ierr = PetscSFGetLeafRanks(sf, &niranks, &iranks, NULL, NULL);CHKERRQ(ierr); 3018 ierr = PetscMalloc1(njranks + niranks + 1, &data->neighbors);CHKERRQ(ierr); 3019 ierr = PetscArraycpy(data->neighbors + 1, jranks, njranks);CHKERRQ(ierr); 3020 ierr = PetscArraycpy(data->neighbors + njranks + 1, iranks, niranks);CHKERRQ(ierr); 3021 n = njranks + niranks; 3022 ierr = PetscSortRemoveDupsMPIInt(&n, data->neighbors + 1);CHKERRQ(ierr); 3023 /* The following cast should never fail: can't have more neighbors than PETSC_MPI_INT_MAX */ 3024 ierr = PetscMPIIntCast(n, data->neighbors);CHKERRQ(ierr); 3025 } 3026 if (nranks) *nranks = data->neighbors[0]; 3027 if (ranks) { 3028 if (data->neighbors[0]) *ranks = data->neighbors + 1; 3029 else *ranks = NULL; 3030 } 3031 PetscFunctionReturn(0); 3032 } 3033 3034 PETSC_INTERN PetscErrorCode DMInterpolateSolution_Plex(DM, DM, Mat, Vec, Vec); 3035 3036 static PetscErrorCode DMInitialize_Plex(DM dm) 3037 { 3038 PetscErrorCode ierr; 3039 3040 PetscFunctionBegin; 3041 dm->ops->view = DMView_Plex; 3042 dm->ops->load = DMLoad_Plex; 3043 dm->ops->setfromoptions = DMSetFromOptions_Plex; 3044 dm->ops->clone = DMClone_Plex; 3045 dm->ops->setup = DMSetUp_Plex; 3046 dm->ops->createlocalsection = DMCreateLocalSection_Plex; 3047 dm->ops->createdefaultconstraints = DMCreateDefaultConstraints_Plex; 3048 dm->ops->createglobalvector = DMCreateGlobalVector_Plex; 3049 dm->ops->createlocalvector = DMCreateLocalVector_Plex; 3050 dm->ops->getlocaltoglobalmapping = NULL; 3051 dm->ops->createfieldis = NULL; 3052 dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex; 3053 dm->ops->createcoordinatefield = DMCreateCoordinateField_Plex; 3054 dm->ops->getcoloring = NULL; 3055 dm->ops->creatematrix = DMCreateMatrix_Plex; 3056 dm->ops->createinterpolation = DMCreateInterpolation_Plex; 3057 dm->ops->createmassmatrix = DMCreateMassMatrix_Plex; 3058 dm->ops->createmassmatrixlumped = DMCreateMassMatrixLumped_Plex; 3059 dm->ops->createinjection = DMCreateInjection_Plex; 3060 dm->ops->refine = DMRefine_Plex; 3061 dm->ops->coarsen = DMCoarsen_Plex; 3062 dm->ops->refinehierarchy = DMRefineHierarchy_Plex; 3063 dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Plex; 3064 dm->ops->extrude = DMExtrude_Plex; 3065 dm->ops->globaltolocalbegin = NULL; 3066 dm->ops->globaltolocalend = NULL; 3067 dm->ops->localtoglobalbegin = NULL; 3068 dm->ops->localtoglobalend = NULL; 3069 dm->ops->destroy = DMDestroy_Plex; 3070 dm->ops->createsubdm = DMCreateSubDM_Plex; 3071 dm->ops->createsuperdm = DMCreateSuperDM_Plex; 3072 dm->ops->getdimpoints = DMGetDimPoints_Plex; 3073 dm->ops->locatepoints = DMLocatePoints_Plex; 3074 dm->ops->projectfunctionlocal = DMProjectFunctionLocal_Plex; 3075 dm->ops->projectfunctionlabellocal = DMProjectFunctionLabelLocal_Plex; 3076 dm->ops->projectfieldlocal = DMProjectFieldLocal_Plex; 3077 dm->ops->projectfieldlabellocal = DMProjectFieldLabelLocal_Plex; 3078 dm->ops->projectbdfieldlabellocal = DMProjectBdFieldLabelLocal_Plex; 3079 dm->ops->computel2diff = DMComputeL2Diff_Plex; 3080 dm->ops->computel2gradientdiff = DMComputeL2GradientDiff_Plex; 3081 dm->ops->computel2fielddiff = DMComputeL2FieldDiff_Plex; 3082 dm->ops->getneighbors = DMGetNeighbors_Plex; 3083 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMPlexInsertBoundaryValues_C",DMPlexInsertBoundaryValues_Plex);CHKERRQ(ierr); 3084 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMPlexInsertTimeDerviativeBoundaryValues_C",DMPlexInsertTimeDerivativeBoundaryValues_Plex);CHKERRQ(ierr); 3085 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMSetUpGLVisViewer_C",DMSetUpGLVisViewer_Plex);CHKERRQ(ierr); 3086 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMCreateNeumannOverlap_C",DMCreateNeumannOverlap_Plex);CHKERRQ(ierr); 3087 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMPlexGetOverlap_C",DMPlexGetOverlap_Plex);CHKERRQ(ierr); 3088 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMPlexDistributeGetDefault_C",DMPlexDistributeGetDefault_Plex);CHKERRQ(ierr); 3089 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMPlexDistributeSetDefault_C",DMPlexDistributeSetDefault_Plex);CHKERRQ(ierr); 3090 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMInterpolateSolution_C",DMInterpolateSolution_Plex);CHKERRQ(ierr); 3091 PetscFunctionReturn(0); 3092 } 3093 3094 PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm) 3095 { 3096 DM_Plex *mesh = (DM_Plex *) dm->data; 3097 PetscErrorCode ierr; 3098 3099 PetscFunctionBegin; 3100 mesh->refct++; 3101 (*newdm)->data = mesh; 3102 ierr = PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);CHKERRQ(ierr); 3103 ierr = DMInitialize_Plex(*newdm);CHKERRQ(ierr); 3104 PetscFunctionReturn(0); 3105 } 3106 3107 /*MC 3108 DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram. 3109 In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is 3110 specified by a PetscSection object. Ownership in the global representation is determined by 3111 ownership of the underlying DMPlex points. This is specified by another PetscSection object. 3112 3113 Options Database Keys: 3114 + -dm_refine_pre - Refine mesh before distribution 3115 + -dm_refine_uniform_pre - Choose uniform or generator-based refinement 3116 + -dm_refine_volume_limit_pre - Cell volume limit after pre-refinement using generator 3117 . -dm_distribute - Distribute mesh across processes 3118 . -dm_distribute_overlap - Number of cells to overlap for distribution 3119 . -dm_refine - Refine mesh after distribution 3120 . -dm_plex_hash_location - Use grid hashing for point location 3121 . -dm_plex_hash_box_faces <n,m,p> - The number of divisions in each direction of the grid hash 3122 . -dm_plex_partition_balance - Attempt to evenly divide points on partition boundary between processes 3123 . -dm_plex_remesh_bd - Allow changes to the boundary on remeshing 3124 . -dm_plex_max_projection_height - Maxmimum mesh point height used to project locally 3125 . -dm_plex_regular_refinement - Use special nested projection algorithm for regular refinement 3126 . -dm_plex_check_all - Perform all shecks below 3127 . -dm_plex_check_symmetry - Check that the adjacency information in the mesh is symmetric 3128 . -dm_plex_check_skeleton <celltype> - Check that each cell has the correct number of vertices 3129 . -dm_plex_check_faces <celltype> - Check that the faces of each cell give a vertex order this is consistent with what we expect from the cell type 3130 . -dm_plex_check_geometry - Check that cells have positive volume 3131 . -dm_view :mesh.tex:ascii_latex - View the mesh in LaTeX/TikZ 3132 . -dm_plex_view_scale <num> - Scale the TikZ 3133 - -dm_plex_print_fem <num> - View FEM assembly information, such as element vectors and matrices 3134 3135 Level: intermediate 3136 3137 .seealso: DMType, DMPlexCreate(), DMCreate(), DMSetType() 3138 M*/ 3139 3140 PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm) 3141 { 3142 DM_Plex *mesh; 3143 PetscInt unit; 3144 PetscErrorCode ierr; 3145 3146 PetscFunctionBegin; 3147 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3148 ierr = PetscNewLog(dm,&mesh);CHKERRQ(ierr); 3149 dm->data = mesh; 3150 3151 mesh->refct = 1; 3152 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);CHKERRQ(ierr); 3153 mesh->maxConeSize = 0; 3154 mesh->cones = NULL; 3155 mesh->coneOrientations = NULL; 3156 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);CHKERRQ(ierr); 3157 mesh->maxSupportSize = 0; 3158 mesh->supports = NULL; 3159 mesh->refinementUniform = PETSC_TRUE; 3160 mesh->refinementLimit = -1.0; 3161 mesh->distDefault = PETSC_TRUE; 3162 mesh->interpolated = DMPLEX_INTERPOLATED_INVALID; 3163 mesh->interpolatedCollective = DMPLEX_INTERPOLATED_INVALID; 3164 3165 mesh->facesTmp = NULL; 3166 3167 mesh->tetgenOpts = NULL; 3168 mesh->triangleOpts = NULL; 3169 ierr = PetscPartitionerCreate(PetscObjectComm((PetscObject)dm), &mesh->partitioner);CHKERRQ(ierr); 3170 mesh->remeshBd = PETSC_FALSE; 3171 3172 mesh->subpointMap = NULL; 3173 3174 for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0; 3175 3176 mesh->regularRefinement = PETSC_FALSE; 3177 mesh->depthState = -1; 3178 mesh->celltypeState = -1; 3179 mesh->globalVertexNumbers = NULL; 3180 mesh->globalCellNumbers = NULL; 3181 mesh->anchorSection = NULL; 3182 mesh->anchorIS = NULL; 3183 mesh->createanchors = NULL; 3184 mesh->computeanchormatrix = NULL; 3185 mesh->parentSection = NULL; 3186 mesh->parents = NULL; 3187 mesh->childIDs = NULL; 3188 mesh->childSection = NULL; 3189 mesh->children = NULL; 3190 mesh->referenceTree = NULL; 3191 mesh->getchildsymmetry = NULL; 3192 mesh->vtkCellHeight = 0; 3193 mesh->useAnchors = PETSC_FALSE; 3194 3195 mesh->maxProjectionHeight = 0; 3196 3197 mesh->neighbors = NULL; 3198 3199 mesh->printSetValues = PETSC_FALSE; 3200 mesh->printFEM = 0; 3201 mesh->printTol = 1.0e-10; 3202 3203 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr); 3204 PetscFunctionReturn(0); 3205 } 3206 3207 /*@ 3208 DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram. 3209 3210 Collective 3211 3212 Input Parameter: 3213 . comm - The communicator for the DMPlex object 3214 3215 Output Parameter: 3216 . mesh - The DMPlex object 3217 3218 Level: beginner 3219 3220 @*/ 3221 PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh) 3222 { 3223 PetscErrorCode ierr; 3224 3225 PetscFunctionBegin; 3226 PetscValidPointer(mesh,2); 3227 ierr = DMCreate(comm, mesh);CHKERRQ(ierr); 3228 ierr = DMSetType(*mesh, DMPLEX);CHKERRQ(ierr); 3229 PetscFunctionReturn(0); 3230 } 3231 3232 /*@C 3233 DMPlexBuildFromCellListParallel - Build distributed DMPLEX topology from a list of vertices for each cell (common mesh generator output) 3234 3235 Input Parameters: 3236 + dm - The DM 3237 . numCells - The number of cells owned by this process 3238 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE 3239 . NVertices - The global number of vertices, or PETSC_DECIDE 3240 . numCorners - The number of vertices for each cell 3241 - cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 3242 3243 Output Parameters: 3244 + vertexSF - (Optional) SF describing complete vertex ownership 3245 - verticesAdjSaved - (Optional) vertex adjacency array 3246 3247 Notes: 3248 Two triangles sharing a face 3249 $ 3250 $ 2 3251 $ / | \ 3252 $ / | \ 3253 $ / | \ 3254 $ 0 0 | 1 3 3255 $ \ | / 3256 $ \ | / 3257 $ \ | / 3258 $ 1 3259 would have input 3260 $ numCells = 2, numVertices = 4 3261 $ cells = [0 1 2 1 3 2] 3262 $ 3263 which would result in the DMPlex 3264 $ 3265 $ 4 3266 $ / | \ 3267 $ / | \ 3268 $ / | \ 3269 $ 2 0 | 1 5 3270 $ \ | / 3271 $ \ | / 3272 $ \ | / 3273 $ 3 3274 3275 Vertices are implicitly numbered consecutively 0,...,NVertices. 3276 Each rank owns a chunk of numVertices consecutive vertices. 3277 If numVertices is PETSC_DECIDE, PETSc will distribute them as evenly as possible using PetscLayout. 3278 If both NVertices and numVertices are PETSC_DECIDE, NVertices is computed by PETSc as the maximum vertex index in cells + 1. 3279 If only NVertices is PETSC_DECIDE, it is computed as the sum of numVertices over all ranks. 3280 3281 The cell distribution is arbitrary non-overlapping, independent of the vertex distribution. 3282 3283 Not currently supported in Fortran. 3284 3285 Level: advanced 3286 3287 .seealso: DMPlexBuildFromCellList(), DMPlexCreateFromCellListParallelPetsc(), DMPlexBuildCoordinatesFromCellListParallel() 3288 @*/ 3289 PetscErrorCode DMPlexBuildFromCellListParallel(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt NVertices, PetscInt numCorners, const PetscInt cells[], PetscSF *vertexSF, PetscInt **verticesAdjSaved) 3290 { 3291 PetscSF sfPoint; 3292 PetscLayout layout; 3293 PetscInt numVerticesAdj, *verticesAdj, *cones, c, p, dim; 3294 PetscMPIInt rank, size; 3295 PetscErrorCode ierr; 3296 3297 PetscFunctionBegin; 3298 PetscValidLogicalCollectiveInt(dm,NVertices,4); 3299 ierr = PetscLogEventBegin(DMPLEX_BuildFromCellList,dm,0,0,0);CHKERRQ(ierr); 3300 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRMPI(ierr); 3301 ierr = MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);CHKERRMPI(ierr); 3302 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3303 /* Get/check global number of vertices */ 3304 { 3305 PetscInt NVerticesInCells, i; 3306 const PetscInt len = numCells * numCorners; 3307 3308 /* NVerticesInCells = max(cells) + 1 */ 3309 NVerticesInCells = PETSC_MIN_INT; 3310 for (i=0; i<len; i++) if (cells[i] > NVerticesInCells) NVerticesInCells = cells[i]; 3311 ++NVerticesInCells; 3312 ierr = MPI_Allreduce(MPI_IN_PLACE, &NVerticesInCells, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject) dm));CHKERRMPI(ierr); 3313 3314 if (numVertices == PETSC_DECIDE && NVertices == PETSC_DECIDE) NVertices = NVerticesInCells; 3315 else PetscCheckFalse(NVertices != PETSC_DECIDE && NVertices < NVerticesInCells,PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Specified global number of vertices %D must be greater than or equal to the number of vertices in cells %D",NVertices,NVerticesInCells); 3316 } 3317 /* Count locally unique vertices */ 3318 { 3319 PetscHSetI vhash; 3320 PetscInt off = 0; 3321 3322 ierr = PetscHSetICreate(&vhash);CHKERRQ(ierr); 3323 for (c = 0; c < numCells; ++c) { 3324 for (p = 0; p < numCorners; ++p) { 3325 ierr = PetscHSetIAdd(vhash, cells[c*numCorners+p]);CHKERRQ(ierr); 3326 } 3327 } 3328 ierr = PetscHSetIGetSize(vhash, &numVerticesAdj);CHKERRQ(ierr); 3329 if (!verticesAdjSaved) { ierr = PetscMalloc1(numVerticesAdj, &verticesAdj);CHKERRQ(ierr); } 3330 else { verticesAdj = *verticesAdjSaved; } 3331 ierr = PetscHSetIGetElems(vhash, &off, verticesAdj);CHKERRQ(ierr); 3332 ierr = PetscHSetIDestroy(&vhash);CHKERRQ(ierr); 3333 PetscCheckFalse(off != numVerticesAdj,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid number of local vertices %D should be %D", off, numVerticesAdj); 3334 } 3335 ierr = PetscSortInt(numVerticesAdj, verticesAdj);CHKERRQ(ierr); 3336 /* Create cones */ 3337 ierr = DMPlexSetChart(dm, 0, numCells+numVerticesAdj);CHKERRQ(ierr); 3338 for (c = 0; c < numCells; ++c) {ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr);} 3339 ierr = DMSetUp(dm);CHKERRQ(ierr); 3340 ierr = DMPlexGetCones(dm,&cones);CHKERRQ(ierr); 3341 for (c = 0; c < numCells; ++c) { 3342 for (p = 0; p < numCorners; ++p) { 3343 const PetscInt gv = cells[c*numCorners+p]; 3344 PetscInt lv; 3345 3346 /* Positions within verticesAdj form 0-based local vertex numbering; 3347 we need to shift it by numCells to get correct DAG points (cells go first) */ 3348 ierr = PetscFindInt(gv, numVerticesAdj, verticesAdj, &lv);CHKERRQ(ierr); 3349 PetscCheckFalse(lv < 0,PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not find global vertex %D in local connectivity", gv); 3350 cones[c*numCorners+p] = lv+numCells; 3351 } 3352 } 3353 /* Build point sf */ 3354 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject)dm), &layout);CHKERRQ(ierr); 3355 ierr = PetscLayoutSetSize(layout, NVertices);CHKERRQ(ierr); 3356 ierr = PetscLayoutSetLocalSize(layout, numVertices);CHKERRQ(ierr); 3357 ierr = PetscLayoutSetBlockSize(layout, 1);CHKERRQ(ierr); 3358 ierr = PetscSFCreateByMatchingIndices(layout, numVerticesAdj, verticesAdj, NULL, numCells, numVerticesAdj, verticesAdj, NULL, numCells, vertexSF, &sfPoint);CHKERRQ(ierr); 3359 ierr = PetscLayoutDestroy(&layout);CHKERRQ(ierr); 3360 if (!verticesAdjSaved) { ierr = PetscFree(verticesAdj);CHKERRQ(ierr); } 3361 ierr = PetscObjectSetName((PetscObject) sfPoint, "point SF");CHKERRQ(ierr); 3362 if (dm->sf) { 3363 const char *prefix; 3364 3365 ierr = PetscObjectGetOptionsPrefix((PetscObject)dm->sf, &prefix);CHKERRQ(ierr); 3366 ierr = PetscObjectSetOptionsPrefix((PetscObject)sfPoint, prefix);CHKERRQ(ierr); 3367 } 3368 ierr = DMSetPointSF(dm, sfPoint);CHKERRQ(ierr); 3369 ierr = PetscSFDestroy(&sfPoint);CHKERRQ(ierr); 3370 if (vertexSF) {ierr = PetscObjectSetName((PetscObject)(*vertexSF), "Vertex Ownership SF");CHKERRQ(ierr);} 3371 /* Fill in the rest of the topology structure */ 3372 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 3373 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 3374 ierr = PetscLogEventEnd(DMPLEX_BuildFromCellList,dm,0,0,0);CHKERRQ(ierr); 3375 PetscFunctionReturn(0); 3376 } 3377 3378 /*@C 3379 DMPlexBuildCoordinatesFromCellListParallel - Build DM coordinates from a list of coordinates for each owned vertex (common mesh generator output) 3380 3381 Input Parameters: 3382 + dm - The DM 3383 . spaceDim - The spatial dimension used for coordinates 3384 . sfVert - SF describing complete vertex ownership 3385 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 3386 3387 Level: advanced 3388 3389 Notes: 3390 Not currently supported in Fortran. 3391 3392 .seealso: DMPlexBuildCoordinatesFromCellList(), DMPlexCreateFromCellListParallelPetsc(), DMPlexBuildFromCellListParallel() 3393 @*/ 3394 PetscErrorCode DMPlexBuildCoordinatesFromCellListParallel(DM dm, PetscInt spaceDim, PetscSF sfVert, const PetscReal vertexCoords[]) 3395 { 3396 PetscSection coordSection; 3397 Vec coordinates; 3398 PetscScalar *coords; 3399 PetscInt numVertices, numVerticesAdj, coordSize, v, vStart, vEnd; 3400 PetscErrorCode ierr; 3401 3402 PetscFunctionBegin; 3403 ierr = PetscLogEventBegin(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0);CHKERRQ(ierr); 3404 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 3405 PetscCheckFalse(vStart < 0 || vEnd < 0,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM is not set up properly. DMPlexBuildFromCellList() should be called first."); 3406 ierr = DMSetCoordinateDim(dm, spaceDim);CHKERRQ(ierr); 3407 ierr = PetscSFGetGraph(sfVert, &numVertices, &numVerticesAdj, NULL, NULL);CHKERRQ(ierr); 3408 PetscCheckFalse(vEnd - vStart != numVerticesAdj,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Supplied sfVert has wrong number of leaves = %D != %D = vEnd - vStart",numVerticesAdj,vEnd - vStart); 3409 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 3410 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 3411 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 3412 ierr = PetscSectionSetChart(coordSection, vStart, vEnd);CHKERRQ(ierr); 3413 for (v = vStart; v < vEnd; ++v) { 3414 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 3415 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 3416 } 3417 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 3418 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 3419 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr); 3420 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 3421 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 3422 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 3423 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 3424 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 3425 { 3426 MPI_Datatype coordtype; 3427 3428 /* Need a temp buffer for coords if we have complex/single */ 3429 ierr = MPI_Type_contiguous(spaceDim, MPIU_SCALAR, &coordtype);CHKERRMPI(ierr); 3430 ierr = MPI_Type_commit(&coordtype);CHKERRMPI(ierr); 3431 #if defined(PETSC_USE_COMPLEX) 3432 { 3433 PetscScalar *svertexCoords; 3434 PetscInt i; 3435 ierr = PetscMalloc1(numVertices*spaceDim,&svertexCoords);CHKERRQ(ierr); 3436 for (i=0; i<numVertices*spaceDim; i++) svertexCoords[i] = vertexCoords[i]; 3437 ierr = PetscSFBcastBegin(sfVert, coordtype, svertexCoords, coords,MPI_REPLACE);CHKERRQ(ierr); 3438 ierr = PetscSFBcastEnd(sfVert, coordtype, svertexCoords, coords,MPI_REPLACE);CHKERRQ(ierr); 3439 ierr = PetscFree(svertexCoords);CHKERRQ(ierr); 3440 } 3441 #else 3442 ierr = PetscSFBcastBegin(sfVert, coordtype, vertexCoords, coords,MPI_REPLACE);CHKERRQ(ierr); 3443 ierr = PetscSFBcastEnd(sfVert, coordtype, vertexCoords, coords,MPI_REPLACE);CHKERRQ(ierr); 3444 #endif 3445 ierr = MPI_Type_free(&coordtype);CHKERRMPI(ierr); 3446 } 3447 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 3448 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 3449 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 3450 ierr = PetscLogEventEnd(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0);CHKERRQ(ierr); 3451 PetscFunctionReturn(0); 3452 } 3453 3454 /*@ 3455 DMPlexCreateFromCellListParallelPetsc - Create distributed DMPLEX from a list of vertices for each cell (common mesh generator output) 3456 3457 Input Parameters: 3458 + comm - The communicator 3459 . dim - The topological dimension of the mesh 3460 . numCells - The number of cells owned by this process 3461 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE 3462 . NVertices - The global number of vertices, or PETSC_DECIDE 3463 . numCorners - The number of vertices for each cell 3464 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 3465 . cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 3466 . spaceDim - The spatial dimension used for coordinates 3467 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 3468 3469 Output Parameters: 3470 + dm - The DM 3471 . vertexSF - (Optional) SF describing complete vertex ownership 3472 - verticesAdjSaved - (Optional) vertex adjacency array 3473 3474 Notes: 3475 This function is just a convenient sequence of DMCreate(), DMSetType(), DMSetDimension(), 3476 DMPlexBuildFromCellListParallel(), DMPlexInterpolate(), DMPlexBuildCoordinatesFromCellListParallel() 3477 3478 See DMPlexBuildFromCellListParallel() for an example and details about the topology-related parameters. 3479 See DMPlexBuildCoordinatesFromCellListParallel() for details about the geometry-related parameters. 3480 3481 Level: intermediate 3482 3483 .seealso: DMPlexCreateFromCellListPetsc(), DMPlexBuildFromCellListParallel(), DMPlexBuildCoordinatesFromCellListParallel(), DMPlexCreateFromDAG(), DMPlexCreate() 3484 @*/ 3485 PetscErrorCode DMPlexCreateFromCellListParallelPetsc(MPI_Comm comm, PetscInt dim, PetscInt numCells, PetscInt numVertices, PetscInt NVertices, PetscInt numCorners, PetscBool interpolate, const PetscInt cells[], PetscInt spaceDim, const PetscReal vertexCoords[], PetscSF *vertexSF, PetscInt **verticesAdj, DM *dm) 3486 { 3487 PetscSF sfVert; 3488 PetscErrorCode ierr; 3489 3490 PetscFunctionBegin; 3491 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 3492 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 3493 PetscValidLogicalCollectiveInt(*dm, dim, 2); 3494 PetscValidLogicalCollectiveInt(*dm, spaceDim, 9); 3495 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 3496 ierr = DMPlexBuildFromCellListParallel(*dm, numCells, numVertices, NVertices, numCorners, cells, &sfVert, verticesAdj);CHKERRQ(ierr); 3497 if (interpolate) { 3498 DM idm; 3499 3500 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 3501 ierr = DMDestroy(dm);CHKERRQ(ierr); 3502 *dm = idm; 3503 } 3504 ierr = DMPlexBuildCoordinatesFromCellListParallel(*dm, spaceDim, sfVert, vertexCoords);CHKERRQ(ierr); 3505 if (vertexSF) *vertexSF = sfVert; 3506 else {ierr = PetscSFDestroy(&sfVert);CHKERRQ(ierr);} 3507 PetscFunctionReturn(0); 3508 } 3509 3510 /*@ 3511 DMPlexCreateFromCellListParallel - Deprecated, use DMPlexCreateFromCellListParallelPetsc() 3512 3513 Level: deprecated 3514 3515 .seealso: DMPlexCreateFromCellListParallelPetsc() 3516 @*/ 3517 PetscErrorCode DMPlexCreateFromCellListParallel(MPI_Comm comm, PetscInt dim, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, PetscBool interpolate, const int cells[], PetscInt spaceDim, const PetscReal vertexCoords[], PetscSF *vertexSF, DM *dm) 3518 { 3519 PetscErrorCode ierr; 3520 PetscInt i; 3521 PetscInt *pintCells; 3522 3523 PetscFunctionBegin; 3524 PetscCheckFalse(sizeof(int) > sizeof(PetscInt),comm, PETSC_ERR_ARG_SIZ, "Size of int %zd greater than size of PetscInt %zd. Reconfigure PETSc --with-64-bit-indices=1", sizeof(int), sizeof(PetscInt)); 3525 if (sizeof(int) == sizeof(PetscInt)) { 3526 pintCells = (PetscInt *) cells; 3527 } else { 3528 ierr = PetscMalloc1(numCells*numCorners, &pintCells);CHKERRQ(ierr); 3529 for (i = 0; i < numCells*numCorners; i++) { 3530 pintCells[i] = (PetscInt) cells[i]; 3531 } 3532 } 3533 ierr = DMPlexCreateFromCellListParallelPetsc(comm, dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, pintCells, spaceDim, vertexCoords, vertexSF, NULL, dm);CHKERRQ(ierr); 3534 if (sizeof(int) != sizeof(PetscInt)) { 3535 ierr = PetscFree(pintCells);CHKERRQ(ierr); 3536 } 3537 PetscFunctionReturn(0); 3538 } 3539 3540 /*@C 3541 DMPlexBuildFromCellList - Build DMPLEX topology from a list of vertices for each cell (common mesh generator output) 3542 3543 Input Parameters: 3544 + dm - The DM 3545 . numCells - The number of cells owned by this process 3546 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE 3547 . numCorners - The number of vertices for each cell 3548 - cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell 3549 3550 Level: advanced 3551 3552 Notes: 3553 Two triangles sharing a face 3554 $ 3555 $ 2 3556 $ / | \ 3557 $ / | \ 3558 $ / | \ 3559 $ 0 0 | 1 3 3560 $ \ | / 3561 $ \ | / 3562 $ \ | / 3563 $ 1 3564 would have input 3565 $ numCells = 2, numVertices = 4 3566 $ cells = [0 1 2 1 3 2] 3567 $ 3568 which would result in the DMPlex 3569 $ 3570 $ 4 3571 $ / | \ 3572 $ / | \ 3573 $ / | \ 3574 $ 2 0 | 1 5 3575 $ \ | / 3576 $ \ | / 3577 $ \ | / 3578 $ 3 3579 3580 If numVertices is PETSC_DECIDE, it is computed by PETSc as the maximum vertex index in cells + 1. 3581 3582 Not currently supported in Fortran. 3583 3584 .seealso: DMPlexBuildFromCellListParallel(), DMPlexBuildCoordinatesFromCellList(), DMPlexCreateFromCellListPetsc() 3585 @*/ 3586 PetscErrorCode DMPlexBuildFromCellList(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const PetscInt cells[]) 3587 { 3588 PetscInt *cones, c, p, dim; 3589 PetscErrorCode ierr; 3590 3591 PetscFunctionBegin; 3592 ierr = PetscLogEventBegin(DMPLEX_BuildFromCellList,dm,0,0,0);CHKERRQ(ierr); 3593 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3594 /* Get/check global number of vertices */ 3595 { 3596 PetscInt NVerticesInCells, i; 3597 const PetscInt len = numCells * numCorners; 3598 3599 /* NVerticesInCells = max(cells) + 1 */ 3600 NVerticesInCells = PETSC_MIN_INT; 3601 for (i=0; i<len; i++) if (cells[i] > NVerticesInCells) NVerticesInCells = cells[i]; 3602 ++NVerticesInCells; 3603 3604 if (numVertices == PETSC_DECIDE) numVertices = NVerticesInCells; 3605 else PetscCheckFalse(numVertices < NVerticesInCells,PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Specified number of vertices %D must be greater than or equal to the number of vertices in cells %D",numVertices,NVerticesInCells); 3606 } 3607 ierr = DMPlexSetChart(dm, 0, numCells+numVertices);CHKERRQ(ierr); 3608 for (c = 0; c < numCells; ++c) { 3609 ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr); 3610 } 3611 ierr = DMSetUp(dm);CHKERRQ(ierr); 3612 ierr = DMPlexGetCones(dm,&cones);CHKERRQ(ierr); 3613 for (c = 0; c < numCells; ++c) { 3614 for (p = 0; p < numCorners; ++p) { 3615 cones[c*numCorners+p] = cells[c*numCorners+p]+numCells; 3616 } 3617 } 3618 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 3619 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 3620 ierr = PetscLogEventEnd(DMPLEX_BuildFromCellList,dm,0,0,0);CHKERRQ(ierr); 3621 PetscFunctionReturn(0); 3622 } 3623 3624 /*@C 3625 DMPlexBuildCoordinatesFromCellList - Build DM coordinates from a list of coordinates for each owned vertex (common mesh generator output) 3626 3627 Input Parameters: 3628 + dm - The DM 3629 . spaceDim - The spatial dimension used for coordinates 3630 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex 3631 3632 Level: advanced 3633 3634 Notes: 3635 Not currently supported in Fortran. 3636 3637 .seealso: DMPlexBuildCoordinatesFromCellListParallel(), DMPlexCreateFromCellListPetsc(), DMPlexBuildFromCellList() 3638 @*/ 3639 PetscErrorCode DMPlexBuildCoordinatesFromCellList(DM dm, PetscInt spaceDim, const PetscReal vertexCoords[]) 3640 { 3641 PetscSection coordSection; 3642 Vec coordinates; 3643 DM cdm; 3644 PetscScalar *coords; 3645 PetscInt v, vStart, vEnd, d; 3646 PetscErrorCode ierr; 3647 3648 PetscFunctionBegin; 3649 ierr = PetscLogEventBegin(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0);CHKERRQ(ierr); 3650 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 3651 PetscCheckFalse(vStart < 0 || vEnd < 0,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM is not set up properly. DMPlexBuildFromCellList() should be called first."); 3652 ierr = DMSetCoordinateDim(dm, spaceDim);CHKERRQ(ierr); 3653 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 3654 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 3655 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr); 3656 ierr = PetscSectionSetChart(coordSection, vStart, vEnd);CHKERRQ(ierr); 3657 for (v = vStart; v < vEnd; ++v) { 3658 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr); 3659 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr); 3660 } 3661 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 3662 3663 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 3664 ierr = DMCreateLocalVector(cdm, &coordinates);CHKERRQ(ierr); 3665 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr); 3666 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 3667 ierr = VecGetArrayWrite(coordinates, &coords);CHKERRQ(ierr); 3668 for (v = 0; v < vEnd-vStart; ++v) { 3669 for (d = 0; d < spaceDim; ++d) { 3670 coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d]; 3671 } 3672 } 3673 ierr = VecRestoreArrayWrite(coordinates, &coords);CHKERRQ(ierr); 3674 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 3675 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 3676 ierr = PetscLogEventEnd(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0);CHKERRQ(ierr); 3677 PetscFunctionReturn(0); 3678 } 3679 3680 /*@ 3681 DMPlexCreateFromCellListPetsc - Create DMPLEX from a list of vertices for each cell (common mesh generator output), but only process 0 takes in the input 3682 3683 Collective on comm 3684 3685 Input Parameters: 3686 + comm - The communicator 3687 . dim - The topological dimension of the mesh 3688 . numCells - The number of cells, only on process 0 3689 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE, only on process 0 3690 . numCorners - The number of vertices for each cell, only on process 0 3691 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically 3692 . cells - An array of numCells*numCorners numbers, the vertices for each cell, only on process 0 3693 . spaceDim - The spatial dimension used for coordinates 3694 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex, only on process 0 3695 3696 Output Parameter: 3697 . dm - The DM, which only has points on process 0 3698 3699 Notes: 3700 This function is just a convenient sequence of DMCreate(), DMSetType(), DMSetDimension(), DMPlexBuildFromCellList(), 3701 DMPlexInterpolate(), DMPlexBuildCoordinatesFromCellList() 3702 3703 See DMPlexBuildFromCellList() for an example and details about the topology-related parameters. 3704 See DMPlexBuildCoordinatesFromCellList() for details about the geometry-related parameters. 3705 See DMPlexCreateFromCellListParallelPetsc() for parallel input 3706 3707 Level: intermediate 3708 3709 .seealso: DMPlexCreateFromCellListParallelPetsc(), DMPlexBuildFromCellList(), DMPlexBuildCoordinatesFromCellList(), DMPlexCreateFromDAG(), DMPlexCreate() 3710 @*/ 3711 PetscErrorCode DMPlexCreateFromCellListPetsc(MPI_Comm comm, PetscInt dim, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, PetscBool interpolate, const PetscInt cells[], PetscInt spaceDim, const PetscReal vertexCoords[], DM *dm) 3712 { 3713 PetscMPIInt rank; 3714 PetscErrorCode ierr; 3715 3716 PetscFunctionBegin; 3717 PetscCheckFalse(!dim,comm, PETSC_ERR_ARG_OUTOFRANGE, "This is not appropriate for 0-dimensional meshes. Consider either creating the DM using DMPlexCreateFromDAG(), by hand, or using DMSwarm."); 3718 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 3719 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 3720 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 3721 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 3722 if (!rank) {ierr = DMPlexBuildFromCellList(*dm, numCells, numVertices, numCorners, cells);CHKERRQ(ierr);} 3723 else {ierr = DMPlexBuildFromCellList(*dm, 0, 0, 0, NULL);CHKERRQ(ierr);} 3724 if (interpolate) { 3725 DM idm; 3726 3727 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 3728 ierr = DMDestroy(dm);CHKERRQ(ierr); 3729 *dm = idm; 3730 } 3731 if (!rank) {ierr = DMPlexBuildCoordinatesFromCellList(*dm, spaceDim, vertexCoords);CHKERRQ(ierr);} 3732 else {ierr = DMPlexBuildCoordinatesFromCellList(*dm, spaceDim, NULL);CHKERRQ(ierr);} 3733 PetscFunctionReturn(0); 3734 } 3735 3736 /*@ 3737 DMPlexCreateFromCellList - Deprecated, use DMPlexCreateFromCellListPetsc() 3738 3739 Level: deprecated 3740 3741 .seealso: DMPlexCreateFromCellListPetsc() 3742 @*/ 3743 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) 3744 { 3745 PetscErrorCode ierr; 3746 PetscInt i; 3747 PetscInt *pintCells; 3748 PetscReal *prealVC; 3749 3750 PetscFunctionBegin; 3751 PetscCheckFalse(sizeof(int) > sizeof(PetscInt),comm, PETSC_ERR_ARG_SIZ, "Size of int %zd greater than size of PetscInt %zd. Reconfigure PETSc --with-64-bit-indices=1", sizeof(int), sizeof(PetscInt)); 3752 if (sizeof(int) == sizeof(PetscInt)) { 3753 pintCells = (PetscInt *) cells; 3754 } else { 3755 ierr = PetscMalloc1(numCells*numCorners, &pintCells);CHKERRQ(ierr); 3756 for (i = 0; i < numCells*numCorners; i++) { 3757 pintCells[i] = (PetscInt) cells[i]; 3758 } 3759 } 3760 PetscCheckFalse(sizeof(double) > sizeof(PetscReal),comm, PETSC_ERR_ARG_SIZ, "Size of double %zd greater than size of PetscReal %zd. Reconfigure PETSc --with-precision=<higher precision>.", sizeof(double), sizeof(PetscReal)); 3761 if (sizeof(double) == sizeof(PetscReal)) { 3762 prealVC = (PetscReal *) vertexCoords; 3763 } else { 3764 ierr = PetscMalloc1(numVertices*spaceDim, &prealVC);CHKERRQ(ierr); 3765 for (i = 0; i < numVertices*spaceDim; i++) { 3766 prealVC[i] = (PetscReal) vertexCoords[i]; 3767 } 3768 } 3769 ierr = DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, pintCells, spaceDim, prealVC, dm);CHKERRQ(ierr); 3770 if (sizeof(int) != sizeof(PetscInt)) { 3771 ierr = PetscFree(pintCells);CHKERRQ(ierr); 3772 } 3773 if (sizeof(double) != sizeof(PetscReal)) { 3774 ierr = PetscFree(prealVC);CHKERRQ(ierr); 3775 } 3776 PetscFunctionReturn(0); 3777 } 3778 3779 /*@ 3780 DMPlexCreateFromDAG - This takes as input the adjacency-list representation of the Directed Acyclic Graph (Hasse Diagram) encoding a mesh, and produces a DM 3781 3782 Input Parameters: 3783 + dm - The empty DM object, usually from DMCreate() and DMSetDimension() 3784 . depth - The depth of the DAG 3785 . numPoints - Array of size depth + 1 containing the number of points at each depth 3786 . coneSize - The cone size of each point 3787 . cones - The concatenation of the cone points for each point, the cone list must be oriented correctly for each point 3788 . coneOrientations - The orientation of each cone point 3789 - vertexCoords - An array of numPoints[0]*spacedim numbers representing the coordinates of each vertex, with spacedim the value set via DMSetCoordinateDim() 3790 3791 Output Parameter: 3792 . dm - The DM 3793 3794 Note: Two triangles sharing a face would have input 3795 $ depth = 1, numPoints = [4 2], coneSize = [3 3 0 0 0 0] 3796 $ cones = [2 3 4 3 5 4], coneOrientations = [0 0 0 0 0 0] 3797 $ vertexCoords = [-1.0 0.0 0.0 -1.0 0.0 1.0 1.0 0.0] 3798 $ 3799 which would result in the DMPlex 3800 $ 3801 $ 4 3802 $ / | \ 3803 $ / | \ 3804 $ / | \ 3805 $ 2 0 | 1 5 3806 $ \ | / 3807 $ \ | / 3808 $ \ | / 3809 $ 3 3810 $ 3811 $ Notice that all points are numbered consecutively, unlike DMPlexCreateFromCellListPetsc() 3812 3813 Level: advanced 3814 3815 .seealso: DMPlexCreateFromCellListPetsc(), DMPlexCreate() 3816 @*/ 3817 PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[]) 3818 { 3819 Vec coordinates; 3820 PetscSection coordSection; 3821 PetscScalar *coords; 3822 PetscInt coordSize, firstVertex = -1, pStart = 0, pEnd = 0, p, v, dim, dimEmbed, d, off; 3823 PetscErrorCode ierr; 3824 3825 PetscFunctionBegin; 3826 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3827 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 3828 PetscCheckFalse(dimEmbed < dim,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Embedding dimension %D cannot be less than intrinsic dimension %d",dimEmbed,dim); 3829 for (d = 0; d <= depth; ++d) pEnd += numPoints[d]; 3830 ierr = DMPlexSetChart(dm, pStart, pEnd);CHKERRQ(ierr); 3831 for (p = pStart; p < pEnd; ++p) { 3832 ierr = DMPlexSetConeSize(dm, p, coneSize[p-pStart]);CHKERRQ(ierr); 3833 if (firstVertex < 0 && !coneSize[p - pStart]) { 3834 firstVertex = p - pStart; 3835 } 3836 } 3837 PetscCheckFalse(firstVertex < 0 && numPoints[0],PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected %D vertices but could not find any", numPoints[0]); 3838 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */ 3839 for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) { 3840 ierr = DMPlexSetCone(dm, p, &cones[off]);CHKERRQ(ierr); 3841 ierr = DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);CHKERRQ(ierr); 3842 } 3843 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr); 3844 ierr = DMPlexStratify(dm);CHKERRQ(ierr); 3845 /* Build coordinates */ 3846 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 3847 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 3848 ierr = PetscSectionSetFieldComponents(coordSection, 0, dimEmbed);CHKERRQ(ierr); 3849 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0]);CHKERRQ(ierr); 3850 for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) { 3851 ierr = PetscSectionSetDof(coordSection, v, dimEmbed);CHKERRQ(ierr); 3852 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed);CHKERRQ(ierr); 3853 } 3854 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 3855 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 3856 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 3857 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 3858 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 3859 ierr = VecSetBlockSize(coordinates, dimEmbed);CHKERRQ(ierr); 3860 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr); 3861 if (vertexCoords) { 3862 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 3863 for (v = 0; v < numPoints[0]; ++v) { 3864 PetscInt off; 3865 3866 ierr = PetscSectionGetOffset(coordSection, v+firstVertex, &off);CHKERRQ(ierr); 3867 for (d = 0; d < dimEmbed; ++d) { 3868 coords[off+d] = vertexCoords[v*dimEmbed+d]; 3869 } 3870 } 3871 } 3872 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 3873 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 3874 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 3875 PetscFunctionReturn(0); 3876 } 3877 3878 /*@C 3879 DMPlexCreateCellVertexFromFile - Create a DMPlex mesh from a simple cell-vertex file. 3880 3881 + comm - The MPI communicator 3882 . filename - Name of the .dat file 3883 - interpolate - Create faces and edges in the mesh 3884 3885 Output Parameter: 3886 . dm - The DM object representing the mesh 3887 3888 Note: The format is the simplest possible: 3889 $ Ne 3890 $ v0 v1 ... vk 3891 $ Nv 3892 $ x y z marker 3893 3894 Level: beginner 3895 3896 .seealso: DMPlexCreateFromFile(), DMPlexCreateMedFromFile(), DMPlexCreateGmsh(), DMPlexCreate() 3897 @*/ 3898 PetscErrorCode DMPlexCreateCellVertexFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 3899 { 3900 DMLabel marker; 3901 PetscViewer viewer; 3902 Vec coordinates; 3903 PetscSection coordSection; 3904 PetscScalar *coords; 3905 char line[PETSC_MAX_PATH_LEN]; 3906 PetscInt dim = 3, cdim = 3, coordSize, v, c, d; 3907 PetscMPIInt rank; 3908 int snum, Nv, Nc, Ncn, Nl; 3909 PetscErrorCode ierr; 3910 3911 PetscFunctionBegin; 3912 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 3913 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 3914 ierr = PetscViewerSetType(viewer, PETSCVIEWERASCII);CHKERRQ(ierr); 3915 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 3916 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 3917 if (rank == 0) { 3918 ierr = PetscViewerRead(viewer, line, 4, NULL, PETSC_STRING);CHKERRQ(ierr); 3919 snum = sscanf(line, "%d %d %d %d", &Nc, &Nv, &Ncn, &Nl); 3920 PetscCheckFalse(snum != 4,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 3921 } else { 3922 Nc = Nv = Ncn = Nl = 0; 3923 } 3924 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 3925 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 3926 ierr = DMPlexSetChart(*dm, 0, Nc+Nv);CHKERRQ(ierr); 3927 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 3928 ierr = DMSetCoordinateDim(*dm, cdim);CHKERRQ(ierr); 3929 /* Read topology */ 3930 if (rank == 0) { 3931 char format[PETSC_MAX_PATH_LEN]; 3932 PetscInt cone[8]; 3933 int vbuf[8], v; 3934 3935 for (c = 0; c < Ncn; ++c) {format[c*3+0] = '%'; format[c*3+1] = 'd'; format[c*3+2] = ' ';} 3936 format[Ncn*3-1] = '\0'; 3937 for (c = 0; c < Nc; ++c) {ierr = DMPlexSetConeSize(*dm, c, Ncn);CHKERRQ(ierr);} 3938 ierr = DMSetUp(*dm);CHKERRQ(ierr); 3939 for (c = 0; c < Nc; ++c) { 3940 ierr = PetscViewerRead(viewer, line, Ncn, NULL, PETSC_STRING);CHKERRQ(ierr); 3941 switch (Ncn) { 3942 case 2: snum = sscanf(line, format, &vbuf[0], &vbuf[1]);break; 3943 case 3: snum = sscanf(line, format, &vbuf[0], &vbuf[1], &vbuf[2]);break; 3944 case 4: snum = sscanf(line, format, &vbuf[0], &vbuf[1], &vbuf[2], &vbuf[3]);break; 3945 case 6: snum = sscanf(line, format, &vbuf[0], &vbuf[1], &vbuf[2], &vbuf[3], &vbuf[4], &vbuf[5]);break; 3946 case 8: snum = sscanf(line, format, &vbuf[0], &vbuf[1], &vbuf[2], &vbuf[3], &vbuf[4], &vbuf[5], &vbuf[6], &vbuf[7]);break; 3947 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No cell shape with %D vertices", Ncn); 3948 } 3949 PetscCheckFalse(snum != Ncn,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 3950 for (v = 0; v < Ncn; ++v) cone[v] = vbuf[v] + Nc; 3951 /* Hexahedra are inverted */ 3952 if (Ncn == 8) { 3953 PetscInt tmp = cone[1]; 3954 cone[1] = cone[3]; 3955 cone[3] = tmp; 3956 } 3957 ierr = DMPlexSetCone(*dm, c, cone);CHKERRQ(ierr); 3958 } 3959 } 3960 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 3961 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 3962 /* Read coordinates */ 3963 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 3964 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 3965 ierr = PetscSectionSetFieldComponents(coordSection, 0, cdim);CHKERRQ(ierr); 3966 ierr = PetscSectionSetChart(coordSection, Nc, Nc + Nv);CHKERRQ(ierr); 3967 for (v = Nc; v < Nc+Nv; ++v) { 3968 ierr = PetscSectionSetDof(coordSection, v, cdim);CHKERRQ(ierr); 3969 ierr = PetscSectionSetFieldDof(coordSection, v, 0, cdim);CHKERRQ(ierr); 3970 } 3971 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 3972 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 3973 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 3974 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 3975 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 3976 ierr = VecSetBlockSize(coordinates, cdim);CHKERRQ(ierr); 3977 ierr = VecSetType(coordinates, VECSTANDARD);CHKERRQ(ierr); 3978 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 3979 if (rank == 0) { 3980 char format[PETSC_MAX_PATH_LEN]; 3981 double x[3]; 3982 int l, val[3]; 3983 3984 if (Nl) { 3985 for (l = 0; l < Nl; ++l) {format[l*3+0] = '%'; format[l*3+1] = 'd'; format[l*3+2] = ' ';} 3986 format[Nl*3-1] = '\0'; 3987 ierr = DMCreateLabel(*dm, "marker");CHKERRQ(ierr); 3988 ierr = DMGetLabel(*dm, "marker", &marker);CHKERRQ(ierr); 3989 } 3990 for (v = 0; v < Nv; ++v) { 3991 ierr = PetscViewerRead(viewer, line, 3+Nl, NULL, PETSC_STRING);CHKERRQ(ierr); 3992 snum = sscanf(line, "%lg %lg %lg", &x[0], &x[1], &x[2]); 3993 PetscCheckFalse(snum != 3,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 3994 switch (Nl) { 3995 case 0: snum = 0;break; 3996 case 1: snum = sscanf(line, format, &val[0]);break; 3997 case 2: snum = sscanf(line, format, &val[0], &val[1]);break; 3998 case 3: snum = sscanf(line, format, &val[0], &val[1], &val[2]);break; 3999 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Request support for %D labels", Nl); 4000 } 4001 PetscCheckFalse(snum != Nl,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line); 4002 for (d = 0; d < cdim; ++d) coords[v*cdim+d] = x[d]; 4003 for (l = 0; l < Nl; ++l) {ierr = DMLabelSetValue(marker, v+Nc, val[l]);CHKERRQ(ierr);} 4004 } 4005 } 4006 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 4007 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 4008 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 4009 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 4010 if (interpolate) { 4011 DM idm; 4012 DMLabel bdlabel; 4013 4014 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 4015 ierr = DMDestroy(dm);CHKERRQ(ierr); 4016 *dm = idm; 4017 4018 if (!Nl) { 4019 ierr = DMCreateLabel(*dm, "marker");CHKERRQ(ierr); 4020 ierr = DMGetLabel(*dm, "marker", &bdlabel);CHKERRQ(ierr); 4021 ierr = DMPlexMarkBoundaryFaces(*dm, PETSC_DETERMINE, bdlabel);CHKERRQ(ierr); 4022 ierr = DMPlexLabelComplete(*dm, bdlabel);CHKERRQ(ierr); 4023 } 4024 } 4025 PetscFunctionReturn(0); 4026 } 4027 4028 /*@C 4029 DMPlexCreateFromFile - This takes a filename and produces a DM 4030 4031 Input Parameters: 4032 + comm - The communicator 4033 . filename - A file name 4034 . plexname - The object name of the resulting DM, also used for intra-datafile lookup by some formats 4035 - interpolate - Flag to create intermediate mesh pieces (edges, faces) 4036 4037 Output Parameter: 4038 . dm - The DM 4039 4040 Options Database Keys: 4041 . -dm_plex_create_from_hdf5_xdmf - use the PETSC_VIEWER_HDF5_XDMF format for reading HDF5 4042 4043 Use -dm_plex_create_ prefix to pass options to the internal PetscViewer, e.g. 4044 $ -dm_plex_create_viewer_hdf5_collective 4045 4046 Notes: 4047 Using PETSCVIEWERHDF5 type with PETSC_VIEWER_HDF5_PETSC format, one can save multiple DMPlex 4048 meshes in a single HDF5 file. This in turn requires one to name the DMPlex object with PetscObjectSetName() 4049 before saving it with DMView() and before loading it with DMLoad() for identification of the mesh object. 4050 The input parameter name is thus used to name the DMPlex object when DMPlexCreateFromFile() internally 4051 calls DMLoad(). Currently, name is ignored for other viewer types and/or formats. 4052 4053 Level: beginner 4054 4055 .seealso: DMPlexCreateFromDAG(), DMPlexCreateFromCellListPetsc(), DMPlexCreate(), PetscObjectSetName(), DMView(), DMLoad() 4056 @*/ 4057 PetscErrorCode DMPlexCreateFromFile(MPI_Comm comm, const char filename[], const char plexname[], PetscBool interpolate, DM *dm) 4058 { 4059 const char *extGmsh = ".msh"; 4060 const char *extGmsh2 = ".msh2"; 4061 const char *extGmsh4 = ".msh4"; 4062 const char *extCGNS = ".cgns"; 4063 const char *extExodus = ".exo"; 4064 const char *extGenesis = ".gen"; 4065 const char *extFluent = ".cas"; 4066 const char *extHDF5 = ".h5"; 4067 const char *extMed = ".med"; 4068 const char *extPLY = ".ply"; 4069 const char *extEGADSLite = ".egadslite"; 4070 const char *extEGADS = ".egads"; 4071 const char *extIGES = ".igs"; 4072 const char *extSTEP = ".stp"; 4073 const char *extCV = ".dat"; 4074 size_t len; 4075 PetscBool isGmsh, isGmsh2, isGmsh4, isCGNS, isExodus, isGenesis, isFluent, isHDF5, isMed, isPLY, isEGADSLite, isEGADS, isIGES, isSTEP, isCV; 4076 PetscMPIInt rank; 4077 PetscErrorCode ierr; 4078 4079 PetscFunctionBegin; 4080 PetscValidCharPointer(filename, 2); 4081 if (plexname) PetscValidCharPointer(plexname, 3); 4082 PetscValidPointer(dm, 5); 4083 ierr = DMInitializePackage();CHKERRQ(ierr); 4084 ierr = PetscLogEventBegin(DMPLEX_CreateFromFile,0,0,0,0);CHKERRQ(ierr); 4085 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 4086 ierr = PetscStrlen(filename, &len);CHKERRQ(ierr); 4087 PetscCheckFalse(!len,comm, PETSC_ERR_ARG_WRONG, "Filename must be a valid path"); 4088 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGmsh, 4, &isGmsh);CHKERRQ(ierr); 4089 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extGmsh2, 5, &isGmsh2);CHKERRQ(ierr); 4090 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extGmsh4, 5, &isGmsh4);CHKERRQ(ierr); 4091 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extCGNS, 5, &isCGNS);CHKERRQ(ierr); 4092 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extExodus, 4, &isExodus);CHKERRQ(ierr); 4093 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGenesis, 4, &isGenesis);CHKERRQ(ierr); 4094 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extFluent, 4, &isFluent);CHKERRQ(ierr); 4095 ierr = PetscStrncmp(&filename[PetscMax(0,len-3)], extHDF5, 3, &isHDF5);CHKERRQ(ierr); 4096 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extMed, 4, &isMed);CHKERRQ(ierr); 4097 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extPLY, 4, &isPLY);CHKERRQ(ierr); 4098 ierr = PetscStrncmp(&filename[PetscMax(0,len-10)], extEGADSLite, 10, &isEGADSLite);CHKERRQ(ierr); 4099 ierr = PetscStrncmp(&filename[PetscMax(0,len-6)], extEGADS, 6, &isEGADS);CHKERRQ(ierr); 4100 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extIGES, 4, &isIGES);CHKERRQ(ierr); 4101 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extSTEP, 4, &isSTEP);CHKERRQ(ierr); 4102 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extCV, 4, &isCV);CHKERRQ(ierr); 4103 if (isGmsh || isGmsh2 || isGmsh4) { 4104 ierr = DMPlexCreateGmshFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 4105 } else if (isCGNS) { 4106 ierr = DMPlexCreateCGNSFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 4107 } else if (isExodus || isGenesis) { 4108 ierr = DMPlexCreateExodusFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 4109 } else if (isFluent) { 4110 ierr = DMPlexCreateFluentFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 4111 } else if (isHDF5) { 4112 PetscBool load_hdf5_xdmf = PETSC_FALSE; 4113 PetscViewer viewer; 4114 4115 /* PETSC_VIEWER_HDF5_XDMF is used if the filename ends with .xdmf.h5, or if -dm_plex_create_from_hdf5_xdmf option is present */ 4116 ierr = PetscStrncmp(&filename[PetscMax(0,len-8)], ".xdmf", 5, &load_hdf5_xdmf);CHKERRQ(ierr); 4117 ierr = PetscOptionsGetBool(NULL, NULL, "-dm_plex_create_from_hdf5_xdmf", &load_hdf5_xdmf, NULL);CHKERRQ(ierr); 4118 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 4119 ierr = PetscViewerSetType(viewer, PETSCVIEWERHDF5);CHKERRQ(ierr); 4120 ierr = PetscViewerSetOptionsPrefix(viewer, "dm_plex_create_");CHKERRQ(ierr); 4121 ierr = PetscViewerSetFromOptions(viewer);CHKERRQ(ierr); 4122 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 4123 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 4124 4125 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 4126 ierr = PetscObjectSetName((PetscObject)(*dm), plexname);CHKERRQ(ierr); 4127 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 4128 if (load_hdf5_xdmf) {ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_XDMF);CHKERRQ(ierr);} 4129 ierr = DMLoad(*dm, viewer);CHKERRQ(ierr); 4130 if (load_hdf5_xdmf) {ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);} 4131 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 4132 4133 if (interpolate) { 4134 DM idm; 4135 4136 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 4137 ierr = DMDestroy(dm);CHKERRQ(ierr); 4138 *dm = idm; 4139 } 4140 } else if (isMed) { 4141 ierr = DMPlexCreateMedFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 4142 } else if (isPLY) { 4143 ierr = DMPlexCreatePLYFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 4144 } else if (isEGADSLite || isEGADS || isIGES || isSTEP) { 4145 if (isEGADSLite) {ierr = DMPlexCreateEGADSLiteFromFile(comm, filename, dm);CHKERRQ(ierr);} 4146 else {ierr = DMPlexCreateEGADSFromFile(comm, filename, dm);CHKERRQ(ierr);} 4147 if (!interpolate) { 4148 DM udm; 4149 4150 ierr = DMPlexUninterpolate(*dm, &udm);CHKERRQ(ierr); 4151 ierr = DMDestroy(dm);CHKERRQ(ierr); 4152 *dm = udm; 4153 } 4154 } else if (isCV) { 4155 ierr = DMPlexCreateCellVertexFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr); 4156 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot load file %s: unrecognized extension", filename); 4157 ierr = PetscStrlen(plexname, &len);CHKERRQ(ierr); 4158 if (len) {ierr = PetscObjectSetName((PetscObject)(*dm), plexname);CHKERRQ(ierr);} 4159 ierr = PetscLogEventEnd(DMPLEX_CreateFromFile,0,0,0,0);CHKERRQ(ierr); 4160 PetscFunctionReturn(0); 4161 } 4162