DMPlexConstructGhostCells

Synopsis

#include "petscdmplex.h"   
PetscErrorCode DMPlexConstructCohesiveCells(DM dm, DMLabel label, DM *dmSplit)

Input Parameters

Output Parameters

Note: If no label exists of that name, one will be created marking all boundary faces

See Also

PetscFunctionBegin; PetscValidHeaderSpecific(dm, DM_CLASSID, 1); PetscValidPointer(numGhostCells, 3); PetscValidPointer(dmGhosted, 4); ierr = DMCreate(PetscObjectComm((PetscObject)dm), &gdm);CHKERRQ(ierr); ierr = DMSetType(gdm, DMPLEX);CHKERRQ(ierr); ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); ierr = DMPlexSetDimension(gdm, dim);CHKERRQ(ierr); ierr = DMPlexGetLabel(dm, name, &label);CHKERRQ(ierr); if (!label) { /* Get label for boundary faces */ ierr = DMPlexCreateLabel(dm, name);CHKERRQ(ierr); ierr = DMPlexGetLabel(dm, name, &label);CHKERRQ(ierr); ierr = DMPlexMarkBoundaryFaces(dm, label);CHKERRQ(ierr); } ierr = DMPlexConstructGhostCells_Internal(dm, label, numGhostCells, gdm);CHKERRQ(ierr); ierr = DMSetFromOptions(gdm);CHKERRQ(ierr); *dmGhosted = gdm; PetscFunctionReturn(0); }

#undef __FUNCT__ #define __FUNCT__ "DMPlexConstructCohesiveCells_Internal" static PetscErrorCode DMPlexConstructCohesiveCells_Internal(DM dm, DMLabel label, DM sdm) { MPI_Comm comm; IS valueIS, *pointIS; const PetscInt *values, **splitPoints; PetscSection coordSection; Vec coordinates; PetscScalar *coords; PetscInt *depthShift, *depthOffset, *pMaxNew, *numSplitPoints, *coneNew, *supportNew; PetscInt shift = 100, depth = 0, dep, dim, d, numSP = 0, sp, maxConeSize, maxSupportSize, numLabels, p, v; PetscErrorCode ierr;

PetscFunctionBegin; ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); /* Count split points and add cohesive cells */ if (label) { ierr = DMLabelGetValueIS(label, &valueIS);CHKERRQ(ierr); ierr = ISGetLocalSize(valueIS, &numSP);CHKERRQ(ierr); ierr = ISGetIndices(valueIS, &values);CHKERRQ(ierr); } ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); ierr = DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);CHKERRQ(ierr); ierr = PetscMalloc5(depth+1,PetscInt,&depthShift,depth+1,PetscInt,&depthOffset,depth+1,PetscInt,&pMaxNew,maxConeSize*3,PetscInt,&coneNew,maxSupportSize,PetscInt,&supportNew);CHKERRQ(ierr); ierr = PetscMalloc3(depth+1,IS,&pointIS,depth+1,PetscInt,&numSplitPoints,depth+1,const PetscInt*,&splitPoints);CHKERRQ(ierr); ierr = PetscMemzero(depthShift, (depth+1) * sizeof(PetscInt));CHKERRQ(ierr); for (d = 0; d <= depth; ++d) { ierr = DMPlexGetDepthStratum(dm, d, NULL, &pMaxNew[d]);CHKERRQ(ierr); numSplitPoints[d] = 0; splitPoints[d] = NULL; pointIS[d] = NULL; } for (sp = 0; sp < numSP; ++sp) { const PetscInt dep = values[sp];

if ((dep < 0) || (dep > depth)) continue; ierr = DMLabelGetStratumSize(label, dep, &depthShift[dep]);CHKERRQ(ierr); ierr = DMLabelGetStratumIS(label, dep, &pointIS[dep]);CHKERRQ(ierr); if (pointIS[dep]) { ierr = ISGetLocalSize(pointIS[dep], &numSplitPoints[dep]);CHKERRQ(ierr); ierr = ISGetIndices(pointIS[dep], &splitPoints[dep]);CHKERRQ(ierr); } } if (depth >= 0) { /* Calculate number of additional points */ depthShift[depth] = depthShift[depth-1]; /* There is a cohesive cell for every split face */ depthShift[1] += depthShift[0]; /* There is a cohesive edge for every split vertex */ /* Calculate hybrid bound for each dimension */ pMaxNew[0] += depthShift[depth]; if (depth > 1) pMaxNew[dim-1] += depthShift[depth] + depthShift[0]; if (depth > 2) pMaxNew[1] += depthShift[depth] + depthShift[0] + depthShift[dim-1];

/* Calculate point offset for each dimension */ depthOffset[depth] = 0; depthOffset[0] = depthOffset[depth] + depthShift[depth]; if (depth > 1) depthOffset[dim-1] = depthOffset[0] + depthShift[0]; if (depth > 2) depthOffset[1] = depthOffset[dim-1] + depthShift[dim-1]; } ierr = DMPlexShiftSizes_Internal(dm, depthShift, sdm);CHKERRQ(ierr); /* Step 3: Set cone/support sizes for new points */ for (dep = 0; dep <= depth; ++dep) { for (p = 0; p < numSplitPoints[dep]; ++p) { const PetscInt oldp = splitPoints[dep][p]; const PetscInt newp = depthOffset[dep] + oldp; const PetscInt splitp = pMaxNew[dep] + p; const PetscInt *support; PetscInt coneSize, supportSize, q, e;

ierr = DMPlexGetConeSize(dm, oldp, &coneSize);CHKERRQ(ierr); ierr = DMPlexSetConeSize(sdm, splitp, coneSize);CHKERRQ(ierr); ierr = DMPlexGetSupportSize(dm, oldp, &supportSize);CHKERRQ(ierr); ierr = DMPlexSetSupportSize(sdm, splitp, supportSize);CHKERRQ(ierr); if (dep == depth-1) { const PetscInt ccell = pMaxNew[depth] + p; /* Add cohesive cells, they are prisms */ ierr = DMPlexSetConeSize(sdm, ccell, 2 + coneSize);CHKERRQ(ierr); } else if (dep == 0) { const PetscInt cedge = pMaxNew[1] + (depthShift[1] - depthShift[0]) + p;

ierr = DMPlexGetSupport(dm, oldp, &support);CHKERRQ(ierr); /* Split old vertex: Edges in old split faces and new cohesive edge */ for (e = 0, q = 0; e < supportSize; ++e) { PetscInt val;

ierr = DMLabelGetValue(label, support[e], &val);CHKERRQ(ierr); if ((val == 1) || (val == (shift + 1))) ++q; } ierr = DMPlexSetSupportSize(sdm, newp, q+1);CHKERRQ(ierr); /* Split new vertex: Edges in new split faces and new cohesive edge */ for (e = 0, q = 0; e < supportSize; ++e) { PetscInt val;

ierr = DMLabelGetValue(label, support[e], &val);CHKERRQ(ierr); if ((val == 1) || (val == -(shift + 1))) ++q; } ierr = DMPlexSetSupportSize(sdm, splitp, q+1);CHKERRQ(ierr); /* Add cohesive edges */ ierr = DMPlexSetConeSize(sdm, cedge, 2);CHKERRQ(ierr); /* Punt for now on support, you loop over closure, extract faces, check which ones are in the label */ } else if (dep == dim-2) { ierr = DMPlexGetSupport(dm, oldp, &support);CHKERRQ(ierr); /* Split old edge: Faces in positive side cells and old split faces */ for (e = 0, q = 0; e < supportSize; ++e) { PetscInt val;

ierr = DMLabelGetValue(label, support[e], &val);CHKERRQ(ierr); if ((val == dim-1) || (val == (shift + dim-1))) ++q; } ierr = DMPlexSetSupportSize(sdm, newp, q);CHKERRQ(ierr); /* Split new edge: Faces in negative side cells and new split faces */ for (e = 0, q = 0; e < supportSize; ++e) { PetscInt val;

ierr = DMLabelGetValue(label, support[e], &val);CHKERRQ(ierr); if ((val == dim-1) || (val == -(shift + dim-1))) ++q; } ierr = DMPlexSetSupportSize(sdm, splitp, q);CHKERRQ(ierr); } } } /* Step 4: Setup split DM */ ierr = DMSetUp(sdm);CHKERRQ(ierr); ierr = DMPlexShiftPoints_Internal(dm, depthShift, sdm);CHKERRQ(ierr); /* Step 6: Set cones and supports for new points */ for (dep = 0; dep <= depth; ++dep) { for (p = 0; p < numSplitPoints[dep]; ++p) { const PetscInt oldp = splitPoints[dep][p]; const PetscInt newp = depthOffset[dep] + oldp; const PetscInt splitp = pMaxNew[dep] + p; const PetscInt *cone, *support, *ornt; PetscInt coneSize, supportSize, q, v, e, s;

ierr = DMPlexGetConeSize(dm, oldp, &coneSize);CHKERRQ(ierr); ierr = DMPlexGetCone(dm, oldp, &cone);CHKERRQ(ierr); ierr = DMPlexGetConeOrientation(dm, oldp, &ornt);CHKERRQ(ierr); ierr = DMPlexGetSupportSize(dm, oldp, &supportSize);CHKERRQ(ierr); ierr = DMPlexGetSupport(dm, oldp, &support);CHKERRQ(ierr); if (dep == depth-1) { const PetscInt ccell = pMaxNew[depth] + p; const PetscInt *supportF;

/* Split face: copy in old face to new face to start */ ierr = DMPlexGetSupport(sdm, newp, &supportF);CHKERRQ(ierr); ierr = DMPlexSetSupport(sdm, splitp, supportF);CHKERRQ(ierr); /* Split old face: old vertices/edges in cone so no change */ /* Split new face: new vertices/edges in cone */ for (q = 0; q < coneSize; ++q) { ierr = PetscFindInt(cone[q], numSplitPoints[dim-2], splitPoints[dim-2], &v);CHKERRQ(ierr);

coneNew[2+q] = pMaxNew[dim-2] + v; } ierr = DMPlexSetCone(sdm, splitp, &coneNew[2]);CHKERRQ(ierr); ierr = DMPlexSetConeOrientation(sdm, splitp, ornt);CHKERRQ(ierr); /* Cohesive cell: Old and new split face, then new cohesive edges */ coneNew[0] = newp; coneNew[1] = splitp; for (q = 0; q < coneSize; ++q) { coneNew[2+q] = (pMaxNew[1] - pMaxNew[dim-2]) + (depthShift[1] - depthShift[0]) + coneNew[2+q]; } ierr = DMPlexSetCone(sdm, ccell, coneNew);CHKERRQ(ierr);

for (s = 0; s < supportSize; ++s) { PetscInt val;

ierr = DMLabelGetValue(label, support[s], &val);CHKERRQ(ierr); if (val < 0) { /* Split old face: Replace negative side cell with cohesive cell */ ierr = DMPlexInsertSupport(sdm, newp, s, ccell);CHKERRQ(ierr); } else { /* Split new face: Replace positive side cell with cohesive cell */ ierr = DMPlexInsertSupport(sdm, splitp, s, ccell);CHKERRQ(ierr); } } } else if (dep == 0) { const PetscInt cedge = pMaxNew[1] + (depthShift[1] - depthShift[0]) + p;

/* Split old vertex: Edges in old split faces and new cohesive edge */ for (e = 0, q = 0; e < supportSize; ++e) { PetscInt val;

ierr = DMLabelGetValue(label, support[e], &val);CHKERRQ(ierr); if ((val == 1) || (val == (shift + 1))) { supportNew[q++] = depthOffset[1] + support[e]; } } supportNew[q] = cedge;

ierr = DMPlexSetSupport(sdm, newp, supportNew);CHKERRQ(ierr); /* Split new vertex: Edges in new split faces and new cohesive edge */ for (e = 0, q = 0; e < supportSize; ++e) { PetscInt val, edge;

ierr = DMLabelGetValue(label, support[e], &val);CHKERRQ(ierr); if (val == 1) { ierr = PetscFindInt(support[e], numSplitPoints[1], splitPoints[1], &edge);CHKERRQ(ierr); if (edge < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Edge %d is not a split edge", support[e]); supportNew[q++] = pMaxNew[1] + edge; } else if (val == -(shift + 1)) { supportNew[q++] = depthOffset[1] + support[e]; } } supportNew[q] = cedge; ierr = DMPlexSetSupport(sdm, splitp, supportNew);CHKERRQ(ierr); /* Cohesive edge: Old and new split vertex, punting on support */ coneNew[0] = newp; coneNew[1] = splitp; ierr = DMPlexSetCone(sdm, cedge, coneNew);CHKERRQ(ierr); } else if (dep == dim-2) { /* Split old edge: old vertices in cone so no change */ /* Split new edge: new vertices in cone */ for (q = 0; q < coneSize; ++q) { ierr = PetscFindInt(cone[q], numSplitPoints[dim-3], splitPoints[dim-3], &v);CHKERRQ(ierr);

coneNew[q] = pMaxNew[dim-3] + v; } ierr = DMPlexSetCone(sdm, splitp, coneNew);CHKERRQ(ierr); /* Split old edge: Faces in positive side cells and old split faces */ for (e = 0, q = 0; e < supportSize; ++e) { PetscInt val;

ierr = DMLabelGetValue(label, support[e], &val);CHKERRQ(ierr); if ((val == dim-1) || (val == (shift + dim-1))) { supportNew[q++] = depthOffset[dim-1] + support[e]; } } ierr = DMPlexSetSupport(sdm, newp, supportNew);CHKERRQ(ierr); /* Split new edge: Faces in negative side cells and new split faces */ for (e = 0, q = 0; e < supportSize; ++e) { PetscInt val, face;

ierr = DMLabelGetValue(label, support[e], &val);CHKERRQ(ierr); if (val == dim-1) { ierr = PetscFindInt(support[e], numSplitPoints[dim-1], splitPoints[dim-1], &face);CHKERRQ(ierr); if (face < 0) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Face %d is not a split face", support[e]); supportNew[q++] = pMaxNew[dim-1] + face; } else if (val == -(shift + dim-1)) { supportNew[q++] = depthOffset[dim-1] + support[e]; } } ierr = DMPlexSetSupport(sdm, splitp, supportNew);CHKERRQ(ierr); } } } /* Step 6b: Replace split points in negative side cones */ for (sp = 0; sp < numSP; ++sp) { PetscInt dep = values[sp]; IS pIS; PetscInt numPoints; const PetscInt *points;

if (dep >= 0) continue; ierr = DMLabelGetStratumIS(label, dep, &pIS);CHKERRQ(ierr); if (!pIS) continue; dep = -dep - shift; ierr = ISGetLocalSize(pIS, &numPoints);CHKERRQ(ierr); ierr = ISGetIndices(pIS, &points);CHKERRQ(ierr); for (p = 0; p < numPoints; ++p) { const PetscInt oldp = points[p]; const PetscInt newp = depthOffset[dep] + oldp; const PetscInt *cone; PetscInt coneSize, c; PetscBool replaced = PETSC_FALSE;

/* Negative edge: replace split vertex */ /* Negative cell: replace split face */ ierr = DMPlexGetConeSize(sdm, newp, &coneSize);CHKERRQ(ierr); ierr = DMPlexGetCone(sdm, newp, &cone);CHKERRQ(ierr); for (c = 0; c < coneSize; ++c) { const PetscInt coldp = cone[c] - depthOffset[dep-1]; PetscInt csplitp, cp, val;

ierr = DMLabelGetValue(label, coldp, &val);CHKERRQ(ierr); if (val == dep-1) { ierr = PetscFindInt(coldp, numSplitPoints[dep-1], splitPoints[dep-1], &cp);CHKERRQ(ierr); if (cp < 0) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "Point %d is not a split point of dimension %d", oldp, dep-1); csplitp = pMaxNew[dep-1] + cp; ierr = DMPlexInsertCone(sdm, newp, c, csplitp);CHKERRQ(ierr); replaced = PETSC_TRUE; } } if (!replaced) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "The cone of point %d does not contain split points", oldp); } ierr = ISRestoreIndices(pIS, &points);CHKERRQ(ierr); ierr = ISDestroy(&pIS);CHKERRQ(ierr); } /* Step 7: Stratify */ ierr = DMPlexStratify(sdm);CHKERRQ(ierr); /* Step 8: Coordinates */ ierr = DMPlexShiftCoordinates_Internal(dm, depthShift, sdm);CHKERRQ(ierr); ierr = DMPlexGetCoordinateSection(sdm, &coordSection);CHKERRQ(ierr); ierr = DMGetCoordinatesLocal(sdm, &coordinates);CHKERRQ(ierr); ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); for (v = 0; v < (numSplitPoints ? numSplitPoints[0] : 0); ++v) { const PetscInt newp = depthOffset[0] + splitPoints[0][v]; const PetscInt splitp = pMaxNew[0] + v; PetscInt dof, off, soff, d;

ierr = PetscSectionGetDof(coordSection, newp, &dof);CHKERRQ(ierr); ierr = PetscSectionGetOffset(coordSection, newp, &off);CHKERRQ(ierr); ierr = PetscSectionGetOffset(coordSection, splitp, &soff);CHKERRQ(ierr); for (d = 0; d < dof; ++d) coords[soff+d] = coords[off+d]; } ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); /* Step 9: SF, if I can figure this out we can split the mesh in parallel */ ierr = DMPlexShiftSF_Internal(dm, depthShift, sdm);CHKERRQ(ierr); /* Step 10: Labels */ ierr = DMPlexShiftLabels_Internal(dm, depthShift, sdm);CHKERRQ(ierr); ierr = DMPlexGetNumLabels(sdm, &numLabels);CHKERRQ(ierr); for (dep = 0; dep <= depth; ++dep) { for (p = 0; p < numSplitPoints[dep]; ++p) { const PetscInt newp = depthOffset[dep] + splitPoints[dep][p]; const PetscInt splitp = pMaxNew[dep] + p; PetscInt l;

for (l = 0; l < numLabels; ++l) { DMLabel mlabel; const char *lname; PetscInt val;

ierr = DMPlexGetLabelName(sdm, l, &lname);CHKERRQ(ierr); ierr = DMPlexGetLabel(sdm, lname, &mlabel);CHKERRQ(ierr); ierr = DMLabelGetValue(mlabel, newp, &val);CHKERRQ(ierr); if (val >= 0) { ierr = DMLabelSetValue(mlabel, splitp, val);CHKERRQ(ierr); if (dep == 0) { const PetscInt cedge = pMaxNew[1] + (depthShift[1] - depthShift[0]) + p; ierr = DMLabelSetValue(mlabel, cedge, val);CHKERRQ(ierr); } } } } } for (sp = 0; sp < numSP; ++sp) { const PetscInt dep = values[sp];

if ((dep < 0) || (dep > depth)) continue; if (pointIS[dep]) {ierr = ISRestoreIndices(pointIS[dep], &splitPoints[dep]);CHKERRQ(ierr);} ierr = ISDestroy(&pointIS[dep]);CHKERRQ(ierr); } if (label) { ierr = ISRestoreIndices(valueIS, &values);CHKERRQ(ierr); ierr = ISDestroy(&valueIS);CHKERRQ(ierr); } ierr = PetscFree5(depthShift, depthOffset, pMaxNew, coneNew, supportNew);CHKERRQ(ierr); ierr = PetscFree3(pointIS, numSplitPoints, splitPoints);CHKERRQ(ierr); PetscFunctionReturn(0); }

#undef __FUNCT__ #define __FUNCT__ "DMPlexConstructCohesiveCells" /*@C DMPlexConstructCohesiveCells - Construct cohesive cells which split the face along an internal interface

Collective on dm

Input Parameters

Output Parameters