zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Merge remote-tracking branch 'origin/master' into feature/procAgglom

Browse Source
This commit is contained in:
mattijs
2013-05-02 12:38:25 +01:00
parent 6be54d77a1 c7dbe491c6
commit c3737a206b
63 changed files with 2471 additions and 1730 deletions
Download Patch File Download Diff File Expand all files Collapse all files

7
src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/AMIMethod/AMIMethod.C
View File

@ -210,7 +210,12 @@ void Foam::AMIMethod<SourcePatch, TargetPatch>::resetTree()
(
new indexedOctree<treeType>
(
treeType(false, tgtPatch_),
treeType
(
false,
tgtPatch_,
indexedOctree<treeType>::perturbTol()
),
bb, // overall search domain
8, // maxLevel
10, // leaf size

1
src/meshTools/Make/files
View File

@ -160,6 +160,7 @@ $(intersectedSurface)/intersectedSurface.C
$(intersectedSurface)/edgeSurface.C
triSurface/triSurfaceSearch/triSurfaceSearch.C
triSurface/triSurfaceSearch/triSurfaceRegionSearch.C
triSurface/triangleFuncs/triangleFuncs.C
triSurface/surfaceFeatures/surfaceFeatures.C
triSurface/triSurfaceTools/triSurfaceTools.C

3
src/meshTools/cellClassification/cellClassification.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -34,6 +34,7 @@ License
#include "ListOps.H"
#include "meshTools.H"
#include "cpuTime.H"
#include "triSurface.H"
#include "globalMeshData.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

4
src/meshTools/cellClassification/cellClassification.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -107,9 +107,7 @@ namespace Foam
class triSurfaceSearch;
class meshSearch;
class polyMesh;
class polyMesh;
class primitiveMesh;
class triSurface;
/*---------------------------------------------------------------------------*\
Class cellClassification Declaration

61
src/meshTools/indexedOctree/treeDataEdge.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -96,6 +96,24 @@ Foam::treeDataEdge::treeDataEdge
}
Foam::treeDataEdge::findNearestOp::findNearestOp
(
const indexedOctree<treeDataEdge>& tree
)
:
tree_(tree)
{}
Foam::treeDataEdge::findIntersectOp::findIntersectOp
(
const indexedOctree<treeDataEdge>& tree
)
:
tree_(tree)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::pointField Foam::treeDataEdge::shapePoints() const
@ -114,13 +132,13 @@ Foam::pointField Foam::treeDataEdge::shapePoints() const
//- Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
// Only makes sense for closed surfaces.
Foam::label Foam::treeDataEdge::getVolumeType
Foam::volumeType Foam::treeDataEdge::getVolumeType
(
const indexedOctree<treeDataEdge>& oc,
const point& sample
) const
{
return indexedOctree<treeDataEdge>::UNKNOWN;
return volumeType::UNKNOWN;
}
@ -165,9 +183,7 @@ bool Foam::treeDataEdge::overlaps
}
// Calculate nearest point to sample. Updates (if any) nearestDistSqr, minIndex,
// nearestPoint.
void Foam::treeDataEdge::findNearest
void Foam::treeDataEdge::findNearestOp::operator()
(
const labelUList& indices,
const point& sample,
@ -177,13 +193,15 @@ void Foam::treeDataEdge::findNearest
point& nearestPoint
) const
{
const treeDataEdge& shape = tree_.shapes();
forAll(indices, i)
{
const label index = indices[i];
const edge& e = edges_[edgeLabels_[index]];
const edge& e = shape.edges()[shape.edgeLabels()[index]];
pointHit nearHit = e.line(points_).nearestDist(sample);
pointHit nearHit = e.line(shape.points()).nearestDist(sample);
scalar distSqr = sqr(nearHit.distance());
@ -197,9 +215,7 @@ void Foam::treeDataEdge::findNearest
}
//- Calculates nearest (to line) point in shape.
// Returns point and distance (squared)
void Foam::treeDataEdge::findNearest
void Foam::treeDataEdge::findNearestOp::operator()
(
const labelUList& indices,
const linePointRef& ln,
@ -210,6 +226,8 @@ void Foam::treeDataEdge::findNearest
point& nearestPoint
) const
{
const treeDataEdge& shape = tree_.shapes();
// Best so far
scalar nearestDistSqr = magSqr(linePoint - nearestPoint);
@ -217,13 +235,13 @@ void Foam::treeDataEdge::findNearest
{
const label index = indices[i];
const edge& e = edges_[edgeLabels_[index]];
const edge& e = shape.edges()[shape.edgeLabels()[index]];
// Note: could do bb test ? Worthwhile?
// Nearest point on line
point ePoint, lnPt;
scalar dist = e.line(points_).nearestDist(ln, ePoint, lnPt);
scalar dist = e.line(shape.points()).nearestDist(ln, ePoint, lnPt);
scalar distSqr = sqr(dist);
if (distSqr < nearestDistSqr)
@ -252,4 +270,21 @@ void Foam::treeDataEdge::findNearest
}
bool Foam::treeDataEdge::findIntersectOp::operator()
(
const label index,
const point& start,
const point& end,
point& result
) const
{
notImplemented
(
"treeDataEdge::intersects(const label, const point&,"
"const point&, point&)"
);
return false;
}
// ************************************************************************* //

111
src/meshTools/indexedOctree/treeDataEdge.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -37,6 +37,7 @@ SourceFiles
#include "treeBoundBoxList.H"
#include "linePointRef.H"
#include "volumeType.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -86,6 +87,58 @@ class treeDataEdge
public:
class findNearestOp
{
const indexedOctree<treeDataEdge>& tree_;
public:
findNearestOp(const indexedOctree<treeDataEdge>& tree);
void operator()
(
const labelUList& indices,
const point& sample,
scalar& nearestDistSqr,
label& minIndex,
point& nearestPoint
) const;
void operator()
(
const labelUList& indices,
const linePointRef& ln,
treeBoundBox& tightest,
label& minIndex,
point& linePoint,
point& nearestPoint
) const;
};
class findIntersectOp
{
const indexedOctree<treeDataEdge>& tree_;
public:
findIntersectOp(const indexedOctree<treeDataEdge>& tree);
//- Calculate intersection of triangle with ray. Sets result
// accordingly
bool operator()
(
const label index,
const point& start,
const point& end,
point& intersectionPoint
) const;
};
// Declare name of the class and its debug switch
ClassName("treeDataEdge");
@ -116,6 +169,16 @@ public:
// Access
const edgeList& edges() const
{
return edges_;
}
const pointField& points() const
{
return points_;
}
const labelList& edgeLabels() const
{
return edgeLabels_;
@ -135,7 +198,7 @@ public:
//- Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
// Only makes sense for closed surfaces.
label getVolumeType
volumeType getVolumeType
(
const indexedOctree<treeDataEdge>&,
const point&
@ -155,50 +218,6 @@ public:
const point& centre,
const scalar radiusSqr
) const;
//- Calculates nearest (to sample) point in shape.
// Returns actual point and distance (squared)
void findNearest
(
const labelUList& indices,
const point& sample,
scalar& nearestDistSqr,
label& nearestIndex,
point& nearestPoint
) const;
//- Calculates nearest (to line) point in shape.
// Returns point and distance (squared)
void findNearest
(
const labelUList& indices,
const linePointRef& ln,
treeBoundBox& tightest,
label& minIndex,
point& linePoint,
point& nearestPoint
) const;
//- Calculate intersection of shape with ray. Sets result
// accordingly
bool intersects
(
const label index,
const point& start,
const point& end,
point& result
) const
{
notImplemented
(
"treeDataEdge::intersects(const label, const point&,"
"const point&, point&)"
);
return false;
}
};

76
src/meshTools/indexedOctree/treeDataFace.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -145,6 +145,24 @@ Foam::treeDataFace::treeDataFace
}
Foam::treeDataFace::findNearestOp::findNearestOp
(
const indexedOctree<treeDataFace>& tree
)
:
tree_(tree)
{}
Foam::treeDataFace::findIntersectOp::findIntersectOp
(
const indexedOctree<treeDataFace>& tree
)
:
tree_(tree)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::pointField Foam::treeDataFace::shapePoints() const
@ -162,7 +180,7 @@ Foam::pointField Foam::treeDataFace::shapePoints() const
//- Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
// Only makes sense for closed surfaces.
Foam::label Foam::treeDataFace::getVolumeType
Foam::volumeType Foam::treeDataFace::getVolumeType
(
const indexedOctree<treeDataFace>& oc,
const point& sample
@ -415,7 +433,7 @@ Foam::label Foam::treeDataFace::getVolumeType
// - tolerances are wrong. (if e.g. face has zero area)
// - or (more likely) surface is not closed.
return indexedOctree<treeDataFace>::UNKNOWN;
return volumeType::UNKNOWN;
}
@ -477,9 +495,7 @@ bool Foam::treeDataFace::overlaps
}
// Calculate nearest point to sample. Updates (if any) nearestDistSqr, minIndex,
// nearestPoint.
void Foam::treeDataFace::findNearest
void Foam::treeDataFace::findNearestOp::operator()
(
const labelUList& indices,
const point& sample,
@ -489,13 +505,15 @@ void Foam::treeDataFace::findNearest
point& nearestPoint
) const
{
const treeDataFace& shape = tree_.shapes();
forAll(indices, i)
{
const label index = indices[i];
const face& f = mesh_.faces()[faceLabels_[index]];
const face& f = shape.mesh().faces()[shape.faceLabels()[index]];
pointHit nearHit = f.nearestPoint(sample, mesh_.points());
pointHit nearHit = f.nearestPoint(sample, shape.mesh().points());
scalar distSqr = sqr(nearHit.distance());
if (distSqr < nearestDistSqr)
@ -508,7 +526,33 @@ void Foam::treeDataFace::findNearest
}
bool Foam::treeDataFace::intersects
void Foam::treeDataFace::findNearestOp::operator()
(
const labelUList& indices,
const linePointRef& ln,
treeBoundBox& tightest,
label& minIndex,
point& linePoint,
point& nearestPoint
) const
{
notImplemented
(
"treeDataFace::findNearestOp::operator()"
"("
" const labelUList&,"
" const linePointRef&,"
" treeBoundBox&,"
" label&,"
" point&,"
" point&"
") const"
);
}
bool Foam::treeDataFace::findIntersectOp::operator()
(
const label index,
const point& start,
@ -516,10 +560,12 @@ bool Foam::treeDataFace::intersects
point& intersectionPoint
) const
{
const treeDataFace& shape = tree_.shapes();
// Do quick rejection test
if (cacheBb_)
if (shape.cacheBb_)
{
const treeBoundBox& faceBb = bbs_[index];
const treeBoundBox& faceBb = shape.bbs_[index];
if ((faceBb.posBits(start) & faceBb.posBits(end)) != 0)
{
@ -528,16 +574,16 @@ bool Foam::treeDataFace::intersects
}
}
const label faceI = faceLabels_[index];
const label faceI = shape.faceLabels_[index];
const vector dir(end - start);
pointHit inter = mesh_.faces()[faceI].intersection
pointHit inter = shape.mesh_.faces()[faceI].intersection
(
start,
dir,
mesh_.faceCentres()[faceI],
mesh_.points(),
shape.mesh_.faceCentres()[faceI],
shape.mesh_.points(),
intersection::HALF_RAY
);

103
src/meshTools/indexedOctree/treeDataFace.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -39,6 +39,8 @@ SourceFiles
#include "indexedOctree.H"
#include "treeBoundBoxList.H"
#include "PackedBoolList.H"
#include "primitiveMesh.H"
#include "volumeType.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -46,7 +48,7 @@ namespace Foam
{
// Forward declaration of classes
class primitiveMesh;
//class primitiveMesh;
//template<class Type> class indexedOctree;
class polyPatch;
@ -90,6 +92,58 @@ class treeDataFace
public:
class findNearestOp
{
const indexedOctree<treeDataFace>& tree_;
public:
findNearestOp(const indexedOctree<treeDataFace>& tree);
void operator()
(
const labelUList& indices,
const point& sample,
scalar& nearestDistSqr,
label& minIndex,
point& nearestPoint
) const;
void operator()
(
const labelUList& indices,
const linePointRef& ln,
treeBoundBox& tightest,
label& minIndex,
point& linePoint,
point& nearestPoint
) const;
};
class findIntersectOp
{
const indexedOctree<treeDataFace>& tree_;
public:
findIntersectOp(const indexedOctree<treeDataFace>& tree);
//- Calculate intersection of triangle with ray. Sets result
// accordingly
bool operator()
(
const label index,
const point& start,
const point& end,
point& intersectionPoint
) const;
};
// Declare name of the class and its debug switch
ClassName("treeDataFace");
@ -147,7 +201,7 @@ public:
//- Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
// Only makes sense for closed surfaces.
label getVolumeType
volumeType getVolumeType
(
const indexedOctree<treeDataFace>&,
const point&
@ -159,49 +213,6 @@ public:
const label index,
const treeBoundBox& sampleBb
) const;
//- Calculates nearest (to sample) point in shape.
// Returns actual point and distance (squared)
void findNearest
(
const labelUList& indices,
const point& sample,
scalar& nearestDistSqr,
label& nearestIndex,
point& nearestPoint
) const;
//- Calculates nearest (to line) point in shape.
// Returns point and distance (squared)
void findNearest
(
const labelUList& indices,
const linePointRef& ln,
treeBoundBox& tightest,
label& minIndex,
point& linePoint,
point& nearestPoint
) const
{
notImplemented
(
"treeDataFace::findNearest"
"(const labelUList&, const linePointRef&, ..)"
);
}
//- Calculate intersection of shape with ray. Sets result
// accordingly
bool intersects
(
const label index,
const point& start,
const point& end,
point& result
) const;
};

69
src/meshTools/indexedOctree/treeDataPoint.C
View File

@ -57,6 +57,24 @@ Foam::treeDataPoint::treeDataPoint
{}
Foam::treeDataPoint::findNearestOp::findNearestOp
(
const indexedOctree<treeDataPoint>& tree
)
:
tree_(tree)
{}
Foam::treeDataPoint::findIntersectOp::findIntersectOp
(
const indexedOctree<treeDataPoint>& tree
)
:
tree_(tree)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::pointField Foam::treeDataPoint::shapePoints() const
@ -74,13 +92,13 @@ Foam::pointField Foam::treeDataPoint::shapePoints() const
//- Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
// Only makes sense for closed surfaces.
Foam::label Foam::treeDataPoint::getVolumeType
Foam::volumeType Foam::treeDataPoint::getVolumeType
(
const indexedOctree<treeDataPoint>& oc,
const point& sample
) const
{
return indexedOctree<treeDataPoint>::UNKNOWN;
return volumeType::UNKNOWN;
}
@ -115,9 +133,7 @@ bool Foam::treeDataPoint::overlaps
}
// Calculate nearest point to sample. Updates (if any) nearestDistSqr, minIndex,
// nearestPoint.
void Foam::treeDataPoint::findNearest
void Foam::treeDataPoint::findNearestOp::operator()
(
const labelUList& indices,
const point& sample,
@ -127,12 +143,19 @@ void Foam::treeDataPoint::findNearest
point& nearestPoint
) const
{
const treeDataPoint& shape = tree_.shapes();
forAll(indices, i)
{
const label index = indices[i];
label pointI = (useSubset_ ? pointLabels_[index] : index);
label pointI =
(
shape.useSubset()
? shape.pointLabels()[index]
: index
);
const point& pt = points_[pointI];
const point& pt = shape.points()[pointI];
scalar distSqr = magSqr(pt - sample);
@ -146,9 +169,7 @@ void Foam::treeDataPoint::findNearest
}
//- Calculates nearest (to line) point in shape.
// Returns point and distance (squared)
void Foam::treeDataPoint::findNearest
void Foam::treeDataPoint::findNearestOp::operator()
(
const labelUList& indices,
const linePointRef& ln,
@ -159,6 +180,8 @@ void Foam::treeDataPoint::findNearest
point& nearestPoint
) const
{
const treeDataPoint& shape = tree_.shapes();
// Best so far
scalar nearestDistSqr = GREAT;
if (minIndex >= 0)
@ -169,9 +192,14 @@ void Foam::treeDataPoint::findNearest
forAll(indices, i)
{
const label index = indices[i];
label pointI = (useSubset_ ? pointLabels_[index] : index);
label pointI =
(
shape.useSubset()
? shape.pointLabels()[index]
: index
);
const point& shapePt = points_[pointI];
const point& shapePt = shape.points()[pointI];
if (tightest.contains(shapePt))
{
@ -206,4 +234,21 @@ void Foam::treeDataPoint::findNearest
}
bool Foam::treeDataPoint::findIntersectOp::operator()
(
const label index,
const point& start,
const point& end,
point& result
) const
{
notImplemented
(
"treeDataPoint::intersects(const label, const point&,"
"const point&, point&)"
);
return false;
}
// ************************************************************************* //

106
src/meshTools/indexedOctree/treeDataPoint.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -43,6 +43,7 @@ SourceFiles
#include "pointField.H"
#include "treeBoundBox.H"
#include "linePointRef.H"
#include "volumeType.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -69,6 +70,58 @@ class treeDataPoint
public:
class findNearestOp
{
const indexedOctree<treeDataPoint>& tree_;
public:
findNearestOp(const indexedOctree<treeDataPoint>& tree);
void operator()
(
const labelUList& indices,
const point& sample,
scalar& nearestDistSqr,
label& minIndex,
point& nearestPoint
) const;
void operator()
(
const labelUList& indices,
const linePointRef& ln,
treeBoundBox& tightest,
label& minIndex,
point& linePoint,
point& nearestPoint
) const;
};
class findIntersectOp
{
const indexedOctree<treeDataPoint>& tree_;
public:
findIntersectOp(const indexedOctree<treeDataPoint>& tree);
//- Calculate intersection of triangle with ray. Sets result
// accordingly
bool operator()
(
const label index,
const point& start,
const point& end,
point& intersectionPoint
) const;
};
// Declare name of the class and its debug switch
ClassName("treeDataPoint");
@ -106,6 +159,11 @@ public:
return points_;
}
bool useSubset() const
{
return useSubset_;
}
//- Get representative point cloud for all shapes inside
// (one point per shape)
pointField shapePoints() const;
@ -115,7 +173,7 @@ public:
//- Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
// Only makes sense for closed surfaces.
label getVolumeType
volumeType getVolumeType
(
const indexedOctree<treeDataPoint>&,
const point&
@ -135,50 +193,6 @@ public:
const point& centre,
const scalar radiusSqr
) const;
//- Calculates nearest (to sample) point in shape.
// Returns actual point and distance (squared)
void findNearest
(
const labelUList& indices,
const point& sample,
scalar& nearestDistSqr,
label& nearestIndex,
point& nearestPoint
) const;
//- Calculates nearest (to line) point in shape.
// Returns point and distance (squared)
void findNearest
(
const labelUList& indices,
const linePointRef& ln,
treeBoundBox& tightest,
label& minIndex,
point& linePoint,
point& nearestPoint
) const;
//- Calculate intersection of shape with ray. Sets result
// accordingly
bool intersects
(
const label index,
const point& start,
const point& end,
point& result
) const
{
notImplemented
(
"treeDataPoint::intersects(const label, const point&,"
"const point&, point&)"
);
return false;
}
};

252
src/meshTools/indexedOctree/treeDataPrimitivePatch.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -26,12 +26,8 @@ License
#include "treeDataPrimitivePatch.H"
#include "indexedOctree.H"
#include "triangleFuncs.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
template<class PatchType>
Foam::scalar Foam::treeDataPrimitivePatch<PatchType>::tolSqr = sqr(1e-6);
#include "triSurfaceTools.H"
#include "triFace.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@ -70,6 +66,75 @@ void Foam::treeDataPrimitivePatch<PatchType>::update()
}
template<class PatchType>
bool Foam::treeDataPrimitivePatch<PatchType>::findIntersection
(
const indexedOctree<treeDataPrimitivePatch<PatchType> >& tree,
const label index,
const point& start,
const point& end,
point& intersectionPoint
)
{
const treeDataPrimitivePatch<PatchType>& shape = tree.shapes();
const PatchType& patch = shape.patch();
const pointField& points = patch.points();
const typename PatchType::FaceType& f = patch[index];
// Do quick rejection test
if (shape.cacheBb_)
{
const treeBoundBox& faceBb = shape.bbs_[index];
if ((faceBb.posBits(start) & faceBb.posBits(end)) != 0)
{
// start and end in same block outside of faceBb.
return false;
}
}
const vector dir(end - start);
pointHit inter;
if (f.size() == 3)
{
inter = triPointRef
(
points[f[0]],
points[f[1]],
points[f[2]]
).intersection(start, dir, intersection::HALF_RAY, shape.planarTol_);
}
else
{
const pointField& faceCentres = patch.faceCentres();
inter = f.intersection
(
start,
dir,
faceCentres[index],
points,
intersection::HALF_RAY
);
}
if (inter.hit() && inter.distance() <= 1)
{
// Note: no extra test on whether intersection is in front of us
// since using half_ray
intersectionPoint = inter.hitPoint();
return true;
}
else
{
return false;
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
@ -77,16 +142,50 @@ template<class PatchType>
Foam::treeDataPrimitivePatch<PatchType>::treeDataPrimitivePatch
(
const bool cacheBb,
const PatchType& patch
const PatchType& patch,
const scalar planarTol
)
:
patch_(patch),
cacheBb_(cacheBb)
cacheBb_(cacheBb),
planarTol_(planarTol)
{
update();
}
template<class PatchType>
Foam::treeDataPrimitivePatch<PatchType>::findNearestOp::findNearestOp
(
const indexedOctree<treeDataPrimitivePatch<PatchType> >& tree
)
:
tree_(tree)
{}
template<class PatchType>
Foam::treeDataPrimitivePatch<PatchType>::findIntersectOp::findIntersectOp
(
const indexedOctree<treeDataPrimitivePatch<PatchType> >& tree
)
:
tree_(tree)
{}
template<class PatchType>
Foam::treeDataPrimitivePatch<PatchType>::findAllIntersectOp::findAllIntersectOp
(
const indexedOctree<treeDataPrimitivePatch<PatchType> >& tree,
DynamicList<label>& shapeMask
)
:
tree_(tree),
shapeMask_(shapeMask)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class PatchType>
@ -106,7 +205,7 @@ Foam::pointField Foam::treeDataPrimitivePatch<PatchType>::shapePoints() const
//- Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
// Only makes sense for closed surfaces.
template<class PatchType>
Foam::label Foam::treeDataPrimitivePatch<PatchType>:: getVolumeType
Foam::volumeType Foam::treeDataPrimitivePatch<PatchType>::getVolumeType
(
const indexedOctree<treeDataPrimitivePatch<PatchType> >& oc,
const point& sample
@ -136,7 +235,6 @@ Foam::label Foam::treeDataPrimitivePatch<PatchType>:: getVolumeType
<< abort(FatalError);
}
// Get actual intersection point on face
label faceI = info.index();
@ -147,7 +245,7 @@ Foam::label Foam::treeDataPrimitivePatch<PatchType>:: getVolumeType
}
const pointField& points = patch_.localPoints();
const face& f = patch_.localFaces()[faceI];
const typename PatchType::FaceType& f = patch_.localFaces()[faceI];
// Retest to classify where on face info is. Note: could be improved. We
// already have point.
@ -188,9 +286,10 @@ Foam::label Foam::treeDataPrimitivePatch<PatchType>:: getVolumeType
const scalar typDimSqr = mag(area) + VSMALL;
forAll(f, fp)
{
if ((magSqr(points[f[fp]] - curPt)/typDimSqr) < tolSqr)
if ((magSqr(points[f[fp]] - curPt)/typDimSqr) < planarTol_)
{
// Face intersection point equals face vertex fp
@ -207,7 +306,7 @@ Foam::label Foam::treeDataPrimitivePatch<PatchType>:: getVolumeType
const point fc(f.centre(points));
if ((magSqr(fc - curPt)/typDimSqr) < tolSqr)
if ((magSqr(fc - curPt)/typDimSqr) < planarTol_)
{
// Face intersection point equals face centre. Normal at face centre
// is already average of face normals
@ -240,7 +339,7 @@ Foam::label Foam::treeDataPrimitivePatch<PatchType>:: getVolumeType
pointHit edgeHit = e.line(points).nearestDist(sample);
if ((magSqr(edgeHit.rawPoint() - curPt)/typDimSqr) < tolSqr)
if ((magSqr(edgeHit.rawPoint() - curPt)/typDimSqr) < planarTol_)
{
// Face intersection point lies on edge e
@ -285,7 +384,7 @@ Foam::label Foam::treeDataPrimitivePatch<PatchType>:: getVolumeType
fc
).nearestDist(sample);
if ((magSqr(edgeHit.rawPoint() - curPt)/typDimSqr) < tolSqr)
if ((magSqr(edgeHit.rawPoint() - curPt)/typDimSqr) < planarTol_)
{
// Face intersection point lies on edge between two face triangles
@ -336,7 +435,7 @@ Foam::label Foam::treeDataPrimitivePatch<PatchType>:: getVolumeType
// - tolerances are wrong. (if e.g. face has zero area)
// - or (more likely) surface is not closed.
return indexedOctree<treeDataPrimitivePatch>::UNKNOWN;
return volumeType::UNKNOWN;
}
@ -368,7 +467,7 @@ bool Foam::treeDataPrimitivePatch<PatchType>::overlaps
// 2. Check if one or more face points inside
const pointField& points = patch_.points();
const face& f = patch_[index];
const typename PatchType::FaceType& f = patch_[index];
if (cubeBb.containsAny(points, f))
{
@ -379,21 +478,35 @@ bool Foam::treeDataPrimitivePatch<PatchType>::overlaps
// go through cube. Use triangle-bounding box intersection.
const point fc = f.centre(points);
forAll(f, fp)
if (f.size() == 3)
{
bool triIntersects = triangleFuncs::intersectBb
return triangleFuncs::intersectBb
(
points[f[fp]],
points[f[f.fcIndex(fp)]],
fc,
points[f[0]],
points[f[1]],
points[f[2]],
cubeBb
);
if (triIntersects)
}
else
{
forAll(f, fp)
{
return true;
bool triIntersects = triangleFuncs::intersectBb
(
points[f[fp]],
points[f[f.fcIndex(fp)]],
fc,
cubeBb
);
if (triIntersects)
{
return true;
}
}
}
return false;
}
@ -439,10 +552,8 @@ bool Foam::treeDataPrimitivePatch<PatchType>::overlaps
}
// Calculate nearest point to sample. Updates (if any) nearestDistSqr, minIndex,
// nearestPoint.
template<class PatchType>
void Foam::treeDataPrimitivePatch<PatchType>::findNearest
void Foam::treeDataPrimitivePatch<PatchType>::findNearestOp::operator()
(
const labelUList& indices,
const point& sample,
@ -452,13 +563,15 @@ void Foam::treeDataPrimitivePatch<PatchType>::findNearest
point& nearestPoint
) const
{
const pointField& points = patch_.points();
const treeDataPrimitivePatch<PatchType>& shape = tree_.shapes();
const PatchType& patch = shape.patch();
const pointField& points = patch.points();
forAll(indices, i)
{
const label index = indices[i];
const face& f = patch_[index];
const typename PatchType::FaceType& f = patch[index];
pointHit nearHit = f.nearestPoint(sample, points);
scalar distSqr = sqr(nearHit.distance());
@ -474,7 +587,34 @@ void Foam::treeDataPrimitivePatch<PatchType>::findNearest
template<class PatchType>
bool Foam::treeDataPrimitivePatch<PatchType>::intersects
void Foam::treeDataPrimitivePatch<PatchType>::findNearestOp::operator()
(
const labelUList& indices,
const linePointRef& ln,
treeBoundBox& tightest,
label& minIndex,
point& linePoint,
point& nearestPoint
) const
{
notImplemented
(
"treeDataPrimitivePatch<PatchType>::findNearestOp::operator()"
"("
" const labelUList&,"
" const linePointRef&,"
" treeBoundBox&,"
" label&,"
" point&,"
" point&"
") const"
);
}
template<class PatchType>
bool Foam::treeDataPrimitivePatch<PatchType>::findIntersectOp::operator()
(
const label index,
const point& start,
@ -482,43 +622,25 @@ bool Foam::treeDataPrimitivePatch<PatchType>::intersects
point& intersectionPoint
) const
{
// Do quick rejection test
if (cacheBb_)
{
const treeBoundBox& faceBb = bbs_[index];
return findIntersection(tree_, index, start, end, intersectionPoint);
}
if ((faceBb.posBits(start) & faceBb.posBits(end)) != 0)
{
// start and end in same block outside of faceBb.
return false;
}
}
const pointField& points = patch_.points();
const face& f = patch_[index];
const point fc = f.centre(points);
const vector dir(end - start);
pointHit inter = patch_[index].intersection
(
start,
dir,
fc,
points,
intersection::HALF_RAY
);
if (inter.hit() && inter.distance() <= 1)
{
// Note: no extra test on whether intersection is in front of us
// since using half_ray
intersectionPoint = inter.hitPoint();
return true;
}
else
template<class PatchType>
bool Foam::treeDataPrimitivePatch<PatchType>::findAllIntersectOp::operator()
(
const label index,
const point& start,
const point& end,
point& intersectionPoint
) const
{
if (!shapeMask_.empty() && findIndex(shapeMask_, index) != -1)
{
return false;
}
return findIntersection(tree_, index, start, end, intersectionPoint);
}

152
src/meshTools/indexedOctree/treeDataPrimitivePatch.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -35,9 +35,8 @@ SourceFiles
#ifndef treeDataPrimitivePatch_H
#define treeDataPrimitivePatch_H
#include "PrimitivePatch.H"
//#include "indexedOctree.H"
#include "treeBoundBoxList.H"
#include "volumeType.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -64,11 +63,6 @@ class treeDataPrimitivePatch
:
public treeDataPrimitivePatchName
{
// Static data
//- tolerance on linear dimensions
static scalar tolSqr;
// Private data
//- Underlying geometry
@ -77,6 +71,9 @@ class treeDataPrimitivePatch
//- Whether to precalculate and store face bounding box
const bool cacheBb_;
//- Tolerance to use for intersection tests
const scalar planarTol_;
//- face bounding boxes (valid only if cacheBb_)
treeBoundBoxList bbs_;
@ -89,15 +86,107 @@ class treeDataPrimitivePatch
//- Initialise all member data
void update();
//- Find intersection of line with shapes
static bool findIntersection
(
const indexedOctree<treeDataPrimitivePatch<PatchType> >& tree,
const label index,
const point& start,
const point& end,
point& intersectionPoint
);
public:
class findNearestOp
{
const indexedOctree<treeDataPrimitivePatch>& tree_;
public:
findNearestOp(const indexedOctree<treeDataPrimitivePatch>& tree);
void operator()
(
const labelUList& indices,
const point& sample,
scalar& nearestDistSqr,
label& minIndex,
point& nearestPoint
) const;
//- Calculates nearest (to line) point in shape.
// Returns point and distance (squared)
void operator()
(
const labelUList& indices,
const linePointRef& ln,
treeBoundBox& tightest,
label& minIndex,
point& linePoint,
point& nearestPoint
) const;
};
class findIntersectOp
{
const indexedOctree<treeDataPrimitivePatch>& tree_;
public:
findIntersectOp(const indexedOctree<treeDataPrimitivePatch>& tree);
//- Calculate intersection of triangle with ray. Sets result
// accordingly
bool operator()
(
const label index,
const point& start,
const point& end,
point& intersectionPoint
) const;
};
class findAllIntersectOp
{
const indexedOctree<treeDataPrimitivePatch>& tree_;
DynamicList<label>& shapeMask_;
public:
findAllIntersectOp
(
const indexedOctree<treeDataPrimitivePatch>& tree,
DynamicList<label>& shapeMask
);
//- Calculate intersection of triangle with ray. Sets result
// accordingly
bool operator()
(
const label index,
const point& start,
const point& end,
point& intersectionPoint
) const;
};
// Constructors
//- Construct from patch.
treeDataPrimitivePatch
(
const bool cacheBb,
const PatchType&
const PatchType&,
const scalar planarTol
);
@ -125,7 +214,7 @@ public:
//- Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
// Only makes sense for closed surfaces.
label getVolumeType
volumeType getVolumeType
(
const indexedOctree<treeDataPrimitivePatch<PatchType> >&,
const point&
@ -145,49 +234,6 @@ public:
const point& centre,
const scalar radiusSqr
) const;
//- Calculates nearest (to sample) point in shape.
// Returns actual point and distance (squared)
void findNearest
(
const labelUList& indices,
const point& sample,
scalar& nearestDistSqr,
label& nearestIndex,
point& nearestPoint
) const;
//- Calculates nearest (to line) point in shape.
// Returns point and distance (squared)
void findNearest
(
const labelUList& indices,
const linePointRef& ln,
treeBoundBox& tightest,
label& minIndex,
point& linePoint,
point& nearestPoint
) const
{
notImplemented
(
"treeDataPrimitivePatch::findNearest"
"(const labelUList&, const linePointRef&, ..)"
);
}
//- Calculate intersection of shape with ray. Sets result
// accordingly
bool intersects
(
const label index,
const point& start,
const point& end,
point& result
) const;
};

477
src/meshTools/indexedOctree/treeDataTriSurface.C
View File

@ -25,491 +25,58 @@ License
#include "treeDataTriSurface.H"
#include "triSurfaceTools.H"
#include "triangleFuncs.H"
#include "indexedOctree.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(treeDataTriSurface, 0);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// // Fast distance to triangle calculation. From
// // "Distance Between Point and Triangle in 3D"
// // David Eberly, Magic Software Inc. Aug. 2003.
// // Works on function Q giving distance to point and tries to minimize this.
// Foam::scalar Foam::treeDataTriSurface::nearestCoords
// (
// const point& base,
// const point& E0,
// const point& E1,
// const scalar a,
// const scalar b,
// const scalar c,
// const point& P,
// scalar& s,
// scalar& t
// )
// {
// // distance vector
// const vector D(base - P);
// // Precalculate distance factors.
// const scalar d = E0 & D;
// const scalar e = E1 & D;
// // Do classification
// const scalar det = a*c - b*b;
// s = b*e - c*d;
// t = b*d - a*e;
// if (s+t < det)
// {
// if (s < 0)
// {
// if (t < 0)
// {
// //region 4
// if (e > 0)
// {
// //min on edge t = 0
// t = 0;
// s = (d >= 0 ? 0 : (-d >= a ? 1 : -d/a));
// }
// else
// {
// //min on edge s=0
// s = 0;
// t = (e >= 0 ? 0 : (-e >= c ? 1 : -e/c));
// }
// }
// else
// {
// //region 3. Min on edge s = 0
// s = 0;
// t = (e >= 0 ? 0 : (-e >= c ? 1 : -e/c));
// }
// }
// else if (t < 0)
// {
// //region 5
// t = 0;
// s = (d >= 0 ? 0 : (-d >= a ? 1 : -d/a));
// }
// else
// {
// //region 0
// const scalar invDet = 1/det;
// s *= invDet;
// t *= invDet;
// }
// }
// else
// {
// if (s < 0)
// {
// //region 2
// const scalar tmp0 = b + d;
// const scalar tmp1 = c + e;
// if (tmp1 > tmp0)
// {
// //min on edge s+t=1
// const scalar numer = tmp1 - tmp0;
// const scalar denom = a-2*b+c;
// s = (numer >= denom ? 1 : numer/denom);
// t = 1 - s;
// }
// else
// {
// //min on edge s=0
// s = 0;
// t = (tmp1 <= 0 ? 1 : (e >= 0 ? 0 : - e/c));
// }
// }
// else if (t < 0)
// {
// //region 6
// const scalar tmp0 = b + d;
// const scalar tmp1 = c + e;
// if (tmp1 > tmp0)
// {
// //min on edge s+t=1
// const scalar numer = tmp1 - tmp0;
// const scalar denom = a-2*b+c;
// s = (numer >= denom ? 1 : numer/denom);
// t = 1 - s;
// }
// else
// {
// //min on edge t=0
// t = 0;
// s = (tmp1 <= 0 ? 1 : (d >= 0 ? 0 : - d/a));
// }
// }
// else
// {
// //region 1
// const scalar numer = c+e-(b+d);
// if (numer <= 0)
// {
// s = 0;
// }
// else
// {
// const scalar denom = a-2*b+c;
// s = (numer >= denom ? 1 : numer/denom);
// }
// }
// t = 1 - s;
// }
// // Calculate distance.
// // Note: abs should not be needed but truncation error causes problems
// // with points very close to one of the triangle vertices.
// // (seen up to -9e-15). Alternatively add some small value.
// const scalar f = D & D;
// return Foam::mag(a*s*s + 2*b*s*t + c*t*t + 2*d*s + 2*e*t + f);
// }
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
Foam::treeDataTriSurface::treeDataTriSurface
template<>
Foam::volumeType Foam::treeDataPrimitivePatch<Foam::triSurface>::getVolumeType
(
const triSurface& surface,
const scalar planarTol
)
:
surface_(surface),
planarTol_(planarTol)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::pointField Foam::treeDataTriSurface::shapePoints() const
{
const pointField& points = surface_.points();
pointField centres(surface_.size());
forAll(surface_, triI)
{
centres[triI] = surface_[triI].centre(points);
}
return centres;
}
//- Get type of sample (inside/outside/mixed) w.r.t. surface.
Foam::label Foam::treeDataTriSurface::getVolumeType
(
const indexedOctree<treeDataTriSurface>& tree,
const indexedOctree<treeDataPrimitivePatch<triSurface> >& oc,
const point& sample
) const
{
// Find nearest point
const treeBoundBox& treeBb = tree.bb();
// Find nearest face to sample
pointIndexHit info = oc.findNearest(sample, sqr(GREAT));
pointIndexHit pHit = tree.findNearest
(
sample,
max
(
Foam::sqr(GREAT),
Foam::magSqr(treeBb.span())
)
);
if (!pHit.hit())
if (info.index() == -1)
{
FatalErrorIn("treeDataTriSurface::getVolumeType(..)")
<< "treeBb:" << treeBb
<< " sample:" << sample
<< " pHit:" << pHit
FatalErrorIn
(
"treeDataPrimitivePatch::getSampleType"
"(indexedOctree<treeDataPrimitivePatch>&, const point&)"
) << "Could not find " << sample << " in octree."
<< abort(FatalError);
}
// Get actual intersection point on face
label faceI = info.index();
triSurfaceTools::sideType t = triSurfaceTools::surfaceSide
(
surface_,
patch_,
sample,
pHit.index()
faceI
);
if (t == triSurfaceTools::UNKNOWN)
{
return indexedOctree<treeDataTriSurface>::UNKNOWN;
return volumeType::UNKNOWN;
}
else if (t == triSurfaceTools::INSIDE)
{
return indexedOctree<treeDataTriSurface>::INSIDE;
return volumeType::INSIDE;
}
else if (t == triSurfaceTools::OUTSIDE)
{
return indexedOctree<treeDataTriSurface>::OUTSIDE;
return volumeType::OUTSIDE;
}
else
{
FatalErrorIn("treeDataTriSurface::getVolumeType(..)")
FatalErrorIn("treeDataPrimitivePatch<PatchType>::getVolumeType(..)")
<< "problem" << abort(FatalError);
return indexedOctree<treeDataTriSurface>::UNKNOWN;
return volumeType::UNKNOWN;
}
}
// Check if any point on triangle is inside cubeBb.
bool Foam::treeDataTriSurface::overlaps
(
const label index,
const treeBoundBox& cubeBb
) const
{
const pointField& points = surface_.points();
const labelledTri& f = surface_[index];
treeBoundBox triBb(points, surface_[index]);
//- For testing: robust one
//return cubeBb.overlaps(triBb);
//- Exact test of triangle intersecting bb
// Quick rejection. If whole bounding box of tri is outside cubeBb then
// there will be no intersection.
if (!cubeBb.overlaps(triBb))
{
return false;
}
// Check if one or more triangle point inside
if (cubeBb.containsAny(points, f))
{
return true;
}
// Triangle points
const point& p0 = points[f[0]];
const point& p1 = points[f[1]];
const point& p2 = points[f[2]];
// Now we have the difficult case: all points are outside but connecting
// edges might go through cube. Use fast intersection of bounding box.
//return triangleFuncs::intersectBbExact(p0, p1, p2, cubeBb);
return triangleFuncs::intersectBb(p0, p1, p2, cubeBb);
}
// Calculate nearest point to sample. Updates (if any) nearestDistSqr, minIndex,
// nearestPoint.
void Foam::treeDataTriSurface::findNearest
(
const labelUList& indices,
const point& sample,
scalar& nearestDistSqr,
label& minIndex,
point& nearestPoint
) const
{
const pointField& points = surface_.points();
forAll(indices, i)
{
label index = indices[i];
const triSurface::FaceType& f = surface_[index];
////- Possible optimization: do quick rejection of triangle if bounding
//// sphere does not intersect triangle bounding box. From simplistic
//// test was not found to speed up things.
//
//// Triangle bounding box.
//point triBbMin = min(p0, min(p1, p2));
//point triBbMax = max(p0, max(p1, p2));
//
//if
//(
// indexedOctree<treeDataTriSurface>::intersects
// (
// triBbMin,
// triBbMax,
// nearestDistSqr,
// sample
// )
//)
{
// // Get spanning vectors of triangle
// vector base(p1);
// vector E0(p0 - p1);
// vector E1(p2 - p1);
// scalar a(E0& E0);
// scalar b(E0& E1);
// scalar c(E1& E1);
// // Get nearest point in s,t coordinates (s is along E0, t
// // is along E1)
// scalar s;
// scalar t;
// scalar distSqr = nearestCoords
// (
// base,
// E0,
// E1,
// a,
// b,
// c,
// sample,
// s,
// t
// );
pointHit pHit = f.nearestPoint(sample, points);
scalar distSqr = sqr(pHit.distance());
if (distSqr < nearestDistSqr)
{
nearestDistSqr = distSqr;
minIndex = index;
nearestPoint = pHit.rawPoint();
}
}
}
}
// Calculate nearest point to line. Updates (if any) nearestDistSqr, minIndex,
// nearestPoint.
void Foam::treeDataTriSurface::findNearest
(
const labelUList& indices,
const linePointRef& ln,
treeBoundBox& tightest,
label& minIndex,
point& linePoint,
point& nearestPoint
) const
{
// Best so far
scalar nearestDistSqr = VGREAT;
if (minIndex >= 0)
{
nearestDistSqr = magSqr(linePoint - nearestPoint);
}
const pointField& points = surface_.points();
forAll(indices, i)
{
label index = indices[i];
const triSurface::FaceType& f = surface_[index];
triPointRef tri(f.tri(points));
pointHit lineInfo(point::max);
pointHit pHit = tri.nearestPoint(ln, lineInfo);
scalar distSqr = sqr(pHit.distance());
if (distSqr < nearestDistSqr)
{
nearestDistSqr = distSqr;
minIndex = index;
linePoint = lineInfo.rawPoint();
nearestPoint = pHit.rawPoint();
{
point& minPt = tightest.min();
minPt = min(ln.start(), ln.end());
minPt.x() -= pHit.distance();
minPt.y() -= pHit.distance();
minPt.z() -= pHit.distance();
}
{
point& maxPt = tightest.max();
maxPt = max(ln.start(), ln.end());
maxPt.x() += pHit.distance();
maxPt.y() += pHit.distance();
maxPt.z() += pHit.distance();
}
}
}
}
bool Foam::treeDataTriSurface::intersects
(
const label index,
const point& start,
const point& end,
point& intersectionPoint
) const
{
const pointField& points = surface_.points();
const triSurface::FaceType& f = surface_[index];
// Do quick rejection test
treeBoundBox triBb(points, f);
const direction startBits(triBb.posBits(start));
const direction endBits(triBb.posBits(end));
if ((startBits & endBits) != 0)
{
// start and end in same block outside of triBb.
return false;
}
const vector dir(end - start);
// Use relative tolerance (from octree) to determine intersection.
pointHit inter = f.intersection
(
start,
dir,
points,
intersection::HALF_RAY,
planarTol_
);
if (inter.hit() && inter.distance() <= 1)
{
// Note: no extra test on whether intersection is in front of us
// since using half_ray.
intersectionPoint = inter.hitPoint();
return true;
}
else
{
return false;
}
//- Using exact intersections
//pointHit info = f.tri(points).intersectionExact(start, end);
//
//if (info.hit())
//{
// intersectionPoint = info.hitPoint();
//}
//return info.hit();
}
// ************************************************************************* //

142
src/meshTools/indexedOctree/treeDataTriSurface.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -25,7 +25,7 @@ Class
Foam::treeDataTriSurface
Description
Encapsulates data for (indexedOc)tree searches on triSurface.
Encapsulates data for (indexedOc)tree searches on a triSurface.
SourceFiles
treeDataTriSurface.C
@ -35,139 +35,27 @@ SourceFiles
#ifndef treeDataTriSurface_H
#define treeDataTriSurface_H
#include "treeDataPrimitivePatch.H"
#include "triSurface.H"
#include "point.H"
#include "indexedOctree.H"
#include "volumeType.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
typedef treeDataPrimitivePatch<triSurface> treeDataTriSurface;
// Forward declaration of classes
class treeBoundBox;
class treeDataTriSurface;
template<class Type> class indexedOctree;
//- Template specialisation of getVolumeType for treeDataTriSurface
template<>
volumeType treeDataPrimitivePatch<triSurface>::getVolumeType
(
const indexedOctree<treeDataPrimitivePatch<triSurface> >& oc,
const point& sample
) const;
}
/*---------------------------------------------------------------------------*\
Class treeDataTriSurface Declaration
\*---------------------------------------------------------------------------*/
class treeDataTriSurface
{
// Private data
//- Reference to triSurface
const triSurface& surface_;
//- Tolerance to use for intersection tests
const scalar planarTol_;
// Private Member Functions
// //- fast triangle nearest point calculation. Returns point in E0, E1
// // coordinate system: base + s*E0 + t*E1
// static scalar nearestCoords
// (
// const point& base,
// const point& E0,
// const point& E1,
// const scalar a,
// const scalar b,
// const scalar c,
// const point& P,
// scalar& s,
// scalar& t
// );
public:
// Declare name of the class and its debug switch
ClassName("treeDataTriSurface");
// Constructors
//- Construct from triSurface and tolerance for intersection
// tests. Holds reference.
treeDataTriSurface(const triSurface&, const scalar planarTol);
// Member Functions
// Access
inline const triSurface& surface() const
{
return surface_;
}
inline label size() const
{
return surface_.size();
}
//- Get representative point cloud for all shapes inside
// (one point per shape)
pointField shapePoints() const;
// Search
//- Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
// Only makes sense for closed surfaces.
label getVolumeType
(
const indexedOctree<treeDataTriSurface>&,
const point&
) const;
//- Does (bb of) shape at index overlap bb
bool overlaps
(
const label index,
const treeBoundBox& sampleBb
) const;
//- Calculates nearest (to sample) point in shape.
// Returns actual point and distance (squared)
void findNearest
(
const labelUList& indices,
const point& sample,
scalar& nearestDistSqr,
label& nearestIndex,
point& nearestPoint
) const;
//- Calculates nearest (to line) point in shape.
// Returns point and distance (squared)
void findNearest
(
const labelUList& indices,
const linePointRef& ln,
treeBoundBox& tightest,
label& minIndex,
point& linePoint,
point& nearestPoint
) const;
//- Calculate intersection of triangle with ray. Sets result
// accordingly
bool intersects
(
const label index,
const point& start,
const point& end,
point& result
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

6
src/meshTools/meshSearch/meshSearch.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -947,9 +947,7 @@ Foam::List<Foam::pointIndexHit> Foam::meshSearch::intersections
bool Foam::meshSearch::isInside(const point& p) const
{
return
boundaryTree().getVolumeType(p)
== indexedOctree<treeDataFace>::INSIDE;
return (boundaryTree().getVolumeType(p) == volumeType::INSIDE);
}

6
src/meshTools/searchableSurface/searchableBox.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -583,7 +583,7 @@ void Foam::searchableBox::getVolumeType
) const
{
volType.setSize(points.size());
volType = INSIDE;
volType = volumeType::INSIDE;
forAll(points, pointI)
{
@ -593,7 +593,7 @@ void Foam::searchableBox::getVolumeType
{
if (pt[dir] < min()[dir] || pt[dir] > max()[dir])
{
volType[pointI] = OUTSIDE;
volType[pointI] = volumeType::OUTSIDE;
break;
}
}

12
src/meshTools/searchableSurface/searchableCylinder.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -673,7 +673,7 @@ void Foam::searchableCylinder::getVolumeType
) const
{
volType.setSize(points.size());
volType = INSIDE;
volType = volumeType::INSIDE;
forAll(points, pointI)
{
@ -687,12 +687,12 @@ void Foam::searchableCylinder::getVolumeType
if (parallel < 0)
{
// left of point1 endcap
volType[pointI] = OUTSIDE;
volType[pointI] = volumeType::OUTSIDE;
}
else if (parallel > magDir_)
{
// right of point2 endcap
volType[pointI] = OUTSIDE;
volType[pointI] = volumeType::OUTSIDE;
}
else
{
@ -701,11 +701,11 @@ void Foam::searchableCylinder::getVolumeType
if (mag(v) > radius_)
{
volType[pointI] = OUTSIDE;
volType[pointI] = volumeType::OUTSIDE;
}
else
{
volType[pointI] = INSIDE;
volType[pointI] = volumeType::INSIDE;
}
}
}

8
src/meshTools/searchableSurface/searchableSphere.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -334,7 +334,7 @@ void Foam::searchableSphere::getVolumeType
) const
{
volType.setSize(points.size());
volType = INSIDE;
volType = volumeType::INSIDE;
forAll(points, pointI)
{
@ -342,11 +342,11 @@ void Foam::searchableSphere::getVolumeType
if (magSqr(pt - centre_) <= sqr(radius_))
{
volType[pointI] = INSIDE;
volType[pointI] = volumeType::INSIDE;
}
else
{
volType[pointI] = OUTSIDE;
volType[pointI] = volumeType::OUTSIDE;
}
}
}

22
src/meshTools/searchableSurface/searchableSurface.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -30,26 +30,12 @@ License
namespace Foam
{
defineTypeNameAndDebug(searchableSurface, 0);
defineRunTimeSelectionTable(searchableSurface, dict);
template<>
const char* Foam::NamedEnum
<
Foam::searchableSurface::volumeType,
4
>::names[] =
{
"unknown",
"mixed",
"inside",
"outside"
};
defineTypeNameAndDebug(searchableSurface, 0);
defineRunTimeSelectionTable(searchableSurface, dict);
}
const Foam::NamedEnum<Foam::searchableSurface::volumeType, 4>
Foam::searchableSurface::volumeTypeNames;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

35
src/meshTools/searchableSurface/searchableSurface.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -48,6 +48,7 @@ SourceFiles
#include "pointIndexHit.H"
#include "linePointRef.H"
#include "objectRegistry.H"
#include "volumeType.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -67,25 +68,6 @@ class searchableSurface
:
public regIOobject
{
public:
// Data types
//- Volume types
enum volumeType
{
UNKNOWN = 0,
MIXED = 1, // not used. only here to maintain consistency with
// indexedOctree volumeType.
INSIDE = 2,
OUTSIDE = 3
};
static const NamedEnum<volumeType, 4> volumeTypeNames;
private:
// Private data
const word name_;
@ -278,6 +260,19 @@ public:
List<pointIndexHit>&
) const = 0;
//- Find the nearest locations for the supplied points to a
// particular region in the searchable surface.
virtual void findNearest
(
const pointField& samples,
const scalarField& nearestDistSqr,
const labelList& regionIndices,
List<pointIndexHit>& info
) const
{
findNearest(samples, nearestDistSqr, info);
}
//- Find first intersection on segment from start to end.
// Note: searchableSurfacesQueries expects no
// intersection to be found if start==end. Is problem?

34
src/meshTools/searchableSurface/searchableSurfaces.C
View File

@ -293,6 +293,18 @@ Foam::label Foam::searchableSurfaces::findSurfaceID
}
Foam::label Foam::searchableSurfaces::findSurfaceRegionID
(
const word& surfaceName,
const word& regionName
) const
{
label surfaceIndex = findSurfaceID(surfaceName);
return findIndex(regionNames()[surfaceIndex], regionName);
}
// Find any intersection
void Foam::searchableSurfaces::findAnyIntersection
(
@ -382,6 +394,28 @@ void Foam::searchableSurfaces::findNearest
}
// Find nearest. Return -1 or nearest point
void Foam::searchableSurfaces::findNearest
(
const pointField& samples,
const scalarField& nearestDistSqr,
const labelList& regionIndices,
labelList& nearestSurfaces,
List<pointIndexHit>& nearestInfo
) const
{
searchableSurfacesQueries::findNearest
(
*this,
allSurfaces_,
samples,
nearestDistSqr,
regionIndices,
nearestSurfaces,
nearestInfo
);
}
//- Calculate bounding box
Foam::boundBox Foam::searchableSurfaces::bounds() const
{

16
src/meshTools/searchableSurface/searchableSurfaces.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -140,6 +140,11 @@ public:
//- Find index of surface. Return -1 if not found.
label findSurfaceID(const word& name) const;
label findSurfaceRegionID
(
const word& surfaceName,
const word& regionName
) const;
// Multiple point queries.
@ -185,6 +190,15 @@ public:
List<pointIndexHit>&
) const;
void findNearest
(
const pointField& samples,
const scalarField& nearestDistSqr,
const labelList& regionIndices,
labelList& nearestSurfaces,
List<pointIndexHit>& nearestInfo
) const;
//- Calculate bounding box
boundBox bounds() const;

114
src/meshTools/searchableSurface/searchableSurfacesQueries.C
View File

@ -324,9 +324,8 @@ bool Foam::searchableSurfacesQueries::morphTet
void Foam::searchableSurfacesQueries::mergeHits
(
const point& start,
const scalar mergeDist,
const label testI, // index of surface
const label testI, // index of surface
const List<pointIndexHit>& surfHits, // hits on surface
labelList& allSurfaces,
@ -338,7 +337,7 @@ void Foam::searchableSurfacesQueries::mergeHits
scalarList surfDistSqr(surfHits.size());
forAll(surfHits, i)
{
surfDistSqr[i] = magSqr(surfHits[i].hitPoint()-start);
surfDistSqr[i] = magSqr(surfHits[i].hitPoint() - start);
}
forAll(surfDistSqr, i)
@ -346,11 +345,7 @@ void Foam::searchableSurfacesQueries::mergeHits
label index = findLower(allDistSqr, surfDistSqr[i]);
// Check if equal to lower.
if
(
index >= 0
&& (mag(allDistSqr[index]-surfDistSqr[i]) < mergeDist)
)
if (index >= 0)
{
// Same. Do not count.
//Pout<< "point:" << surfHits[i].hitPoint()
@ -361,12 +356,9 @@ void Foam::searchableSurfacesQueries::mergeHits
else
{
// Check if equal to higher
label next = index+1;
if
(
next < allDistSqr.size()
&& (mag(allDistSqr[next]-surfDistSqr[i]) < mergeDist)
)
label next = index + 1;
if (next < allDistSqr.size())
{
//Pout<< "point:" << surfHits[i].hitPoint()
// << " considered same as:" << allInfo[next].hitPoint()
@ -510,21 +502,6 @@ void Foam::searchableSurfacesQueries::findAllIntersections
if (surfacesToTest.size() > 1)
{
// Small vector to increment start vector by to find next intersection
// along line. Constant factor added to make sure that start==end still
// ends iteration in findAllIntersections. Also SMALL is just slightly
// too small.
const vectorField smallVec
(
1E2*SMALL*(end-start)
+ vector(ROOTVSMALL,ROOTVSMALL,ROOTVSMALL)
);
// Tolerance used to check whether points are equal. Note: used to
// compare distance^2. Note that we use the maximum possible tolerance
// (reached at intersections close to the end point)
const scalarField mergeDist(2*mag(smallVec)*mag(end-start));
// Test the other surfaces and merge (according to distance from start).
for (label testI = 1; testI < surfacesToTest.size(); testI++)
{
@ -541,7 +518,6 @@ void Foam::searchableSurfacesQueries::findAllIntersections
mergeHits
(
start[pointI], // Current segment
mergeDist[pointI],
testI, // Surface and its hits
surfHits[pointI],
@ -691,13 +667,75 @@ void Foam::searchableSurfacesQueries::findNearest
}
// Find nearest. Return -1 or nearest point
void Foam::searchableSurfacesQueries::findNearest
(
const PtrList<searchableSurface>& allSurfaces,
const labelList& surfacesToTest,
const pointField& samples,
const scalarField& nearestDistSqr,
const labelList& regionIndices,
labelList& nearestSurfaces,
List<pointIndexHit>& nearestInfo
)
{
if (regionIndices.empty())
{
findNearest
(
allSurfaces,
surfacesToTest,
samples,
nearestDistSqr,
nearestSurfaces,
nearestInfo
);
}
// Initialise
nearestSurfaces.setSize(samples.size());
nearestSurfaces = -1;
nearestInfo.setSize(samples.size());
// Work arrays
scalarField minDistSqr(nearestDistSqr);
List<pointIndexHit> hitInfo(samples.size());
forAll(surfacesToTest, testI)
{
allSurfaces[surfacesToTest[testI]].findNearest
(
samples,
minDistSqr,
regionIndices,
hitInfo
);
// Update minDistSqr and arguments
forAll(hitInfo, pointI)
{
if (hitInfo[pointI].hit())
{
minDistSqr[pointI] = magSqr
(
hitInfo[pointI].hitPoint()
- samples[pointI]
);
nearestInfo[pointI] = hitInfo[pointI];
nearestSurfaces[pointI] = testI;
}
}
}
}
void Foam::searchableSurfacesQueries::signedDistance
(
const PtrList<searchableSurface>& allSurfaces,
const labelList& surfacesToTest,
const pointField& samples,
const scalarField& nearestDistSqr,
const searchableSurface::volumeType illegalHandling,
const volumeType illegalHandling,
labelList& nearestSurfaces,
scalarField& distance
)
@ -737,7 +775,7 @@ void Foam::searchableSurfacesQueries::signedDistance
}
// Calculate sideness of these surface points
List<searchableSurface::volumeType> volType;
List<volumeType> volType;
allSurfaces[surfacesToTest[testI]].getVolumeType(surfPoints, volType);
// Push back to original
@ -746,13 +784,13 @@ void Foam::searchableSurfacesQueries::signedDistance
label pointI = surfIndices[i];
scalar dist = mag(samples[pointI] - nearestInfo[pointI].hitPoint());
searchableSurface::volumeType vT = volType[i];
volumeType vT = volType[i];
if (vT == searchableSurface::OUTSIDE)
if (vT == volumeType::OUTSIDE)
{
distance[pointI] = dist;
}
else if (vT == searchableSurface::INSIDE)
else if (vT == volumeType::INSIDE)
{
distance[i] = -dist;
}
@ -760,12 +798,12 @@ void Foam::searchableSurfacesQueries::signedDistance
{
switch (illegalHandling)
{
case searchableSurface::OUTSIDE:
case volumeType::OUTSIDE:
{
distance[pointI] = dist;
break;
}
case searchableSurface::INSIDE:
case volumeType::INSIDE:
{
distance[pointI] = -dist;
break;
@ -779,7 +817,7 @@ void Foam::searchableSurfacesQueries::signedDistance
<< " surface:"
<< allSurfaces[surfacesToTest[testI]].name()
<< " volType:"
<< searchableSurface::volumeTypeNames[vT]
<< volumeType::names[vT]
<< exit(FatalError);
break;
}

17
src/meshTools/searchableSurface/searchableSurfacesQueries.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -116,7 +116,6 @@ class searchableSurfacesQueries
static void mergeHits
(
const point& start,
const scalar mergeDist,
const label surfI,
const List<pointIndexHit>& surfHits,
@ -184,6 +183,18 @@ public:
List<pointIndexHit>&
);
//- Find nearest points to a specific region of the surface
static void findNearest
(
const PtrList<searchableSurface>& allSurfaces,
const labelList& surfacesToTest,
const pointField& samples,
const scalarField& nearestDistSqr,
const labelList& regionIndices,
labelList& nearestSurfaces,
List<pointIndexHit>& nearestInfo
);
//- Find signed distance to nearest surface. Outside is positive.
// illegalHandling: how to handle non-inside or outside
// OUTSIDE : treat as outside
@ -195,7 +206,7 @@ public:
const labelList& surfacesToTest,
const pointField& samples,
const scalarField& nearestDistSqr,
const searchableSurface::volumeType illegalHandling,
const volumeType illegalHandling,
labelList& nearestSurfaces,
scalarField& distance
);

319
src/meshTools/searchableSurface/triSurfaceMesh.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -29,7 +29,7 @@ License
#include "EdgeMap.H"
#include "triSurfaceFields.H"
#include "Time.H"
#include "PackedBoolList.H"
#include "PatchTools.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -218,99 +218,6 @@ bool Foam::triSurfaceMesh::isSurfaceClosed() const
}
// Gets all intersections after initial one. Adds smallVec and starts tracking
// from there.
void Foam::triSurfaceMesh::getNextIntersections
(
const indexedOctree<treeDataTriSurface>& octree,
const point& start,
const point& end,
const vector& smallVec,
DynamicList<pointIndexHit, 1, 1>& hits
)
{
const vector dirVec(end-start);
const scalar magSqrDirVec(magSqr(dirVec));
// Initial perturbation amount
vector perturbVec(smallVec);
while (true)
{
// Start tracking from last hit.
point pt = hits.last().hitPoint() + perturbVec;
if (((pt-start)&dirVec) > magSqrDirVec)
{
return;
}
// See if any intersection between pt and end
pointIndexHit inter = octree.findLine(pt, end);
if (!inter.hit())
{
return;
}
// Check if already found this intersection
bool duplicateHit = false;
forAllReverse(hits, i)
{
if (hits[i].index() == inter.index())
{
duplicateHit = true;
break;
}
}
if (duplicateHit)
{
// Hit same triangle again. Increase perturbVec and try again.
perturbVec *= 2;
}
else
{
// Proper hit
hits.append(inter);
// Restore perturbVec
perturbVec = smallVec;
}
}
}
void Foam::triSurfaceMesh::calcBounds(boundBox& bb, label& nPoints) const
{
// Unfortunately nPoints constructs meshPoints() so do compact version
// ourselves
const triSurface& s = static_cast<const triSurface&>(*this);
PackedBoolList pointIsUsed(points()().size());
nPoints = 0;
bb = boundBox::invertedBox;
forAll(s, faceI)
{
const triSurface::FaceType& f = s[faceI];
forAll(f, fp)
{
label pointI = f[fp];
if (pointIsUsed.set(pointI, 1u))
{
bb.min() = ::Foam::min(bb.min(), points()()[pointI]);
bb.max() = ::Foam::max(bb.max(), points()()[pointI]);
nPoints++;
}
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io, const triSurface& s)
@ -330,9 +237,8 @@ Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io, const triSurface& s)
)
),
triSurface(s),
tolerance_(indexedOctree<treeDataTriSurface>::perturbTol()),
triSurfaceRegionSearch(s),
minQuality_(-1),
maxTreeDepth_(10),
surfaceClosed_(-1)
{
bounds() = boundBox(points());
@ -377,9 +283,8 @@ Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io)
searchableSurface::objectPath()
)
),
tolerance_(indexedOctree<treeDataTriSurface>::perturbTol()),
triSurfaceRegionSearch(static_cast<const triSurface&>(*this)),
minQuality_(-1),
maxTreeDepth_(10),
surfaceClosed_(-1)
{
bounds() = boundBox(points());
@ -427,9 +332,8 @@ Foam::triSurfaceMesh::triSurfaceMesh
searchableSurface::objectPath()
)
),
tolerance_(indexedOctree<treeDataTriSurface>::perturbTol()),
triSurfaceRegionSearch(static_cast<const triSurface&>(*this), dict),
minQuality_(-1),
maxTreeDepth_(10),
surfaceClosed_(-1)
{
scalar scaleFactor = 0;
@ -445,13 +349,6 @@ Foam::triSurfaceMesh::triSurfaceMesh
bounds() = boundBox(points());
// Have optional non-standard search tolerance for gappy surfaces.
if (dict.readIfPresent("tolerance", tolerance_) && tolerance_ > 0)
{
Info<< searchableSurface::name() << " : using intersection tolerance "
<< tolerance_ << endl;
}
// Have optional minimum quality for normal calculation
if (dict.readIfPresent("minQuality", minQuality_) && minQuality_ > 0)
{
@ -459,13 +356,6 @@ Foam::triSurfaceMesh::triSurfaceMesh
<< " : ignoring triangles with quality < "
<< minQuality_ << " for normals calculation." << endl;
}
// Have optional non-standard tree-depth to limit storage.
if (dict.readIfPresent("maxTreeDepth", maxTreeDepth_) && maxTreeDepth_ > 0)
{
Info<< searchableSurface::name() << " : using maximum tree depth "
<< maxTreeDepth_ << endl;
}
}
@ -479,7 +369,7 @@ Foam::triSurfaceMesh::~triSurfaceMesh()
void Foam::triSurfaceMesh::clearOut()
{
tree_.clear();
triSurfaceRegionSearch::clearOut();
edgeTree_.clear();
triSurface::clearOut();
}
@ -521,64 +411,12 @@ bool Foam::triSurfaceMesh::overlaps(const boundBox& bb) const
void Foam::triSurfaceMesh::movePoints(const pointField& newPoints)
{
tree_.clear();
triSurfaceRegionSearch::clearOut();
edgeTree_.clear();
triSurface::movePoints(newPoints);
}
const Foam::indexedOctree<Foam::treeDataTriSurface>&
Foam::triSurfaceMesh::tree() const
{
if (tree_.empty())
{
// Calculate bb without constructing local point numbering.
treeBoundBox bb;
label nPoints;
calcBounds(bb, nPoints);
if (nPoints != points()().size())
{
WarningIn("triSurfaceMesh::tree() const")
<< "Surface " << searchableSurface::name()
<< " does not have compact point numbering."
<< " Of " << points()().size() << " only " << nPoints
<< " are used. This might give problems in some routines."
<< endl;
}
// Random number generator. Bit dodgy since not exactly random ;-)
Random rndGen(65431);
// Slightly extended bb. Slightly off-centred just so on symmetric
// geometry there are less face/edge aligned items.
bb = bb.extend(rndGen, 1e-4);
bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
tree_.reset
(
new indexedOctree<treeDataTriSurface>
(
treeDataTriSurface(*this, tolerance_),
bb,
maxTreeDepth_, // maxLevel
10, // leafsize
3.0 // duplicity
)
);
indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
}
return tree_();
}
const Foam::indexedOctree<Foam::treeDataEdge>&
Foam::triSurfaceMesh::edgeTree() const
{
@ -597,7 +435,12 @@ Foam::triSurfaceMesh::edgeTree() const
treeBoundBox bb;
label nPoints;
calcBounds(bb, nPoints);
PatchTools::calcBounds
(
static_cast<const triSurface&>(*this),
bb,
nPoints
);
// Random number generator. Bit dodgy since not exactly random ;-)
Random rndGen(65431);
@ -609,8 +452,8 @@ Foam::triSurfaceMesh::edgeTree() const
bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
indexedOctree<treeDataEdge>::perturbTol() = tolerance_;
scalar oldTol = indexedOctree<treeDataEdge>::perturbTol();
indexedOctree<treeDataEdge>::perturbTol() = tolerance();
edgeTree_.reset
(
@ -624,7 +467,7 @@ Foam::triSurfaceMesh::edgeTree() const
bEdges // selected edges
),
bb, // bb
maxTreeDepth_, // maxLevel
maxTreeDepth(), // maxLevel
10, // leafsize
3.0 // duplicity
)
@ -636,12 +479,6 @@ Foam::triSurfaceMesh::edgeTree() const
}
Foam::scalar Foam::triSurfaceMesh::tolerance() const
{
return tolerance_;
}
const Foam::wordList& Foam::triSurfaceMesh::regions() const
{
if (regions_.empty())
@ -682,23 +519,25 @@ void Foam::triSurfaceMesh::findNearest
List<pointIndexHit>& info
) const
{
const indexedOctree<treeDataTriSurface>& octree = tree();
triSurfaceSearch::findNearest(samples, nearestDistSqr, info);
}
info.setSize(samples.size());
scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
forAll(samples, i)
{
static_cast<pointIndexHit&>(info[i]) = octree.findNearest
(
samples[i],
nearestDistSqr[i]
);
}
indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
void Foam::triSurfaceMesh::findNearest
(
const pointField& samples,
const scalarField& nearestDistSqr,
const labelList& regionIndices,
List<pointIndexHit>& info
) const
{
triSurfaceRegionSearch::findNearest
(
samples,
nearestDistSqr,
regionIndices,
info
);
}
@ -709,23 +548,7 @@ void Foam::triSurfaceMesh::findLine
List<pointIndexHit>& info
) const
{
const indexedOctree<treeDataTriSurface>& octree = tree();
info.setSize(start.size());
scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
forAll(start, i)
{
static_cast<pointIndexHit&>(info[i]) = octree.findLine
(
start[i],
end[i]
);
}
indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
triSurfaceSearch::findLine(start, end, info);
}
@ -736,23 +559,7 @@ void Foam::triSurfaceMesh::findLineAny
List<pointIndexHit>& info
) const
{
const indexedOctree<treeDataTriSurface>& octree = tree();
info.setSize(start.size());
scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
forAll(start, i)
{
static_cast<pointIndexHit&>(info[i]) = octree.findLineAny
(
start[i],
end[i]
);
}
indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
triSurfaceSearch::findLineAny(start, end, info);
}
@ -763,57 +570,7 @@ void Foam::triSurfaceMesh::findLineAll
List<List<pointIndexHit> >& info
) const
{
const indexedOctree<treeDataTriSurface>& octree = tree();
info.setSize(start.size());
scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
// Work array
DynamicList<pointIndexHit, 1, 1> hits;
// Tolerances:
// To find all intersections we add a small vector to the last intersection
// This is chosen such that
// - it is significant (SMALL is smallest representative relative tolerance;
// we need something bigger since we're doing calculations)
// - if the start-end vector is zero we still progress
const vectorField dirVec(end-start);
const vectorField smallVec
(
indexedOctree<treeDataTriSurface>::perturbTol()*dirVec
+ vector(ROOTVSMALL,ROOTVSMALL,ROOTVSMALL)
);
forAll(start, pointI)
{
// See if any intersection between pt and end
pointIndexHit inter = octree.findLine(start[pointI], end[pointI]);
if (inter.hit())
{
hits.clear();
hits.append(inter);
getNextIntersections
(
octree,
start[pointI],
end[pointI],
smallVec[pointI],
hits
);
info[pointI].transfer(hits);
}
else
{
info[pointI].clear();
}
}
indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
triSurfaceSearch::findLineAll(start, end, info);
}
@ -975,7 +732,7 @@ void Foam::triSurfaceMesh::getVolumeType
volType.setSize(points.size());
scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
indexedOctree<treeDataTriSurface>::perturbTol() = tolerance();
forAll(points, pointI)
{
@ -983,7 +740,7 @@ void Foam::triSurfaceMesh::getVolumeType
// - use cached volume type per each tree node
// - cheat conversion since same values
volType[pointI] = static_cast<volumeType>(tree().getVolumeType(pt));
volType[pointI] = tree().getVolumeType(pt);
}
indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;

34
src/meshTools/searchableSurface/triSurfaceMesh.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -46,8 +46,11 @@ SourceFiles
#include "objectRegistry.H"
#include "indexedOctree.H"
#include "treeDataTriSurface.H"
#include "treeDataPrimitivePatch.H"
#include "treeDataEdge.H"
#include "EdgeMap.H"
#include "triSurface.H"
#include "triSurfaceRegionSearch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -62,25 +65,17 @@ class triSurfaceMesh
:
public searchableSurface,
public objectRegistry, // so we can store fields
public triSurface
public triSurface,
public triSurfaceRegionSearch
{
private:
// Private member data
//- Optional tolerance to use in searches
scalar tolerance_;
//- Optional min triangle quality. Triangles below this get
// ignored for normal calculation
scalar minQuality_;
//- Optional max tree depth of octree
label maxTreeDepth_;
//- Search tree (triangles)
mutable autoPtr<indexedOctree<treeDataTriSurface> > tree_;
//- Search tree for boundary edges.
mutable autoPtr<indexedOctree<treeDataEdge> > edgeTree_;
@ -136,11 +131,6 @@ private:
void operator=(const triSurfaceMesh&);
protected:
//- Calculate (number of)used points and their bounding box
void calcBounds(boundBox& bb, label& nPoints) const;
public:
//- Runtime type information
@ -176,13 +166,9 @@ public:
//- Move points
virtual void movePoints(const pointField&);
//- Demand driven contruction of octree
const indexedOctree<treeDataTriSurface>& tree() const;
//- Demand driven contruction of octree for boundary edges
const indexedOctree<treeDataEdge>& edgeTree() const;
scalar tolerance() const;
// searchableSurface implementation
@ -214,6 +200,14 @@ public:
List<pointIndexHit>&
) const;
virtual void findNearest
(
const pointField& sample,
const scalarField& nearestDistSqr,
const labelList& regionIndices,
List<pointIndexHit>&
) const;
virtual void findLine
(
const pointField& start,

3
src/meshTools/sets/pointSources/surfaceToPoint/surfaceToPoint.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -26,6 +26,7 @@ License
#include "surfaceToPoint.H"
#include "polyMesh.H"
#include "triSurfaceSearch.H"
#include "triSurface.H"
#include "cpuTime.H"
#include "addToRunTimeSelectionTable.H"

3
src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -694,6 +694,7 @@ void Foam::surfaceIntersection::doCutEdges
}
while (doTrack);
}
intersection::setPlanarTol(oldTol);
}

14
src/meshTools/triSurface/orientedSurface/orientedSurface.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -258,7 +258,8 @@ void Foam::orientedSurface::findZoneSide
zoneFaceI = -1;
isOutside = false;
List<pointIndexHit> hits;
pointField start(1, outsidePoint);
List<List<pointIndexHit> > hits(1, List<pointIndexHit>());
forAll(faceZone, faceI)
{
@ -273,7 +274,7 @@ void Foam::orientedSurface::findZoneSide
// Check if normal different enough to decide upon
if (magD > SMALL && (mag(n & d/magD) > 1e-6))
{
point end = fc + d;
pointField end(1, fc + d);
//Info<< "Zone " << zoneI << " : Shooting ray"
// << " from " << outsidePoint
@ -281,12 +282,13 @@ void Foam::orientedSurface::findZoneSide
// << " to pierce triangle " << faceI
// << " with centre " << fc << endl;
surfSearches.findLineAll(outsidePoint, end, hits);
surfSearches.findLineAll(start, end, hits);
label zoneIndex = -1;
forAll(hits, i)
forAll(hits[0], i)
{
if (hits[i].index() == faceI)
if (hits[0][i].index() == faceI)
{
zoneIndex = i;
break;

246
src/meshTools/triSurface/triSurfaceSearch/triSurfaceRegionSearch.C Normal file
View File

@ -0,0 +1,246 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ 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 3 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, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "triSurfaceRegionSearch.H"
#include "indexedOctree.H"
#include "triSurface.H"
#include "PatchTools.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::triSurfaceRegionSearch::triSurfaceRegionSearch(const triSurface& surface)
:
triSurfaceSearch(surface),
indirectRegionPatches_(),
treeByRegion_()
{}
Foam::triSurfaceRegionSearch::triSurfaceRegionSearch
(
const triSurface& surface,
const dictionary& dict
)
:
triSurfaceSearch(surface, dict),
indirectRegionPatches_(),
treeByRegion_()
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::triSurfaceRegionSearch::~triSurfaceRegionSearch()
{
clearOut();
}
void Foam::triSurfaceRegionSearch::clearOut()
{
triSurfaceSearch::clearOut();
treeByRegion_.clear();
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
const Foam::PtrList<Foam::triSurfaceRegionSearch::treeType>&
Foam::triSurfaceRegionSearch::treeByRegion() const
{
if (treeByRegion_.empty())
{
Map<label> regionSizes;
forAll(surface(), fI)
{
const label regionI = surface()[fI].region();
regionSizes(regionI)++;
}
label nRegions = regionSizes.size();
indirectRegionPatches_.setSize(nRegions);
treeByRegion_.setSize(nRegions);
labelListList regionsAddressing(nRegions);
forAll(regionsAddressing, regionI)
{
regionsAddressing[regionI] = labelList(regionSizes[regionI], -1);
}
labelList nFacesInRegions(nRegions, 0);
forAll(surface(), fI)
{
const label regionI = surface()[fI].region();
regionsAddressing[regionI][nFacesInRegions[regionI]++] = fI;
}
forAll(regionsAddressing, regionI)
{
scalar oldTol = treeType::perturbTol();
treeType::perturbTol() = tolerance();
indirectRegionPatches_.set
(
regionI,
new indirectTriSurface
(
IndirectList<labelledTri>
(
surface(),
regionsAddressing[regionI]
),
surface().points()
)
);
// Calculate bb without constructing local point numbering.
treeBoundBox bb;
label nPoints;
PatchTools::calcBounds
(
indirectRegionPatches_[regionI],
bb,
nPoints
);
// if (nPoints != surface().points().size())
// {
// WarningIn("triSurfaceRegionSearch::treeByRegion() const")
// << "Surface does not have compact point numbering. "
// << "Of " << surface().points().size()
// << " only " << nPoints
// << " are used. This might give problems in some routines."
// << endl;
// }
// Random number generator. Bit dodgy since not exactly random ;-)
Random rndGen(65431);
// Slightly extended bb. Slightly off-centred just so on symmetric
// geometry there are fewer face/edge aligned items.
bb = bb.extend(rndGen, 1e-4);
bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
treeByRegion_.set
(
regionI,
new treeType
(
treeDataIndirectTriSurface
(
true,
indirectRegionPatches_[regionI],
tolerance()
),
bb,
maxTreeDepth(), // maxLevel
10, // leafsize
3.0 // duplicity
)
);
treeType::perturbTol() = oldTol;
}
}
return treeByRegion_;
}
void Foam::triSurfaceRegionSearch::findNearest
(
const pointField& samples,
const scalarField& nearestDistSqr,
const labelList& regionIndices,
List<pointIndexHit>& info
) const
{
if (regionIndices.empty())
{
triSurfaceSearch::findNearest(samples, nearestDistSqr, info);
}
else
{
scalar oldTol = treeType::perturbTol();
treeType::perturbTol() = tolerance();
const PtrList<treeType>& octrees = treeByRegion();
info.setSize(samples.size());
forAll(octrees, treeI)
{
if (findIndex(regionIndices, treeI) == -1)
{
continue;
}
const treeType& octree = octrees[treeI];
forAll(samples, i)
{
// if (!octree.bb().contains(samples[i]))
// {
// continue;
// }
pointIndexHit currentRegionHit = octree.findNearest
(
samples[i],
nearestDistSqr[i],
treeDataIndirectTriSurface::findNearestOp(octree)
);
if
(
currentRegionHit.hit()
&&
(
!info[i].hit()
||
(
magSqr(currentRegionHit.hitPoint() - samples[i])
< magSqr(info[i].hitPoint() - samples[i])
)
)
)
{
info[i] = currentRegionHit;
}
}
}
treeType::perturbTol() = oldTol;
}
}
// ************************************************************************* //

144
src/meshTools/triSurface/triSurfaceSearch/triSurfaceRegionSearch.H Normal file
View File

@ -0,0 +1,144 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ 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 3 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, see <http://www.gnu.org/licenses/>.
Class
Foam::triSurfaceRegionSearch
Description
Helper class to search on triSurface. Creates an octree for each region of
the surface and only searches on the specified regions.
SourceFiles
triSurfaceRegionSearch.C
\*---------------------------------------------------------------------------*/
#ifndef triSurfaceRegionSearch_H
#define triSurfaceRegionSearch_H
#include "pointField.H"
#include "pointIndexHit.H"
#include "triSurfaceSearch.H"
#include "labelledTri.H"
#include "IndirectList.H"
#include "PtrList.H"
#include "indexedOctree.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class triSurfaceRegionSearch Declaration
\*---------------------------------------------------------------------------*/
class triSurfaceRegionSearch
:
public triSurfaceSearch
{
// Private typedefs
typedef PrimitivePatch
<
labelledTri,
IndirectList,
const pointField&
> indirectTriSurface;
typedef treeDataPrimitivePatch<indirectTriSurface>
treeDataIndirectTriSurface;
typedef indexedOctree<treeDataIndirectTriSurface> treeType;
// Private data
//- Surface is split into patches by region
mutable PtrList<indirectTriSurface> indirectRegionPatches_;
//- Search tree for each region
mutable PtrList<treeType> treeByRegion_;
// Private Member Functions
//- Disallow default bitwise copy construct
triSurfaceRegionSearch(const triSurfaceRegionSearch&);
//- Disallow default bitwise assignment
void operator=(const triSurfaceRegionSearch&);
public:
// Constructors
//- Construct from surface. Holds reference to surface!
explicit triSurfaceRegionSearch(const triSurface&);
//- Construct from surface and dictionary. Holds reference to surface!
triSurfaceRegionSearch(const triSurface&, const dictionary& dict);
//- Destructor
~triSurfaceRegionSearch();
//- Clear storage
void clearOut();
// Member Functions
// Access
//- Demand driven construction of octree for each region.
// @todo Currently creates a tree for each region; could optimise
// by only constructing trees when they are in regionIndices
const PtrList<treeType>& treeByRegion() const;
// Query
//- Find the nearest point on the surface out of the regions
// supplied in the list regionIndices. Ignores regions that are
// not specified
void findNearest
(
const pointField& samples,
const scalarField& nearestDistSqr,
const labelList& regionIndices,
List<pointIndexHit>& info
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

443
src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -24,62 +24,203 @@ License
\*---------------------------------------------------------------------------*/
#include "triSurfaceSearch.H"
#include "indexedOctree.H"
#include "boolList.H"
#include "treeDataTriSurface.H"
#include "triSurface.H"
#include "line.H"
#include "cpuTime.H"
#include "PatchTools.H"
#include "volumeType.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
const Foam::point Foam::triSurfaceSearch::greatPoint(GREAT, GREAT, GREAT);
bool Foam::triSurfaceSearch::checkUniqueHit
(
const pointIndexHit& currHit,
const DynamicList<pointIndexHit, 1, 1>& hits,
const vector& lineVec
) const
{
// Classify the hit
label nearType = -1;
label nearLabel = -1;
const labelledTri& f = surface()[currHit.index()];
f.nearestPointClassify
(
currHit.hitPoint(),
surface().points(),
nearType,
nearLabel
);
if (nearType == 1)
{
// near point
}
else if (nearType == 2)
{
// near edge
// check if the other face of the edge is already hit
const labelList& fEdges = surface().faceEdges()[currHit.index()];
const label edgeI = fEdges[nearLabel];
const labelList& edgeFaces = surface().edgeFaces()[edgeI];
forAll(edgeFaces, fI)
{
const label edgeFaceI = edgeFaces[fI];
if (edgeFaceI != currHit.index())
{
forAll(hits, hI)
{
const pointIndexHit& hit = hits[hI];
if (hit.index() == edgeFaceI)
{
// Check normals
const vector currHitNormal =
surface().faceNormals()[currHit.index()];
const vector existingHitNormal =
surface().faceNormals()[edgeFaceI];
const label signCurrHit =
pos(currHitNormal & lineVec);
const label signExistingHit =
pos(existingHitNormal & lineVec);
if (signCurrHit == signExistingHit)
{
return false;
}
}
}
}
}
}
return true;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from surface. Holds reference!
Foam::triSurfaceSearch::triSurfaceSearch(const triSurface& surface)
:
surface_(surface),
tolerance_(indexedOctree<treeDataTriSurface>::perturbTol()),
maxTreeDepth_(10),
treePtr_(NULL)
{}
Foam::triSurfaceSearch::triSurfaceSearch
(
const triSurface& surface,
const dictionary& dict
)
:
surface_(surface),
tolerance_(indexedOctree<treeDataTriSurface>::perturbTol()),
maxTreeDepth_(10),
treePtr_(NULL)
{
// Random number generator. Bit dodgy since not exactly random ;-)
Random rndGen(65431);
// Have optional non-standard search tolerance for gappy surfaces.
if (dict.readIfPresent("tolerance", tolerance_) && tolerance_ > 0)
{
Info<< " using intersection tolerance " << tolerance_ << endl;
}
// Slightly extended bb. Slightly off-centred just so on symmetric
// geometry there are less face/edge aligned items.
treeBoundBox treeBb
(
treeBoundBox(surface_.points(), surface_.meshPoints()).extend
(
rndGen,
1e-4
)
);
treeBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
treeBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
// Have optional non-standard tree-depth to limit storage.
if (dict.readIfPresent("maxTreeDepth", maxTreeDepth_) && maxTreeDepth_ > 0)
{
Info<< " using maximum tree depth " << maxTreeDepth_ << endl;
}
}
treePtr_.reset
(
new indexedOctree<treeDataTriSurface>
(
treeDataTriSurface
(
surface_,
indexedOctree<treeDataTriSurface>::perturbTol()
),
treeBb,
8, // maxLevel
10, // leafsize
3.0 // duplicity
)
);
Foam::triSurfaceSearch::triSurfaceSearch
(
const triSurface& surface,
const scalar tolerance,
const label maxTreeDepth
)
:
surface_(surface),
tolerance_(tolerance),
maxTreeDepth_(maxTreeDepth),
treePtr_(NULL)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::triSurfaceSearch::~triSurfaceSearch()
{
clearOut();
}
void Foam::triSurfaceSearch::clearOut()
{
treePtr_.clear();
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
const Foam::indexedOctree<Foam::treeDataTriSurface>&
Foam::triSurfaceSearch::tree() const
{
if (treePtr_.empty())
{
// Calculate bb without constructing local point numbering.
treeBoundBox bb;
label nPoints;
PatchTools::calcBounds(surface(), bb, nPoints);
if (nPoints != surface().points().size())
{
WarningIn("triSurfaceSearch::tree() const")
<< "Surface does not have compact point numbering."
<< " Of " << surface().points().size() << " only " << nPoints
<< " are used. This might give problems in some routines."
<< endl;
}
// Random number generator. Bit dodgy since not exactly random ;-)
Random rndGen(65431);
// Slightly extended bb. Slightly off-centred just so on symmetric
// geometry there are less face/edge aligned items.
bb = bb.extend(rndGen, 1e-4);
bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
treePtr_.reset
(
new indexedOctree<treeDataTriSurface>
(
treeDataTriSurface(true, surface_, tolerance_),
bb,
maxTreeDepth_, // maxLevel
10, // leafsize
3.0 // duplicity
)
);
indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
}
return treePtr_();
}
// Determine inside/outside for samples
Foam::boolList Foam::triSurfaceSearch::calcInside
(
@ -96,11 +237,7 @@ Foam::boolList Foam::triSurfaceSearch::calcInside
{
inside[sampleI] = false;
}
else if
(
tree().getVolumeType(sample)
== indexedOctree<treeDataTriSurface>::INSIDE
)
else if (tree().getVolumeType(sample) == volumeType::INSIDE)
{
inside[sampleI] = true;
}
@ -113,65 +250,31 @@ Foam::boolList Foam::triSurfaceSearch::calcInside
}
Foam::labelList Foam::triSurfaceSearch::calcNearestTri
void Foam::triSurfaceSearch::findNearest
(
const pointField& samples,
const vector& span
const scalarField& nearestDistSqr,
List<pointIndexHit>& info
) const
{
labelList nearest(samples.size());
scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
indexedOctree<treeDataTriSurface>::perturbTol() = tolerance();
const scalar nearestDistSqr = 0.25*magSqr(span);
const indexedOctree<treeDataTriSurface>& octree = tree();
pointIndexHit hitInfo;
info.setSize(samples.size());
forAll(samples, sampleI)
forAll(samples, i)
{
hitInfo = tree().findNearest(samples[sampleI], nearestDistSqr);
if (hitInfo.hit())
{
nearest[sampleI] = hitInfo.index();
}
else
{
nearest[sampleI] = -1;
}
static_cast<pointIndexHit&>(info[i]) = octree.findNearest
(
samples[i],
nearestDistSqr[i],
treeDataTriSurface::findNearestOp(octree)
);
}
return nearest;
}
// Nearest point on surface
Foam::tmp<Foam::pointField> Foam::triSurfaceSearch::calcNearest
(
const pointField& samples,
const vector& span
) const
{
const scalar nearestDistSqr = 0.25*magSqr(span);
tmp<pointField> tnearest(new pointField(samples.size()));
pointField& nearest = tnearest();
pointIndexHit hitInfo;
forAll(samples, sampleI)
{
hitInfo = tree().findNearest(samples[sampleI], nearestDistSqr);
if (hitInfo.hit())
{
nearest[sampleI] = hitInfo.hitPoint();
}
else
{
nearest[sampleI] = greatPoint;
}
}
return tnearest;
indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
}
@ -188,84 +291,126 @@ const
}
void Foam::triSurfaceSearch::findLine
(
const pointField& start,
const pointField& end,
List<pointIndexHit>& info
) const
{
const indexedOctree<treeDataTriSurface>& octree = tree();
info.setSize(start.size());
scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
indexedOctree<treeDataTriSurface>::perturbTol() = tolerance();
forAll(start, i)
{
static_cast<pointIndexHit&>(info[i]) = octree.findLine
(
start[i],
end[i]
);
}
indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
}
void Foam::triSurfaceSearch::findLineAny
(
const pointField& start,
const pointField& end,
List<pointIndexHit>& info
) const
{
const indexedOctree<treeDataTriSurface>& octree = tree();
info.setSize(start.size());
scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
indexedOctree<treeDataTriSurface>::perturbTol() = tolerance();
forAll(start, i)
{
static_cast<pointIndexHit&>(info[i]) = octree.findLineAny
(
start[i],
end[i]
);
}
indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
}
void Foam::triSurfaceSearch::findLineAll
(
const point& start,
const point& end,
List<pointIndexHit>& hits
)
const
const pointField& start,
const pointField& end,
List<List<pointIndexHit> >& info
) const
{
// See if any intersection between pt and end
pointIndexHit inter = tree().findLine(start, end);
const indexedOctree<treeDataTriSurface>& octree = tree();
if (inter.hit())
info.setSize(start.size());
scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
indexedOctree<treeDataTriSurface>::perturbTol() = tolerance();
// Work array
DynamicList<pointIndexHit, 1, 1> hits;
DynamicList<label> shapeMask;
treeDataTriSurface::findAllIntersectOp allIntersectOp(octree, shapeMask);
forAll(start, pointI)
{
hits.setSize(1);
hits[0] = inter;
const vector dirVec(end-start);
const scalar magSqrDirVec(magSqr(dirVec));
const vector smallVec
(
indexedOctree<treeDataTriSurface>::perturbTol()*dirVec
+ vector(ROOTVSMALL,ROOTVSMALL,ROOTVSMALL)
);
// Initial perturbation amount
vector perturbVec(smallVec);
hits.clear();
shapeMask.clear();
while (true)
{
// Start tracking from last hit.
point pt = hits.last().hitPoint() + perturbVec;
if (((pt-start)&dirVec) > magSqrDirVec)
{
return;
}
// See if any intersection between pt and end
pointIndexHit inter = tree().findLine(pt, end);
pointIndexHit inter = octree.findLine
(
start[pointI],
end[pointI],
allIntersectOp
);
if (!inter.hit())
if (inter.hit())
{
return;
}
vector lineVec = end[pointI] - start[pointI];
lineVec /= mag(lineVec) + VSMALL;
// Check if already found this intersection
bool duplicateHit = false;
forAllReverse(hits, i)
{
if (hits[i].index() == inter.index())
if
(
checkUniqueHit
(
inter,
hits,
lineVec
)
)
{
duplicateHit = true;
break;
hits.append(inter);
}
}
if (duplicateHit)
{
// Hit same triangle again. Increase perturbVec and try again.
perturbVec *= 2;
shapeMask.append(inter.index());
}
else
{
// Proper hit
label sz = hits.size();
hits.setSize(sz+1);
hits[sz] = inter;
// Restore perturbVec
perturbVec = smallVec;
break;
}
}
info[pointI].transfer(hits);
}
else
{
hits.clear();
}
indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
}

104
src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -48,11 +48,9 @@ namespace Foam
// Forward declaration of classes
class triSurface;
class treeDataTriSurface;
template<class Type> class indexedOctree;
/*---------------------------------------------------------------------------*\
Class triSurfaceSearch Declaration
Class triSurfaceSearch Declaration
\*---------------------------------------------------------------------------*/
class triSurfaceSearch
@ -62,12 +60,26 @@ class triSurfaceSearch
//- Reference to surface to work on
const triSurface& surface_;
//- Optional tolerance to use in searches
scalar tolerance_;
//- Optional max tree depth of octree
label maxTreeDepth_;
//- Octree for searches
autoPtr<indexedOctree<treeDataTriSurface> > treePtr_;
mutable autoPtr<indexedOctree<treeDataTriSurface> > treePtr_;
// Private Member Functions
//- Check whether the current hit on the surface which lies on lineVec
// is unique.
bool checkUniqueHit
(
const pointIndexHit& currHit,
const DynamicList<pointIndexHit, 1, 1>& hits,
const vector& lineVec
) const;
//- Disallow default bitwise copy construct
triSurfaceSearch(const triSurfaceSearch&);
@ -75,51 +87,64 @@ class triSurfaceSearch
//- Disallow default bitwise assignment
void operator=(const triSurfaceSearch&);
public:
// Static data members
//- Point far away; used for illegal finds
static const point greatPoint;
// Constructors
//- Construct from surface. Holds reference to surface!
triSurfaceSearch(const triSurface&);
explicit triSurfaceSearch(const triSurface&);
//- Construct from surface and dictionary.
triSurfaceSearch(const triSurface&, const dictionary& dict);
//- Construct from components
triSurfaceSearch
(
const triSurface& surface,
const scalar tolerance,
const label maxTreeDepth
);
//- Destructor
~triSurfaceSearch();
//- Clear storage
void clearOut();
// Member Functions
const indexedOctree<treeDataTriSurface>& tree() const
{
return treePtr_();
}
//- Demand driven construction of the octree
const indexedOctree<treeDataTriSurface>& tree() const;
//- Return reference to the surface.
const triSurface& surface() const
{
return surface_;
}
//- Return tolerance to use in searches
scalar tolerance() const
{
return tolerance_;
}
//- Return max tree depth of octree
label maxTreeDepth() const
{
return maxTreeDepth_;
}
//- Calculate for each searchPoint inside/outside status.
boolList calcInside(const pointField& searchPoints) const;
//- Calculate index of nearest triangle (or -1) for each sample.
// Looks only in box of size 2*span around sample.
labelList calcNearestTri
void findNearest
(
const pointField& samples,
const vector& span
) const;
//- Calculate nearest points (to searchPoints) on surface.
// Looks only in box of size 2*span around sample. Returns greatPoint
// if not found.
tmp<pointField> calcNearest
(
const pointField& samples,
const vector& span
const scalarField& nearestDistSqr,
List<pointIndexHit>& info
) const;
//- Calculate nearest point on surface for single searchPoint. Returns
@ -129,14 +154,27 @@ public:
// - index() : surface triangle label
pointIndexHit nearest(const point&, const vector& span) const;
void findLine
(
const pointField& start,
const pointField& end,
List<pointIndexHit>& info
) const;
void findLineAny
(
const pointField& start,
const pointField& end,
List<pointIndexHit>& info
) const;
//- Calculate all intersections from start to end
void findLineAll
(
const point& start,
const point& end,
List<pointIndexHit>&
const pointField& start,
const pointField& end,
List<List<pointIndexHit> >& info
) const;
};

4
src/meshTools/twoDPointCorrector/twoDPointCorrector.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -31,7 +31,7 @@ Description
and thus introduce a third dimension into a 2-D problem.
The operation is performed by looping through all edges approximately
normal to the plane and enforcing their orthoginality onto the plane by
normal to the plane and enforcing their orthogonality onto the plane by
adjusting points on their ends.
SourceFiles