Commit 21a6de1d authored by franjo's avatar franjo

Added functionality for calculating intersections between the volume

mesh and the surface mesh


git-svn-id: https://pl5.projectlocker.com/igui/meshGeneration/svn@4 fdcce57e-7e00-11e2-b579-49867b4cea03
parent 9492038b
......@@ -22,6 +22,9 @@ renameBoundaryPatches = utilities/surfaceTools/renameBoundaryPatches
hexHelpers = utilities/hexahedra/hexHelpers
intersectionTools = utilities/intersectionTools
findCellsIntersectingSurface = $(intersectionTools)/findCellsIntersectingSurface
meshOptimizer = utilities/smoothers/geometry/meshOptimizer
tetMeshOptimisation = $(meshOptimizer)/tetMeshOptimisation
simplexSmoother = $(tetMeshOptimisation)/advancedSmoothers/simplexSmoother
......@@ -182,6 +185,8 @@ $(triSurf)/triSurf.C
$(hexHelpers)/hexHelpers.C
$(findCellsIntersectingSurface)/findCellsIntersectingSurface.C
$(tetrahedra)/tessellationElement.C
$(tetrahedra)/tetTessellation.C
$(tetrahedra)/tetTessellationFunctions.C
......
......@@ -39,6 +39,9 @@ SourceFiles
#include "plane.H"
#include "face.H"
#include "triangle.H"
#include "tetrahedron.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
class boolList;
......@@ -48,7 +51,7 @@ class triSurface;
namespace Foam
{
/*---------------------------------------------------------------------------*\
Namespace help functions Declaration
\*---------------------------------------------------------------------------*/
......@@ -79,7 +82,7 @@ namespace help
const edge& e,
const pointField& ep
);
//- check if the point p belongs to the plane
inline bool vertexInPlane(const point& p, const plane& pl);
......@@ -90,7 +93,7 @@ namespace help
const point& edgePoint1,
const point& p
);
//- find the vertex on the edge nearest to the point p
inline point nearestPointOnTheEdgeExact
(
......@@ -98,7 +101,7 @@ namespace help
const point& edgePoint1,
const point& p
);
//- find and return the distance between the edge and the point p
inline scalar distanceOfPointFromTheEdge
(
......@@ -106,7 +109,7 @@ namespace help
const point& edgePoint1,
const point& p
);
//- find the nearest points on the edge and the line
inline bool nearestEdgePointToTheLine
(
......@@ -117,7 +120,7 @@ namespace help
point& nearestOnEdge,
point& nearestOnLine
);
//- check if the edge intersects the plane
inline bool planeIntersectsEdge
(
......@@ -126,7 +129,23 @@ namespace help
const plane& pl,
point& intersection
);
//- check if a vertex lies inside the tetrahedron
inline bool pointInTetrahedron
(
const point& p,
const tetrahedron<point, point>& tet
);
//- check if a line intersects the triangle, and return the intersection
inline bool triLineIntersection
(
const triangle<point, point>& tria,
const point& lineStart,
const point& lineEnd,
point& intersection
);
//- check if a line intersects the triangle and return the intersection
inline bool triLineIntersection
(
......@@ -136,7 +155,7 @@ namespace help
const point&,
point&
);
//- check if the line intersects the bounding box
inline bool boundBoxLineIntersection
(
......@@ -144,7 +163,7 @@ namespace help
const point&,
const boundBox&
);
//- check if the surface triangle and the face intersect
inline bool doFaceAndTriangleIntersect
(
......@@ -153,7 +172,7 @@ namespace help
const face& f,
const pointField& facePoints
);
//- find the nearest vertex on the surface triangle to the given point
inline point nearestPointOnTheTriangle
(
......@@ -161,7 +180,7 @@ namespace help
const triSurface&,
const point&
);
//- check if the point is inside or outside the face
inline bool pointInsideFace
(
......@@ -170,7 +189,7 @@ namespace help
const vector& n,
const pointField& fp
);
//- check if the point is inside or outside the face
inline bool pointInsideFace
(
......@@ -178,10 +197,10 @@ namespace help
const face& f,
const pointField& fp
);
//- check if the vertex is on the positive side of the face plane
inline bool isVertexVisible(const point& p, const plane& pl);
//- find number of face groups within a given range
inline label numberOfFaceGroups
(
......@@ -190,7 +209,7 @@ namespace help
const scalar range,
const triSurface& surface
);
//- find the number of edge groups within the given range
inline label numberOfEdgeGroups
(
......
......@@ -233,23 +233,71 @@ inline bool planeIntersectsEdge
return false;
}
inline bool pointInTetrahedron
(
const point& p,
const tetrahedron<point, point>& tet
)
{
const vector v0 = tet.a() - tet.d();
const vector v1 = tet.b() - tet.d();
const vector v2 = tet.c() - tet.d();
const vector sp = p - tet.d();
matrix3D mat;
FixedList<scalar, 3> source;
for(label i=0;i<3;++i)
{
mat[i][0] = v0[i];
mat[i][1] = v1[i];
mat[i][2] = v2[i];
source[i] = sp[i];
}
//- check the determinant of the transformation
const scalar det = mat.determinant();
if( mag(det) < VSMALL )
return false;
//- get the coordinates of the point in the barycentric corrdinate system
const scalar u0 = mat.solveFirst(source);
if( (u0 < -SMALL) || (u0 > (1.0+SMALL)) )
return false;
const scalar u1 = mat.solveSecond(source);
if( (u1 < -SMALL) || ((u0+u1) > (1.0+SMALL)) )
return false;
const scalar u2 = mat.solveThird(source);
if( (u2 < -SMALL) || (u2 > (1.0+SMALL)) )
return false;
const scalar u3 = 1.0 - u0 - u1 - u2;
if( (u3 < -SMALL) || (u3 > (1.0+SMALL)) )
return false;
return true;
}
inline bool triLineIntersection
(
const triSurface& surface,
const label tI,
const point& s,
const point& e,
const triangle<point, point>& tria,
const point& lineStart,
const point& lineEnd,
point& intersection
)
{
const pointField& points = surface.points();
const labelledTri& ltri = surface.localFaces()[tI];
const point& p0 = points[ltri[0]];
const vector v(s - e);
const vector v0 = points[ltri[1]] - p0;
const vector v1 = points[ltri[2]] - p0;
const vector sp = s - p0;
const point& p0 = tria.a();
const vector v(lineStart - lineEnd);
const vector v0 = tria.b() - p0;
const vector v1 = tria.c() - p0;
const vector sp = lineStart - p0;
matrix3D mat;
FixedList<scalar, 3> source;
......@@ -281,10 +329,32 @@ inline bool triLineIntersection
if( (u1 < -SMALL) || ((u0+u1) > (1.0+SMALL)) )
return false;
intersection = s - t * v;
intersection = lineStart - t * v;
return true;
}
inline bool triLineIntersection
(
const triSurface& surface,
const label tI,
const point& s,
const point& e,
point& intersection
)
{
const pointField& pts = surface.localPoints();
const labelledTri& tri = surface[tI];
triangle<point, point> tria
(
pts[tri[0]],
pts[tri[1]],
pts[tri[2]]
);
return triLineIntersection(tria, s, e, intersection);
}
inline bool boundBoxLineIntersection
(
const point& s,
......@@ -361,16 +431,28 @@ inline bool doFaceAndTriangleIntersect
forAll(triEdges, eI)
{
const edge& e = triEdges[eI];
help::planeIntersectsEdge
(
triPoints[e.start()],
triPoints[e.end()],
pl,
intersection
);
if( help::pointInsideFace(intersection, f, n, facePoints) )
return true;
forAll(f, pI)
{
triangle<point, point> tria
(
facePoints[f[pI]],
facePoints[f.nextLabel(pI)],
centre
);
const bool currIntersection =
help::triLineIntersection
(
tria,
triPoints[e.start()],
triPoints[e.end()],
intersection
);
if( currIntersection )
return true;
}
}
//- check if any face edges intersect the triangle
......
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2005-2007 Franjo Juretic
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
Description
\*---------------------------------------------------------------------------*/
#include "findCellsIntersectingSurface.H"
#include "polyMeshGen.H"
#include "polyMeshGenAddressing.H"
#include "triSurf.H"
#include "boundBox.H"
#include "helperFunctions.H"
#include "meshOctree.H"
#include "meshOctreeCreator.H"
#include "HashSet.H"
#include <omp.h>
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void findCellsIntersectingSurface::generateOctree(const triSurf& surf)
{
octreePtr_ = new meshOctree(surf);
meshOctreeCreator(*octreePtr_).createOctreeWithRefinedBoundary(15, 15);
}
void findCellsIntersectingSurface::findIntersectedCells()
{
const pointFieldPMG& points = mesh_.points();
const faceListPMG& faces = mesh_.faces();
const cellListPMG& cells = mesh_.cells();
const labelList& owner = mesh_.owner();
const vectorField& faceCentres = mesh_.addressingData().faceCentres();
const vectorField& cellCentres = mesh_.addressingData().cellCentres();
meshOctreeModifier octreeModifier(*octreePtr_);
const LongList<meshOctreeCube*>& leaves = octreeModifier.leavesAccess();
intersectedCells_.setSize(cells.size());
const triSurf& surf = octreePtr_->surface();
# pragma omp parallel for schedule(dynamic, 40)
forAll(cells, cellI)
{
bool intersected(false);
const cell& c = cells[cellI];
//- find the bounding box of the cell
boundBox bb(cellCentres[cellI], cellCentres[cellI]);
forAll(c, fI)
{
const face& f = faces[c[fI]];
forAll(f, pI)
{
bb.min() = Foam::min(bb.min(), points[f[pI]]);
bb.max() = Foam::max(bb.max(), points[f[pI]]);
}
}
//- find surface triangles within the bounding box
DynList<label> leavesInBox;
octreePtr_->findLeavesContainedInBox(bb, leavesInBox);
labelHashSet triangles;
forAll(leavesInBox, boxI)
{
const meshOctreeCube& oc = *leaves[leavesInBox[boxI]];
if( oc.hasContainedElements() )
{
const meshOctreeSlot& slot = *oc.slotPtr();
const label ceI = oc.containedElements();
forAllRow(slot.containedTriangles_, ceI, tI)
triangles.insert(slot.containedTriangles_(ceI, tI));
}
}
//- check if any triangle in the surface mesh
//- intersects any of the cell's faces
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;
break;
}
}
if( intersected )
break;
}
//- check if any of the surface vertices is contained within the cell
if( !intersected )
{
labelHashSet nodes;
forAllConstIter(labelHashSet, triangles, tIter)
{
const labelledTri& tri = surf[tIter.key()];
for(label i=0;i<3;++i)
nodes.insert(tri[i]);
}
forAll(c, fI)
{
const face& f = faces[c[fI]];
if( owner[c[fI]] == cellI )
{
forAll(f, pI)
{
tetrahedron<point, point> tet
(
points[f[pI]],
points[f.prevLabel(pI)],
faceCentres[c[fI]],
cellCentres[cellI]
);
forAllConstIter(labelHashSet, nodes, nIter)
{
const point& p = surf.points()[nIter.key()];
if( help::pointInTetrahedron(p, tet) )
{
intersected = true;
break;
}
}
if( intersected )
break;
}
}
else
{
forAll(f, pI)
{
tetrahedron<point, point> tet
(
points[f[pI]],
points[f.nextLabel(pI)],
faceCentres[c[fI]],
cellCentres[cellI]
);
forAllConstIter(labelHashSet, nodes, nIter)
{
const point& p = surf.points()[nIter.key()];
if( help::pointInTetrahedron(p, tet) )
{
intersected = true;
break;
}
}
if( intersected )
break;
}
}
if( intersected )
break;
}
}
intersectedCells_[cellI] = intersected;
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
findCellsIntersectingSurface::findCellsIntersectingSurface
(
polyMeshGen& mesh,
const meshOctree& octree
)
:
mesh_(mesh),
octreePtr_(const_cast<meshOctree*>(&octree)),
octreeGenerated_(false)
{
findIntersectedCells();
}
findCellsIntersectingSurface::findCellsIntersectingSurface
(
polyMeshGen& mesh,
const triSurf& surface
)
:
mesh_(mesh),
octreePtr_(NULL),
octreeGenerated_(true)
{
generateOctree(surface);
findIntersectedCells();
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
findCellsIntersectingSurface::~findCellsIntersectingSurface()
{
if( octreeGenerated_ )
deleteDemandDrivenData(octreePtr_);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
const boolList& findCellsIntersectingSurface::intersectedCells() const
{
return intersectedCells_;
}
void findCellsIntersectingSurface::addIntersectedCellsToSubset
(
const word subsetName
)
{
const label setId = mesh_.addCellSubset(subsetName);
forAll(intersectedCells_, cellI)
if( intersectedCells_[cellI] )
mesh_.addCellToSubset(setId, cellI);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //
/*---------------------------------------------------------------------------*\
========= |