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 branch 'master' of ssh://noisy/home/noisy3/OpenFOAM/OpenFOAM-dev

Browse Source
This commit is contained in:
andy
2010-10-07 11:22:38 +01:00
parent 590e2c3a52 95a2b3effd
commit 1b7ad868bb
70 changed files with 1878 additions and 20751 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
src/OpenFOAM/Make/files
View File

@ -296,6 +296,7 @@ primitiveShapes = meshes/primitiveShapes
$(primitiveShapes)/line/line.C
$(primitiveShapes)/plane/plane.C
$(primitiveShapes)/triangle/intersection.C
$(primitiveShapes)/objectHit/pointIndexHitIOList.C
meshShapes = meshes/meshShapes
$(meshShapes)/edge/edge.C

120
src/OpenFOAM/containers/Lists/IndirectList/IndirectList.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2010-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -25,7 +25,11 @@ Class
Foam::IndirectList
Description
A List with indirect addressing
A List with indirect addressing.
SeeAlso
Foam::UIndirectList for a version without any allocation for the
addressing.
SourceFiles
IndirectListI.H
@ -35,24 +39,79 @@ SourceFiles
#ifndef IndirectList_H
#define IndirectList_H
#include "List.H"
#include "UIndirectList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class IndirectListAddressing Declaration
\*---------------------------------------------------------------------------*/
//- A helper class for storing addresses.
class IndirectListAddressing
{
// Private data
//- Storage for the list addressing
List<label> addressing_;
// Private Member Functions
//- Disallow default bitwise copy construct
IndirectListAddressing(const IndirectListAddressing&);
//- Disallow default bitwise assignment
void operator=(const IndirectListAddressing&);
protected:
// Constructors
//- Construct by copying the addressing array
explicit inline IndirectListAddressing(const UList<label>& addr);
//- Construct by transferring addressing array
explicit inline IndirectListAddressing(const Xfer<List<label> >& addr);
// Member Functions
// Access
//- Return the list addressing
inline const List<label>& addressing() const;
// Edit
//- Reset addressing
inline void resetAddressing(const UList<label>&);
inline void resetAddressing(const Xfer<List<label> >&);
};
/*---------------------------------------------------------------------------*\
Class IndirectList Declaration
\*---------------------------------------------------------------------------*/
template<class T>
class IndirectList
:
private IndirectListAddressing,
public UIndirectList<T>
{
// Private data
// Private Member Functions
UList<T>& completeList_;
List<label> addressing_;
//- Disable default assignment operator
void operator=(const IndirectList<T>&);
//- Disable assignment from UIndirectList
void operator=(const UIndirectList<T>&);
public:
@ -65,59 +124,32 @@ public:
//- Construct given the complete list and by transferring addressing
inline IndirectList(const UList<T>&, const Xfer<List<label> >&);
//- Copy constructor
inline IndirectList(const IndirectList<T>&);
//- Construct from UIndirectList
explicit inline IndirectList(const UIndirectList<T>&);
// Member Functions
// Access
//- Return the number of elements in the list
inline label size() const;
//- Return true if the list is empty (ie, size() is zero).
inline bool empty() const;
//- Return the first element of the list.
inline T& first();
//- Return first element of the list.
inline const T& first() const;
//- Return the last element of the list.
inline T& last();
//- Return the last element of the list.
inline const T& last() const;
//- Return the complete list
inline const UList<T>& completeList() const;
//- Return the list addressing
inline const List<label>& addressing() const;
using UIndirectList<T>::addressing;
// Edit
//- Reset addressing
inline void resetAddressing(const UList<label>&);
inline void resetAddressing(const Xfer<List<label> >&);
using IndirectListAddressing::resetAddressing;
// Member Operators
//- Return the addressed elements as a List
inline List<T> operator()() const;
//- Return non-const access to an element
inline T& operator[](const label);
//- Return const access to an element
inline const T& operator[](const label) const;
//- Assignment from UList of addressed elements
inline void operator=(const UList<T>&);
//- Assignment of all entries to the given value
inline void operator=(const T&);
//- Assignment operator
using UIndirectList<T>::operator=;
};

185
src/OpenFOAM/containers/Lists/IndirectList/IndirectListI.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2010-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -25,6 +25,25 @@ License
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
inline Foam::IndirectListAddressing::IndirectListAddressing
(
const UList<label>& addr
)
:
addressing_(addr)
{}
inline Foam::IndirectListAddressing::IndirectListAddressing
(
const Xfer<List<label> >& addr
)
:
addressing_(addr)
{}
template<class T>
inline Foam::IndirectList<T>::IndirectList
(
@ -32,8 +51,12 @@ inline Foam::IndirectList<T>::IndirectList
const UList<label>& addr
)
:
completeList_(const_cast<UList<T>&>(completeList)),
addressing_(addr)
IndirectListAddressing(addr),
UIndirectList<T>
(
completeList,
IndirectListAddressing::addressing()
)
{}
@ -44,71 +67,55 @@ inline Foam::IndirectList<T>::IndirectList
const Xfer<List<label> >& addr
)
:
completeList_(const_cast<UList<T>&>(completeList)),
addressing_(addr)
IndirectListAddressing(addr),
UIndirectList<T>
(
completeList,
IndirectListAddressing::addressing()
)
{}
template<class T>
inline Foam::IndirectList<T>::IndirectList
(
const IndirectList<T>& lst
)
:
IndirectListAddressing(lst.addressing()),
UIndirectList<T>
(
lst.completeList(),
IndirectListAddressing::addressing()
)
{}
template<class T>
inline Foam::IndirectList<T>::IndirectList
(
const UIndirectList<T>& lst
)
:
IndirectListAddressing(lst.addressing()),
UIndirectList<T>
(
lst.completeList(),
IndirectListAddressing::addressing()
)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class T>
inline Foam::label Foam::IndirectList<T>::size() const
{
return addressing_.size();
}
template<class T>
inline bool Foam::IndirectList<T>::empty() const
{
return addressing_.empty();
}
template<class T>
inline T& Foam::IndirectList<T>::first()
{
return completeList_[addressing_.first()];
}
template<class T>
inline const T& Foam::IndirectList<T>::first() const
{
return completeList_[addressing_.first()];
}
template<class T>
inline T& Foam::IndirectList<T>::last()
{
return completeList_[addressing_.last()];
}
template<class T>
inline const T& Foam::IndirectList<T>::last() const
{
return completeList_[addressing_.last()];
}
template<class T>
inline const Foam::UList<T>& Foam::IndirectList<T>::completeList() const
{
return completeList_;
}
template<class T>
inline const Foam::List<Foam::label>& Foam::IndirectList<T>::addressing() const
inline const Foam::List<Foam::label>&
Foam::IndirectListAddressing::addressing() const
{
return addressing_;
}
template<class T>
inline void Foam::IndirectList<T>::resetAddressing
inline void Foam::IndirectListAddressing::resetAddressing
(
const UList<label>& addr
)
@ -117,8 +124,7 @@ inline void Foam::IndirectList<T>::resetAddressing
}
template<class T>
inline void Foam::IndirectList<T>::resetAddressing
inline void Foam::IndirectListAddressing::resetAddressing
(
const Xfer<List<label> >& addr
)
@ -127,63 +133,4 @@ inline void Foam::IndirectList<T>::resetAddressing
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
template<class T>
inline Foam::List<T> Foam::IndirectList<T>::operator()() const
{
List<T> result(size());
forAll(*this, i)
{
result[i] = operator[](i);
}
return result;
}
template<class T>
inline T& Foam::IndirectList<T>::operator[](const label i)
{
return completeList_[addressing_[i]];
}
template<class T>
inline const T& Foam::IndirectList<T>::operator[](const label i) const
{
return completeList_[addressing_[i]];
}
template<class T>
inline void Foam::IndirectList<T>::operator=(const UList<T>& ae)
{
if (addressing_.size() != ae.size())
{
FatalErrorIn("IndirectList<T>::operator=(const UList<T>&)")
<< "Addressing and list of addressed elements "
"have different sizes: "
<< addressing_.size() << " " << ae.size()
<< abort(FatalError);
}
forAll(addressing_, i)
{
completeList_[addressing_[i]] = ae[i];
}
}
template<class T>
inline void Foam::IndirectList<T>::operator=(const T& t)
{
forAll(addressing_, i)
{
completeList_[addressing_[i]] = t;
}
}
// ************************************************************************* //

37
src/OpenFOAM/containers/Lists/List/List.C
View File

@ -266,24 +266,6 @@ Foam::List<T>::List(const SLList<T>& lst)
}
// Construct as copy of IndirectList<T>
template<class T>
Foam::List<T>::List(const IndirectList<T>& lst)
:
UList<T>(NULL, lst.size())
{
if (this->size_)
{
this->v_ = new T[this->size_];
forAll(*this, i)
{
this->operator[](i) = lst[i];
}
}
}
// Construct as copy of UIndirectList<T>
template<class T>
Foam::List<T>::List(const UIndirectList<T>& lst)
@ -517,25 +499,6 @@ void Foam::List<T>::operator=(const SLList<T>& lst)
}
// Assignment operator. Takes linear time.
template<class T>
void Foam::List<T>::operator=(const IndirectList<T>& lst)
{
if (lst.size() != this->size_)
{
if (this->v_) delete[] this->v_;
this->v_ = 0;
this->size_ = lst.size();
if (this->size_) this->v_ = new T[this->size_];
}
forAll(*this, i)
{
this->operator[](i) = lst[i];
}
}
// Assignment operator. Takes linear time.
template<class T>
void Foam::List<T>::operator=(const UIndirectList<T>& lst)

6
src/OpenFOAM/containers/Lists/List/List.H
View File

@ -131,9 +131,6 @@ public:
//- Construct as copy of SLList<T>
explicit List(const SLList<T>&);
//- Construct as copy of IndirectList<T>
explicit List(const IndirectList<T>&);
//- Construct as copy of UIndirectList<T>
explicit List(const UIndirectList<T>&);
@ -219,9 +216,6 @@ public:
//- Assignment from SLList operator. Takes linear time.
void operator=(const SLList<T>&);
//- Assignment from IndirectList operator. Takes linear time.
void operator=(const IndirectList<T>&);
//- Assignment from UIndirectList operator. Takes linear time.
void operator=(const UIndirectList<T>&);

22
src/OpenFOAM/meshes/meshShapes/face/face.H
View File

@ -126,6 +126,14 @@ class face
public:
//- Return types for classify
enum proxType
{
NONE,
POINT, // Close to point
EDGE // Close to edge
};
// Static data members
static const char* const typeName;
@ -249,6 +257,20 @@ public:
const pointField& meshPoints
) const;
//- Return nearest point to face and classify it:
// + near point (nearType=POINT, nearLabel=0, 1, 2)
// + near edge (nearType=EDGE, nearLabel=0, 1, 2)
// Note: edges are counted from starting vertex so
// e.g. edge n is from f[n] to f[0], where the face has n + 1
// points
pointHit nearestPointClassify
(
const point& p,
const pointField& meshPoints,
label& nearType,
label& nearLabel
) const;
//- Return contact sphere diameter
scalar contactSphereDiameter
(

55
src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C
View File

@ -181,6 +181,22 @@ Foam::pointHit Foam::face::nearestPoint
const point& p,
const pointField& meshPoints
) const
{
// Dummy labels
label nearType = -1;
label nearLabel = -1;
return nearestPointClassify(p, meshPoints, nearType, nearLabel);
}
Foam::pointHit Foam::face::nearestPointClassify
(
const point& p,
const pointField& meshPoints,
label& nearType,
label& nearLabel
) const
{
const face& f = *this;
point ctr = centre(meshPoints);
@ -188,6 +204,9 @@ Foam::pointHit Foam::face::nearestPoint
// Initialize to miss, distance=GREAT
pointHit nearest(p);
nearType = -1;
nearLabel = -1;
label nPoints = f.size();
point nextPoint = ctr;
@ -196,8 +215,10 @@ Foam::pointHit Foam::face::nearestPoint
{
nextPoint = meshPoints[f[fcIndex(pI)]];
label tmpNearType = -1;
label tmpNearLabel = -1;
// Note: for best accuracy, centre point always comes last
//
triPointRef tri
(
meshPoints[f[pI]],
@ -205,12 +226,42 @@ Foam::pointHit Foam::face::nearestPoint
ctr
);
pointHit curHit = tri.nearestPoint(p);
pointHit curHit = tri.nearestPointClassify
(
p,
tmpNearType,
tmpNearLabel
);
if (Foam::mag(curHit.distance()) < Foam::mag(nearest.distance()))
{
nearest.setDistance(curHit.distance());
// Assume at first that the near type is NONE on the
// triangle (i.e. on the face of the triangle) then it is
// therefore also for the face.
nearType = NONE;
if (tmpNearType == triPointRef::EDGE && tmpNearLabel == 0)
{
// If the triangle edge label is 0, then this is also
// an edge of the face, if not, it is on the face
nearType = EDGE;
nearLabel = pI;
}
else if (tmpNearType == triPointRef::POINT && tmpNearLabel < 2)
{
// If the triangle point label is 0 or 1, then this is
// also a point of the face, if not, it is on the face
nearType = POINT;
nearLabel = pI + tmpNearLabel;
}
if (curHit.hit())
{
nearest.setHit();

0
src/meshTools/octree/PointIndexHit.H → src/OpenFOAM/meshes/primitiveShapes/objectHit/PointIndexHit.H
View File

0
src/meshTools/octree/pointHitSort.H → src/OpenFOAM/meshes/primitiveShapes/objectHit/pointHitSort.H
View File

0
src/meshTools/octree/pointIndexHit.H → src/OpenFOAM/meshes/primitiveShapes/objectHit/pointIndexHit.H
View File

0
src/meshTools/octree/pointIndexHitIOList.C → src/OpenFOAM/meshes/primitiveShapes/objectHit/pointIndexHitIOList.C
View File

0
src/meshTools/octree/pointIndexHitIOList.H → src/OpenFOAM/meshes/primitiveShapes/objectHit/pointIndexHitIOList.H
View File

44
src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H
View File

@ -79,21 +79,6 @@ class triangle
PointRef a_, b_, c_;
// Private Member Functions
//- Fast distance to triangle calculation. From
// "Distance Between Point and Trangle in 3D"
// David Eberly, Magic Software Inc. Aug. 2002.
// Works on function Q giving distance to point and tries to
// minimize this.
static pointHit nearestPoint
(
const Point& baseVertex,
const vector& E0,
const vector& E1,
const point& P
);
public:
@ -202,24 +187,27 @@ public:
const scalar tol = 0.0
) const;
//- Return nearest point to p on triangle
inline pointHit nearestPoint
(
const point& p
) const;
//- Classify point in triangle plane w.r.t. triangle edges.
// - inside (true returned)/outside (false returned)
// - near point (nearType=POINT, nearLabel=0, 1, 2)
// - near edge (nearType=EDGE, nearLabel=0, 1, 2)
//- Find the nearest point to p on the triangle and classify it:
// + near point (nearType=POINT, nearLabel=0, 1, 2)
// + near edge (nearType=EDGE, nearLabel=0, 1, 2)
// Note: edges are counted from starting
// vertex so e.g. edge 2 is from f[2] to f[0]
// tol is fraction to account for truncation error. Is only used
// when comparing normalized (0..1) numbers.
pointHit nearestPointClassify
(
const point& p,
label& nearType,
label& nearLabel
) const;
//- Return nearest point to p on triangle
inline pointHit nearestPoint(const point& p) const;
//- Classify nearest point to p in triangle plane
// w.r.t. triangle edges and points. Returns inside
// (true)/outside (false).
bool classify
(
const point& p,
const scalar tol,
label& nearType,
label& nearLabel
) const;

426
src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H
View File

@ -32,158 +32,6 @@ License
namespace Foam
{
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class Point, class PointRef>
pointHit triangle<Point, PointRef>::nearestPoint
(
const Point& baseVertex,
const vector& E0,
const vector& E1,
const point& P
)
{
// Distance vector
const vector D(baseVertex - P);
// Some geometrical factors
const scalar a = E0 & E0;
const scalar b = E0 & E1;
const scalar c = E1 & E1;
// Precalculate distance factors
const scalar d = E0 & D;
const scalar e = E1 & D;
const scalar f = D & D;
// Do classification
const scalar det = a*c - b*b;
scalar s = b*e - c*d;
scalar t = b*d - a*e;
bool inside = false;
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;
inside = true;
}
}
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: Foam::mag used since truncation error causes negative distances
// with points very close to one of the triangle vertices.
// (Up to -2.77556e-14 seen). Could use +SMALL but that not large enough.
return pointHit
(
inside,
baseVertex + s*E0 + t*E1,
Foam::sqrt
(
Foam::mag(a*s*s + 2*b*s*t + c*t*t + 2*d*s + 2*e*t + f)
),
!inside
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Point, class PointRef>
@ -247,7 +95,7 @@ inline Point triangle<Point, PointRef>::centre() const
template<class Point, class PointRef>
inline scalar triangle<Point, PointRef>::mag() const
{
return ::Foam::mag(normal());
return Foam::mag(normal());
}
@ -536,7 +384,7 @@ inline pointHit triangle<Point, PointRef>::ray
inter.setMiss(eligible);
// The miss point is the nearest point on the triangle
inter.setPoint(nearestPoint(a_, E0, E1, p).rawPoint());
inter.setPoint(nearestPoint(p).rawPoint());
// The distance to the miss is the distance between the
// original point and plane of intersection
@ -633,18 +481,130 @@ inline pointHit triangle<Point, PointRef>::intersection
}
template<class Point, class PointRef>
pointHit triangle<Point, PointRef>::nearestPointClassify
(
const point& p,
label& nearType,
label& nearLabel
) const
{
// Adapted from:
// Real-time collision detection, Christer Ericson, 2005, 136-142
// Check if P in vertex region outside A
vector ab = b_ - a_;
vector ac = c_ - a_;
vector ap = p - a_;
scalar d1 = ab & ap;
scalar d2 = ac & ap;
if (d1 <= 0.0 && d2 <= 0.0)
{
// barycentric coordinates (1, 0, 0)
nearType = POINT;
nearLabel = 0;
return pointHit(false, a_, Foam::mag(a_ - p), true);
}
// Check if P in vertex region outside B
vector bp = p - b_;
scalar d3 = ab & bp;
scalar d4 = ac & bp;
if (d3 >= 0.0 && d4 <= d3)
{
// barycentric coordinates (0, 1, 0)
nearType = POINT;
nearLabel = 1;
return pointHit(false, b_, Foam::mag(b_ - p), true);
}
// Check if P in edge region of AB, if so return projection of P onto AB
scalar vc = d1*d4 - d3*d2;
if (vc <= 0.0 && d1 >= 0.0 && d3 <= 0.0)
{
// barycentric coordinates (1-v, v, 0)
scalar v = d1/(d1 - d3);
point nearPt = a_ + v*ab;
nearType = EDGE;
nearLabel = 0;
return pointHit(false, nearPt, Foam::mag(nearPt - p), true);
}
// Check if P in vertex region outside C
vector cp = p - c_;
scalar d5 = ab & cp;
scalar d6 = ac & cp;
if (d6 >= 0.0 && d5 <= d6)
{
// barycentric coordinates (0, 0, 1)
nearType = POINT;
nearLabel = 2;
return pointHit(false, c_, Foam::mag(c_ - p), true);
}
// Check if P in edge region of AC, if so return projection of P onto AC
scalar vb = d5*d2 - d1*d6;
if (vb <= 0.0 && d2 >= 0.0 && d6 <= 0.0)
{
// barycentric coordinates (1-w, 0, w)
scalar w = d2/(d2 - d6);
point nearPt = a_ + w*ac;
nearType = EDGE;
nearLabel = 2;
return pointHit(false, nearPt, Foam::mag(nearPt - p), true);
}
// Check if P in edge region of BC, if so return projection of P onto BC
scalar va = d3*d6 - d5*d4;
if (va <= 0.0 && (d4 - d3) >= 0.0 && (d5 - d6) >= 0.0)
{
// barycentric coordinates (0, 1-w, w)
scalar w = (d4 - d3)/((d4 - d3) + (d5 - d6));
point nearPt = b_ + w*(c_ - b_);
nearType = EDGE;
nearLabel = 1;
return pointHit(false, nearPt, Foam::mag(nearPt - p), true);
}
// P inside face region. Compute Q through its barycentric
// coordinates (u, v, w)
scalar denom = 1.0/(va + vb + vc);
scalar v = vb * denom;
scalar w = vc * denom;
// = u*a + v*b + w*c, u = va*denom = 1.0 - v - w
point nearPt = a_ + ab*v + ac*w;
nearType = NONE,
nearLabel = -1;
return pointHit(true, nearPt, Foam::mag(nearPt - p), false);
}
template<class Point, class PointRef>
inline pointHit triangle<Point, PointRef>::nearestPoint
(
const point& p
) const
{
// Express triangle in terms of baseVertex (point a_) and
// two edge vectors
vector E0 = b_ - a_;
vector E1 = c_ - a_;
// Dummy labels
label nearType = -1;
label nearLabel = -1;
return nearestPoint(a_, E0, E1, p);
return nearestPointClassify(p, nearType, nearLabel);
}
@ -652,160 +612,14 @@ template<class Point, class PointRef>
inline bool triangle<Point, PointRef>::classify
(
const point& p,
const scalar tol,
label& nearType,
label& nearLabel
) const
{
const vector E0 = b_ - a_;
const vector E1 = c_ - a_;
const vector n = 0.5*(E0 ^ E1);
// Get largest component of normal
scalar magX = Foam::mag(n.x());
scalar magY = Foam::mag(n.y());
scalar magZ = Foam::mag(n.z());
label i0 = -1;
if ((magX >= magY) && (magX >= magZ))
{
i0 = 0;
}
else if ((magY >= magX) && (magY >= magZ))
{
i0 = 1;
}
else
{
i0 = 2;
}
// Get other components
label i1 = (i0 + 1) % 3;
label i2 = (i1 + 1) % 3;
scalar u1 = E0[i1];
scalar v1 = E0[i2];
scalar u2 = E1[i1];
scalar v2 = E1[i2];
scalar det = v2*u1 - u2*v1;
scalar u0 = p[i1] - a_[i1];
scalar v0 = p[i2] - a_[i2];
scalar alpha = 0;
scalar beta = 0;
bool hit = false;
if (Foam::mag(u1) < ROOTVSMALL)
{
beta = u0/u2;
alpha = (v0 - beta*v2)/v1;
hit = ((alpha >= 0) && ((alpha + beta) <= 1));
}
else
{
beta = (v0*u1 - u0*v1)/det;
alpha = (u0 - beta*u2)/u1;
hit = ((alpha >= 0) && ((alpha + beta) <= 1));
}
//
// Now alpha, beta are the coordinates in the triangle local coordinate
// system E0, E1
//
//Pout<< "alpha:" << alpha << endl;
//Pout<< "beta:" << beta << endl;
//Pout<< "hit:" << hit << endl;
//Pout<< "tol:" << tol << endl;
if (hit)
{
// alpha,beta might get negative due to precision errors
alpha = max(0.0, min(1.0, alpha));
beta = max(0.0, min(1.0, beta));
}
nearType = NONE;
nearLabel = -1;
if (Foam::mag(alpha+beta-1) <= tol)
{
// On edge between vert 1-2 (=edge 1)
if (Foam::mag(alpha) <= tol)
{
nearType = POINT;
nearLabel = 2;
}
else if (Foam::mag(beta) <= tol)
{
nearType = POINT;
nearLabel = 1;
}
else if ((alpha >= 0) && (alpha <= 1) && (beta >= 0) && (beta <= 1))
{
nearType = EDGE;
nearLabel = 1;
}
}
else if (Foam::mag(alpha) <= tol)
{
// On edge between vert 2-0 (=edge 2)
if (Foam::mag(beta) <= tol)
{
nearType = POINT;
nearLabel = 0;
}
else if (Foam::mag(beta-1) <= tol)
{
nearType = POINT;
nearLabel = 2;
}
else if ((beta >= 0) && (beta <= 1))
{
nearType = EDGE;
nearLabel = 2;
}
}
else if (Foam::mag(beta) <= tol)
{
// On edge between vert 0-1 (= edge 0)
if (Foam::mag(alpha) <= tol)
{
nearType = POINT;
nearLabel = 0;
}
else if (Foam::mag(alpha-1) <= tol)
{
nearType = POINT;
nearLabel = 1;
}
else if ((alpha >= 0) && (alpha <= 1))
{
nearType = EDGE;
nearLabel = 0;
}
}
return hit;
return nearestPointClassify(p, nearType, nearLabel).hit();
}
// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
template<class point, class pointRef>

7
src/finiteVolume/fields/fvPatchFields/derived/directMappedVelocityFluxFixedValue/directMappedVelocityFluxFixedValueFvPatchField.C
View File

@ -112,13 +112,6 @@ directMappedVelocityFluxFixedValueFvPatchField
<< " in file " << dimensionedInternalField().objectPath()
<< exit(FatalError);
}
// Force calculation of schedule (uses parallel comms)
const directMappedPolyPatch& mpp = refCast<const directMappedPolyPatch>
(
this->patch().patch()
);
(void)mpp.map().schedule();
}

34
src/finiteVolume/fvMesh/singleCellFvMesh/singleCellFvMesh.C
View File

@ -42,24 +42,26 @@ void Foam::singleCellFvMesh::agglomerateMesh
const polyBoundaryMesh& oldPatches = mesh.boundaryMesh();
// Check agglomeration within patch face range and continuous
labelList nAgglom(oldPatches.size());
labelList nAgglom(oldPatches.size(), 0);
forAll(oldPatches, patchI)
{
const polyPatch& pp = oldPatches[patchI];
nAgglom[patchI] = max(agglom[patchI])+1;
forAll(pp, i)
if (pp.size() > 0)
{
if (agglom[patchI][i] < 0 || agglom[patchI][i] >= pp.size())
nAgglom[patchI] = max(agglom[patchI])+1;
forAll(pp, i)
{
FatalErrorIn
(
"singleCellFvMesh::agglomerateMesh(..)"
) << "agglomeration on patch " << patchI
<< " is out of range 0.." << pp.size()-1
<< exit(FatalError);
if (agglom[patchI][i] < 0 || agglom[patchI][i] >= pp.size())
{
FatalErrorIn
(
"singleCellFvMesh::agglomerateMesh(..)"
) << "agglomeration on patch " << patchI
<< " is out of range 0.." << pp.size()-1
<< exit(FatalError);
}
}
}
}
@ -155,6 +157,8 @@ void Foam::singleCellFvMesh::agglomerateMesh
forAll(oldPatches, patchI)
{
patchStarts[patchI] = coarseI;
const polyPatch& pp = oldPatches[patchI];
if (pp.size() > 0)
@ -170,8 +174,6 @@ void Foam::singleCellFvMesh::agglomerateMesh
// From agglomeration to compact patch face
labelList agglomToFace(nAgglom[patchI], -1);
patchStarts[patchI] = coarseI;
forAll(pp, i)
{
label myAgglom = agglom[patchI][i];
@ -223,9 +225,9 @@ void Foam::singleCellFvMesh::agglomerateMesh
);
}
}
patchSizes[patchI] = coarseI-patchStarts[patchI];
}
patchSizes[patchI] = coarseI-patchStarts[patchI];
}
//Pout<< "patchStarts:" << patchStarts << endl;

1
src/meshTools/Make/files
View File

@ -44,7 +44,6 @@ octree/octreeDataFace.C
octree/treeBoundBox.C
octree/treeNodeName.C
octree/treeLeafName.C
octree/pointIndexHitIOList.C
indexedOctree/indexedOctreeName.C
indexedOctree/treeDataCell.C

198
src/meshTools/indexedOctree/indexedOctree.C
View File

@ -2056,6 +2056,180 @@ void Foam::indexedOctree<Type>::findBox
}
template <class Type>
template <class CompareOp>
void Foam::indexedOctree<Type>::findNear
(
const scalar nearDist,
const bool okOrder,
const indexedOctree<Type>& tree1,
const labelBits index1,
const treeBoundBox& bb1,
const indexedOctree<Type>& tree2,
const labelBits index2,
const treeBoundBox& bb2,
CompareOp& cop
)
{
const vector nearSpan(nearDist, nearDist, nearDist);
if (tree1.isNode(index1))
{
const node& nod1 = tree1.nodes()[tree1.getNode(index1)];
const treeBoundBox searchBox
(
bb1.min()-nearSpan,
bb1.max()+nearSpan
);
if (tree2.isNode(index2))
{
if (bb2.overlaps(searchBox))
{
const node& nod2 = tree2.nodes()[tree2.getNode(index2)];
for (direction i2 = 0; i2 < nod2.subNodes_.size(); i2++)
{
labelBits subIndex2 = nod2.subNodes_[i2];
const treeBoundBox subBb2
(
tree2.isNode(subIndex2)
? tree2.nodes()[tree2.getNode(subIndex2)].bb_
: bb2.subBbox(i2)
);
findNear
(
nearDist,
!okOrder,
tree2,
subIndex2,
subBb2,
tree1,
index1,
bb1,
cop
);
}
}
}
else if (tree2.isContent(index2))
{
// index2 is leaf, index1 not yet.
for (direction i1 = 0; i1 < nod1.subNodes_.size(); i1++)
{
labelBits subIndex1 = nod1.subNodes_[i1];
const treeBoundBox subBb1
(
tree1.isNode(subIndex1)
? tree1.nodes()[tree1.getNode(subIndex1)].bb_
: bb1.subBbox(i1)
);
findNear
(
nearDist,
!okOrder,
tree2,
index2,
bb2,
tree1,
subIndex1,
subBb1,
cop
);
}
}
}
else if (tree1.isContent(index1))
{
const treeBoundBox searchBox
(
bb1.min()-nearSpan,
bb1.max()+nearSpan
);
if (tree2.isNode(index2))
{
const node& nod2 =
tree2.nodes()[tree2.getNode(index2)];
if (bb2.overlaps(searchBox))
{
for (direction i2 = 0; i2 < nod2.subNodes_.size(); i2++)
{
labelBits subIndex2 = nod2.subNodes_[i2];
const treeBoundBox subBb2
(
tree2.isNode(subIndex2)
? tree2.nodes()[tree2.getNode(subIndex2)].bb_
: bb2.subBbox(i2)
);
findNear
(
nearDist,
!okOrder,
tree2,
subIndex2,
subBb2,
tree1,
index1,
bb1,
cop
);
}
}
}
else if (tree2.isContent(index2))
{
// Both are leaves. Check n^2.
const labelList& indices1 =
tree1.contents()[tree1.getContent(index1)];
const labelList& indices2 =
tree2.contents()[tree2.getContent(index2)];
forAll(indices1, i)
{
label shape1 = indices1[i];
forAll(indices2, j)
{
label shape2 = indices2[j];
if ((&tree1 != &tree2) || (shape1 != shape2))
{
if (okOrder)
{
cop
(
nearDist,
tree1.shapes(),
shape1,
tree2.shapes(),
shape2
);
}
else
{
cop
(
nearDist,
tree2.shapes(),
shape2,
tree1.shapes(),
shape1
);
}
}
}
}
}
}
}
// Number of elements in node.
template <class Type>
Foam::label Foam::indexedOctree<Type>::countElements
@ -2620,6 +2794,30 @@ Foam::indexedOctree<Type>::getVolumeType
}
template <class Type>
template <class CompareOp>
void Foam::indexedOctree<Type>::findNear
(
const scalar nearDist,
const indexedOctree<Type>& tree2,
CompareOp& cop
) const
{
findNear
(
nearDist,
true,
*this,
nodePlusOctant(0, 0),
bb(),
tree2,
nodePlusOctant(0, 0),
tree2.bb(),
cop
);
}
// Print contents of nodeI
template <class Type>
void Foam::indexedOctree<Type>::print

26
src/meshTools/indexedOctree/indexedOctree.H
View File

@ -336,6 +336,22 @@ private:
labelHashSet& elements
) const;
template <class CompareOp>
static void findNear
(
const scalar nearDist,
const bool okOrder,
const indexedOctree<Type>& tree1,
const labelBits index1,
const treeBoundBox& bb1,
const indexedOctree<Type>& tree2,
const labelBits index2,
const treeBoundBox& bb2,
CompareOp& cop
);
// Other
//- Count number of elements on this and sublevels
@ -581,6 +597,16 @@ public:
const point& sample
);
//- Find near pairs and apply CompareOp to them.
// tree2 can be *this or different tree.
template <class CompareOp>
void findNear
(
const scalar nearDist,
const indexedOctree<Type>& tree2,
CompareOp& cop
) const;
// Write

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

@ -35,145 +35,145 @@ defineTypeNameAndDebug(Foam::treeDataTriSurface, 0);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// Fast distance to triangle calculation. From
// "Distance Between Point and Trangle 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);
// // 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;
// // Precalculate distance factors.
// const scalar d = E0 & D;
// const scalar e = E1 & D;
// Do classification
const scalar det = a*c - b*b;
// // Do classification
// const scalar det = a*c - b*b;
s = b*e - c*d;
t = b*d - a*e;
// 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;
}
// 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.
// // 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);
}
// 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 * * * * * * * * * * * * * * //
@ -234,9 +234,7 @@ Foam::label Foam::treeDataTriSurface::getVolumeType
(
surface_,
sample,
pHit.index(),
pHit.hitPoint(),
indexedOctree<treeDataTriSurface>::perturbTol()
pHit.index()
);
if (t == triSurfaceTools::UNKNOWN)
@ -353,39 +351,43 @@ void Foam::treeDataTriSurface::findNearest
// )
//)
{
// Get spanning vectors of triangle
vector base(p1);
vector E0(p0 - p1);
vector E1(p2 - p1);
// // 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);
// 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;
// // 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,
// scalar distSqr = nearestCoords
// (
// base,
// E0,
// E1,
// a,
// b,
// c,
// sample,
s,
t
);
// s,
// t
// );
pointHit pHit = triPointRef(p0, p1, p2).nearestPoint(sample);
scalar distSqr = sqr(pHit.distance());
if (distSqr < nearestDistSqr)
{
nearestDistSqr = distSqr;
minIndex = index;
nearestPoint = base + s*E0 + t*E1;
nearestPoint = pHit.rawPoint();
}
}
}

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

@ -60,20 +60,20 @@ class treeDataTriSurface
// 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
);
// //- 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:

5
src/meshTools/triSurface/booleanOps/surfaceIntersection/edgeIntersections.C
View File

@ -430,10 +430,7 @@ bool Foam::edgeIntersections::offsetPerturb
point ctr = tri.centre();
// Get measure for tolerance.
scalar tolDim = 0.001*mag(tri.a() - ctr);
tri.classify(pHit.hitPoint(), tolDim, nearType, nearLabel);
tri.classify(pHit.hitPoint(), nearType, nearLabel);
if (nearType == triPointRef::POINT || nearType == triPointRef::EDGE)
{

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

@ -315,7 +315,7 @@ void Foam::surfaceIntersection::classifyHit
surf2Pts[f2[0]],
surf2Pts[f2[1]],
surf2Pts[f2[2]]
).classify(pHit.hitPoint(), tolDim, nearType, nearLabel);
).classify(pHit.hitPoint(), nearType, nearLabel);
// Classify points on edge of surface1
label edgeEnd =

2
src/meshTools/triSurface/octreeData/octreeDataTriSurface.C
View File

@ -43,7 +43,7 @@ Foam::scalar Foam::octreeDataTriSurface::tol(1E-6);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// Fast distance to triangle calculation. From
// "Distance Between Point and Trangle in 3D"
// "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.
void Foam::octreeDataTriSurface::nearestCoords

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

@ -234,9 +234,7 @@ void Foam::orientedSurface::propagateOrientation
(
s,
samplePoint,
nearestFaceI,
nearestPt,
10*SMALL
nearestFaceI
);
if (side == triSurfaceTools::UNKNOWN)

82
src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C
View File

@ -2121,12 +2121,13 @@ Foam::vector Foam::triSurfaceTools::surfaceNormal
label nearType;
label nearLabel;
triPointRef
(
points[f[0]],
points[f[1]],
points[f[2]]
).classify(nearestPt, 1E-6, nearType, nearLabel);
).classify(nearestPt, nearType, nearLabel);
if (nearType == triPointRef::NONE)
{
@ -2199,28 +2200,61 @@ Foam::triSurfaceTools::sideType Foam::triSurfaceTools::surfaceSide
(
const triSurface& surf,
const point& sample,
const label nearestFaceI, // nearest face
const point& nearestPoint, // nearest point on nearest face
const scalar tol
const label nearestFaceI
)
{
const labelledTri& f = surf[nearestFaceI];
const pointField& points = surf.points();
// Find where point is on triangle. Note tolerance needed. Is relative
// one (used in comparing normalized [0..1] triangle coordinates).
// Find where point is on triangle.
label nearType, nearLabel;
triPointRef
pointHit pHit = triPointRef
(
points[f[0]],
points[f[1]],
points[f[2]]
).classify(nearestPoint, tol, nearType, nearLabel);
).nearestPointClassify(sample, nearType, nearLabel);
const point& nearestPoint(pHit.rawPoint());
if (nearType == triPointRef::NONE)
{
vector sampleNearestVec = (sample - nearestPoint);
// Nearest to face interior. Use faceNormal to determine side
scalar c = (sample - nearestPoint) & surf.faceNormals()[nearestFaceI];
scalar c = sampleNearestVec & surf.faceNormals()[nearestFaceI];
// // If the sample is essentially on the face, do not check for
// // it being perpendicular.
// scalar magSampleNearestVec = mag(sampleNearestVec);
// if (magSampleNearestVec > SMALL)
// {
// c /= magSampleNearestVec*mag(surf.faceNormals()[nearestFaceI]);
// if (mag(c) < 0.99)
// {
// FatalErrorIn("triSurfaceTools::surfaceSide")
// << "nearestPoint identified as being on triangle face "
// << "but vector from nearestPoint to sample is not "
// << "perpendicular to the normal." << nl
// << "sample: " << sample << nl
// << "nearestPoint: " << nearestPoint << nl
// << "sample - nearestPoint: "
// << sample - nearestPoint << nl
// << "normal: " << surf.faceNormals()[nearestFaceI] << nl
// << "mag(sample - nearestPoint): "
// << mag(sample - nearestPoint) << nl
// << "normalised dot product: " << c << nl
// << "triangle vertices: " << nl
// << " " << points[f[0]] << nl
// << " " << points[f[1]] << nl
// << " " << points[f[2]] << nl
// << abort(FatalError);
// }
// }
if (c > 0)
{
@ -2239,13 +2273,13 @@ Foam::triSurfaceTools::sideType Foam::triSurfaceTools::surfaceSide
// Get the edge. Assume order of faceEdges same as face vertices.
label edgeI = surf.faceEdges()[nearestFaceI][nearLabel];
//if (debug)
//{
// if (debug)
// {
// // Check order of faceEdges same as face vertices.
// const edge& e = surf.edges()[edgeI];
// const labelList& meshPoints = surf.meshPoints();
// const edge meshEdge(meshPoints[e[0]], meshPoints[e[1]]);
//
// if
// (
// meshEdge
@ -2259,7 +2293,7 @@ Foam::triSurfaceTools::sideType Foam::triSurfaceTools::surfaceSide
// << " in face " << f
// << abort(FatalError);
// }
//}
// }
return edgeSide(surf, sample, nearestPoint, edgeI);
}
@ -2717,7 +2751,14 @@ void Foam::triSurfaceTools::calcInterpolationWeights
triPointRef tri(f.tri(points));
pointHit nearest = tri.nearestPoint(samplePt);
label nearType, nearLabel;
pointHit nearest = tri.nearestPointClassify
(
samplePt,
nearType,
nearLabel
);
if (nearest.hit())
{
@ -2741,14 +2782,6 @@ void Foam::triSurfaceTools::calcInterpolationWeights
minDistance = nearest.distance();
// Outside triangle. Store nearest.
label nearType, nearLabel;
tri.classify
(
nearest.rawPoint(),
1E-6, // relative tolerance
nearType,
nearLabel
);
if (nearType == triPointRef::POINT)
{
@ -2779,12 +2812,12 @@ void Foam::triSurfaceTools::calcInterpolationWeights
max
(
0,
mag(nearest.rawPoint()-p0)/mag(p1-p0)
mag(nearest.rawPoint() - p0)/mag(p1 - p0)
)
);
// Interpolate
weights[0] = 1-s;
weights[0] = 1 - s;
weights[1] = s;
weights[2] = -GREAT;
@ -2830,7 +2863,6 @@ Foam::surfaceLocation Foam::triSurfaceTools::classify
triPointRef(s[triI].tri(s.points())).classify
(
trianglePoint,
1E-6,
elemType,
index
);

4
src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H
View File

@ -458,9 +458,7 @@ public:
(
const triSurface& surf,
const point& sample,
const label nearestFaceI, // nearest face
const point& nearestPt, // nearest point on nearest face
const scalar tol // tolerance for nearness test.
const label nearestFaceI
);
// Triangulation of faces

7
src/postProcessing/functionObjects/field/streamLine/streamLineParticle.C
View File

@ -168,7 +168,6 @@ bool Foam::streamLineParticle::move(streamLineParticle::trackData& td)
td.keepParticle
&& !td.switchProcessor
&& lifeTime_ > 0
&& tEnd > ROOTVSMALL
)
{
// TBD: implement subcycling so step through cells in more than
@ -191,6 +190,12 @@ bool Foam::streamLineParticle::move(streamLineParticle::trackData& td)
tEnd -= dt;
stepFraction() = 1.0 - tEnd/deltaT;
if (tEnd <= ROOTVSMALL)
{
// Force removal
lifeTime_ = 0;
}
}
if (!td.keepParticle || lifeTime_ == 0)

72
src/sampling/sampledSurface/distanceSurface/distanceSurface.C
View File

@ -89,20 +89,29 @@ void Foam::distanceSurface::createGeometry()
if (signed_)
{
vectorField normal;
surfPtr_().getNormal(nearest, normal);
List<searchableSurface::volumeType> volType;
forAll(nearest, i)
surfPtr_().getVolumeType(cc, volType);
forAll(volType, i)
{
vector d(cc[i]-nearest[i].hitPoint());
searchableSurface::volumeType vT = volType[i];
if ((d&normal[i]) > 0)
if (vT == searchableSurface::OUTSIDE)
{
fld[i] = Foam::mag(d);
fld[i] = Foam::mag(cc[i] - nearest[i].hitPoint());
}
else if (vT == searchableSurface::INSIDE)
{
fld[i] = -Foam::mag(cc[i] - nearest[i].hitPoint());
}
else
{
fld[i] = -Foam::mag(d);
FatalErrorIn
(
"void Foam::distanceSurface::createGeometry()"
) << "getVolumeType failure, neither INSIDE or OUTSIDE"
<< exit(FatalError);
}
}
}
@ -132,20 +141,30 @@ void Foam::distanceSurface::createGeometry()
if (signed_)
{
vectorField normal;
surfPtr_().getNormal(nearest, normal);
List<searchableSurface::volumeType> volType;
forAll(nearest, i)
surfPtr_().getVolumeType(cc, volType);
forAll(volType, i)
{
vector d(cc[i]-nearest[i].hitPoint());
searchableSurface::volumeType vT = volType[i];
if ((d&normal[i]) > 0)
if (vT == searchableSurface::OUTSIDE)
{
fld[i] = Foam::mag(d);
fld[i] = Foam::mag(cc[i] - nearest[i].hitPoint());
}
else if (vT == searchableSurface::INSIDE)
{
fld[i] = -Foam::mag(cc[i] - nearest[i].hitPoint());
}
else
{
fld[i] = -Foam::mag(d);
FatalErrorIn
(
"void Foam::distanceSurface::createGeometry()"
) << "getVolumeType failure, "
<< "neither INSIDE or OUTSIDE"
<< exit(FatalError);
}
}
}
@ -179,20 +198,31 @@ void Foam::distanceSurface::createGeometry()
if (signed_)
{
vectorField normal;
surfPtr_().getNormal(nearest, normal);
List<searchableSurface::volumeType> volType;
forAll(nearest, i)
surfPtr_().getVolumeType(pts, volType);
forAll(volType, i)
{
vector d(pts[i]-nearest[i].hitPoint());
searchableSurface::volumeType vT = volType[i];
if ((d&normal[i]) > 0)
if (vT == searchableSurface::OUTSIDE)
{
pointDistance_[i] = Foam::mag(d);
pointDistance_[i] =
Foam::mag(pts[i] - nearest[i].hitPoint());
}
else if (vT == searchableSurface::INSIDE)
{
pointDistance_[i] =
-Foam::mag(pts[i] - nearest[i].hitPoint());
}
else
{
pointDistance_[i] = -Foam::mag(d);
FatalErrorIn
(
"void Foam::distanceSurface::createGeometry()"
) << "getVolumeType failure, neither INSIDE or OUTSIDE"
<< exit(FatalError);
}
}
}