mirror of
https://develop.openfoam.com/Development/openfoam.git
synced 2025-11-28 03:28:01 +00:00
789 lines
24 KiB
C
789 lines
24 KiB
C
/*---------------------------------------------------------------------------*\
|
|
========= |
|
|
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
|
\\ / O peration |
|
|
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
|
|
\\/ M anipulation |
|
|
-------------------------------------------------------------------------------
|
|
License
|
|
This file is part of OpenFOAM.
|
|
|
|
OpenFOAM is free software: you can redistribute it and/or modify it
|
|
under the terms of the GNU General Public License as published by
|
|
the Free Software Foundation, either version 3 of the License, or
|
|
(at your option) any later version.
|
|
|
|
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
|
|
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
|
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
|
for more details.
|
|
|
|
You should have received a copy of the GNU General Public License
|
|
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
|
|
|
|
\*---------------------------------------------------------------------------*/
|
|
|
|
#include "refinementFeatures.H"
|
|
#include "Time.H"
|
|
#include "Tuple2.H"
|
|
#include "DynamicField.H"
|
|
|
|
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
|
|
|
|
void Foam::refinementFeatures::read
|
|
(
|
|
const objectRegistry& io,
|
|
const PtrList<dictionary>& featDicts
|
|
)
|
|
{
|
|
forAll(featDicts, featI)
|
|
{
|
|
const dictionary& dict = featDicts[featI];
|
|
|
|
fileName featFileName(dict.lookup("file"));
|
|
|
|
|
|
// Try reading extendedEdgeMesh first
|
|
|
|
IOobject extFeatObj
|
|
(
|
|
featFileName, // name
|
|
io.time().constant(), // instance
|
|
"extendedFeatureEdgeMesh", // local
|
|
io.time(), // registry
|
|
IOobject::MUST_READ,
|
|
IOobject::NO_WRITE,
|
|
false
|
|
);
|
|
|
|
const fileName fName(extFeatObj.filePath());
|
|
|
|
if (!fName.empty() && extendedEdgeMesh::canRead(fName))
|
|
{
|
|
autoPtr<extendedEdgeMesh> eMeshPtr = extendedEdgeMesh::New
|
|
(
|
|
fName
|
|
);
|
|
|
|
Info<< "Read extendedFeatureEdgeMesh " << extFeatObj.name()
|
|
<< nl << incrIndent;
|
|
eMeshPtr().writeStats(Info);
|
|
Info<< decrIndent << endl;
|
|
|
|
set(featI, new extendedFeatureEdgeMesh(extFeatObj, eMeshPtr()));
|
|
}
|
|
else
|
|
{
|
|
// Try reading edgeMesh
|
|
|
|
IOobject featObj
|
|
(
|
|
featFileName, // name
|
|
io.time().constant(), // instance
|
|
"triSurface", // local
|
|
io.time(), // registry
|
|
IOobject::MUST_READ,
|
|
IOobject::NO_WRITE,
|
|
false
|
|
);
|
|
|
|
const fileName fName(featObj.filePath());
|
|
|
|
if (fName.empty())
|
|
{
|
|
FatalIOErrorIn
|
|
(
|
|
"refinementFeatures::read"
|
|
"(const objectRegistry&"
|
|
", const PtrList<dictionary>&)",
|
|
dict
|
|
) << "Could not open " << featObj.objectPath()
|
|
<< exit(FatalIOError);
|
|
}
|
|
|
|
|
|
// Read as edgeMesh
|
|
autoPtr<edgeMesh> eMeshPtr = edgeMesh::New(fName);
|
|
const edgeMesh& eMesh = eMeshPtr();
|
|
|
|
Info<< "Read edgeMesh " << featObj.name() << nl
|
|
<< incrIndent;
|
|
eMesh.writeStats(Info);
|
|
Info<< decrIndent << endl;
|
|
|
|
|
|
// Analyse for feature points. These are all classified as mixed
|
|
// points for lack of anything better
|
|
const labelListList& pointEdges = eMesh.pointEdges();
|
|
|
|
labelList oldToNew(eMesh.points().size(), -1);
|
|
DynamicField<point> newPoints(eMesh.points().size());
|
|
forAll(pointEdges, pointI)
|
|
{
|
|
if (pointEdges[pointI].size() > 2)
|
|
{
|
|
oldToNew[pointI] = newPoints.size();
|
|
newPoints.append(eMesh.points()[pointI]);
|
|
}
|
|
//else if (pointEdges[pointI].size() == 2)
|
|
//MEJ: do something based on a feature angle?
|
|
}
|
|
label nFeatures = newPoints.size();
|
|
forAll(oldToNew, pointI)
|
|
{
|
|
if (oldToNew[pointI] == -1)
|
|
{
|
|
oldToNew[pointI] = newPoints.size();
|
|
newPoints.append(eMesh.points()[pointI]);
|
|
}
|
|
}
|
|
|
|
|
|
const edgeList& edges = eMesh.edges();
|
|
edgeList newEdges(edges.size());
|
|
forAll(edges, edgeI)
|
|
{
|
|
const edge& e = edges[edgeI];
|
|
newEdges[edgeI] = edge
|
|
(
|
|
oldToNew[e[0]],
|
|
oldToNew[e[1]]
|
|
);
|
|
}
|
|
|
|
// Construct an extendedEdgeMesh with
|
|
// - all points on more than 2 edges : mixed feature points
|
|
// - all edges : external edges
|
|
|
|
extendedEdgeMesh eeMesh
|
|
(
|
|
newPoints, // pts
|
|
newEdges, // eds
|
|
0, // (point) concaveStart
|
|
0, // (point) mixedStart
|
|
nFeatures, // (point) nonFeatureStart
|
|
edges.size(), // (edge) internalStart
|
|
edges.size(), // (edge) flatStart
|
|
edges.size(), // (edge) openStart
|
|
edges.size(), // (edge) multipleStart
|
|
vectorField(0), // normals
|
|
List<extendedEdgeMesh::sideVolumeType>(0),// normalVolumeTypes
|
|
vectorField(0), // edgeDirections
|
|
labelListList(0), // normalDirections
|
|
labelListList(0), // edgeNormals
|
|
labelListList(0), // featurePointNormals
|
|
labelListList(0), // featurePointEdges
|
|
identity(newEdges.size()) // regionEdges
|
|
);
|
|
|
|
//Info<< "Constructed extendedFeatureEdgeMesh " << featObj.name()
|
|
// << nl << incrIndent;
|
|
//eeMesh.writeStats(Info);
|
|
//Info<< decrIndent << endl;
|
|
|
|
set(featI, new extendedFeatureEdgeMesh(featObj, eeMesh));
|
|
}
|
|
|
|
const extendedEdgeMesh& eMesh = operator[](featI);
|
|
|
|
if (dict.found("levels"))
|
|
{
|
|
List<Tuple2<scalar, label> > distLevels(dict["levels"]);
|
|
|
|
if (dict.size() < 1)
|
|
{
|
|
FatalErrorIn
|
|
(
|
|
"refinementFeatures::read"
|
|
"(const objectRegistry&"
|
|
", const PtrList<dictionary>&)"
|
|
) << " : levels should be at least size 1" << endl
|
|
<< "levels : " << dict["levels"]
|
|
<< exit(FatalError);
|
|
}
|
|
|
|
distances_[featI].setSize(distLevels.size());
|
|
levels_[featI].setSize(distLevels.size());
|
|
|
|
forAll(distLevels, j)
|
|
{
|
|
distances_[featI][j] = distLevels[j].first();
|
|
levels_[featI][j] = distLevels[j].second();
|
|
|
|
// Check in incremental order
|
|
if (j > 0)
|
|
{
|
|
if
|
|
(
|
|
(distances_[featI][j] <= distances_[featI][j-1])
|
|
|| (levels_[featI][j] > levels_[featI][j-1])
|
|
)
|
|
{
|
|
FatalErrorIn
|
|
(
|
|
"refinementFeatures::read"
|
|
"(const objectRegistry&"
|
|
", const PtrList<dictionary>&)"
|
|
) << " : Refinement should be specified in order"
|
|
<< " of increasing distance"
|
|
<< " (and decreasing refinement level)." << endl
|
|
<< "Distance:" << distances_[featI][j]
|
|
<< " refinementLevel:" << levels_[featI][j]
|
|
<< exit(FatalError);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
// Look up 'level' for single level
|
|
levels_[featI] = labelList(1, readLabel(dict.lookup("level")));
|
|
distances_[featI] = scalarField(1, 0.0);
|
|
}
|
|
|
|
Info<< "Refinement level according to distance to "
|
|
<< featFileName << " (" << eMesh.points().size() << " points, "
|
|
<< eMesh.edges().size() << " edges)." << endl;
|
|
forAll(levels_[featI], j)
|
|
{
|
|
Info<< " level " << levels_[featI][j]
|
|
<< " for all cells within " << distances_[featI][j]
|
|
<< " metre." << endl;
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
void Foam::refinementFeatures::buildTrees(const label featI)
|
|
{
|
|
const extendedEdgeMesh& eMesh = operator[](featI);
|
|
const pointField& points = eMesh.points();
|
|
const edgeList& edges = eMesh.edges();
|
|
|
|
// Calculate bb of all points
|
|
treeBoundBox bb(points);
|
|
|
|
// Random number generator. Bit dodgy since not exactly random ;-)
|
|
Random rndGen(65431);
|
|
|
|
// Slightly extended bb. Slightly off-centred just so on symmetric
|
|
// geometry there are less face/edge aligned items.
|
|
bb = bb.extend(rndGen, 1e-4);
|
|
bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
|
|
bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
|
|
|
|
edgeTrees_.set
|
|
(
|
|
featI,
|
|
new indexedOctree<treeDataEdge>
|
|
(
|
|
treeDataEdge
|
|
(
|
|
false, // do not cache bb
|
|
edges,
|
|
points,
|
|
identity(edges.size())
|
|
),
|
|
bb, // overall search domain
|
|
8, // maxLevel
|
|
10, // leafsize
|
|
3.0 // duplicity
|
|
)
|
|
);
|
|
|
|
|
|
labelList featurePoints(identity(eMesh.nonFeatureStart()));
|
|
|
|
pointTrees_.set
|
|
(
|
|
featI,
|
|
new indexedOctree<treeDataPoint>
|
|
(
|
|
treeDataPoint(points, featurePoints),
|
|
bb, // overall search domain
|
|
8, // maxLevel
|
|
10, // leafsize
|
|
3.0 // duplicity
|
|
)
|
|
);
|
|
}
|
|
|
|
|
|
const Foam::PtrList<Foam::indexedOctree<Foam::treeDataEdge> >&
|
|
Foam::refinementFeatures::regionEdgeTrees() const
|
|
{
|
|
if (!regionEdgeTreesPtr_.valid())
|
|
{
|
|
regionEdgeTreesPtr_.reset
|
|
(
|
|
new PtrList<indexedOctree<treeDataEdge> >(size())
|
|
);
|
|
PtrList<indexedOctree<treeDataEdge> >& trees = regionEdgeTreesPtr_();
|
|
|
|
forAll(*this, featI)
|
|
{
|
|
const extendedEdgeMesh& eMesh = operator[](featI);
|
|
const pointField& points = eMesh.points();
|
|
const edgeList& edges = eMesh.edges();
|
|
|
|
// Calculate bb of all points
|
|
treeBoundBox bb(points);
|
|
|
|
// Random number generator. Bit dodgy since not exactly random ;-)
|
|
Random rndGen(65431);
|
|
|
|
// Slightly extended bb. Slightly off-centred just so on symmetric
|
|
// geometry there are less face/edge aligned items.
|
|
bb = bb.extend(rndGen, 1e-4);
|
|
bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
|
|
bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
|
|
|
|
trees.set
|
|
(
|
|
featI,
|
|
new indexedOctree<treeDataEdge>
|
|
(
|
|
treeDataEdge
|
|
(
|
|
false, // do not cache bb
|
|
edges,
|
|
points,
|
|
eMesh.regionEdges()
|
|
),
|
|
bb, // overall search domain
|
|
8, // maxLevel
|
|
10, // leafsize
|
|
3.0 // duplicity
|
|
)
|
|
);
|
|
}
|
|
}
|
|
return regionEdgeTreesPtr_();
|
|
}
|
|
|
|
|
|
// Find maximum level of a shell.
|
|
void Foam::refinementFeatures::findHigherLevel
|
|
(
|
|
const pointField& pt,
|
|
const label featI,
|
|
labelList& maxLevel
|
|
) const
|
|
{
|
|
const labelList& levels = levels_[featI];
|
|
|
|
const scalarField& distances = distances_[featI];
|
|
|
|
// Collect all those points that have a current maxLevel less than
|
|
// (any of) the shell. Also collect the furthest distance allowable
|
|
// to any shell with a higher level.
|
|
|
|
pointField candidates(pt.size());
|
|
labelList candidateMap(pt.size());
|
|
scalarField candidateDistSqr(pt.size());
|
|
label candidateI = 0;
|
|
|
|
forAll(maxLevel, pointI)
|
|
{
|
|
forAllReverse(levels, levelI)
|
|
{
|
|
if (levels[levelI] > maxLevel[pointI])
|
|
{
|
|
candidates[candidateI] = pt[pointI];
|
|
candidateMap[candidateI] = pointI;
|
|
candidateDistSqr[candidateI] = sqr(distances[levelI]);
|
|
candidateI++;
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
candidates.setSize(candidateI);
|
|
candidateMap.setSize(candidateI);
|
|
candidateDistSqr.setSize(candidateI);
|
|
|
|
// Do the expensive nearest test only for the candidate points.
|
|
const indexedOctree<treeDataEdge>& tree = edgeTrees_[featI];
|
|
|
|
List<pointIndexHit> nearInfo(candidates.size());
|
|
forAll(candidates, candidateI)
|
|
{
|
|
nearInfo[candidateI] = tree.findNearest
|
|
(
|
|
candidates[candidateI],
|
|
candidateDistSqr[candidateI]
|
|
);
|
|
}
|
|
|
|
// Update maxLevel
|
|
forAll(nearInfo, candidateI)
|
|
{
|
|
if (nearInfo[candidateI].hit())
|
|
{
|
|
// Check which level it actually is in.
|
|
label minDistI = findLower
|
|
(
|
|
distances,
|
|
mag(nearInfo[candidateI].hitPoint()-candidates[candidateI])
|
|
);
|
|
|
|
label pointI = candidateMap[candidateI];
|
|
|
|
// pt is inbetween shell[minDistI] and shell[minDistI+1]
|
|
maxLevel[pointI] = levels[minDistI+1];
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
|
|
|
|
Foam::refinementFeatures::refinementFeatures
|
|
(
|
|
const objectRegistry& io,
|
|
const PtrList<dictionary>& featDicts
|
|
)
|
|
:
|
|
PtrList<extendedFeatureEdgeMesh>(featDicts.size()),
|
|
distances_(featDicts.size()),
|
|
levels_(featDicts.size()),
|
|
edgeTrees_(featDicts.size()),
|
|
pointTrees_(featDicts.size())
|
|
{
|
|
// Read features
|
|
read(io, featDicts);
|
|
|
|
// Search engines
|
|
forAll(*this, i)
|
|
{
|
|
buildTrees(i);
|
|
}
|
|
}
|
|
|
|
|
|
//Foam::refinementFeatures::refinementFeatures
|
|
//(
|
|
// const objectRegistry& io,
|
|
// const PtrList<dictionary>& featDicts,
|
|
// const scalar minCos
|
|
//)
|
|
//:
|
|
// PtrList<extendedFeatureEdgeMesh>(featDicts.size()),
|
|
// distances_(featDicts.size()),
|
|
// levels_(featDicts.size()),
|
|
// edgeTrees_(featDicts.size()),
|
|
// pointTrees_(featDicts.size())
|
|
//{
|
|
// // Read features
|
|
// read(io, featDicts);
|
|
//
|
|
// // Search engines
|
|
// forAll(*this, i)
|
|
// {
|
|
// const edgeMesh& eMesh = operator[](i);
|
|
// const pointField& points = eMesh.points();
|
|
// const edgeList& edges = eMesh.edges();
|
|
// const labelListList& pointEdges = eMesh.pointEdges();
|
|
//
|
|
// DynamicList<label> featurePoints;
|
|
// forAll(pointEdges, pointI)
|
|
// {
|
|
// const labelList& pEdges = pointEdges[pointI];
|
|
// if (pEdges.size() > 2)
|
|
// {
|
|
// featurePoints.append(pointI);
|
|
// }
|
|
// else if (pEdges.size() == 2)
|
|
// {
|
|
// // Check the angle
|
|
// const edge& e0 = edges[pEdges[0]];
|
|
// const edge& e1 = edges[pEdges[1]];
|
|
//
|
|
// const point& p = points[pointI];
|
|
// const point& p0 = points[e0.otherVertex(pointI)];
|
|
// const point& p1 = points[e1.otherVertex(pointI)];
|
|
//
|
|
// vector v0 = p-p0;
|
|
// scalar v0Mag = mag(v0);
|
|
//
|
|
// vector v1 = p1-p;
|
|
// scalar v1Mag = mag(v1);
|
|
//
|
|
// if
|
|
// (
|
|
// v0Mag > SMALL
|
|
// && v1Mag > SMALL
|
|
// && ((v0/v0Mag & v1/v1Mag) < minCos)
|
|
// )
|
|
// {
|
|
// featurePoints.append(pointI);
|
|
// }
|
|
// }
|
|
// }
|
|
//
|
|
// Info<< "Detected " << featurePoints.size()
|
|
// << " featurePoints out of " << points.size()
|
|
// << " points on feature " << i //eMesh.name()
|
|
// << " when using feature cos " << minCos << endl;
|
|
//
|
|
// buildTrees(i, featurePoints);
|
|
// }
|
|
//}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
|
|
|
void Foam::refinementFeatures::findNearestEdge
|
|
(
|
|
const pointField& samples,
|
|
const scalarField& nearestDistSqr,
|
|
labelList& nearFeature,
|
|
List<pointIndexHit>& nearInfo,
|
|
vectorField& nearNormal
|
|
) const
|
|
{
|
|
nearFeature.setSize(samples.size());
|
|
nearFeature = -1;
|
|
nearInfo.setSize(samples.size());
|
|
nearInfo = pointIndexHit();
|
|
nearNormal.setSize(samples.size());
|
|
nearNormal = vector::zero;
|
|
|
|
forAll(edgeTrees_, featI)
|
|
{
|
|
const indexedOctree<treeDataEdge>& tree = edgeTrees_[featI];
|
|
|
|
if (tree.shapes().size() > 0)
|
|
{
|
|
forAll(samples, sampleI)
|
|
{
|
|
const point& sample = samples[sampleI];
|
|
|
|
scalar distSqr;
|
|
if (nearInfo[sampleI].hit())
|
|
{
|
|
distSqr = magSqr(nearInfo[sampleI].hitPoint()-sample);
|
|
}
|
|
else
|
|
{
|
|
distSqr = nearestDistSqr[sampleI];
|
|
}
|
|
|
|
pointIndexHit info = tree.findNearest(sample, distSqr);
|
|
|
|
if (info.hit())
|
|
{
|
|
nearFeature[sampleI] = featI;
|
|
nearInfo[sampleI] = pointIndexHit
|
|
(
|
|
info.hit(),
|
|
info.hitPoint(),
|
|
tree.shapes().edgeLabels()[info.index()]
|
|
);
|
|
|
|
const treeDataEdge& td = tree.shapes();
|
|
const edge& e = td.edges()[nearInfo[sampleI].index()];
|
|
nearNormal[sampleI] = e.vec(td.points());
|
|
nearNormal[sampleI] /= mag(nearNormal[sampleI])+VSMALL;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
void Foam::refinementFeatures::findNearestRegionEdge
|
|
(
|
|
const pointField& samples,
|
|
const scalarField& nearestDistSqr,
|
|
labelList& nearFeature,
|
|
List<pointIndexHit>& nearInfo,
|
|
vectorField& nearNormal
|
|
) const
|
|
{
|
|
nearFeature.setSize(samples.size());
|
|
nearFeature = -1;
|
|
nearInfo.setSize(samples.size());
|
|
nearInfo = pointIndexHit();
|
|
nearNormal.setSize(samples.size());
|
|
nearNormal = vector::zero;
|
|
|
|
|
|
const PtrList<indexedOctree<treeDataEdge> >& regionTrees =
|
|
regionEdgeTrees();
|
|
|
|
forAll(regionTrees, featI)
|
|
{
|
|
const indexedOctree<treeDataEdge>& regionTree = regionTrees[featI];
|
|
|
|
forAll(samples, sampleI)
|
|
{
|
|
const point& sample = samples[sampleI];
|
|
|
|
scalar distSqr;
|
|
if (nearInfo[sampleI].hit())
|
|
{
|
|
distSqr = magSqr(nearInfo[sampleI].hitPoint()-sample);
|
|
}
|
|
else
|
|
{
|
|
distSqr = nearestDistSqr[sampleI];
|
|
}
|
|
|
|
// Find anything closer than current best
|
|
pointIndexHit info = regionTree.findNearest(sample, distSqr);
|
|
|
|
if (info.hit())
|
|
{
|
|
const treeDataEdge& td = regionTree.shapes();
|
|
|
|
nearFeature[sampleI] = featI;
|
|
nearInfo[sampleI] = pointIndexHit
|
|
(
|
|
info.hit(),
|
|
info.hitPoint(),
|
|
regionTree.shapes().edgeLabels()[info.index()]
|
|
);
|
|
|
|
const edge& e = td.edges()[nearInfo[sampleI].index()];
|
|
nearNormal[sampleI] = e.vec(td.points());
|
|
nearNormal[sampleI] /= mag(nearNormal[sampleI])+VSMALL;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
//void Foam::refinementFeatures::findNearestPoint
|
|
//(
|
|
// const pointField& samples,
|
|
// const scalarField& nearestDistSqr,
|
|
// labelList& nearFeature,
|
|
// labelList& nearIndex
|
|
//) const
|
|
//{
|
|
// nearFeature.setSize(samples.size());
|
|
// nearFeature = -1;
|
|
// nearIndex.setSize(samples.size());
|
|
// nearIndex = -1;
|
|
//
|
|
// forAll(pointTrees_, featI)
|
|
// {
|
|
// const indexedOctree<treeDataPoint>& tree = pointTrees_[featI];
|
|
//
|
|
// if (tree.shapes().pointLabels().size() > 0)
|
|
// {
|
|
// forAll(samples, sampleI)
|
|
// {
|
|
// const point& sample = samples[sampleI];
|
|
//
|
|
// scalar distSqr;
|
|
// if (nearFeature[sampleI] != -1)
|
|
// {
|
|
// label nearFeatI = nearFeature[sampleI];
|
|
// const indexedOctree<treeDataPoint>& nearTree =
|
|
// pointTrees_[nearFeatI];
|
|
// label featPointI =
|
|
// nearTree.shapes().pointLabels()[nearIndex[sampleI]];
|
|
// const point& featPt =
|
|
// operator[](nearFeatI).points()[featPointI];
|
|
// distSqr = magSqr(featPt-sample);
|
|
// }
|
|
// else
|
|
// {
|
|
// distSqr = nearestDistSqr[sampleI];
|
|
// }
|
|
//
|
|
// pointIndexHit info = tree.findNearest(sample, distSqr);
|
|
//
|
|
// if (info.hit())
|
|
// {
|
|
// nearFeature[sampleI] = featI;
|
|
// nearIndex[sampleI] = info.index();
|
|
// }
|
|
// }
|
|
// }
|
|
// }
|
|
//}
|
|
|
|
|
|
void Foam::refinementFeatures::findNearestPoint
|
|
(
|
|
const pointField& samples,
|
|
const scalarField& nearestDistSqr,
|
|
labelList& nearFeature,
|
|
List<pointIndexHit>& nearInfo
|
|
) const
|
|
{
|
|
nearFeature.setSize(samples.size());
|
|
nearFeature = -1;
|
|
nearInfo.setSize(samples.size());
|
|
nearInfo = pointIndexHit();
|
|
|
|
forAll(pointTrees_, featI)
|
|
{
|
|
const indexedOctree<treeDataPoint>& tree = pointTrees_[featI];
|
|
|
|
if (tree.shapes().pointLabels().size() > 0)
|
|
{
|
|
forAll(samples, sampleI)
|
|
{
|
|
const point& sample = samples[sampleI];
|
|
|
|
scalar distSqr;
|
|
if (nearFeature[sampleI] != -1)
|
|
{
|
|
distSqr = magSqr(nearInfo[sampleI].hitPoint()-sample);
|
|
}
|
|
else
|
|
{
|
|
distSqr = nearestDistSqr[sampleI];
|
|
}
|
|
|
|
pointIndexHit info = tree.findNearest(sample, distSqr);
|
|
|
|
if (info.hit())
|
|
{
|
|
nearFeature[sampleI] = featI;
|
|
nearInfo[sampleI] = pointIndexHit
|
|
(
|
|
info.hit(),
|
|
info.hitPoint(),
|
|
tree.shapes().pointLabels()[info.index()]
|
|
);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
void Foam::refinementFeatures::findHigherLevel
|
|
(
|
|
const pointField& pt,
|
|
const labelList& ptLevel,
|
|
labelList& maxLevel
|
|
) const
|
|
{
|
|
// Maximum level of any shell. Start off with level of point.
|
|
maxLevel = ptLevel;
|
|
|
|
forAll(*this, featI)
|
|
{
|
|
findHigherLevel(pt, featI, maxLevel);
|
|
}
|
|
}
|
|
|
|
|
|
Foam::scalar Foam::refinementFeatures::maxDistance() const
|
|
{
|
|
scalar overallMax = -GREAT;
|
|
forAll(distances_, featI)
|
|
{
|
|
overallMax = max(overallMax, max(distances_[featI]));
|
|
}
|
|
return overallMax;
|
|
}
|
|
|
|
|
|
// ************************************************************************* //
|