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

removed octreeDataFaceList; created generic PrimitivePatch search one

Browse Source
This commit is contained in:
mattijs
2009-11-11 17:19:52 +00:00
parent dd746d4e4b
commit d7e8802896
8 changed files with 888 additions and 855 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
src/dynamicMesh/Make/files
View File

@ -49,7 +49,6 @@ perfectInterface/perfectInterface.C
boundaryMesh/boundaryMesh.C
boundaryPatch/boundaryPatch.C
boundaryMesh/octreeDataFaceList.C
setUpdater/setUpdater.C
meshModifiers = meshCut/meshModifiers

138
src/dynamicMesh/boundaryMesh/boundaryMesh.C
View File

@ -29,11 +29,12 @@ License
#include "polyMesh.H"
#include "repatchPolyTopoChanger.H"
#include "faceList.H"
#include "octree.H"
#include "octreeDataFaceList.H"
#include "indexedOctree.H"
#include "treeDataPrimitivePatch.H"
#include "triSurface.H"
#include "SortableList.H"
#include "OFstream.H"
#include "uindirectPrimitivePatch.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -892,6 +893,17 @@ Foam::labelList Foam::boundaryMesh::getNearest
<< endl;
}
uindirectPrimitivePatch leftPatch
(
UIndirectList<face>(mesh(), leftFaces),
mesh().points()
);
uindirectPrimitivePatch rightPatch
(
UIndirectList<face>(mesh(), rightFaces),
mesh().points()
);
// Overall bb
treeBoundBox overallBb(mesh().localPoints());
@ -911,29 +923,35 @@ Foam::labelList Foam::boundaryMesh::getNearest
bbMax.z() += 2*tol;
// Create the octrees
octree<octreeDataFaceList> leftTree
indexedOctree
<
treeDataPrimitivePatch<face, UIndirectList, const pointField&>
> leftTree
(
overallBb,
octreeDataFaceList
treeDataPrimitivePatch<face, UIndirectList, const pointField&>
(
mesh(),
leftFaces
false, // cacheBb
leftPatch
),
1,
20,
10
overallBb,
10, // maxLevel
10, // leafSize
3.0 // duplicity
);
octree<octreeDataFaceList> rightTree
indexedOctree
<
treeDataPrimitivePatch<face, UIndirectList, const pointField&>
> rightTree
(
overallBb,
octreeDataFaceList
treeDataPrimitivePatch<face, UIndirectList, const pointField&>
(
mesh(),
rightFaces
false, // cacheBb
rightPatch
),
1,
20,
10
overallBb,
10, // maxLevel
10, // leafSize
3.0 // duplicity
);
if (debug)
@ -953,7 +971,7 @@ Foam::labelList Foam::boundaryMesh::getNearest
treeBoundBox tightest;
const scalar searchDim = mag(searchSpan);
const scalar searchDimSqr = magSqr(searchSpan);
forAll(nearestBFaceI, patchFaceI)
{
@ -982,50 +1000,25 @@ Foam::labelList Foam::boundaryMesh::getNearest
}
// Search right tree
tightest.min() = ctr - searchSpan;
tightest.max() = ctr + searchSpan;
scalar rightDist = searchDim;
label rightI = rightTree.findNearest(ctr, tightest, rightDist);
pointIndexHit rightInfo = rightTree.findNearest(ctr, searchDimSqr);
// Search left tree. Note: could start from rightDist bounding box
// instead of starting from top.
tightest.min() = ctr - searchSpan;
tightest.max() = ctr + searchSpan;
scalar leftDist = searchDim;
label leftI = leftTree.findNearest(ctr, tightest, leftDist);
pointIndexHit leftInfo = leftTree.findNearest(ctr, searchDimSqr);
if (rightI == -1)
if (rightInfo.hit())
{
// No face found in right tree.
if (leftI == -1)
{
// No face found in left tree.
nearestBFaceI[patchFaceI] = -1;
}
else
{
// Found in left but not in right. Choose left regardless if
// correct sign. Note: do we want this?
nearestBFaceI[patchFaceI] = leftFaces[leftI];
}
}
else
{
if (leftI == -1)
{
// Found in right but not in left. Choose right regardless if
// correct sign. Note: do we want this?
nearestBFaceI[patchFaceI] = rightFaces[rightI];
}
else
if (leftInfo.hit())
{
// Found in both trees. Compare normals.
label rightFaceI = rightFaces[rightInfo.index()];
label leftFaceI = leftFaces[leftInfo.index()];
scalar rightSign = n & ns[rightFaces[rightI]];
scalar leftSign = n & ns[leftFaces[leftI]];
label rightDist = mag(rightInfo.hitPoint()-ctr);
label leftDist = mag(leftInfo.hitPoint()-ctr);
scalar rightSign = n & ns[rightFaceI];
scalar leftSign = n & ns[leftFaceI];
if
(
@ -1036,11 +1029,11 @@ Foam::labelList Foam::boundaryMesh::getNearest
// Both same sign. Choose nearest.
if (rightDist < leftDist)
{
nearestBFaceI[patchFaceI] = rightFaces[rightI];
nearestBFaceI[patchFaceI] = rightFaceI;
}
else
{
nearestBFaceI[patchFaceI] = leftFaces[leftI];
nearestBFaceI[patchFaceI] = leftFaceI;
}
}
else
@ -1059,11 +1052,11 @@ Foam::labelList Foam::boundaryMesh::getNearest
// Different sign and nearby. Choosing matching normal
if (rightSign > 0)
{
nearestBFaceI[patchFaceI] = rightFaces[rightI];
nearestBFaceI[patchFaceI] = rightFaceI;
}
else
{
nearestBFaceI[patchFaceI] = leftFaces[leftI];
nearestBFaceI[patchFaceI] = leftFaceI;
}
}
else
@ -1071,15 +1064,38 @@ Foam::labelList Foam::boundaryMesh::getNearest
// Different sign but faraway. Choosing nearest.
if (rightDist < leftDist)
{
nearestBFaceI[patchFaceI] = rightFaces[rightI];
nearestBFaceI[patchFaceI] = rightFaceI;
}
else
{
nearestBFaceI[patchFaceI] = leftFaces[leftI];
nearestBFaceI[patchFaceI] = leftFaceI;
}
}
}
}
else
{
// Found in right but not in left. Choose right regardless if
// correct sign. Note: do we want this?
label rightFaceI = rightFaces[rightInfo.index()];
nearestBFaceI[patchFaceI] = rightFaceI;
}
}
else
{
// No face found in right tree.
if (leftInfo.hit())
{
// Found in left but not in right. Choose left regardless if
// correct sign. Note: do we want this?
nearestBFaceI[patchFaceI] = leftFaces[leftInfo.index()];
}
else
{
// No face found in left tree.
nearestBFaceI[patchFaceI] = -1;
}
}
}

573
src/dynamicMesh/boundaryMesh/octreeDataFaceList.C
View File

@ -1,573 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Description
\*---------------------------------------------------------------------------*/
#include "octreeDataFaceList.H"
#include "octree.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(Foam::octreeDataFaceList, 0);
Foam::scalar Foam::octreeDataFaceList::tol = 1E-6;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::octreeDataFaceList::calcBb()
{
allBb_.setSize(faceLabels_.size());
allBb_ = treeBoundBox
(
vector(GREAT, GREAT, GREAT),
vector(-GREAT, -GREAT, -GREAT)
);
forAll (faceLabels_, faceLabelI)
{
label faceI = faceLabels_[faceLabelI];
// Update bb of face
treeBoundBox& myBb = allBb_[faceLabelI];
const face& f = mesh_.localFaces()[faceI];
forAll(f, fp)
{
const point& coord = mesh_.localPoints()[f[fp]];
myBb.min() = min(myBb.min(), coord);
myBb.max() = max(myBb.max(), coord);
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from all faces in bMesh
Foam::octreeDataFaceList::octreeDataFaceList(const bMesh& mesh)
:
mesh_(mesh),
faceLabels_(mesh_.size()),
allBb_(mesh_.size())
{
forAll(faceLabels_, faceI)
{
faceLabels_[faceI] = faceI;
}
// Generate tight fitting bounding box
calcBb();
}
// Construct from selected faces in bMesh
Foam::octreeDataFaceList::octreeDataFaceList
(
const bMesh& mesh,
const labelList& faceLabels
)
:
mesh_(mesh),
faceLabels_(faceLabels),
allBb_(faceLabels.size())
{
// Generate tight fitting bounding box
calcBb();
}
// Construct as copy
Foam::octreeDataFaceList::octreeDataFaceList(const octreeDataFaceList& shapes)
:
mesh_(shapes.mesh()),
faceLabels_(shapes.faceLabels()),
allBb_(shapes.allBb())
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::octreeDataFaceList::~octreeDataFaceList()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::label Foam::octreeDataFaceList::getSampleType
(
const octree<octreeDataFaceList>& oc,
const point& sample
) const
{
// Need to determine whether sample is 'inside' or 'outside'
// Done by finding nearest face. This gives back a face which is
// guaranteed to contain nearest point. This point can be
// - in interior of face: compare to face normal
// - on edge of face: compare to edge normal
// - on point of face: compare to point normal
// Unfortunately the octree does not give us back the intersection point
// or where on the face it has hit so we have to recreate all that
// information.
// Find nearest face to sample
treeBoundBox tightest(treeBoundBox::greatBox);
scalar tightestDist = GREAT;
label index = oc.findNearest(sample, tightest, tightestDist);
if (index == -1)
{
FatalErrorIn
(
"octreeDataFaceList::getSampleType"
"(octree<octreeDataFaceList>&, const point&)"
) << "Could not find " << sample << " in octree."
<< abort(FatalError);
}
label faceI = faceLabels_[index];
// Get actual intersection point on face
if (debug & 2)
{
Pout<< "getSampleType : sample:" << sample
<< " nearest face:" << faceI;
}
const face& f = mesh_.localFaces()[faceI];
const pointField& points = mesh_.localPoints();
pointHit curHit = f.nearestPoint(sample, points);
//
// 1] Check whether sample is above face
//
if (curHit.hit())
{
// Simple case. Compare to face normal.
if (debug & 2)
{
Pout<< " -> face hit:" << curHit.hitPoint()
<< " comparing to face normal " << mesh_.faceNormals()[faceI]
<< endl;
}
return octree<octreeDataFaceList>::getVolType
(
mesh_.faceNormals()[faceI],
sample - curHit.hitPoint()
);
}
if (debug & 2)
{
Pout<< " -> face miss:" << curHit.missPoint();
}
//
// 2] Check whether intersection is on one of the face vertices or
// face centre
//
// typical dimension as sqrt of face area.
scalar typDim = sqrt(mag(f.normal(points))) + VSMALL;
forAll(f, fp)
{
if ((mag(points[f[fp]] - curHit.missPoint())/typDim) < tol)
{
// Face intersection point equals face vertex fp
if (debug & 2)
{
Pout<< " -> face point hit :" << points[f[fp]]
<< " point normal:" << mesh_.pointNormals()[f[fp]]
<< " distance:"
<< mag(points[f[fp]] - curHit.missPoint())/typDim
<< endl;
}
return octree<octreeDataFaceList>::getVolType
(
mesh_.pointNormals()[f[fp]],
sample - curHit.missPoint()
);
}
}
// Get face centre
point ctr(f.centre(points));
if ((mag(ctr - curHit.missPoint())/typDim) < tol)
{
// Face intersection point equals face centre. Normal at face centre
// is already average of face normals
if (debug & 2)
{
Pout<< " -> centre hit:" << ctr
<< " distance:"
<< mag(ctr - curHit.missPoint())/typDim
<< endl;
}
return octree<octreeDataFaceList>::getVolType
(
mesh_.faceNormals()[faceI],
sample - curHit.missPoint()
);
}
//
// 3] Get the 'real' edge the face intersection is on
//
const labelList& myEdges = mesh_.faceEdges()[faceI];
forAll(myEdges, myEdgeI)
{
const edge& e = mesh_.edges()[myEdges[myEdgeI]];
pointHit edgeHit =
line<point, const point&>
(
points[e.start()],
points[e.end()]
).nearestDist(sample);
point edgePoint;
if (edgeHit.hit())
{
edgePoint = edgeHit.hitPoint();
}
else
{
edgePoint = edgeHit.missPoint();
}
if ((mag(edgePoint - curHit.missPoint())/typDim) < tol)
{
// Face intersection point lies on edge e
// Calculate edge normal (wrong: uses face normals instead of
// triangle normals)
const labelList& myFaces = mesh_.edgeFaces()[myEdges[myEdgeI]];
vector edgeNormal(vector::zero);
forAll(myFaces, myFaceI)
{
edgeNormal += mesh_.faceNormals()[myFaces[myFaceI]];
}
if (debug & 2)
{
Pout<< " -> real edge hit point:" << edgePoint
<< " comparing to edge normal:" << edgeNormal
<< endl;
}
// Found face intersection point on this edge. Compare to edge
// normal
return octree<octreeDataFaceList>::getVolType
(
edgeNormal,
sample - curHit.missPoint()
);
}
}
//
// 4] Get the internal edge (vertex - faceCentre) the face intersection
// is on
//
forAll(f, fp)
{
pointHit edgeHit =
line<point, const point&>
(
points[f[fp]],
ctr
).nearestDist(sample);
point edgePoint;
if (edgeHit.hit())
{
edgePoint = edgeHit.hitPoint();
}
else
{
edgePoint = edgeHit.missPoint();
}
if ((mag(edgePoint - curHit.missPoint())/typDim) < tol)
{
// Face intersection point lies on edge between two face triangles
// Calculate edge normal as average of the two triangle normals
label fpPrev = f.rcIndex(fp);
label fpNext = f.fcIndex(fp);
vector e = points[f[fp]] - ctr;
vector ePrev = points[f[fpPrev]] - ctr;
vector eNext = points[f[fpNext]] - ctr;
vector nLeft = ePrev ^ e;
nLeft /= mag(nLeft) + VSMALL;
vector nRight = e ^ eNext;
nRight /= mag(nRight) + VSMALL;
if (debug & 2)
{
Pout<< " -> internal edge hit point:" << edgePoint
<< " comparing to edge normal "
<< 0.5*(nLeft + nRight)
<< endl;
}
// Found face intersection point on this edge. Compare to edge
// normal
return octree<octreeDataFaceList>::getVolType
(
0.5*(nLeft + nRight),
sample - curHit.missPoint()
);
}
}
if (debug & 2)
{
Pout<< "Did not find sample " << sample
<< " anywhere related to nearest face " << faceI << endl
<< "Face:";
forAll(f, fp)
{
Pout<< " vertex:" << f[fp] << " coord:" << points[f[fp]]
<< endl;
}
}
// Can't determine status of sample with respect to nearest face.
// Either
// - tolerances are wrong. (if e.g. face has zero area)
// - or (more likely) surface is not closed.
return octree<octreeDataFaceList>::UNKNOWN;
}
bool Foam::octreeDataFaceList::overlaps
(
const label index,
const treeBoundBox& sampleBb
) const
{
return sampleBb.overlaps(allBb_[index]);
}
bool Foam::octreeDataFaceList::contains
(
const label,
const point&
) const
{
notImplemented
(
"octreeDataFaceList::contains(const label, const point&)"
);
return false;
}
bool Foam::octreeDataFaceList::intersects
(
const label index,
const point& start,
const point& end,
point& intersectionPoint
) const
{
label faceI = faceLabels_[index];
const face& f = mesh_.localFaces()[faceI];
const vector dir(end - start);
// Disable picking up intersections behind us.
scalar oldTol = intersection::setPlanarTol(0.0);
pointHit inter =
f.ray
(
start,
dir,
mesh_.localPoints(),
intersection::HALF_RAY,
intersection::VECTOR
);
intersection::setPlanarTol(oldTol);
if (inter.hit() && inter.distance() <= mag(dir))
{
intersectionPoint = inter.hitPoint();
return true;
}
else
{
return false;
}
}
bool Foam::octreeDataFaceList::findTightest
(
const label index,
const point& sample,
treeBoundBox& tightest
) const
{
// Get nearest and furthest away vertex
point myNear, myFar;
allBb_[index].calcExtremities(sample, myNear, myFar);
const point dist = myFar - sample;
scalar myFarDist = mag(dist);
point tightestNear, tightestFar;
tightest.calcExtremities(sample, tightestNear, tightestFar);
scalar tightestFarDist = mag(tightestFar - sample);
if (tightestFarDist < myFarDist)
{
// Keep current tightest.
return false;
}
else
{
// Construct bb around sample and myFar
const point dist2(fabs(dist.x()), fabs(dist.y()), fabs(dist.z()));
tightest.min() = sample - dist2;
tightest.max() = sample + dist2;
return true;
}
}
// Determine numerical value of sign of sample compared to shape at index
Foam::scalar Foam::octreeDataFaceList::calcSign
(
const label index,
const point& sample,
vector&
) const
{
label faceI = faceLabels_[index];
const face& f = mesh_.localFaces()[faceI];
point ctr = f.centre(mesh_.localPoints());
vector vec = sample - ctr;
vec /= mag(vec) + VSMALL;
return mesh_.faceNormals()[faceI] & vec;
}
// Calculate nearest point on/in shapei
Foam::scalar Foam::octreeDataFaceList::calcNearest
(
const label index,
const point& sample,
point& nearest
) const
{
label faceI = faceLabels_[index];
const face& f = mesh_.localFaces()[faceI];
pointHit nearHit = f.nearestPoint(sample, mesh_.localPoints());
if (nearHit.hit())
{
nearest = nearHit.hitPoint();
}
else
{
nearest = nearHit.missPoint();
}
if (debug & 1)
{
point ctr = f.centre(mesh_.localPoints());
scalar sign = mesh_.faceNormals()[faceI] & (sample - nearest);
Pout<< "octreeDataFaceList::calcNearest : "
<< "sample:" << sample
<< " faceI:" << faceI
<< " ctr:" << ctr
<< " sign:" << sign
<< " nearest point:" << nearest
<< " distance to face:" << nearHit.distance()
<< endl;
}
return nearHit.distance();
}
void Foam::octreeDataFaceList::write
(
Ostream& os,
const label index
) const
{
os << faceLabels_[index] << " " << allBb_[index];
}
// ************************************************************************* //

220
src/dynamicMesh/boundaryMesh/octreeDataFaceList.H
View File

@ -1,220 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::octreeDataFaceList
Description
Holds data for octree to work on list of faces on a bMesh
(= PrimitivePatch which holds faces, not references them)
Same as octreeDataFace except for that.
SourceFiles
octreeDataFaceList.C
\*---------------------------------------------------------------------------*/
#ifndef octreeDataFaceList_H
#define octreeDataFaceList_H
#include "treeBoundBoxList.H"
#include "faceList.H"
#include "point.H"
#include "className.H"
#include "bMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
template<class Type> class octree;
/*---------------------------------------------------------------------------*\
Class octreeDataFaceList Declaration
\*---------------------------------------------------------------------------*/
class octreeDataFaceList
{
// Static data
//- tolerance on linear dimensions
static scalar tol;
// Static function
static inline label nexti(label max, label i)
{
return (i + 1) % max;
}
// Private data
//- the mesh
const bMesh& mesh_;
//- labels (in mesh indexing) of faces
labelList faceLabels_;
//- bbs for all above faces
treeBoundBoxList allBb_;
// Private Member Functions
//- Set allBb to tight fitting bounding box
void calcBb();
public:
// Declare name of the class and its debug switch
ClassName("octreeDataFaceList");
// Constructors
//- Construct from all faces in bMesh.
octreeDataFaceList(const bMesh& mesh);
//- Construct from selected faces in bMesh.
octreeDataFaceList(const bMesh& mesh, const labelList& faceLabels);
//- Construct as copy
octreeDataFaceList(const octreeDataFaceList&);
// Destructor
~octreeDataFaceList();
// Member Functions
// Access
const bMesh& mesh() const
{
return mesh_;
}
const labelList& faceLabels() const
{
return faceLabels_;
}
const treeBoundBoxList& allBb() const
{
return allBb_;
}
label size() const
{
return allBb_.size();
}
// Search
//- Get type of sample
label getSampleType
(
const octree<octreeDataFaceList>&,
const point&
) const;
//- Does (bb of) shape at index overlap bb
bool overlaps
(
const label index,
const treeBoundBox& sampleBb
) const;
//- Does shape at index contain sample
bool contains
(
const label index,
const point& sample
) const;
//- Segment (from start to end) intersection with shape
// at index. If intersects returns true and sets intersectionPoint
bool intersects
(
const label index,
const point& start,
const point& end,
point& intersectionPoint
) const;
//- Sets newTightest to bounding box (and returns true) if
// nearer to sample than tightest bounding box. Otherwise
// returns false.
bool findTightest
(
const label index,
const point& sample,
treeBoundBox& tightest
) const;
//- Given index get unit normal and calculate (numerical) sign
// of sample.
// Used to determine accuracy of calcNearest or inside/outside.
scalar calcSign
(
const label index,
const point& sample,
vector& n
) const;
//- Calculates nearest (to sample) point in shape.
// Returns point and mag(nearest - sample). Returns GREAT if
// sample does not project onto (triangle decomposition) of face.
scalar calcNearest
(
const label index,
const point& sample,
point& nearest
) const;
// Edit
// Write
//- Write shape at index
void write(Ostream& os, const label index) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

1
src/meshTools/Make/files
View File

@ -50,6 +50,7 @@ indexedOctree/treeDataCell.C
indexedOctree/treeDataEdge.C
indexedOctree/treeDataFace.C
indexedOctree/treeDataPoint.C
indexedOctree/treeDataPrimitivePatchName.C
indexedOctree/treeDataTriSurface.C
searchableSurface = searchableSurface

567
src/meshTools/indexedOctree/treeDataPrimitivePatch.C Normal file
View File

@ -0,0 +1,567 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "treeDataPrimitivePatch.H"
#include "indexedOctree.H"
#include "triangleFuncs.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
template
<
class Face,
template<class> class FaceList,
class PointField,
class PointType
>
Foam::scalar
Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
tolSqr = sqr(1E-6);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template
<
class Face,
template<class> class FaceList,
class PointField,
class PointType
>
Foam::treeBoundBox
Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
calcBb
(
const pointField& points,
const face& f
)
{
treeBoundBox bb(points[f[0]], points[f[0]]);
for (label fp = 1; fp < f.size(); fp++)
{
const point& p = points[f[fp]];
bb.min() = min(bb.min(), p);
bb.max() = max(bb.max(), p);
}
return bb;
}
template
<
class Face,
template<class> class FaceList,
class PointField,
class PointType
>
void Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
update()
{
if (cacheBb_)
{
bbs_.setSize(patch_.size());
forAll(patch_, i)
{
bbs_[i] = calcBb(patch_.points(), patch_[i]);
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
template
<
class Face,
template<class> class FaceList,
class PointField,
class PointType
>
Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
treeDataPrimitivePatch
(
const bool cacheBb,
const PrimitivePatch<Face, FaceList, PointField, PointType>& patch
)
:
patch_(patch),
cacheBb_(cacheBb)
{
update();
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template
<
class Face,
template<class> class FaceList,
class PointField,
class PointType
>
Foam::pointField
Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
points() const
{
pointField cc(patch_.size());
forAll(patch_, i)
{
cc[i] = patch_[i].centre(patch_.points());
}
return cc;
}
//- Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
// Only makes sense for closed surfaces.
template
<
class Face,
template<class> class FaceList,
class PointField,
class PointType
>
Foam::label
Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
getVolumeType
(
const indexedOctree
<
treeDataPrimitivePatch
<
Face,
FaceList,
PointField,
PointType
>
>& oc,
const point& sample
) const
{
// Need to determine whether sample is 'inside' or 'outside'
// Done by finding nearest face. This gives back a face which is
// guaranteed to contain nearest point. This point can be
// - in interior of face: compare to face normal
// - on edge of face: compare to edge normal
// - on point of face: compare to point normal
// Unfortunately the octree does not give us back the intersection point
// or where on the face it has hit so we have to recreate all that
// information.
// Find nearest face to sample
pointIndexHit info = oc.findNearest(sample, sqr(GREAT));
if (info.index() == -1)
{
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();
if (debug & 2)
{
Pout<< "getSampleType : sample:" << sample
<< " nearest face:" << faceI;
}
const pointField& points = patch_.localPoints();
const face& f = patch_.localFaces()[faceI];
// Retest to classify where on face info is. Note: could be improved. We
// already have point.
pointHit curHit = f.nearestPoint(sample, points);
const vector area = f.normal(points);
const point& curPt = curHit.rawPoint();
//
// 1] Check whether sample is above face
//
if (curHit.hit())
{
// Nearest point inside face. Compare to face normal.
if (debug & 2)
{
Pout<< " -> face hit:" << curPt
<< " comparing to face normal " << area << endl;
}
return indexedOctree<treeDataPrimitivePatch>::getSide
(
area,
sample - curPt
);
}
if (debug & 2)
{
Pout<< " -> face miss:" << curPt;
}
//
// 2] Check whether intersection is on one of the face vertices or
// face centre
//
const scalar typDimSqr = mag(area) + VSMALL;
forAll(f, fp)
{
if ((magSqr(points[f[fp]] - curPt)/typDimSqr) < tolSqr)
{
// Face intersection point equals face vertex fp
// Calculate point normal (wrong: uses face normals instead of
// triangle normals)
return indexedOctree<treeDataPrimitivePatch>::getSide
(
patch_.pointNormals()[f[fp]],
sample - curPt
);
}
}
const point fc(f.centre(points));
if ((magSqr(fc - curPt)/typDimSqr) < tolSqr)
{
// Face intersection point equals face centre. Normal at face centre
// is already average of face normals
if (debug & 2)
{
Pout<< " -> centre hit:" << fc
<< " distance:" << magSqr(fc - curPt)/typDimSqr << endl;
}
return indexedOctree<treeDataPrimitivePatch>::getSide
(
area,
sample - curPt
);
}
//
// 3] Get the 'real' edge the face intersection is on
//
const labelList& fEdges = patch_.faceEdges()[faceI];
forAll(fEdges, fEdgeI)
{
label edgeI = fEdges[fEdgeI];
const edge& e = patch_.edges()[edgeI];
pointHit edgeHit = e.line(points).nearestDist(sample);
if ((magSqr(edgeHit.rawPoint() - curPt)/typDimSqr) < tolSqr)
{
// Face intersection point lies on edge e
// Calculate edge normal (wrong: uses face normals instead of
// triangle normals)
const labelList& eFaces = patch_.edgeFaces()[edgeI];
vector edgeNormal(vector::zero);
forAll(eFaces, i)
{
edgeNormal += patch_.faceNormal()[eFaces[i]];
}
if (debug & 2)
{
Pout<< " -> real edge hit point:" << edgeHit.rawPoint()
<< " comparing to edge normal:" << edgeNormal
<< endl;
}
// Found face intersection point on this edge. Compare to edge
// normal
return indexedOctree<treeDataPrimitivePatch>::getSide
(
edgeNormal,
sample - curPt
);
}
}
//
// 4] Get the internal edge the face intersection is on
//
forAll(f, fp)
{
pointHit edgeHit = linePointRef
(
points[f[fp]],
fc
).nearestDist(sample);
if ((magSqr(edgeHit.rawPoint() - curPt)/typDimSqr) < tolSqr)
{
// Face intersection point lies on edge between two face triangles
// Calculate edge normal as average of the two triangle normals
vector e = points[f[fp]] - fc;
vector ePrev = points[f[f.rcIndex(fp)]] - fc;
vector eNext = points[f[f.fcIndex(fp)]] - fc;
vector nLeft = ePrev ^ e;
nLeft /= mag(nLeft) + VSMALL;
vector nRight = e ^ eNext;
nRight /= mag(nRight) + VSMALL;
if (debug & 2)
{
Pout<< " -> internal edge hit point:" << edgeHit.rawPoint()
<< " comparing to edge normal "
<< 0.5*(nLeft + nRight)
<< endl;
}
// Found face intersection point on this edge. Compare to edge
// normal
return indexedOctree<treeDataPrimitivePatch>::getSide
(
0.5*(nLeft + nRight),
sample - curPt
);
}
}
if (debug & 2)
{
Pout<< "Did not find sample " << sample
<< " anywhere related to nearest face " << faceI << endl
<< "Face:";
forAll(f, fp)
{
Pout<< " vertex:" << f[fp] << " coord:" << points[f[fp]]
<< endl;
}
}
// Can't determine status of sample with respect to nearest face.
// Either
// - tolerances are wrong. (if e.g. face has zero area)
// - or (more likely) surface is not closed.
return indexedOctree<treeDataPrimitivePatch>::UNKNOWN;
}
// Check if any point on shape is inside cubeBb.
template
<
class Face,
template<class> class FaceList,
class PointField,
class PointType
>
bool
Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
overlaps
(
const label index,
const treeBoundBox& cubeBb
) const
{
// 1. Quick rejection: bb does not intersect face bb at all
if (cacheBb_)
{
if (!cubeBb.overlaps(bbs_[index]))
{
return false;
}
}
else
{
if (!cubeBb.overlaps(calcBb(patch_.points(), patch_[index])))
{
return false;
}
}
// 2. Check if one or more face points inside
const pointField& points = patch_.points();
const face& f = patch_[index];
forAll(f, fp)
{
if (cubeBb.contains(points[f[fp]]))
{
return true;
}
}
// 3. Difficult case: all points are outside but connecting edges might
// go through cube. Use triangle-bounding box intersection.
const point fc = f.centre(points);
forAll(f, fp)
{
bool triIntersects = triangleFuncs::intersectBb
(
points[f[fp]],
points[f[f.fcIndex(fp)]],
fc,
cubeBb
);
if (triIntersects)
{
return true;
}
}
return false;
}
// Calculate nearest point to sample. Updates (if any) nearestDistSqr, minIndex,
// nearestPoint.
template
<
class Face,
template<class> class FaceList,
class PointField,
class PointType
>
void
Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
findNearest
(
const labelList& indices,
const point& sample,
scalar& nearestDistSqr,
label& minIndex,
point& nearestPoint
) const
{
const pointField& points = patch_.points();
forAll(indices, i)
{
label index = indices[i];
const face& f = patch_[index];
pointHit nearHit = f.nearestPoint(sample, points);
scalar distSqr = sqr(nearHit.distance());
if (distSqr < nearestDistSqr)
{
nearestDistSqr = distSqr;
minIndex = index;
nearestPoint = nearHit.rawPoint();
}
}
}
template
<
class Face,
template<class> class FaceList,
class PointField,
class PointType
>
bool
Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
intersects
(
const label index,
const point& start,
const point& end,
point& intersectionPoint
) const
{
// Do quick rejection test
if (cacheBb_)
{
const treeBoundBox& faceBb = bbs_[index];
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
{
return false;
}
}
// ************************************************************************* //

210
src/meshTools/indexedOctree/treeDataPrimitivePatch.H Normal file
View File

@ -0,0 +1,210 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::treeDataPrimitivePatch
Description
Encapsulation of data needed to search on PrimitivePatches
SourceFiles
treeDataPrimitivePatch.C
\*---------------------------------------------------------------------------*/
#ifndef treeDataPrimitivePatch_H
#define treeDataPrimitivePatch_H
#include "PrimitivePatch.H"
//#include "indexedOctree.H"
#include "treeBoundBoxList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
template<class Type> class indexedOctree;
/*---------------------------------------------------------------------------*\
Class treeDataPrimitivePatchName Declaration
\*---------------------------------------------------------------------------*/
TemplateName(treeDataPrimitivePatch);
/*---------------------------------------------------------------------------*\
Class treeDataPrimitivePatch Declaration
\*---------------------------------------------------------------------------*/
template
<
class Face,
template<class> class FaceList,
class PointField,
class PointType=point
>
class treeDataPrimitivePatch
:
public treeDataPrimitivePatchName
{
// Static data
//- tolerance on linear dimensions
static scalar tolSqr;
// Private data
//- Underlying geometry
const PrimitivePatch<Face, FaceList, PointField, PointType>& patch_;
//- Whether to precalculate and store face bounding box
const bool cacheBb_;
//- face bounding boxes (valid only if cacheBb_)
treeBoundBoxList bbs_;
// Private Member Functions
//- Calculate face bounding box
static treeBoundBox calcBb(const pointField&, const face&);
//- Initialise all member data
void update();
public:
// Constructors
//- Construct from patch.
treeDataPrimitivePatch
(
const bool cacheBb,
const PrimitivePatch<Face, FaceList, PointField, PointType>&
);
// Member Functions
// Access
label size() const
{
return patch_.size();
}
//- Get representative point cloud for all shapes inside
// (one point per shape)
pointField points() const;
// Search
//- Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
// Only makes sense for closed surfaces.
label getVolumeType
(
const indexedOctree
<
treeDataPrimitivePatch
<
Face,
FaceList,
PointField,
PointType
>
>&,
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 labelList& 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 labelList& indices,
const linePointRef& ln,
treeBoundBox& tightest,
label& minIndex,
point& linePoint,
point& nearestPoint
) const
{
notImplemented
(
"treeDataPrimitivePatch::findNearest"
"(const labelList&, 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;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "treeDataPrimitivePatch.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

33
src/meshTools/indexedOctree/treeDataPrimitivePatchName.C Normal file
View File

@ -0,0 +1,33 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "treeDataPrimitivePatch.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(Foam::treeDataPrimitivePatchName, 0);
// ************************************************************************* //