Skip to content
GitLab
Commit c8944ce2 authored by mattijs's avatar mattijs
Browse files

iso surface correction

parent 5a30dd1b
Hide whitespace changes
Inline Side-by-side
src/sampling/sampledSurface/isoSurface/isoSurface.C
View file @ c8944ce2
......@@ -85,9 +85,74 @@ bool Foam::isoSurface::isEdgeOfFaceCut
}
// Get neighbour value and position.
void Foam::isoSurface::getNeighbour
(
const labelList& boundaryRegion,
const volScalarField& cVals,
const label cellI,
const label faceI,
scalar& nbrValue,
point& nbrPoint
) const
{
const labelList& own = mesh_.faceOwner();
const labelList& nei = mesh_.faceNeighbour();
const surfaceScalarField& weights = mesh_.weights();
if (mesh_.isInternalFace(faceI))
{
label nbr = (own[faceI] == cellI ? nei[faceI] : own[faceI]);
nbrValue = cVals[nbr];
nbrPoint = mesh_.cellCentres()[nbr];
}
else
{
label bFaceI = faceI-mesh_.nInternalFaces();
label patchI = boundaryRegion[bFaceI];
const polyPatch& pp = mesh_.boundaryMesh()[patchI];
label patchFaceI = faceI-pp.start();
if (isA<emptyPolyPatch>(pp))
{
// Assume zero gradient
nbrValue = cVals[own[faceI]];
nbrPoint = mesh_.faceCentres()[faceI];
}
else if (pp.coupled())
{
if (!refCast<const coupledPolyPatch>(pp).separated())
{
// Behave as internal face:
// other side value
nbrValue = cVals.boundaryField()[patchI][patchFaceI];
// other side cell centre
nbrPoint = mesh_.C().boundaryField()[patchI][patchFaceI];
}
else
{
// Do some interpolation for now
const scalarField& w = weights.boundaryField()[patchI];
nbrPoint = mesh_.faceCentres()[faceI];
nbrValue =
(1-w[patchFaceI])*cVals[own[faceI]]
+ w[patchFaceI]*cVals.boundaryField()[patchI][patchFaceI];
}
}
else
{
nbrValue = cVals.boundaryField()[patchI][patchFaceI];
nbrPoint = mesh_.faceCentres()[faceI];
}
}
}
// Determine for every face/cell whether it (possibly) generates triangles.
void Foam::isoSurface::calcCutTypes
(
const labelList& boundaryRegion,
const volScalarField& cVals,
const scalarField& pVals
)
......@@ -103,7 +168,20 @@ void Foam::isoSurface::calcCutTypes
{
// CC edge.
bool ownLower = (cVals[own[faceI]] < iso_);
bool neiLower = (cVals[nei[faceI]] < iso_);
scalar nbrValue;
point nbrPoint;
getNeighbour
(
boundaryRegion,
cVals,
own[faceI],
faceI,
nbrValue,
nbrPoint
);
bool neiLower = (nbrValue < iso_);
if (ownLower != neiLower)
{
......@@ -129,15 +207,29 @@ void Foam::isoSurface::calcCutTypes
if (isA<emptyPolyPatch>(pp))
{
// Assume zero gradient so owner and neighbour/boundary value equal
// Still needs special treatment?
forAll(pp, i)
{
bool ownLower = (cVals[own[faceI]] < iso_);
scalar nbrValue;
point nbrPoint;
getNeighbour
(
boundaryRegion,
cVals,
own[faceI],
faceI,
nbrValue,
nbrPoint
);
bool neiLower = (nbrValue < iso_);
const face f = mesh_.faces()[faceI];
if (isEdgeOfFaceCut(pVals, f, ownLower, ownLower))
if (isEdgeOfFaceCut(pVals, f, ownLower, neiLower))
{
faceCutType_[faceI] = CUT;
}
......@@ -150,7 +242,20 @@ void Foam::isoSurface::calcCutTypes
forAll(pp, i)
{
bool ownLower = (cVals[own[faceI]] < iso_);
bool neiLower = (cVals.boundaryField()[patchI][i] < iso_);
scalar nbrValue;
point nbrPoint;
getNeighbour
(
boundaryRegion,
cVals,
own[faceI],
faceI,
nbrValue,
nbrPoint
);
bool neiLower = (nbrValue < iso_);
if (ownLower != neiLower)
{
......@@ -166,6 +271,7 @@ void Foam::isoSurface::calcCutTypes
faceCutType_[faceI] = CUT;
}
}
faceI++;
}
}
......@@ -331,70 +437,6 @@ Foam::pointIndexHit Foam::isoSurface::collapseSurface
}
// Get neighbour value and position.
void Foam::isoSurface::getNeighbour
(
const labelList& boundaryRegion,
const volScalarField& cVals,
const label cellI,
const label faceI,
scalar& nbrValue,
point& nbrPoint
) const
{
const labelList& own = mesh_.faceOwner();
const labelList& nei = mesh_.faceNeighbour();
const surfaceScalarField& weights = mesh_.weights();
if (mesh_.isInternalFace(faceI))
{
label nbr = (own[faceI] == cellI ? nei[faceI] : own[faceI]);
nbrValue = cVals[nbr];
nbrPoint = mesh_.cellCentres()[nbr];
}
else
{
label bFaceI = faceI-mesh_.nInternalFaces();
label patchI = boundaryRegion[bFaceI];
const polyPatch& pp = mesh_.boundaryMesh()[patchI];
label patchFaceI = faceI-pp.start();
if (isA<emptyPolyPatch>(pp))
{
// Assume zero gradient
nbrValue = cVals[own[faceI]];
nbrPoint = mesh_.faceCentres()[faceI];
}
else if (pp.coupled())
{
if (!refCast<const coupledPolyPatch>(pp).separated())
{
// Behave as internal face:
// other side value
nbrValue = cVals.boundaryField()[patchI][patchFaceI];
// other side cell centre
nbrPoint = mesh_.C().boundaryField()[patchI][patchFaceI];
}
else
{
// Do some interpolation for now
const scalarField& w = weights.boundaryField()[patchI];
nbrPoint = mesh_.faceCentres()[faceI];
nbrValue =
(1-w[patchFaceI])*cVals[own[faceI]]
+ w[patchFaceI]*cVals.boundaryField()[patchI][patchFaceI];
}
}
else
{
nbrValue = cVals.boundaryField()[patchI][patchFaceI];
nbrPoint = mesh_.faceCentres()[faceI];
}
}
}
// Determine per cell centre whether all the intersections get collapsed
// to a single point
void Foam::isoSurface::calcSnappedCc
......@@ -618,13 +660,14 @@ void Foam::isoSurface::calcSnappedPoint
forAll(pFaces, pFaceI)
{
// Create points for all intersections close to point
// (i.e. from pyramid edges)
label faceI = pFaces[pFaceI];
const face& f = mesh_.faces()[faceI];
label own = mesh_.faceOwner()[faceI];
// Create points for all intersections close to point
// (i.e. from pyramid edges)
// Get neighbour value
scalar nbrValue;
point nbrPoint;
getNeighbour
......@@ -843,6 +886,7 @@ Foam::triSurface Foam::isoSurface::stitchTriPoints
}
// Determine 'flat hole' situation (see RMT paper).
// Two unconnected triangles get connected because (some of) the edges
// separating them get collapsed. Below only checks for duplicate triangles,
......@@ -870,6 +914,9 @@ Foam::triSurface Foam::isoSurface::stitchTriPoints
if (nbrTriI > triI && (tris[nbrTriI] == tri))
{
//Pout<< "Duplicate : " << triI << " verts:" << tri
// << " to " << nbrTriI << " verts:" << tris[nbrTriI]
// << endl;
dupTriI = nbrTriI;
break;
}
......@@ -1526,11 +1573,17 @@ Foam::isoSurface::isoSurface
{
if (debug)
{
Pout<< "isoSurface :" << nl
<< " isoField : " << cVals.name() << nl
<< " isoValue : " << iso << nl
<< " regularise : " << regularise << nl
<< " mergeTol : " << mergeTol << nl
Pout<< "isoSurface:" << nl
<< " isoField : " << cVals.name() << nl
<< " cell min/max : "
<< min(cVals.internalField()) << " / "
<< max(cVals.internalField()) << nl
<< " point min/max : "
<< min(pVals) << " / "
<< max(pVals) << nl
<< " isoValue : " << iso << nl
<< " regularise : " << regularise << nl
<< " mergeTol : " << mergeTol << nl
<< endl;
}
......@@ -1556,15 +1609,7 @@ Foam::isoSurface::isoSurface
}
}
// Determine if any cut through face/cell
calcCutTypes(cVals, pVals);
// Determine if point is on boundary. Points on boundaries are never
// snapped. Coupled boundaries are handled explicitly so not marked here.
PackedBoolList isBoundaryPoint(mesh_.nPoints());
// Pre-calculate patch-per-face to avoid whichPatch call.
labelList boundaryRegion(mesh_.nFaces()-mesh_.nInternalFaces());
forAll(patches, patchI)
......@@ -1578,29 +1623,11 @@ Foam::isoSurface::isoSurface
boundaryRegion[faceI-mesh_.nInternalFaces()] = patchI;
faceI++;
}
}
// Mark all boundary points that are not physically coupled (so anything
// but collocated coupled patches)
if
(
!pp.coupled()
|| refCast<const coupledPolyPatch>(pp).separated()
)
{
label faceI = pp.start();
forAll(pp, i)
{
const face& f = mesh_.faces()[faceI];
forAll(f, fp)
{
isBoundaryPoint.set(f[fp], 1);
}
faceI++;
}
}
}
// Determine if any cut through face/cell
calcCutTypes(boundaryRegion, cVals, pVals);
DynamicList<point> snappedPoints(nCutCells_);
......@@ -1635,6 +1662,39 @@ Foam::isoSurface::isoSurface
label nCellSnaps = snappedPoints.size();
// Determine if point is on boundary. Points on boundaries are never
// snapped. Coupled boundaries are handled explicitly so not marked here.
PackedBoolList isBoundaryPoint(mesh_.nPoints());
forAll(patches, patchI)
{
// Mark all boundary points that are not physically coupled (so anything
// but collocated coupled patches)
const polyPatch& pp = patches[patchI];
if
(
!pp.coupled()
|| refCast<const coupledPolyPatch>(pp).separated()
)
{
label faceI = pp.start();
forAll(pp, i)
{
const face& f = mesh_.faces()[faceI];
forAll(f, fp)
{
isBoundaryPoint.set(f[fp], 1);
}
faceI++;
}
}
}
// Per point -1 or a point inside snappedPoints.
labelList snappedPoint;
if (regularise)
......@@ -1727,7 +1787,8 @@ Foam::isoSurface::isoSurface
if (debug)
{
Pout<< "isoSurface : generated " << triMeshCells.size()
<< " unmerged triangles." << endl;
<< " unmerged triangles from " << triPoints.size()
<< " unmerged points." << endl;
}
......
src/sampling/sampledSurface/isoSurface/isoSurface.H
View file @ c8944ce2
......@@ -135,9 +135,20 @@ class isoSurface
const bool neiLower
) const;
void getNeighbour
(
const labelList& boundaryRegion,
const volScalarField& cVals,
const label cellI,
const label faceI,
scalar& nbrValue,
point& nbrPoint
) const;
//- Set faceCutType,cellCutType.
void calcCutTypes
(
const labelList& boundaryRegion,
const volScalarField& cVals,
const scalarField& pVals
);
......@@ -156,16 +167,6 @@ class isoSurface
DynamicList<labelledTri, 64>& localTris
);
void getNeighbour
(
const labelList& boundaryRegion,
const volScalarField& cVals,
const label cellI,
const label faceI,
scalar& nbrValue,
point& nbrPoint
) const;
//- Determine per cc whether all near cuts can be snapped to single
// point.
void calcSnappedCc
......@@ -193,37 +194,39 @@ class isoSurface
template<class Type>
Type generatePoint
(
const DynamicList<Type>& snappedPoints,
const scalar s0,
const Type& p0,
const label p0Index,
const bool hasSnap0,
const Type& snapP0,
const scalar s1,
const Type& p1,
const label p1Index
const bool hasSnap1,
const Type& snapP1
) const;
template<class Type>
void generateTriPoints
(
const DynamicList<Type>& snapped,
const scalar s0,
const Type& p0,
const label p0Index,
const bool hasSnap0,
const Type& snapP0,
const scalar s1,
const Type& p1,
const label p1Index,
const bool hasSnap1,
const Type& snapP1,
const scalar s2,
const Type& p2,
const label p2Index,
const bool hasSnap2,
const Type& snapP2,
const scalar s3,
const Type& p3,
const label p3Index,
const bool hasSnap3,
const Type& snapP3,
DynamicList<Type>& points
) const;
......@@ -244,7 +247,8 @@ class isoSurface
const scalar neiVal,
const Type& neiPt,
const label neiSnap,
const bool hasNeiSnap,
const Type& neiSnapPt,
DynamicList<Type>& triPoints,
DynamicList<label>& triMeshCells
......
src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
View file @ c8944ce2
......@@ -33,15 +33,15 @@ License
template<class Type>
Type Foam::isoSurface::generatePoint
(
const DynamicList<Type>& snappedPoints,
const scalar s0,
const Type& p0,
const label p0Index,
const bool hasSnap0,
const Type& snapP0,
const scalar s1,
const Type& p1,
const label p1Index
const bool hasSnap1,
const Type& snapP1
) const
{
scalar d = s1-s0;
......@@ -50,13 +50,13 @@ Type Foam::isoSurface::generatePoint
{
scalar s = (iso_-s0)/d;
if (s >= 0.5 && s <= 1 && p1Index != -1)
if (hasSnap1 && s >= 0.5 && s <= 1)
{
return snappedPoints[p1Index];
return snapP1;
}
else if (s >= 0.0 && s <= 0.5 && p0Index != -1)
else if (hasSnap0 && s >= 0.0 && s <= 0.5)
{
return snappedPoints[p0Index];
return snapP0;
}
else
{
......@@ -75,23 +75,25 @@ Type Foam::isoSurface::generatePoint
template<class Type>
void Foam::isoSurface::generateTriPoints
(
const DynamicList<Type>& snapped,
const scalar s0,
const Type& p0,
const label p0Index,
const bool hasSnap0,
const Type& snapP0,
const scalar s1,
const Type& p1,
const label p1Index,
const bool hasSnap1,
const Type& snapP1,
const scalar s2,
const Type& p2,
const label p2Index,
const bool hasSnap2,
const Type& snapP2,
const scalar s3,
const Type& p3,
const label p3Index,
const bool hasSnap3,
const Type& snapP3,
DynamicList<Type>& points
) const
......@@ -123,29 +125,55 @@ void Foam::isoSurface::generateTriPoints
case 0x0E:
case 0x01:
points.append(generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index));
points.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
points.append
(
generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1)
);
points.append
(
generatePoint(s0,p0,hasSnap0,snapP0,s2,p2,hasSnap2,snapP2)
);
points.append
(
generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3)
);
break;
case 0x0D:
case 0x02:
points.append(generatePoint(snapped,s1,p1,p1Index,s0,p0,p0Index));
points.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
points.append
(
generatePoint(s1,p1,hasSnap1,snapP1,s0,p0,hasSnap0,snapP0)
);
points.append
(
generatePoint(s1,p1,hasSnap1,snapP1,s3,p3,hasSnap3,snapP3)
);
points.append
(
generatePoint(s1,p1,hasSnap1,snapP1,s2,p2,hasSnap2,snapP2)
);
break;
case 0x0C:
case 0x03:
{
Type tp1 = generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index);
Type tp2 = generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index);
Type tp1 =
generatePoint(s0,p0,hasSnap0,snapP0,s2,p2,hasSnap2,snapP2);
Type tp2 =
generatePoint(s1,p1,hasSnap1,snapP1,s3,p3,hasSnap3,snapP3);
points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
points.append
(
generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3)
);
points.append(tp1);
points.append(tp2);
points.append(tp2);
points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));