Skip to content
Permalink
Browse files

meshSearch: Prevent hang in calculation of line-boundary intersections

This fix changes how the intersections loop ignores previously
intersected faces. It now marks them by their index so that subsequent
iterations ignore them.

Before this change, after an intersection was found the start point was
advanced by a small amount to move the past the intersection. The
problem with this was if multiple boundary faces or the end point were
in close proximity to the intersection then the move forward might span
them. This could lead to intersections being missed or counted multiple
times, in some cases indefinitely.

Based on a patch contributed by Mattijs Janssens
Resolves bug report https://bugs.openfoam.org/view.php?id=1147
  • Loading branch information...
Will Bainbridge
Will Bainbridge committed Sep 28, 2018
1 parent 80df542 commit 39ad49bc344b4b4187a2cfac40adf26a87fa1485
@@ -38,6 +38,85 @@ namespace Foam
defineTypeNameAndDebug(meshSearch, 0);

scalar meshSearch::tol_ = 1e-3;

// Intersection operation that checks previous successful hits so that they
// are not duplicated
class findUniqueIntersectOp
:
public treeDataFace::findIntersectOp
{
public:

const indexedOctree<treeDataFace>& tree_;

const List<pointIndexHit>& hits_;

public:

//- Construct from components
findUniqueIntersectOp
(
const indexedOctree<treeDataFace>& tree,
const List<pointIndexHit>& hits
)
:
treeDataFace::findIntersectOp(tree),
tree_(tree),
hits_(hits)
{}

//- Calculate intersection of triangle with ray. Sets result
// accordingly
bool operator()
(
const label index,
const point& start,
const point& end,
point& intersectionPoint
) const
{
const primitiveMesh& mesh = tree_.shapes().mesh();

// Check whether this hit has already happened. If the new face
// index is the same as an existing hit then return no new hit. If
// the new face shares a point with an existing hit face and the
// line passes through both faces in the same direction, then this
// is also assumed to be a duplicate hit.
const label newFacei = tree_.shapes().faceLabels()[index];
const face& newFace = mesh.faces()[newFacei];
const scalar newDot = mesh.faceAreas()[newFacei] & (end - start);
forAll(hits_, hiti)
{
const label oldFacei = hits_[hiti].index();
const face& oldFace = mesh.faces()[oldFacei];
const scalar oldDot =
mesh.faceAreas()[oldFacei] & (end - start);

if
(
hits_[hiti].index() == newFacei
|| (
newDot*oldDot > 0
&& (labelHashSet(newFace) & labelHashSet(oldFace)).size()
)
)
{
return false;
}
}

const bool hit =
treeDataFace::findIntersectOp::operator()
(
index,
start,
end,
intersectionPoint
);

return hit;
}
};
}


@@ -466,25 +545,6 @@ Foam::label Foam::meshSearch::findNearestBoundaryFaceWalk
}


Foam::vector Foam::meshSearch::offset
(
const point& bPoint,
const label bFacei,
const vector& dir
) const
{
// Get the neighbouring cell
label ownerCelli = mesh_.faceOwner()[bFacei];

const point& c = mesh_.cellCentres()[ownerCelli];

// Typical dimension: distance from point on face to cell centre
scalar typDim = mag(c - bPoint);

return tol_*typDim*dir;
}


// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //

Foam::meshSearch::meshSearch
@@ -804,40 +864,21 @@ Foam::List<Foam::pointIndexHit> Foam::meshSearch::intersections
) const
{
DynamicList<pointIndexHit> hits;
pointIndexHit curHit;

vector edgeVec = pEnd - pStart;
edgeVec /= mag(edgeVec);

point pt = pStart;
findUniqueIntersectOp iop(boundaryTree(), hits);

pointIndexHit bHit;
do
while (true)
{
bHit = intersection(pt, pEnd);

if (bHit.hit())
{
hits.append(bHit);
// Get the next hit, or quit
curHit = boundaryTree().findLine(pStart, pEnd, iop);
if (!curHit.hit()) break;

const vector& area = mesh_.faceAreas()[bHit.index()];

scalar typDim = Foam::sqrt(mag(area));

if ((mag(bHit.hitPoint() - pEnd)/typDim) < small)
{
break;
}

// Restart from hitPoint shifted a little bit in the direction
// of the destination

pt =
bHit.hitPoint()
+ offset(bHit.hitPoint(), bHit.index(), edgeVec);
}

} while (bHit.hit());
// Change index into octreeData into face label
curHit.setIndex(boundaryTree().shapes().faceLabels()[curHit.index()]);

hits.append(curHit);
}

hits.shrink();

@@ -135,15 +135,6 @@ class meshSearch
const label seedFacei
) const;

//- Calculate offset vector in direction dir with as length a
// fraction of the cell size (of the cell straddling boundary face)
vector offset
(
const point& bPoint,
const label bFacei,
const vector& dir
) const;


//- Disallow default bitwise copy construct
meshSearch(const meshSearch&);
@@ -56,65 +56,38 @@ void Foam::sampledSets::lineFace::calcSamples
DynamicList<scalar>& samplingCurveDist
)
{
// Test for an inward pointing first hit. If this exists, then the start
// point is outside, and the sampling should begin from the hit.
label bHitI = -1;
// Look for a starting location within the mesh
vector startPt = start;
label startFaceI = -1, startCellI = searchEngine.findCell(start);

// Get all the hits of the line with the boundaries
const List<pointIndexHit> bHits =
searchEngine.intersections(start, end);
if (bHits.size())
{
const label bFaceI = bHits[0].index();
const vector bNormal = normalised(mesh.faceAreas()[bFaceI]);
const scalar bDot = bNormal & (end - start);
if (bDot < 0)
{
bHitI = 0;
}
}

// If an initial inward pointing hit was not found then initialise the
// within the cell containing the start point
vector startPt = vector::max;
label startFaceI = -1, startCellI = -1;
if (bHitI == -1)
{
startPt = start;
startFaceI = -1;
startCellI = searchEngine.findCell(start);
}

// If neither hits nor a containing cell for the start point were found then
// the line is completely outside the mesh. Return without generating any
// sampling points.
if (bHits.size() == 0 && startCellI == -1)
{
return;
}

// If nothing is set by now then something has gone wrong
if (bHitI == -1 && startCellI == -1)
// Loop over the hits, starting new segments each time
label bHitI = startCellI == -1 ? 0 : -1;
label sampleSegmentI = 0;
for (; bHitI < bHits.size(); ++ bHitI)
{
FatalErrorInFunction
<< "The initial location for the line " << start << " to " << end
<< "could not be found" << exit(FatalError);
}

// Loop over the hits, starting a new segment at every second hit
for
(
label sampleSegmentI = 0;
bHitI < bHits.size();
bHitI += 2, sampleSegmentI += 1
)
{
// Set the start point and topology, unless starting within a cell
if (bHitI != -1)
// If a boundary start point, then initialise to the current hit
if (bHitI >= 0)
{
startPt = bHits[bHitI].hitPoint();
startFaceI = bHits[bHitI].index();
startCellI = mesh.faceOwner()[startFaceI];
}

// If the hit points outward, move on to the next one
if (startFaceI != -1)
{
const vector bNormal = normalised(mesh.faceAreas()[startFaceI]);
const scalar bDot = bNormal & (end - start);
if (bDot > 0)
{
continue;
}
}

// Create a particle. If we are starting on a boundary face, track
// backwards into it so that the particle has the correct topology.
passiveParticle sampleParticle(mesh, startPt, startCellI);
@@ -159,6 +132,8 @@ void Foam::sampledSets::lineFace::calcSamples
break;
}
}

++ sampleSegmentI;
}
}

0 comments on commit 39ad49b

Please sign in to comment.
You can’t perform that action at this time.