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