Commit 55165eaf authored by franjo_j@hotmail.com's avatar franjo_j@hotmail.com
Browse files

Intermediate version of removing domains


git-svn-id: https://pl5.projectlocker.com/igui/meshGeneration/svn@25 fdcce57e-7e00-11e2-b579-49867b4cea03
parent bbe6d835
......@@ -286,6 +286,185 @@ inline bool pointInTetrahedron
return true;
}
inline bool nearestEdgePointToTheLine
(
const point& edgePoint0,
const point& edgePoint1,
const point& lp0,
const point& lp1,
point& nearestOnEdge,
point& nearestOnLine
)
{
const vector v = lp1 - lp0;
const vector d = lp0 - edgePoint0;
const vector e = edgePoint1 - edgePoint0;
const scalar vMag = mag(v);
if( vMag < VSMALL )
return false;
const scalar eMag = mag(e);
if( eMag < VSMALL )
{
nearestOnEdge = edgePoint0;
nearestOnLine = nearestPointOnTheEdge(lp0, lp1, nearestOnEdge);
return true;
}
if( mag((v/vMag) & (e/eMag)) > (1.0 - SMALL) )
return false;
tensor mat(tensor::zero);
mat.xx() = (v&v);
mat.xy() = mat.yx() = -1.0 * (v&e);
mat.yy() = (e&e);
mat.zz() = SMALL;
vector source(vector::zero);
source[0] = -1.0 * (d&v);
source[1] = (d&e);
const vector sol = (inv(mat) & source);
nearestOnLine = lp0 + v * sol[0];
if( sol[1] > 1.0 )
{
nearestOnEdge = edgePoint1;
}
else if( sol[1] < 0.0 )
{
nearestOnEdge = edgePoint0;
}
else
{
nearestOnEdge = edgePoint0 + e * sol[1];
}
return true;
}
inline point nearestPointOnTheTriangle
(
const label tI,
const triSurf& surface,
const point& p
)
{
const labelledTri& ltri = surface[tI];
const pointField& points = surface.points();
const vector v0 = points[ltri[1]] - points[ltri[0]];
const vector v1 = points[ltri[2]] - points[ltri[0]];
const vector v2 = p - points[ltri[0]];
const scalar dot00 = (v0 & v0);
const scalar dot01 = (v0 & v1);
const scalar dot02 = (v0 & v2);
const scalar dot11 = (v1 & v1);
const scalar dot12 = (v1 & v2);
// Compute barycentric coordinates
const scalar det = dot00 * dot11 - dot01 * dot01;
if( mag(det) < VSMALL )
{
point nearest(p);
scalar dist(VGREAT);
const edgeList edges = ltri.edges();
forAll(edges, eI)
{
const edge& e = edges[eI];
const point np =
nearestPointOnTheEdge
(
points[e.start()],
points[e.end()],
p
);
if( magSqr(p - np) < dist )
{
nearest = np;
dist = magSqr(p - np);
}
}
return nearest;
}
const scalar u = (dot11 * dot02 - dot01 * dot12) / det;
const scalar v = (dot00 * dot12 - dot01 * dot02) / det;
const point pProj = points[ltri[0]] + u * v0 + v * v1;
//- Check if point is in triangle
if( (u >= -SMALL) && (v >= -SMALL) && ((u + v) <= (1.0+SMALL)) )
{
return pProj;
}
else
{
if( u < -SMALL )
{
const edge e(ltri[0], ltri[2]);
const scalar ed = ((pProj - points[e.start()]) & v1) / magSqr(v1);
if( ed > 1.0 )
{
return points[e.end()];
}
else if( ed < 0.0 )
{
return points[e.start()];
}
else
{
return (points[e.start()] + v1 * ed);
}
}
else if( v < -SMALL )
{
const edge e(ltri[0], ltri[1]);
const scalar ed = ((pProj - points[e.start()]) & v0) / magSqr(v0);
if( ed > 1.0 )
{
return points[e.end()];
}
else if( ed < 0.0 )
{
return points[e.start()];
}
else
{
return (points[e.start()] + v0 * ed);
}
}
else
{
const edge e(ltri[2], ltri[1]);
const vector ev = e.vec(points);
const scalar ed = ((pProj - points[e.start()]) & ev) / magSqr(ev);
if( ed > 1.0 )
{
return points[e.end()];
}
else if( ed < 0.0 )
{
return points[e.start()];
}
else
{
return (points[e.start()] + ev * ed);
}
}
}
return pProj;
}
inline bool triLineIntersection
(
const triangle<point, point>& tria,
......@@ -312,7 +491,7 @@ inline bool triLineIntersection
const scalar det = mat.determinant();
if( mag(det) < VSMALL )
if( mag(det) < SMALL )
return false;
const scalar t = mat.solveThird(source);
......@@ -346,7 +525,38 @@ inline bool triLineIntersection
const pointField& pts = surface.points();
const labelledTri& tri = surface[tI];
triangle<point, point> tria
scalar tolSq = SMALL * tri.mag(pts);
point nearest = nearestPointOnTheTriangle(tI, surface, s);
if( magSqr(nearest - s) < tolSq )
{
intersection = nearest;
return true;
}
nearest = nearestPointOnTheTriangle(tI, surface, e);
if( magSqr(nearest - e) < tolSq )
{
intersection = nearest;
return true;
}
forAll(tri, eI)
{
const point& ep = pts[tri[eI]];
const point& ee = pts[tri[(eI+1)%3]];
point inter1, inter2;
nearestEdgePointToTheLine(s, e, ep, ee, inter1, inter2);
if( magSqr(inter1 - inter2) < tolSq )
{
intersection = inter1;
return true;
}
}
const triangle<point, point> tria
(
pts[tri[0]],
pts[tri[1]],
......@@ -421,21 +631,20 @@ inline bool doFaceAndTriangleIntersect
{
const pointField& triPoints = surface.points();
vector n = f.normal(facePoints);
n /= mag(n);
const point centre = f.centre(facePoints);
plane pl(centre, n);
point intersection(vector::zero);
point intersection;
//- check if any triangle edge intersects the face
const edgeList triEdges = surface[triI].edges();
forAll(triEdges, eI)
const labelledTri& tri = surface[triI];
forAll(tri, eI)
{
const edge& e = triEdges[eI];
const point& s = triPoints[tri[eI]];
const point& e = triPoints[tri[(eI+1)%3]];
forAll(f, pI)
{
triangle<point, point> tria
const triangle<point, point> tria
(
facePoints[f[pI]],
facePoints[f.nextLabel(pI)],
......@@ -446,8 +655,8 @@ inline bool doFaceAndTriangleIntersect
help::triLineIntersection
(
tria,
triPoints[e.start()],
triPoints[e.end()],
s,
e,
intersection
);
......@@ -457,148 +666,28 @@ inline bool doFaceAndTriangleIntersect
}
//- check if any face edges intersect the triangle
const edgeList fEdges = f.edges();
forAll(fEdges, feI)
forAll(f, pI)
{
const edge& e = fEdges[feI];
const point& s = facePoints[f[pI]];
const point& e = facePoints[f.nextLabel(pI)];
if(
const bool intersected =
help::triLineIntersection
(
surface,
triI,
facePoints[e.start()],
facePoints[e.end()],
s,
e,
intersection
)
)
);
if( intersected )
return true;
}
return false;
}
inline point nearestPointOnTheTriangle
(
const label tI,
const triSurf& surface,
const point& p
)
{
const labelledTri& ltri = surface[tI];
const pointField& points = surface.points();
const vector v0 = points[ltri[1]] - points[ltri[0]];
const vector v1 = points[ltri[2]] - points[ltri[0]];
const vector v2 = p - points[ltri[0]];
const scalar dot00 = (v0 & v0);
const scalar dot01 = (v0 & v1);
const scalar dot02 = (v0 & v2);
const scalar dot11 = (v1 & v1);
const scalar dot12 = (v1 & v2);
// Compute barycentric coordinates
const scalar det = dot00 * dot11 - dot01 * dot01;
if( mag(det) < VSMALL )
{
point nearest(p);
scalar dist(VGREAT);
const edgeList edges = ltri.edges();
forAll(edges, eI)
{
const edge& e = edges[eI];
const point np =
nearestPointOnTheEdge
(
points[e.start()],
points[e.end()],
p
);
if( magSqr(p - np) < dist )
{
nearest = np;
dist = magSqr(p - np);
}
}
return nearest;
}
const scalar u = (dot11 * dot02 - dot01 * dot12) / det;
const scalar v = (dot00 * dot12 - dot01 * dot02) / det;
const point pProj = points[ltri[0]] + u * v0 + v * v1;
//- Check if point is in triangle
if( (u >= -SMALL) && (v >= -SMALL) && ((u + v) <= (1.0+SMALL)) )
{
return pProj;
}
else
{
if( u < -SMALL )
{
const edge e(ltri[0], ltri[2]);
const scalar ed = ((pProj - points[e.start()]) & v1) / magSqr(v1);
if( ed > 1.0 )
{
return points[e.end()];
}
else if( ed < 0.0 )
{
return points[e.start()];
}
else
{
return (points[e.start()] + v1 * ed);
}
}
else if( v < -SMALL )
{
const edge e(ltri[0], ltri[1]);
const scalar ed = ((pProj - points[e.start()]) & v0) / magSqr(v0);
if( ed > 1.0 )
{
return points[e.end()];
}
else if( ed < 0.0 )
{
return points[e.start()];
}
else
{
return (points[e.start()] + v0 * ed);
}
}
else
{
const edge e(ltri[2], ltri[1]);
const vector ev = e.vec(points);
const scalar ed = ((pProj - points[e.start()]) & ev) / magSqr(ev);
if( ed > 1.0 )
{
return points[e.end()];
}
else if( ed < 0.0 )
{
return points[e.start()];
}
else
{
return (points[e.start()] + ev * ed);
}
}
}
return pProj;
}
inline bool pointInsideFace
(
const point& p,
......@@ -715,64 +804,6 @@ inline scalar distanceOfPointFromTheEdge
return mag(nearestPointOnTheEdge(edgePoint0, edgePoint1, p) - p);
}
inline bool nearestEdgePointToTheLine
(
const point& edgePoint0,
const point& edgePoint1,
const point& lp0,
const point& lp1,
point& nearestOnEdge,
point& nearestOnLine
)
{
const vector v = lp1 - lp0;
const vector d = lp0 - edgePoint0;
const vector e = edgePoint1 - edgePoint0;
const scalar vMag = mag(v);
if( vMag < VSMALL )
return false;
const scalar eMag = mag(e);
if( eMag < VSMALL )
{
nearestOnEdge = edgePoint0;
nearestOnLine = nearestPointOnTheEdge(lp0, lp1, nearestOnEdge);
return true;
}
if( mag((v/vMag) & (e/eMag)) > (1.0 - SMALL) )
return false;
tensor mat(tensor::zero);
mat.xx() = (v&v);
mat.xy() = mat.yx() = -1.0 * (v&e);
mat.yy() = (e&e);
mat.zz() = SMALL;
vector source(vector::zero);
source[0] = -1.0 * (d&v);
source[1] = (d&e);
const vector sol = (inv(mat) & source);
nearestOnLine = lp0 + v * sol[0];
if( sol[1] > 1.0 )
{
nearestOnEdge = edgePoint1;
}
else if( sol[1] < 0.0 )
{
nearestOnEdge = edgePoint0;
}
else
{
nearestOnEdge = edgePoint0 + e * sol[1];
}
return true;
}
inline label numberOfFaceGroups
(
const labelHashSet& containedElements,
......
......@@ -69,6 +69,8 @@ void findCellsIntersectingSurface::findIntersectedCells()
facetsIntersectingCell_.setSize(cells.size());
const triSurf& surf = octreePtr_->surface();
const VRWGraph& pointFacets = surf.pointFacets();
const pointField& sp = surf.points();
# ifdef USE_OMP
# pragma omp parallel for schedule(dynamic, 40)
......@@ -93,6 +95,10 @@ void findCellsIntersectingSurface::findIntersectedCells()
}
}
const vector spanCorr = 0.01 * bb.span();
bb.max() += spanCorr;
bb.min() -= spanCorr;
//- find surface triangles within the bounding box
DynList<label> leavesInBox;
octreePtr_->findLeavesContainedInBox(bb, leavesInBox);
......@@ -113,7 +119,7 @@ void findCellsIntersectingSurface::findIntersectedCells()
//- remove triangles which do not intersect the bounding box
labelHashSet reasonableCandidates;
const pointField& sp = surf.points();
forAllConstIter(labelHashSet, triangles, tIter)
{
const labelledTri& tri = surf[tIter.key()];
......@@ -126,8 +132,9 @@ void findCellsIntersectingSurface::findIntersectedCells()
obb.max() = Foam::max(obb.max(), v);
}
obb.min() -= vector(VSMALL, VSMALL, VSMALL);
obb.max() += vector(VSMALL, VSMALL, VSMALL);
const vector spanTriCorr = SMALL * obb.span();
obb.min() -= spanTriCorr;
obb.max() += spanTriCorr;
if( obb.overlaps(bb) )
reasonableCandidates.insert(tIter.key());
......@@ -135,49 +142,23 @@ void findCellsIntersectingSurface::findIntersectedCells()
triangles.transfer(reasonableCandidates);
//- check if any triangle in the surface mesh
//- intersects any of the cell's faces
labelHashSet facetsInCell;
forAllConstIter(labelHashSet, triangles, tIter)
{
forAll(c, fI)
{
const face& f = faces[c[fI]];
const bool intersect =
help::doFaceAndTriangleIntersect
(
surf,
tIter.key(),
f,
points
);
if( intersect )
{
intersected = true;
facetsInCell.insert(tIter.key());
break;
}
}
}
//- check if any of the surface vertices is contained within the cell
labelHashSet nodes;
labelHashSet nodes, facetsInCell;
forAllConstIter(labelHashSet, triangles, tIter)
{
if( facetsInCell.found(tIter.key()) )
continue;
const labelledTri& tri = surf[tIter.key()];
for(label i=0;i<3;++i)
forAll(tri, i)
nodes.insert(tri[i]);