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' into sixDofPatch

Browse Source
This commit is contained in:
graham
2010-01-19 13:31:36 +00:00
parent a820021a0a 583e1b1ed9
commit 129f796211
1102 changed files with 503632 additions and 60483 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
src/meshTools/Make/options
View File

@ -6,4 +6,5 @@ EXE_INC = \
LIB_LIBS = \
-ltriSurface \
-ldecompositionMethods \
-lmetisDecompositionMethod \
-llagrangian

16
src/meshTools/PointEdgeWave/PointEdgeWave.C
View File

@ -53,7 +53,7 @@ void Foam::PointEdgeWave<Type>::offset(const label val, labelList& elems)
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// Gets point-point correspondence. Is
// Gets point-point correspondence. Is
// - list of halfA points (in cyclic patch points)
// - list of halfB points (can overlap with A!)
// - for every patchPoint its corresponding point
@ -111,7 +111,7 @@ void Foam::PointEdgeWave<Type>::leaveDomain
) const
{
const labelList& meshPoints = patch.meshPoints();
forAll(patchPointLabels, i)
{
label patchPointI = patchPointLabels[i];
@ -134,7 +134,7 @@ void Foam::PointEdgeWave<Type>::enterDomain
) const
{
const labelList& meshPoints = patch.meshPoints();
forAll(patchPointLabels, i)
{
label patchPointI = patchPointLabels[i];
@ -169,7 +169,7 @@ void Foam::PointEdgeWave<Type>::transform
(
"PointEdgeWave<Type>::transform(const tensorField&, List<Type>&)"
) << "Parallel cyclics not supported" << abort(FatalError);
forAll(pointInfo, i)
{
pointInfo[i].transform(rotTensor[i]);
@ -470,7 +470,7 @@ void Foam::PointEdgeWave<Type>::handleProcPatches()
{
OPstream toNeighbour
(
(
Pstream::blocking,
procPatch.neighbProcNo()
);
@ -505,7 +505,7 @@ void Foam::PointEdgeWave<Type>::handleProcPatches()
);
fromNeighbour >> owner >> ownerIndex >> patchInfo;
}
}
if (debug)
{
@ -642,7 +642,7 @@ void Foam::PointEdgeWave<Type>::handleCyclicPatches()
if (debug)
{
Pout<< "Cyclic patch " << patchI << ' ' << patch.name()
Pout<< "Cyclic patch " << patchI << ' ' << patch.name()
<< " Changed on first half : " << halfAInfo.size()
<< " Changed on second half : " << halfBInfo.size()
<< endl;
@ -679,7 +679,7 @@ void Foam::PointEdgeWave<Type>::handleCyclicPatches()
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Iterate, propagating changedPointsInfo across mesh, until no change (or
// Iterate, propagating changedPointsInfo across mesh, until no change (or
// maxIter reached). Initial point values specified.
template <class Type>
Foam::PointEdgeWave<Type>::PointEdgeWave

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

@ -26,7 +26,7 @@ Class
Foam::cellInfo
Description
Holds information regarding type of cell. Used in inside/outside
Holds information regarding type of cell. Used in inside/outside
determination in cellClassification.
SourceFiles
@ -58,7 +58,7 @@ class cellInfo
// Private data
label type_;
// Private Member Functions
//- Update current cell/face type with neighbouring

8
src/meshTools/cellClassification/cellInfoI.H
View File

@ -87,7 +87,7 @@ inline bool Foam::cellInfo::update
return false;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
@ -182,7 +182,7 @@ inline bool Foam::cellInfo::updateCell
neighbourFaceI,
-1
);
}
}
// Update this with neighbour information
@ -203,7 +203,7 @@ inline bool Foam::cellInfo::updateFace
-1,
neighbourCellI
);
}
}
// Update this with neighbour information
inline bool Foam::cellInfo::updateFace
@ -222,7 +222,7 @@ inline bool Foam::cellInfo::updateFace
-1,
-1
);
}
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //

2
src/meshTools/cellDist/cellDistFuncs.H
View File

@ -131,7 +131,7 @@ public:
//- Size of largest patch (out of supplied subset of patches)
label maxPatchSize(const labelHashSet& patchIDs) const;
//- Sum of patch sizes (out of supplied subset of patches).
//- Sum of patch sizes (out of supplied subset of patches).
// Used in sizing arrays.
label sumPatchSize(const labelHashSet& patchIDs) const;

10
src/meshTools/cellDist/wallPoint/wallPointDataI.H
View File

@ -53,7 +53,7 @@ inline bool wallPointData<Type>::update
data_ = w2.data();
return true;
}
}
scalar diff = distSqr() - dist2;
@ -78,7 +78,7 @@ inline bool wallPointData<Type>::update
return true;
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
@ -140,7 +140,7 @@ inline bool wallPointData<Type>::updateCell
neighbourWallInfo,
tol
);
}
}
// Update this with w2 if w2 nearer to pt.
@ -162,7 +162,7 @@ inline bool wallPointData<Type>::updateFace
neighbourWallInfo,
tol
);
}
}
// Update this with w2 if w2 nearer to pt.
@ -183,7 +183,7 @@ inline bool wallPointData<Type>::updateFace
neighbourWallInfo,
tol
);
}
}
// ************************************************************************* //

12
src/meshTools/cellDist/wallPoint/wallPointI.H
View File

@ -53,7 +53,7 @@ inline bool Foam::wallPoint::update
origin_ = w2.origin();
return true;
}
}
scalar diff = distSqr_ - dist2;
@ -77,7 +77,7 @@ inline bool Foam::wallPoint::update
return true;
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
@ -211,14 +211,14 @@ inline bool Foam::wallPoint::updateCell
const scalar tol
)
{
return
return
update
(
mesh.cellCentres()[thisCellI],
neighbourWallInfo,
tol
);
}
}
// Update this with w2 if w2 nearer to pt.
@ -238,7 +238,7 @@ inline bool Foam::wallPoint::updateFace
neighbourWallInfo,
tol
);
}
}
// Update this with w2 if w2 nearer to pt.
inline bool Foam::wallPoint::updateFace
@ -256,7 +256,7 @@ inline bool Foam::wallPoint::updateFace
neighbourWallInfo,
tol
);
}
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //

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

@ -705,7 +705,7 @@ Foam::treeBoundBox Foam::indexedOctree<Type>::subBbox
// Takes a bb and a point on/close to the edge of the bb and pushes the point
// inside by a small fraction.
// inside by a small fraction.
template <class Type>
Foam::point Foam::indexedOctree<Type>::pushPoint
(
@ -776,7 +776,7 @@ Foam::point Foam::indexedOctree<Type>::pushPoint
// Takes a bb and a point on the edge of the bb and pushes the point
// outside by a small fraction.
// outside by a small fraction.
template <class Type>
Foam::point Foam::indexedOctree<Type>::pushPoint
(
@ -957,6 +957,7 @@ Foam::point Foam::indexedOctree<Type>::pushPointIntoFace
}
else if (nFaces == 1)
{
// Point is on a single face
keepFaceID = faceIndices[0];
}
else
@ -1042,7 +1043,7 @@ Foam::point Foam::indexedOctree<Type>::pushPointIntoFace
}
return facePoint;
return facePoint;
}
@ -1150,7 +1151,7 @@ Foam::point Foam::indexedOctree<Type>::pushPointIntoFace
// // | |
// // +----+
// // Shift point down (away from top):
// //
// //
// // a+----+
// // ----| |
// // | |
@ -1364,7 +1365,7 @@ bool Foam::indexedOctree<Type>::walkToNeighbour
// | a+-+-+
// | |\| |
// +---+-+-+
// \
// \
//
// e.g. ray is at (a) in octant 0(or 4) with faceIDs : LEFTBIT+TOPBIT.
// If we would be in octant 1(or 5) we could go to the correct octant
@ -1782,16 +1783,6 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findLine
label i = 0;
for (; i < 100000; i++)
{
if (verbose)
{
Pout<< "iter:" << i
<< " at startPoint:" << hitInfo.rawPoint() << endl
<< " node:" << nodeI
<< " octant:" << octant
<< " bb:" << subBbox(nodeI, octant) << endl;
}
// Ray-trace to end of current node. Updates point (either on triangle
// in case of hit or on node bounding box in case of miss)
@ -1808,6 +1799,19 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findLine
)
);
if (verbose)
{
Pout<< "iter:" << i
<< " at current:" << hitInfo.rawPoint()
<< " (perturbed:" << startPoint << ")" << endl
<< " node:" << nodeI
<< " octant:" << octant
<< " bb:" << subBbox(nodeI, octant) << endl;
}
// Faces of current bounding box current point is on
direction hitFaceID = 0;
@ -1833,12 +1837,23 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findLine
break;
}
if (hitFaceID == 0)
if (hitFaceID == 0 || hitInfo.rawPoint() == treeEnd)
{
// endpoint inside the tree. Return miss.
break;
}
// Create a point on other side of face.
point perturbedPoint
(
pushPoint
(
octantBb,
hitFaceID,
hitInfo.rawPoint(),
false // push outside of octantBb
)
);
if (verbose)
{
@ -1848,14 +1863,7 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findLine
<< " node:" << nodeI
<< " octant:" << octant
<< " bb:" << subBbox(nodeI, octant) << nl
<< " walking to neighbour containing:"
<< pushPoint
(
octantBb,
hitFaceID,
hitInfo.rawPoint(),
false // push outside of octantBb
)
<< " walking to neighbour containing:" << perturbedPoint
<< endl;
}
@ -1866,13 +1874,7 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findLine
bool ok = walkToNeighbour
(
pushPoint
(
octantBb,
hitFaceID,
hitInfo.rawPoint(),
false // push outside of octantBb
),
perturbedPoint,
hitFaceID, // face(s) that hitInfo is on
nodeI,

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

@ -103,7 +103,7 @@ Foam::label Foam::meshSearch::findNearestCellTree(const point& location) const
const indexedOctree<treeDataPoint>& tree = cellCentreTree();
scalar span = tree.bb().mag();
pointIndexHit info = tree.findNearest(location, Foam::sqr(span));
if (!info.hit())
@ -463,7 +463,7 @@ const Foam::indexedOctree<Foam::treeDataFace>& Foam::meshSearch::boundaryTree()
treeBoundBox overallBb(mesh_.points());
Random rndGen(123456);
overallBb.extend(rndGen, 1E-4);
overallBb = overallBb.extend(rndGen, 1E-4);
overallBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
@ -497,7 +497,7 @@ const Foam::indexedOctree<Foam::treeDataCell>& Foam::meshSearch::cellTree()
treeBoundBox overallBb(mesh_.points());
Random rndGen(123456);
overallBb.extend(rndGen, 1E-4);
overallBb = overallBb.extend(rndGen, 1E-4);
overallBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
@ -516,7 +516,7 @@ const Foam::indexedOctree<Foam::treeDataCell>& Foam::meshSearch::cellTree()
}
return *cellTreePtr_;
}
@ -531,7 +531,7 @@ const Foam::indexedOctree<Foam::treeDataPoint>&
treeBoundBox overallBb(mesh_.cellCentres());
Random rndGen(123456);
overallBb.extend(rndGen, 1E-4);
overallBb = overallBb.extend(rndGen, 1E-4);
overallBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
@ -787,7 +787,7 @@ Foam::label Foam::meshSearch::findCell
{
nearCellI = mesh_.faceOwner()[curHit.index()];
curPoint =
curHit.hitPoint()
curHit.hitPoint()
+ offset(curHit.hitPoint(), curHit.index(), edgeVec);
}
}

2
src/meshTools/meshSearch/meshSearch.H
View File

@ -203,7 +203,7 @@ public:
//- Find nearest cell in terms of cell centre.
// - use octree
// - use linear search
// - if seed is provided walk. (uses findNearestCellWalk;
// - if seed is provided walk. (uses findNearestCellWalk;
// does not handle holes in domain)
label findNearestCell
(

2
src/meshTools/octree/PointIndexHit.H
View File

@ -27,7 +27,7 @@ Class
Description
This class describes the interaction of (usually) a face and a point.
It carries the info of a successful hit and (if successful),
It carries the info of a successful hit and (if successful),
returns the interaction point.
like pointHit but carries face (or cell, edge etc.) index

6
src/meshTools/octree/octreeDataCell.C
View File

@ -80,7 +80,7 @@ Foam::octreeDataCell::octreeDataCell
forAll(pointsi, pointi)
{
const point& p = points[pointsi[pointi]];
bbs_[celli].min() = min(bbs_[celli].min(), p);
bbs_[celli].max() = max(bbs_[celli].max(), p);
}
@ -169,7 +169,7 @@ bool Foam::octreeDataCell::findTightest
else
{
// Construct bb around sample and myFar
const point dist2(fabs(dist.x()), fabs(dist.y()), fabs(dist.z()));
const point dist2(fabs(dist.x()), fabs(dist.y()), fabs(dist.z()));
tightest.min() = sample - dist2;
tightest.max() = sample + dist2;
@ -223,7 +223,7 @@ Foam::scalar Foam::octreeDataCell::calcNearest
);
return GREAT;
}
void Foam::octreeDataCell::write
(

4
src/meshTools/octree/octreeDataEdges.C
View File

@ -169,7 +169,7 @@ bool Foam::octreeDataEdges::findTightest
else
{
// Construct bb around sample and myFar
const point dist2(fabs(dist.x()), fabs(dist.y()), fabs(dist.z()));
const point dist2(fabs(dist.x()), fabs(dist.y()), fabs(dist.z()));
tightest.min() = sample - dist2;
tightest.max() = sample + dist2;
@ -226,7 +226,7 @@ Foam::scalar Foam::octreeDataEdges::calcNearest
return edgeLine.nearestDist(sampleLine, shapePt, sampleLinePt);
}
void Foam::octreeDataEdges::write
(

2
src/meshTools/octree/octreeDataPointTreeLeaf.C
View File

@ -81,7 +81,7 @@ bool Foam::treeLeaf<Foam::octreeDataPoint>::findNearest
if (minIndex != -1)
{
tightestDist = sqrt(minDist2);
// New nearer. Update 'tightest' bounding box
tMin.x() = sample.x() - tightestDist;
tMin.y() = sample.y() - tightestDist;

2
src/meshTools/octree/treeBoundBox.H
View File

@ -254,7 +254,7 @@ public:
// second the octant with boundary nearest to the point etc.
inline void searchOrder
(
const point& pt,
const point& pt,
FixedList<direction, 8>& octantOrder
) const;

68
src/meshTools/octree/treeBoundBoxI.H
View File

@ -242,7 +242,7 @@ inline void Foam::treeBoundBox::searchOrder
if (dist.z() < 0)
{
octant |= treeBoundBox::FRONTHALF;
dist.z() *= -1;
dist.z() *= -1;
}
direction min = 0;
@ -252,44 +252,44 @@ inline void Foam::treeBoundBox::searchOrder
if (dist.x() < dist.y())
{
if (dist.y() < dist.z())
{
min = treeBoundBox::RIGHTHALF;
mid = treeBoundBox::TOPHALF;
max = treeBoundBox::FRONTHALF;
}
else if (dist.z() < dist.x())
{
min = treeBoundBox::FRONTHALF;
{
min = treeBoundBox::RIGHTHALF;
mid = treeBoundBox::TOPHALF;
max = treeBoundBox::FRONTHALF;
}
else if (dist.z() < dist.x())
{
min = treeBoundBox::FRONTHALF;
mid = treeBoundBox::RIGHTHALF;
max = treeBoundBox::TOPHALF;
}
else
{
min = treeBoundBox::RIGHTHALF;
mid = treeBoundBox::FRONTHALF;
max = treeBoundBox::TOPHALF;
}
max = treeBoundBox::TOPHALF;
}
else
{
min = treeBoundBox::RIGHTHALF;
mid = treeBoundBox::FRONTHALF;
max = treeBoundBox::TOPHALF;
}
}
else
{
if (dist.z() < dist.y())
{
min = treeBoundBox::FRONTHALF;
mid = treeBoundBox::TOPHALF;
max = treeBoundBox::RIGHTHALF;
}
else if (dist.x() < dist.z())
{
min = treeBoundBox::TOPHALF;
mid = treeBoundBox::RIGHTHALF;
max = treeBoundBox::FRONTHALF;
}
else
{
min = treeBoundBox::TOPHALF;
mid = treeBoundBox::FRONTHALF;
max = treeBoundBox::RIGHTHALF;
}
{
min = treeBoundBox::FRONTHALF;
mid = treeBoundBox::TOPHALF;
max = treeBoundBox::RIGHTHALF;
}
else if (dist.x() < dist.z())
{
min = treeBoundBox::TOPHALF;
mid = treeBoundBox::RIGHTHALF;
max = treeBoundBox::FRONTHALF;
}
else
{
min = treeBoundBox::TOPHALF;
mid = treeBoundBox::FRONTHALF;
max = treeBoundBox::RIGHTHALF;
}
}
// Primary subOctant

165
src/meshTools/searchableSurface/distributedTriSurfaceMesh.C
View File

@ -36,15 +36,19 @@ License
#include "IFstream.H"
#include "decompositionMethod.H"
#include "vectorList.H"
#include "PackedBoolList.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(distributedTriSurfaceMesh, 0);
addToRunTimeSelectionTable(searchableSurface, distributedTriSurfaceMesh, dict);
defineTypeNameAndDebug(distributedTriSurfaceMesh, 0);
addToRunTimeSelectionTable
(
searchableSurface,
distributedTriSurfaceMesh,
dict
);
}
@ -873,7 +877,7 @@ Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
bbs[procI].setSize(1);
//bbs[procI][0] = boundBox::invertedBox;
bbs[procI][0].min() = point( VGREAT, VGREAT, VGREAT);
bbs[procI][0].max() = point(-VGREAT, -VGREAT, -VGREAT);
bbs[procI][0].max() = point(-VGREAT, -VGREAT, -VGREAT);
}
forAll (s, triI)
@ -908,42 +912,6 @@ Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
}
void Foam::distributedTriSurfaceMesh::calcBounds
(
boundBox& bb,
label& nPoints
) const
{
// Unfortunately nPoints constructs meshPoints() so do compact version
// ourselves
PackedList<1> pointIsUsed(points().size());
pointIsUsed = 0U;
nPoints = 0;
bb.min() = point(VGREAT, VGREAT, VGREAT);
bb.max() = point(-VGREAT, -VGREAT, -VGREAT);
const triSurface& s = static_cast<const triSurface&>(*this);
forAll(s, triI)
{
const labelledTri& f = s[triI];
forAll(f, fp)
{
label pointI = f[fp];
if (pointIsUsed.set(pointI, 1))
{
bb.min() = ::Foam::min(bb.min(), points()[pointI]);
bb.max() = ::Foam::max(bb.max(), points()[pointI]);
nPoints++;
}
}
}
}
// Does any part of triangle overlap bb.
bool Foam::distributedTriSurfaceMesh::overlaps
(
@ -1932,20 +1900,7 @@ void Foam::distributedTriSurfaceMesh::getNormal
{
if (!Pstream::parRun())
{
normal.setSize(info.size());
forAll(info, i)
{
if (info[i].hit())
{
normal[i] = faceNormals()[info[i].index()];
}
else
{
// Set to what?
normal[i] = vector::zero;
}
}
triSurfaceMesh::getNormal(info, normal);
return;
}
@ -2003,70 +1958,64 @@ void Foam::distributedTriSurfaceMesh::getNormal
void Foam::distributedTriSurfaceMesh::getField
(
const word& fieldName,
const List<pointIndexHit>& info,
labelList& values
) const
{
const triSurfaceLabelField& fld = lookupObject<triSurfaceLabelField>
(
fieldName
);
if (!Pstream::parRun())
{
values.setSize(info.size());
forAll(info, i)
{
if (info[i].hit())
{
values[i] = fld[info[i].index()];
}
}
triSurfaceMesh::getField(info, values);
return;
}
// Get query data (= local index of triangle)
// ~~~~~~~~~~~~~~
labelList triangleIndex(info.size());
autoPtr<mapDistribute> mapPtr
(
calcLocalQueries
(
info,
triangleIndex
)
);
const mapDistribute& map = mapPtr();
// Do my tests
// ~~~~~~~~~~~
values.setSize(triangleIndex.size());
forAll(triangleIndex, i)
if (foundObject<triSurfaceLabelField>("values"))
{
label triI = triangleIndex[i];
values[i] = fld[triI];
const triSurfaceLabelField& fld = lookupObject<triSurfaceLabelField>
(
"values"
);
// Get query data (= local index of triangle)
// ~~~~~~~~~~~~~~
labelList triangleIndex(info.size());
autoPtr<mapDistribute> mapPtr
(
calcLocalQueries
(
info,
triangleIndex
)
);
const mapDistribute& map = mapPtr();
// Do my tests
// ~~~~~~~~~~~
values.setSize(triangleIndex.size());
forAll(triangleIndex, i)
{
label triI = triangleIndex[i];
values[i] = fld[triI];
}
// Send back results
// ~~~~~~~~~~~~~~~~~
map.distribute
(
Pstream::nonBlocking,
List<labelPair>(0),
info.size(),
map.constructMap(), // what to send
map.subMap(), // what to receive
values
);
}
// Send back results
// ~~~~~~~~~~~~~~~~~
map.distribute
(
Pstream::nonBlocking,
List<labelPair>(0),
info.size(),
map.constructMap(), // what to send
map.subMap(), // what to receive
values
);
}

14
src/meshTools/searchableSurface/distributedTriSurfaceMesh.H
View File

@ -217,9 +217,6 @@ private:
const triSurface&
);
//- Calculate surface bounding box
void calcBounds(boundBox& bb, label& nPoints) const;
//- Does any part of triangle overlap bb.
static bool overlaps
(
@ -418,7 +415,7 @@ public:
// Should really be split into a routine to determine decomposition
// and one that does actual distribution but determining
// decomposition with duplicate triangle merging requires
// same amoun as work as actual distribution.
// same amount as work as actual distribution.
virtual void distribute
(
const List<treeBoundBox>&,
@ -430,14 +427,9 @@ public:
// Other
//- Specific to triSurfaceMesh: from a set of hits (points and
//- WIP. From a set of hits (points and
// indices) get the specified field. Misses do not get set.
virtual void getField
(
const word& fieldName,
const List<pointIndexHit>&,
labelList& values
) const;
virtual void getField(const List<pointIndexHit>&, labelList&) const;
//- Subset the part of surface that is overlapping bounds.
static triSurface overlappingSurface

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

@ -229,6 +229,21 @@ const Foam::wordList& Foam::searchableBox::regions() const
}
Foam::pointField Foam::searchableBox::coordinates() const
{
pointField ctrs(6);
const pointField pts = treeBoundBox::points();
const faceList& fcs = treeBoundBox::faces;
forAll(fcs, i)
{
ctrs[i] = fcs[i].centre(pts);
}
return ctrs;
}
Foam::pointIndexHit Foam::searchableBox::findNearest
(
const point& sample,

3
src/meshTools/searchableSurface/searchableBox.H
View File

@ -127,6 +127,9 @@ public:
return 6;
}
//- Get representative set of element coordinates
// Usually the element centres (should be of length size()).
virtual pointField coordinates() const;
// Single point queries.

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

@ -40,6 +40,14 @@ addToRunTimeSelectionTable(searchableSurface, searchableCylinder, dict);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::pointField Foam::searchableCylinder::coordinates() const
{
pointField ctrs(1, 0.5*(point1_ + point2_));
return ctrs;
}
Foam::pointIndexHit Foam::searchableCylinder::findNearest
(
const point& sample,

4
src/meshTools/searchableSurface/searchableCylinder.H
View File

@ -150,6 +150,10 @@ public:
return 1;
}
//- Get representative set of element coordinates
// Usually the element centres (should be of length size()).
virtual pointField coordinates() const;
// Multiple point queries.

10
src/meshTools/searchableSurface/searchablePlane.H
View File

@ -26,7 +26,7 @@ Class
Foam::searchablePlane
Description
Searching on plane. See plane.H
Searching on (infinite) plane. See plane.H
SourceFiles
searchablePlane.C
@ -122,6 +122,14 @@ public:
return 1;
}
//- Get representative set of element coordinates
// Usually the element centres (should be of length size()).
virtual pointField coordinates() const
{
//notImplemented("searchablePlane::coordinates()")
return pointField(1, refPoint());
}
// Multiple point queries.

7
src/meshTools/searchableSurface/searchablePlate.H
View File

@ -145,6 +145,13 @@ public:
return 1;
}
//- Get representative set of element coordinates
// Usually the element centres (should be of length size()).
virtual pointField coordinates() const
{
return pointField(1, origin_);
}
// Multiple point queries.

7
src/meshTools/searchableSurface/searchableSphere.H
View File

@ -133,6 +133,13 @@ public:
return 1;
}
//- Get representative set of element coordinates
// Usually the element centres (should be of length size()).
virtual pointField coordinates() const
{
return pointField(1, centre_);
}
// Multiple point queries.

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

@ -190,6 +190,10 @@ public:
return size();
}
//- Get representative set of element coordinates
// Usually the element centres (should be of length size()).
virtual pointField coordinates() const = 0;
// Single point queries.
@ -319,6 +323,18 @@ public:
)
{}
//- WIP. Store element-wise field.
virtual void setField(const labelList& values)
{}
//- WIP. From a set of hits (points and
// indices) get the specified field. Misses do not get set. Return
// empty field if not supported.
virtual void getField(const List<pointIndexHit>&, labelList& values)
const
{
values.clear();
}
};

355
src/meshTools/searchableSurface/searchableSurfaceCollection.C
View File

@ -101,7 +101,11 @@ void Foam::searchableSurfaceCollection::findNearest
minDistSqr[pointI] = distSqr;
nearestInfo[pointI].setPoint(globalPt);
nearestInfo[pointI].setHit();
nearestInfo[pointI].setIndex(hitInfo[pointI].index());
nearestInfo[pointI].setIndex
(
hitInfo[pointI].index()
+ indexOffset_[surfI]
);
nearestSurf[pointI] = surfI;
}
}
@ -110,6 +114,62 @@ void Foam::searchableSurfaceCollection::findNearest
}
// Sort hits into per-surface bins. Misses are rejected. Maintains map back
// to position
void Foam::searchableSurfaceCollection::sortHits
(
const List<pointIndexHit>& info,
List<List<pointIndexHit> >& surfInfo,
labelListList& infoMap
) const
{
// Count hits per surface.
labelList nHits(subGeom_.size(), 0);
forAll(info, pointI)
{
if (info[pointI].hit())
{
label index = info[pointI].index();
label surfI = findLower(indexOffset_, index+1);
nHits[surfI]++;
}
}
// Per surface the hit
surfInfo.setSize(subGeom_.size());
// Per surface the original position
infoMap.setSize(subGeom_.size());
forAll(surfInfo, surfI)
{
surfInfo[surfI].setSize(nHits[surfI]);
infoMap[surfI].setSize(nHits[surfI]);
}
nHits = 0;
forAll(info, pointI)
{
if (info[pointI].hit())
{
label index = info[pointI].index();
label surfI = findLower(indexOffset_, index+1);
// Store for correct surface and adapt indices back to local
// ones
label localI = nHits[surfI]++;
surfInfo[surfI][localI] = pointIndexHit
(
info[pointI].hit(),
info[pointI].rawPoint(),
index-indexOffset_[surfI]
);
infoMap[surfI][localI] = pointI;
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::searchableSurfaceCollection::searchableSurfaceCollection
@ -123,11 +183,13 @@ Foam::searchableSurfaceCollection::searchableSurfaceCollection
scale_(dict.size()),
transform_(dict.size()),
subGeom_(dict.size()),
mergeSubRegions_(dict.lookup("mergeSubRegions"))
mergeSubRegions_(dict.lookup("mergeSubRegions")),
indexOffset_(dict.size()+1)
{
Info<< "SearchableCollection : " << name() << endl;
label surfI = 0;
label startIndex = 0;
forAllConstIter(dictionary, dict, iter)
{
if (dict.isDict(iter().keyword()))
@ -153,8 +215,24 @@ Foam::searchableSurfaceCollection::searchableSurfaceCollection
const searchableSurface& s =
io.db().lookupObject<searchableSurface>(subGeomName);
// I don't know yet how to handle the globalSize combined with
// regionOffset. Would cause non-consecutive indices locally
// if all indices offset by globalSize() of the local region...
if (s.size() != s.globalSize())
{
FatalErrorIn
(
"searchableSurfaceCollection::searchableSurfaceCollection"
"(const IOobject&, const dictionary&)"
) << "Cannot use a distributed surface in a collection."
<< exit(FatalError);
}
subGeom_.set(surfI, &const_cast<searchableSurface&>(s));
indexOffset_[surfI] = startIndex;
startIndex += subGeom_[surfI].size();
Info<< " instance : " << instance_[surfI] << endl;
Info<< " surface : " << s.name() << endl;
Info<< " scale : " << scale_[surfI] << endl;
@ -163,10 +241,13 @@ Foam::searchableSurfaceCollection::searchableSurfaceCollection
surfI++;
}
}
indexOffset_[surfI] = startIndex;
instance_.setSize(surfI);
scale_.setSize(surfI);
transform_.setSize(surfI);
subGeom_.setSize(surfI);
indexOffset_.setSize(surfI+1);
}
@ -212,12 +293,36 @@ const Foam::wordList& Foam::searchableSurfaceCollection::regions() const
Foam::label Foam::searchableSurfaceCollection::size() const
{
label n = 0;
return indexOffset_[indexOffset_.size()-1];
}
Foam::pointField Foam::searchableSurfaceCollection::coordinates() const
{
// Get overall size
pointField coords(size());
// Append individual coordinates
label coordI = 0;
forAll(subGeom_, surfI)
{
n += subGeom_[surfI].size();
const pointField subCoords = subGeom_[surfI].coordinates();
forAll(subCoords, i)
{
coords[coordI++] = transform_[surfI].globalPosition
(
cmptMultiply
(
subCoords[i],
scale_[surfI]
)
);
}
}
return n;
return coords;
}
@ -296,6 +401,11 @@ void Foam::searchableSurfaceCollection::findLine
);
info[pointI] = hitInfo[pointI];
info[pointI].rawPoint() = nearest[pointI];
info[pointI].setIndex
(
hitInfo[pointI].index()
+ indexOffset_[surfI]
);
}
}
}
@ -397,82 +507,42 @@ void Foam::searchableSurfaceCollection::getRegion
}
else
{
// Multiple surfaces. Sort by surface.
// Per surface the hit
List<List<pointIndexHit> > surfInfo;
// Per surface the original position
List<List<label> > infoMap;
sortHits(info, surfInfo, infoMap);
region.setSize(info.size());
region = -1;
// Which region did point come from. Retest for now to see which
// surface it originates from - crap solution! Should use global indices
// in index inside pointIndexHit to do this better.
// Do region tests
pointField samples(info.size());
forAll(info, pointI)
if (mergeSubRegions_)
{
if (info[pointI].hit())
// Actually no need for surfInfo. Just take region for surface.
forAll(infoMap, surfI)
{
samples[pointI] = info[pointI].hitPoint();
}
else
{
samples[pointI] = vector::zero;
}
}
//scalarField minDistSqr(info.size(), SMALL);
scalarField minDistSqr(info.size(), GREAT);
labelList nearestSurf;
List<pointIndexHit> nearestInfo;
findNearest
(
samples,
minDistSqr,
nearestInfo,
nearestSurf
);
// Check
{
forAll(info, pointI)
{
if (info[pointI].hit() && nearestSurf[pointI] == -1)
const labelList& map = infoMap[surfI];
forAll(map, i)
{
FatalErrorIn
(
"searchableSurfaceCollection::getRegion(..)"
) << "pointI:" << pointI
<< " sample:" << samples[pointI]
<< " nearest:" << nearestInfo[pointI]
<< " nearestsurf:" << nearestSurf[pointI]
<< abort(FatalError);
region[map[i]] = regionOffset_[surfI];
}
}
}
forAll(subGeom_, surfI)
else
{
// Collect points from my surface
labelList indices(findIndices(nearestSurf, surfI));
if (mergeSubRegions_)
{
forAll(indices, i)
{
region[indices[i]] = regionOffset_[surfI];
}
}
else
forAll(infoMap, surfI)
{
labelList surfRegion;
subGeom_[surfI].getRegion
(
List<pointIndexHit>
(
UIndirectList<pointIndexHit>(info, indices)
),
surfRegion
);
forAll(indices, i)
subGeom_[surfI].getRegion(surfInfo[surfI], surfRegion);
const labelList& map = infoMap[surfI];
forAll(map, i)
{
region[indices[i]] = regionOffset_[surfI] + surfRegion[i];
region[map[i]] = regionOffset_[surfI] + surfRegion[i];
}
}
}
@ -494,52 +564,26 @@ void Foam::searchableSurfaceCollection::getNormal
}
else
{
// Multiple surfaces. Sort by surface.
// Per surface the hit
List<List<pointIndexHit> > surfInfo;
// Per surface the original position
List<List<label> > infoMap;
sortHits(info, surfInfo, infoMap);
normal.setSize(info.size());
// See above - crap retest to find surface point originates from.
pointField samples(info.size());
forAll(info, pointI)
// Do region tests
forAll(surfInfo, surfI)
{
if (info[pointI].hit())
{
samples[pointI] = info[pointI].hitPoint();
}
else
{
samples[pointI] = vector::zero;
}
}
//scalarField minDistSqr(info.size(), SMALL);
scalarField minDistSqr(info.size(), GREAT);
labelList nearestSurf;
List<pointIndexHit> nearestInfo;
findNearest
(
samples,
minDistSqr,
nearestInfo,
nearestSurf
);
forAll(subGeom_, surfI)
{
// Collect points from my surface
labelList indices(findIndices(nearestSurf, surfI));
vectorField surfNormal;
subGeom_[surfI].getNormal
(
List<pointIndexHit>
(
UIndirectList<pointIndexHit>(info, indices)
),
surfNormal
);
forAll(indices, i)
subGeom_[surfI].getNormal(surfInfo[surfI], surfNormal);
const labelList& map = infoMap[surfI];
forAll(map, i)
{
normal[indices[i]] = surfNormal[i];
normal[map[i]] = surfNormal[i];
}
}
}
@ -561,4 +605,99 @@ void Foam::searchableSurfaceCollection::getVolumeType
}
void Foam::searchableSurfaceCollection::distribute
(
const List<treeBoundBox>& bbs,
const bool keepNonLocal,
autoPtr<mapDistribute>& faceMap,
autoPtr<mapDistribute>& pointMap
)
{
forAll(subGeom_, surfI)
{
// Note:Tranform the bounding boxes? Something like
// pointField bbPoints =
// cmptDivide
// (
// transform_[surfI].localPosition
// (
// bbs[i].points()
// ),
// scale_[surfI]
// );
// treeBoundBox newBb(bbPoints);
// Note: what to do with faceMap, pointMap from multiple surfaces?
subGeom_[surfI].distribute
(
bbs,
keepNonLocal,
faceMap,
pointMap
);
}
}
void Foam::searchableSurfaceCollection::setField(const labelList& values)
{
forAll(subGeom_, surfI)
{
subGeom_[surfI].setField
(
static_cast<const labelList&>
(
SubList<label>
(
values,
subGeom_[surfI].size(),
indexOffset_[surfI]
)
)
);
}
}
void Foam::searchableSurfaceCollection::getField
(
const List<pointIndexHit>& info,
labelList& values
) const
{
if (subGeom_.size() == 0)
{}
else if (subGeom_.size() == 1)
{
subGeom_[0].getField(info, values);
}
else
{
// Multiple surfaces. Sort by surface.
// Per surface the hit
List<List<pointIndexHit> > surfInfo;
// Per surface the original position
List<List<label> > infoMap;
sortHits(info, surfInfo, infoMap);
values.setSize(info.size());
//?Misses do not get set? values = 0;
// Do surface tests
forAll(surfInfo, surfI)
{
labelList surfValues;
subGeom_[surfI].getField(surfInfo[surfI], surfValues);
const labelList& map = infoMap[surfI];
forAll(map, i)
{
values[map[i]] = surfValues[i];
}
}
}
}
// ************************************************************************* //

38
src/meshTools/searchableSurface/searchableSurfaceCollection.H
View File

@ -77,6 +77,10 @@ private:
Switch mergeSubRegions_;
//- offsets for indices coming from different surfaces
// (sized with size() of each surface)
labelList indexOffset_;
//- Region names
mutable wordList regions_;
//- From individual regions to collection regions
@ -95,6 +99,15 @@ private:
labelList& nearestSurf
) const;
//- Sort hits into per-surface bins. Misses are rejected.
// Maintains map back to position
void sortHits
(
const List<pointIndexHit>& info,
List<List<pointIndexHit> >& surfInfo,
labelListList& infoMap
) const;
//- Disallow default bitwise copy construct
searchableSurfaceCollection(const searchableSurfaceCollection&);
@ -161,6 +174,10 @@ public:
//- Range of local indices that can be returned.
virtual label size() const;
//- Get representative set of element coordinates
// Usually the element centres (should be of length size()).
virtual pointField coordinates() const;
// Multiple point queries.
@ -215,6 +232,27 @@ public:
List<volumeType>&
) const;
// Other
//- Set bounds of surface. Bounds currently set as list of
// bounding boxes. The bounds are hints to the surface as for
// the range of queries it can expect. faceMap/pointMap can be
// set if the surface has done any redistribution.
virtual void distribute
(
const List<treeBoundBox>&,
const bool keepNonLocal,
autoPtr<mapDistribute>& faceMap,
autoPtr<mapDistribute>& pointMap
);
//- WIP. Store element-wise field.
virtual void setField(const labelList& values);
//- WIP. From a set of hits (points and
// indices) get the specified field. Misses do not get set. Return
// empty field if not supported.
virtual void getField(const List<pointIndexHit>&, labelList&) const;
// regIOobject implementation

2
src/meshTools/searchableSurface/searchableSurfaceWithGaps.C
View File

@ -84,7 +84,7 @@ Foam::Pair<Foam::vector> Foam::searchableSurfaceWithGaps::offsetVecs
offsets[1] *= gap_;
}
return offsets;
return offsets;
}

42
src/meshTools/searchableSurface/searchableSurfaceWithGaps.H
View File

@ -151,6 +151,13 @@ public:
return surface().size();
}
//- Get representative set of element coordinates
// Usually the element centres (should be of length size()).
virtual pointField coordinates() const
{
return surface().coordinates();
}
// Multiple point queries.
@ -225,6 +232,41 @@ public:
}
// Other
//- Set bounds of surface. Bounds currently set as list of
// bounding boxes. The bounds are hints to the surface as for
// the range of queries it can expect. faceMap/pointMap can be
// set if the surface has done any redistribution.
virtual void distribute
(
const List<treeBoundBox>& bbs,
const bool keepNonLocal,
autoPtr<mapDistribute>& faceMap,
autoPtr<mapDistribute>& pointMap
)
{
subGeom_[0].distribute(bbs, keepNonLocal, faceMap, pointMap);
}
//- WIP. Store element-wise field.
virtual void setField(const labelList& values)
{
subGeom_[0].setField(values);
}
//- WIP. From a set of hits (points and
// indices) get the specified field. Misses do not get set. Return
// empty field if not supported.
virtual void getField
(
const List<pointIndexHit>& info,
labelList& values
) const
{
surface().getField(info, values);
}
// regIOobject implementation
bool writeData(Ostream& os) const

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

@ -30,6 +30,7 @@ License
#include "EdgeMap.H"
#include "triSurfaceFields.H"
#include "Time.H"
#include "PackedBoolList.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -218,6 +219,99 @@ 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[hits.size()-1].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, triI)
{
const labelledTri& f = s[triI];
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)
@ -238,6 +332,7 @@ Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io, const triSurface& s)
),
triSurface(s),
tolerance_(indexedOctree<treeDataTriSurface>::perturbTol()),
maxTreeDepth_(10),
surfaceClosed_(-1)
{}
@ -281,6 +376,7 @@ Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io)
)
),
tolerance_(indexedOctree<treeDataTriSurface>::perturbTol()),
maxTreeDepth_(10),
surfaceClosed_(-1)
{}
@ -327,6 +423,7 @@ Foam::triSurfaceMesh::triSurfaceMesh
)
),
tolerance_(indexedOctree<treeDataTriSurface>::perturbTol()),
maxTreeDepth_(10),
surfaceClosed_(-1)
{
scalar scaleFactor = 0;
@ -347,6 +444,14 @@ Foam::triSurfaceMesh::triSurfaceMesh
Info<< searchableSurface::name() << " : using intersection tolerance "
<< tolerance_ << 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;
}
}
@ -368,6 +473,17 @@ void Foam::triSurfaceMesh::clearOut()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::pointField Foam::triSurfaceMesh::coordinates() const
{
// Use copy to calculate face centres so they don't get stored
return PrimitivePatch<labelledTri, SubList, const pointField&>
(
SubList<labelledTri>(*this, triSurface::size()),
triSurface::points()
).faceCentres();
}
void Foam::triSurfaceMesh::movePoints(const pointField& newPoints)
{
tree_.clear();
@ -381,15 +497,28 @@ const Foam::indexedOctree<Foam::treeDataTriSurface>&
{
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.
treeBoundBox bb
(
treeBoundBox(points(), meshPoints()).extend(rndGen, 1E-4)
);
bb = bb.extend(rndGen, 1E-4);
bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
@ -402,9 +531,9 @@ const Foam::indexedOctree<Foam::treeDataTriSurface>&
(
treeDataTriSurface(*this),
bb,
10, // maxLevel
10, // leafsize
3.0 // duplicity
maxTreeDepth_, // maxLevel
10, // leafsize
3.0 // duplicity
)
);
@ -431,15 +560,17 @@ const Foam::indexedOctree<Foam::treeDataEdge>&
+ nInternalEdges()
);
treeBoundBox bb;
label nPoints;
calcBounds(bb, nPoints);
// 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.
treeBoundBox bb
(
treeBoundBox(points(), meshPoints()).extend(rndGen, 1E-4)
);
bb = bb.extend(rndGen, 1E-4);
bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
@ -454,10 +585,10 @@ const Foam::indexedOctree<Foam::treeDataEdge>&
localPoints(), // points
bEdges // selected edges
),
bb, // bb
8, // maxLevel
10, // leafsize
3.0 // duplicity
bb, // bb
maxTreeDepth_, // maxLevel
10, // leafsize
3.0 // duplicity
)
);
}
@ -620,27 +751,14 @@ void Foam::triSurfaceMesh::findLineAll
hits.clear();
hits.append(inter);
point pt = inter.hitPoint() + smallVec[pointI];
while (((pt-start[pointI])&dirVec[pointI]) <= magSqrDirVec[pointI])
{
// See if any intersection between pt and end
pointIndexHit inter = octree.findLine(pt, end[pointI]);
// Check for not hit or hit same triangle as before (can happen
// if vector along surface of triangle)
if
(
!inter.hit()
|| (inter.index() == hits.last().index())
)
{
break;
}
hits.append(inter);
pt = inter.hitPoint() + smallVec[pointI];
}
getNextIntersections
(
octree,
start[pointI],
end[pointI],
smallVec[pointI],
hits
);
info[pointI].transfer(hits);
}
@ -704,24 +822,53 @@ void Foam::triSurfaceMesh::getNormal
}
void Foam::triSurfaceMesh::setField(const labelList& values)
{
autoPtr<triSurfaceLabelField> fldPtr
(
new triSurfaceLabelField
(
IOobject
(
"values",
objectRegistry::time().timeName(), // instance
"triSurface", // local
*this,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
*this,
dimless,
labelField(values)
)
);
// Store field on triMesh
fldPtr.ptr()->store();
}
void Foam::triSurfaceMesh::getField
(
const word& fieldName,
const List<pointIndexHit>& info,
labelList& values
) const
{
const triSurfaceLabelField& fld = lookupObject<triSurfaceLabelField>
(
fieldName
);
values.setSize(info.size());
forAll(info, i)
if (foundObject<triSurfaceLabelField>("values"))
{
if (info[i].hit())
values.setSize(info.size());
const triSurfaceLabelField& fld = lookupObject<triSurfaceLabelField>
(
"values"
);
forAll(info, i)
{
values[i] = fld[info[i].index()];
if (info[i].hit())
{
values[i] = fld[info[i].index()];
}
}
}
}

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

@ -71,6 +71,9 @@ private:
//- Optional tolerance to use in searches
scalar tolerance_;
//- Optional max tree depth of octree
label maxTreeDepth_;
//- Search tree (triangles)
mutable autoPtr<indexedOctree<treeDataTriSurface> > tree_;
@ -111,6 +114,17 @@ private:
// addressing.
bool isSurfaceClosed() const;
//- Steps to next intersection. Adds smallVec and starts tracking
// from there.
static void getNextIntersections
(
const indexedOctree<treeDataTriSurface>& octree,
const point& start,
const point& end,
const vector& smallVec,
DynamicList<pointIndexHit, 1, 1>& hits
);
//- Disallow default bitwise copy construct
triSurfaceMesh(const triSurfaceMesh&);
@ -118,6 +132,11 @@ 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
@ -174,6 +193,10 @@ public:
return triSurface::size();
}
//- Get representative set of element coordinates
// Usually the element centres (should be of length size()).
virtual pointField coordinates() const;
virtual void findNearest
(
const pointField& sample,
@ -225,6 +248,8 @@ public:
List<volumeType>&
) const;
// Other
//- Set bounds of surface. Bounds currently set as list of
// bounding boxes. The bounds are hints to the surface as for
// the range of queries it can expect. faceMap/pointMap can be
@ -238,16 +263,12 @@ public:
)
{}
// Other
//- WIP. Store element-wise field.
virtual void setField(const labelList& values);
//- Specific to triSurfaceMesh: from a set of hits (points and
//- WIP. From a set of hits (points and
// indices) get the specified field. Misses do not get set.
virtual void getField
(
const word& fieldName,
const List<pointIndexHit>&,
labelList& values
) const;
virtual void getField(const List<pointIndexHit>&, labelList&) const;
// regIOobject implementation

2
src/meshTools/sets/cellSources/cellToCell/cellToCell.C
View File

@ -70,7 +70,7 @@ Foam::cellToCell::cellToCell
Foam::cellToCell::cellToCell
(
const polyMesh& mesh,
const dictionary& dict
const dictionary& dict
)
:
topoSetSource(mesh),

2
src/meshTools/sets/cellSources/faceToCell/faceToCell.C
View File

@ -187,7 +187,7 @@ Foam::faceToCell::faceToCell
Foam::faceToCell::faceToCell
(
const polyMesh& mesh,
const dictionary& dict
const dictionary& dict
)
:
topoSetSource(mesh),

2
src/meshTools/sets/cellSources/nbrToCell/nbrToCell.C
View File

@ -61,7 +61,7 @@ void Foam::nbrToCell::combine(topoSet& set, const bool add) const
boolList isCoupled(mesh_.nFaces()-mesh_.nInternalFaces(), false);
forAll(patches, patchI)
{
{
const polyPatch& pp = patches[patchI];
if (pp.coupled())

4
src/meshTools/sets/cellSources/rotatedBoxToCell/rotatedBoxToCell.C
View File

@ -75,7 +75,7 @@ void Foam::rotatedBoxToCell::combine(topoSet& set, const bool add) const
}
const cellModel& hex = *(cellModeller::lookup("hex"));
// Get outwards pointing faces.
faceList boxFaces(cellShape(hex, boxVerts).faces());
@ -87,7 +87,7 @@ void Foam::rotatedBoxToCell::combine(topoSet& set, const bool add) const
Pout<< "Face:" << i << " position:" << boxFaces[i].centre(boxPoints)
<< " normal:" << boxFaceNormals[i] << endl;
}
}
// Check whether cell centre is inside all faces of box.

4
src/meshTools/sets/cellSources/surfaceToCell/surfaceToCell.C
View File

@ -107,7 +107,7 @@ bool Foam::surfaceToCell::differingPointNormals
const vector& span, // current search span
const label cellI,
const label cellTriI, // nearest (to cell centre) surface triangle
const label cellTriI, // nearest (to cell centre) surface triangle
Map<label>& pointToNearest // cache for nearest triangle to point
) const
@ -394,7 +394,7 @@ Foam::surfaceToCell::surfaceToCell
Foam::surfaceToCell::surfaceToCell
(
const polyMesh& mesh,
const dictionary& dict
const dictionary& dict
)
:
topoSetSource(mesh),

2
src/meshTools/sets/cellSources/surfaceToCell/surfaceToCell.H
View File

@ -131,7 +131,7 @@ class surfaceToCell
//- Depending on surface add to or delete from cellSet.
void combine(topoSet& set, const bool add) const;
//- Check values at construction time.
//- Check values at construction time.
void checkSettings() const;
const triSurfaceSearch& querySurf() const

4
src/meshTools/sets/faceSources/cellToFace/cellToFace.C
View File

@ -74,7 +74,7 @@ void Foam::cellToFace::combine(topoSet& set, const bool add) const
SeriousError<< "Cannot load set "
<< setName_ << endl;
}
cellSet loadedSet(mesh_, setName_);
if (option_ == ALL)
@ -179,7 +179,7 @@ Foam::cellToFace::cellToFace
Foam::cellToFace::cellToFace
(
const polyMesh& mesh,
const dictionary& dict
const dictionary& dict
)
:
topoSetSource(mesh),

2
src/meshTools/sets/pointSources/cellToPoint/cellToPoint.C
View File

@ -126,7 +126,7 @@ Foam::cellToPoint::cellToPoint
Foam::cellToPoint::cellToPoint
(
const polyMesh& mesh,
Istream& is
Istream& is
)
:
topoSetSource(mesh),

2
src/meshTools/sets/pointSources/faceToPoint/faceToPoint.C
View File

@ -106,7 +106,7 @@ Foam::faceToPoint::faceToPoint
Foam::faceToPoint::faceToPoint
(
const polyMesh& mesh,
const dictionary& dict
const dictionary& dict
)
:
topoSetSource(mesh),

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

@ -83,7 +83,7 @@ void Foam::surfaceToPoint::combine(topoSet& set, const bool add) const
if ((isInside && includeInside_) || (!isInside && includeOutside_))
{
addOrDelete(set, pointI, add);
}
}
}
}
@ -149,7 +149,7 @@ Foam::surfaceToPoint::surfaceToPoint
Foam::surfaceToPoint::surfaceToPoint
(
const polyMesh& mesh,
const dictionary& dict
const dictionary& dict
)
:
topoSetSource(mesh),

2
src/meshTools/sets/topoSets/cellSet.H
View File

@ -154,7 +154,7 @@ public:
//- Update any stored data for new labels
virtual void updateMesh(const mapPolyMesh& morphMap);
//- Write maxLen items with label and coordinates.
//- Write maxLen items with label and coordinates.
virtual void writeDebug
(
Ostream& os,

2
src/meshTools/sets/topoSets/cellZoneSet.H
View File

@ -131,7 +131,7 @@ public:
//- Sync cellZoneSet across coupled patches.
virtual void sync(const polyMesh& mesh);
//- Write maxLen items with label and coordinates.
//- Write maxLen items with label and coordinates.
virtual void writeDebug
(
Ostream& os,

2
src/meshTools/sets/topoSets/faceSet.C
View File

@ -153,7 +153,7 @@ void faceSet::sync(const polyMesh& mesh)
}
}
// Receive
// Receive
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];

2
src/meshTools/sets/topoSets/faceSet.H
View File

@ -118,7 +118,7 @@ public:
//- Update any stored data for new labels
virtual void updateMesh(const mapPolyMesh& morphMap);
//- Write maxLen items with label and coordinates.
//- Write maxLen items with label and coordinates.
virtual void writeDebug
(
Ostream& os,

2
src/meshTools/sets/topoSets/faceZoneSet.H
View File

@ -148,7 +148,7 @@ public:
//- Sync faceZoneSet across coupled patches.
virtual void sync(const polyMesh& mesh);
//- Write maxLen items with label and coordinates.
//- Write maxLen items with label and coordinates.
virtual void writeDebug
(
Ostream& os,

2
src/meshTools/sets/topoSets/pointSet.H
View File

@ -121,7 +121,7 @@ public:
//- Update any stored data for new labels
//virtual void updateMesh(const polyTopoChange& meshMod);
//- Write maxLen items with label and coordinates.
//- Write maxLen items with label and coordinates.
virtual void writeDebug
(
Ostream& os,

2
src/meshTools/sets/topoSets/pointZoneSet.H
View File

@ -134,7 +134,7 @@ public:
//- Sync pointZoneSet across coupled patches.
virtual void sync(const polyMesh& mesh);
//- Write maxLen items with label and coordinates.
//- Write maxLen items with label and coordinates.
virtual void writeDebug
(
Ostream& os,

10
src/meshTools/sets/topoSets/topoSet.C
View File

@ -236,7 +236,7 @@ void topoSet::writeDebug
topoSet::const_iterator& iter,
label& elemI
) const
{
{
label n = 0;
for (; (iter != end()) && (n < maxElem); ++iter)
@ -262,7 +262,7 @@ void topoSet::writeDebug
topoSet::const_iterator& iter,
label& elemI
) const
{
{
label n = 0;
for (; (iter != end()) && (n < maxElem); ++iter)
@ -329,7 +329,7 @@ topoSet::topoSet(const IOobject& obj, const word& wantedType)
:
regIOobject(obj)
{
if
if
(
readOpt() == IOobject::MUST_READ
|| (
@ -370,7 +370,7 @@ topoSet::topoSet
)
)
{
if
if
(
readOpt() == IOobject::MUST_READ
|| (
@ -475,7 +475,7 @@ void topoSet::invert(const label maxLen)
insert(cellI);
}
}
}

14
src/meshTools/surfaceSets/surfaceSets.C
View File

@ -65,7 +65,7 @@ License
//}
//
//
//// Returns true if cell uses at least one selected point
//// Returns true if cell uses at least one selected point
//bool Foam::surfaceSets::usesPoint
//(
// const primitiveMesh& mesh,
@ -147,7 +147,7 @@ License
// newPoints[pointI] = tnearest()[0];
// }
// }
// }
// }
// }
//
//
@ -161,7 +161,7 @@ License
// )
// {
// label cellI = iter.key();
//
//
// const cell& cll = cells[cellI];
//
// scalar newVol = cll.mag(newPoints, faces);
@ -284,7 +284,7 @@ void Foam::surfaceSets::getSurfaceSets
cellType.trimCutCells
(
nCutLayers,
cellClassification::OUTSIDE,
cellClassification::OUTSIDE,
cellClassification::INSIDE
);
}
@ -503,7 +503,7 @@ Foam::labelHashSet Foam::surfaceSets::getHangingCells
// mesh,
// surfName,
// surf,
// querySurf,
// querySurf,
// outsidePts,
// false, // includeCut
// false, // includeInside
@ -626,7 +626,7 @@ Foam::labelHashSet Foam::surfaceSets::getHangingCells
// << endl;
// }
// while (nRemoved != 0);
//
//
//
// //
// // Write
@ -642,7 +642,7 @@ Foam::labelHashSet Foam::surfaceSets::getHangingCells
// Pout<< "Writing outside cells (" << outside.size() << ") to cellSet "
// << outside.instance()/outside.local()/outside.name()
// << endl << endl;
//
//
// outside.write();
//}

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

@ -79,7 +79,7 @@ class surfaceSets
////- Min length of connected mesh edges
//static scalar minEdgeLen(const primitiveMesh& mesh, const label pointI);
//
////- Returns true if cell uses at least one selected point
////- Returns true if cell uses at least one selected point
//static bool usesPoint
//(
// const primitiveMesh& mesh,
@ -93,7 +93,7 @@ class surfaceSets
//// removed cells.
//static label removeHangingCells
//(
// const primitiveMesh&,
// const primitiveMesh&,
// const triSurfaceSearch& querySurf,
// labelHashSet& internalCells
//);

2
src/meshTools/triSurface/booleanOps/intersectedSurface/intersectedSurface.H
View File

@ -36,7 +36,7 @@ Description
- split the edges of the surface
- construct a new edgeList with (in this order) the edges from the
intersection ('cuts', i.e. the edges shared with the other surface)
and the (split) edges from the original triangles (from 0 ..
and the (split) edges from the original triangles (from 0 ..
nSurfaceEdges)
- construct face-edge addressing for above edges
- for each face do a right-handed walk to reconstruct faces (splitFace)

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

@ -181,7 +181,7 @@ public:
treeBoundBox& tightest
) const;
//- Given index get unit normal and calculate (numerical) sign
//- Given index get unit normal and calculate (numerical) sign
// of sample.
// Used to determine accuracy of calcNearest or inside/outside.
scalar calcSign

8
src/meshTools/triSurface/surfaceLocation/surfaceLocation.H
View File

@ -85,7 +85,7 @@ public:
//- Construct from components
surfaceLocation
(
(
const pointIndexHit& pih,
const triPointRef::proxType elementType,
const label triangle
@ -111,7 +111,7 @@ public:
{
return elementType_;
}
triPointRef::proxType elementType() const
{
return elementType_;
@ -121,7 +121,7 @@ public:
{
return triangle_;
}
label triangle() const
{
return triangle_;
@ -150,7 +150,7 @@ public:
friend Ostream& operator<<
(
Ostream&,
const InfoProxy<surfaceLocation>&
const InfoProxy<surfaceLocation>&
);
};

20
src/meshTools/triSurface/triSurfaceTools/geompack/geompack.H
View File

@ -5,13 +5,13 @@ void d2vec_part_quick_a ( int n, double a[], int *l, int *r );
void d2vec_permute ( int n, double a[], int p[] );
int *d2vec_sort_heap_index_a ( int n, double a[] );
void d2vec_sort_quick_a ( int n, double a[] );
int diaedg ( double x0, double y0, double x1, double y1, double x2, double y2,
int diaedg ( double x0, double y0, double x1, double y1, double x2, double y2,
double x3, double y3 );
void dmat_transpose_print ( int m, int n, double a[], const char *title );
void dmat_transpose_print_some ( int m, int n, double a[], int ilo, int jlo,
void dmat_transpose_print_some ( int m, int n, double a[], int ilo, int jlo,
int ihi, int jhi, const char *title );
void dmat_uniform ( int m, int n, double b, double c, int *seed, double r[] );
int dtris2 ( int point_num, double point_xy[], int *tri_num,
int dtris2 ( int point_num, double point_xy[], int *tri_num,
int tri_vert[], int tri_nabe[] );
bool dvec_eq ( int n, double a1[], double a2[] );
bool dvec_gt ( int n, double a1[], double a2[] );
@ -24,28 +24,28 @@ int i_modp ( int i, int j );
int i_sign ( int i );
int i_wrap ( int ival, int ilo, int ihi );
void imat_transpose_print ( int m, int n, int a[], const char *title );
void imat_transpose_print_some ( int m, int n, int a[], int ilo, int jlo,
void imat_transpose_print_some ( int m, int n, int a[], int ilo, int jlo,
int ihi, int jhi, const char *title );
void ivec_heap_d ( int n, int a[] );
int *ivec_indicator ( int n );
void ivec_sort_heap_a ( int n, int a[] );
void ivec_sorted_unique ( int n, int a[], int *nuniq );
int lrline ( double xu, double yu, double xv1, double yv1, double xv2, double yv2,
int lrline ( double xu, double yu, double xv1, double yv1, double xv2, double yv2,
double dv );
bool perm_check ( int n, int p[] );
void perm_inv ( int n, int p[] );
int *points_delaunay_naive_2d ( int n, double p[], int *ntri );
int s_len_trim ( const char *s );
int swapec ( int i, int *top, int *btri, int *bedg, int point_num,
double point_xy[], int tri_num, int tri_vert[], int tri_nabe[],
int swapec ( int i, int *top, int *btri, int *bedg, int point_num,
double point_xy[], int tri_num, int tri_vert[], int tri_nabe[],
int stack[] );
void timestamp ( void );
char *timestring ( void );
double *triangle_circumcenter_2d ( double t[] );
bool triangulation_plot_eps ( const char *file_out_name, int g_num, double g_xy[],
bool triangulation_plot_eps ( const char *file_out_name, int g_num, double g_xy[],
int tri_num, int nod_tri[] );
void triangulation_print ( int point_num, double xc[], int tri_num,
void triangulation_print ( int point_num, double xc[], int tri_num,
int tri_vert[], int tri_nabe[] );
void vbedg ( double x, double y, int point_num, double point_xy[], int tri_num,
void vbedg ( double x, double y, int point_num, double point_xy[], int tri_num,
int tri_vert[], int tri_nabe[], int *ltri, int *ledg, int *rtri, int *redg );

2
src/meshTools/triSurface/triangleFuncs/triangleFuncs.C
View File

@ -297,7 +297,7 @@ bool Foam::triangleFuncs::intersectBb
// return true;
// }
// if (tri.intersectionExact(cube3, cube0).hit())
// {
// {
// return true;
// }
//

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

@ -27,7 +27,7 @@ Class
Description
Class applies a two-dimensional correction to mesh motion point field.
The correction guarantees that the mesh does not get twisted during motion
and thus introduce a third dimension into a 2-D problem.