BUG: autoSnapDriver: inplace re-ordering!

This commit is contained in:
mattijs
2013-12-19 11:24:02 +00:00
parent 64b6f921e8
commit 60e3467165
4 changed files with 354 additions and 267 deletions

View File

@ -310,14 +310,13 @@ class autoSnapDriver
//- Find point on nearest feature edge (within searchDist).
// Return point and feature
// and store feature-edge to mesh-point and vice versa
pointIndexHit findNearFeatureEdge
Tuple2<label, pointIndexHit> findNearFeatureEdge
(
const indirectPrimitivePatch& pp,
const scalarField& snapDist,
const label pointI,
const point& estimatedPt,
label& featI,
List<List<DynamicList<point> > >&,
List<List<DynamicList<pointConstraint> > >&,
vectorField&,
@ -330,7 +329,7 @@ class autoSnapDriver
// If another mesh point already referring to this feature
// point and further away, reset that one to a near feature
// edge (using findNearFeatureEdge above)
labelPair findNearFeaturePoint
Tuple2<label, pointIndexHit> findNearFeaturePoint
(
const indirectPrimitivePatch& pp,
const scalarField& snapDist,

View File

@ -170,7 +170,7 @@ void Foam::autoSnapDriver::smoothAndConstrain
if (isMasterEdge[meshEdges[edgeI]])
{
label nbrPointI = edges[pEdges[i]].otherVertex(pointI);
label nbrPointI = edges[edgeI].otherVertex(pointI);
if (constraints[nbrPointI].first() >= nConstraints)
{
dispSum[pointI] += disp[nbrPointI];
@ -213,80 +213,6 @@ void Foam::autoSnapDriver::smoothAndConstrain
}
}
}
//XXXXXX
//TODO: make proper parallel so coupled edges don't have double influence
//void Foam::autoSnapDriver::smoothAndConstrain2
//(
// const bool applyConstraints,
// const indirectPrimitivePatch& pp,
// const List<pointConstraint>& constraints,
// vectorField& disp
//) const
//{
// const fvMesh& mesh = meshRefiner_.mesh();
//
// for (label avgIter = 0; avgIter < 20; avgIter++)
// {
// vectorField dispSum(pp.nPoints(), vector::zero);
// labelList dispCount(pp.nPoints(), 0);
//
// const labelListList& pointEdges = pp.pointEdges();
// const edgeList& edges = pp.edges();
//
// forAll(pointEdges, pointI)
// {
// const labelList& pEdges = pointEdges[pointI];
//
// forAll(pEdges, i)
// {
// label nbrPointI = edges[pEdges[i]].otherVertex(pointI);
// dispSum[pointI] += disp[nbrPointI];
// dispCount[pointI]++;
// }
// }
//
// syncTools::syncPointList
// (
// mesh,
// pp.meshPoints(),
// dispSum,
// plusEqOp<point>(),
// vector::zero,
// mapDistribute::transform()
// );
// syncTools::syncPointList
// (
// mesh,
// pp.meshPoints(),
// dispCount,
// plusEqOp<label>(),
// 0,
// mapDistribute::transform()
// );
//
// // Constraints
// forAll(constraints, pointI)
// {
// if (dispCount[pointI] > 0)// && constraints[pointI].first() <= 1)
// {
// // Mix my displacement with neighbours' displacement
// disp[pointI] =
// 0.5
// *(disp[pointI] + dispSum[pointI]/dispCount[pointI]);
//
// if (applyConstraints)
// {
// disp[pointI] = transform
// (
// constraints[pointI].constraintTransformation(),
// disp[pointI]
// );
// }
// }
// }
// }
//}
//XXXXXX
void Foam::autoSnapDriver::calcNearestFace
@ -716,7 +642,7 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
pNormals = List<point>(pNormals, visitOrder);
pDisp = List<point>(pDisp, visitOrder);
pFc = List<point>(pFc, visitOrder);
pFid = UIndirectList<label>(pFid, visitOrder);
pFid = UIndirectList<label>(pFid, visitOrder)();
}
}
@ -1412,14 +1338,14 @@ Foam::labelPair Foam::autoSnapDriver::findDiagonalAttraction
}
Foam::pointIndexHit Foam::autoSnapDriver::findNearFeatureEdge
Foam::Tuple2<Foam::label, Foam::pointIndexHit>
Foam::autoSnapDriver::findNearFeatureEdge
(
const indirectPrimitivePatch& pp,
const scalarField& snapDist,
const label pointI,
const point& estimatedPt,
label& featI,
List<List<DynamicList<point> > >& edgeAttractors,
List<List<DynamicList<pointConstraint> > >& edgeConstraints,
vectorField& patchAttraction,
@ -1430,16 +1356,18 @@ Foam::pointIndexHit Foam::autoSnapDriver::findNearFeatureEdge
labelList nearEdgeFeat;
List<pointIndexHit> nearEdgeInfo;
vectorField nearNormal;
features.findNearestEdge
(
pointField(1, estimatedPt),
scalarField(1, sqr(snapDist[pointI])),
nearEdgeFeat,
nearEdgeInfo
nearEdgeInfo,
nearNormal
);
const pointIndexHit& nearInfo = nearEdgeInfo[0];
featI = nearEdgeFeat[0];
label featI = nearEdgeFeat[0];
if (nearInfo.hit())
{
@ -1449,12 +1377,7 @@ Foam::pointIndexHit Foam::autoSnapDriver::findNearFeatureEdge
(
nearInfo.hitPoint()
);
pointConstraint c;
const edge e = features[featI].edges()[nearInfo.index()];
vector eVec = e.vec(features[featI].points());
eVec /= mag(eVec)+VSMALL;
c.first() = 2;
c.second() = eVec;
pointConstraint c(Tuple2<label, vector>(2, nearNormal[0]));
edgeConstraints[featI][nearInfo.index()].append(c);
// Store for later use
@ -1462,11 +1385,12 @@ Foam::pointIndexHit Foam::autoSnapDriver::findNearFeatureEdge
nearInfo.hitPoint()-pp.localPoints()[pointI];
patchConstraints[pointI] = c;
}
return nearInfo;
return Tuple2<label, pointIndexHit>(featI, nearInfo);
}
Foam::labelPair Foam::autoSnapDriver::findNearFeaturePoint
Foam::Tuple2<Foam::label, Foam::pointIndexHit>
Foam::autoSnapDriver::findNearFeaturePoint
(
const indirectPrimitivePatch& pp,
const scalarField& snapDist,
@ -1487,26 +1411,23 @@ Foam::labelPair Foam::autoSnapDriver::findNearFeaturePoint
const refinementFeatures& features = meshRefiner_.features();
labelList nearFeat;
labelList nearIndex;
List<pointIndexHit> nearInfo;
features.findNearestPoint
(
pointField(1, estimatedPt),
scalarField(1, sqr(snapDist[pointI])),
nearFeat,
nearIndex
nearInfo
);
label featI = nearFeat[0];
label featPointI = -1;
if (featI != -1)
{
const point& pt = pp.localPoints()[pointI];
const treeDataPoint& shapes =
features.pointTrees()[featI].shapes();
featPointI = shapes.pointLabels()[nearIndex[0]];
const point& featPt = shapes.points()[featPointI];
label featPointI = nearInfo[0].index();
const point& featPt = nearInfo[0].hitPoint();
scalar distSqr = magSqr(featPt-pt);
// Check if already attracted
@ -1537,7 +1458,6 @@ Foam::labelPair Foam::autoSnapDriver::findNearFeaturePoint
patchAttraction[oldPointI] = vector::zero;
patchConstraints[oldPointI] = pointConstraint();
label edgeFeatI;
findNearFeatureEdge
(
pp,
@ -1545,7 +1465,6 @@ Foam::labelPair Foam::autoSnapDriver::findNearFeaturePoint
oldPointI,
pp.localPoints()[oldPointI],
edgeFeatI,
edgeAttractors,
edgeConstraints,
patchAttraction,
@ -1566,7 +1485,7 @@ Foam::labelPair Foam::autoSnapDriver::findNearFeaturePoint
}
}
return labelPair(featI, featPointI);
return Tuple2<label, pointIndexHit>(featI, nearInfo[0]);
}
@ -1637,7 +1556,6 @@ void Foam::autoSnapDriver::determineFeatures
<< featurePointStr().name() << endl;
}
const refinementFeatures& features = meshRefiner_.features();
forAll(pp.localPoints(), pointI)
{
@ -1664,6 +1582,7 @@ void Foam::autoSnapDriver::determineFeatures
constraint
);
if
(
(constraint.first() > patchConstraints[pointI].first())
@ -1696,15 +1615,14 @@ void Foam::autoSnapDriver::determineFeatures
// Behave like when having two surface normals so
// attract to nearest feature edge (with a guess for
// the multipatch point as starting point)
label featI = -1;
pointIndexHit nearInfo = findNearFeatureEdge
Tuple2<label, pointIndexHit> nearInfo =
findNearFeatureEdge
(
pp,
snapDist,
pointI,
multiPatchPt.hitPoint(), //estimatedPt
featI,
edgeAttractors,
edgeConstraints,
@ -1712,14 +1630,15 @@ void Foam::autoSnapDriver::determineFeatures
patchConstraints
);
if (nearInfo.hit())
const pointIndexHit& info = nearInfo.second();
if (info.hit())
{
// Dump
if (featureEdgeStr.valid())
{
featureEdgeStr().write
(
linePointRef(pt, nearInfo.hitPoint())
linePointRef(pt, info.hitPoint())
);
}
}
@ -1729,7 +1648,7 @@ void Foam::autoSnapDriver::determineFeatures
{
missedEdgeStr().write
(
linePointRef(pt, nearInfo.missPoint())
linePointRef(pt, info.missPoint())
);
}
}
@ -1746,15 +1665,13 @@ void Foam::autoSnapDriver::determineFeatures
// Determine nearest point on feature edge. Store constraint
// (calculated from feature edge, alternative would be to
// use constraint calculated from both surfaceNormals)
label featI = -1;
pointIndexHit nearInfo = findNearFeatureEdge
Tuple2<label, pointIndexHit> nearInfo = findNearFeatureEdge
(
pp,
snapDist,
pointI,
estimatedPt,
featI,
edgeAttractors,
edgeConstraints,
@ -1762,14 +1679,15 @@ void Foam::autoSnapDriver::determineFeatures
patchConstraints
);
if (nearInfo.hit())
// Dump to obj
const pointIndexHit& info = nearInfo.second();
if (info.hit())
{
// Dump
if (featureEdgeStr.valid())
{
featureEdgeStr().write
(
linePointRef(pt, nearInfo.hitPoint())
linePointRef(pt, info.hitPoint())
);
}
}
@ -1779,7 +1697,7 @@ void Foam::autoSnapDriver::determineFeatures
{
missedEdgeStr().write
(
linePointRef(pt, nearInfo.missPoint())
linePointRef(pt, info.missPoint())
);
}
}
@ -1789,7 +1707,7 @@ void Foam::autoSnapDriver::determineFeatures
// Mark point on the nearest feature point.
const point estimatedPt(pt + patchAttraction[pointI]);
labelPair nearInfo = findNearFeaturePoint
Tuple2<label, pointIndexHit> nearInfo = findNearFeaturePoint
(
pp,
snapDist,
@ -1807,18 +1725,10 @@ void Foam::autoSnapDriver::determineFeatures
patchConstraints
);
if (nearInfo.first() != -1)
const pointIndexHit& info = nearInfo.second();
if (info.hit() && featurePointStr.valid())
{
// Dump
if (featurePointStr.valid())
{
const treeDataPoint& shapes =
features.pointTrees()[nearInfo.first()].shapes();
const point& featPt =
shapes.points()[nearInfo.second()];
featurePointStr().write(linePointRef(pt, featPt));
}
featurePointStr().write(linePointRef(pt, info.hitPoint()));
}
}
}
@ -2029,22 +1939,20 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
if (pointStatus[pointI] == 0) // baffle edge
{
label featI;
const pointIndexHit nearInfo = findNearFeatureEdge
Tuple2<label, pointIndexHit> nearInfo = findNearFeatureEdge
(
pp,
snapDist,
pointI,
pt,
featI,
edgeAttractors,
edgeConstraints,
rawPatchAttraction,
rawPatchConstraints
);
if (!nearInfo.hit())
if (!nearInfo.second().hit())
{
//Pout<< "*** Failed to find close edge to point " << pt
// << endl;
@ -2053,23 +1961,21 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
else if (pointStatus[pointI] == 1) // baffle point
{
labelList nearFeat;
labelList nearIndex;
List<pointIndexHit> nearInfo;
features.findNearestPoint
(
pointField(1, pt),
scalarField(1, sqr(snapDist[pointI])),
nearFeat,
nearIndex
nearInfo
);
label featI = nearFeat[0];
if (featI != -1)
{
const treeDataPoint& shapes =
features.pointTrees()[featI].shapes();
label featPointI = shapes.pointLabels()[nearIndex[0]];
const point& featPt = shapes.points()[featPointI];
label featPointI = nearInfo[0].index();
const point& featPt = nearInfo[0].hitPoint();
scalar distSqr = magSqr(featPt-pt);
// Check if already attracted
@ -2099,7 +2005,6 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
// The current point is closer so wins. Reset
// the old point to attract to nearest edge
// instead.
label edgeFeatI;
findNearFeatureEdge
(
pp,
@ -2107,7 +2012,6 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
oldPointI,
pp.localPoints()[oldPointI],
edgeFeatI,
edgeAttractors,
edgeConstraints,
rawPatchAttraction,
@ -2130,7 +2034,6 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
// << " for baffle-feature-point " << pt
// << endl;
label featI;
findNearFeatureEdge
(
pp,
@ -2138,7 +2041,6 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
pointI,
pt, // starting point
featI,
edgeAttractors,
edgeConstraints,
rawPatchAttraction,

View File

@ -26,6 +26,7 @@ License
#include "refinementFeatures.H"
#include "Time.H"
#include "Tuple2.H"
#include "DynamicField.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@ -41,7 +42,40 @@ void Foam::refinementFeatures::read
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
@ -53,21 +87,106 @@ void Foam::refinementFeatures::read
false
);
autoPtr<edgeMesh> eMeshPtr = edgeMesh::New(featObj.filePath());
const fileName fName(featObj.filePath());
set
(
featI,
new featureEdgeMesh
if (fName.empty())
{
FatalIOErrorIn
(
featObj,
eMeshPtr->points(),
eMeshPtr->edges()
)
"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
labelList(0) // regionEdges
);
Info<< "Constructed extendedFeatureEdgeMesh " << featObj.name()
<< nl << incrIndent;
eeMesh.writeStats(Info);
Info<< decrIndent << endl;
set(featI, new extendedFeatureEdgeMesh(featObj, eeMesh));
}
const featureEdgeMesh& eMesh = operator[](featI);
const edgeMesh& eMesh = operator[](featI);
//eMesh.mergePoints(meshRefiner_.mergeDistance());
@ -145,7 +264,7 @@ void Foam::refinementFeatures::buildTrees
const labelList& featurePoints
)
{
const featureEdgeMesh& eMesh = operator[](featI);
const edgeMesh& eMesh = operator[](featI);
const pointField& points = eMesh.points();
const edgeList& edges = eMesh.edges();
@ -277,7 +396,7 @@ Foam::refinementFeatures::refinementFeatures
const PtrList<dictionary>& featDicts
)
:
PtrList<featureEdgeMesh>(featDicts.size()),
PtrList<extendedFeatureEdgeMesh>(featDicts.size()),
distances_(featDicts.size()),
levels_(featDicts.size()),
edgeTrees_(featDicts.size()),
@ -289,7 +408,7 @@ Foam::refinementFeatures::refinementFeatures
// Search engines
forAll(*this, i)
{
const featureEdgeMesh& eMesh = operator[](i);
const extendedEdgeMesh& eMesh = operator[](i);
const labelListList& pointEdges = eMesh.pointEdges();
DynamicList<label> featurePoints;
@ -303,81 +422,82 @@ Foam::refinementFeatures::refinementFeatures
Info<< "Detected " << featurePoints.size()
<< " featurePoints out of " << pointEdges.size()
<< " points on feature " << eMesh.name() << endl;
<< " points on feature " << operator[](i).name()
<< endl;
buildTrees(i, featurePoints);
}
}
Foam::refinementFeatures::refinementFeatures
(
const objectRegistry& io,
const PtrList<dictionary>& featDicts,
const scalar minCos
)
:
PtrList<featureEdgeMesh>(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 featureEdgeMesh& 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 " << eMesh.name()
<< " when using feature cos " << minCos << endl;
buildTrees(i, featurePoints);
}
}
//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 * * * * * * * * * * * * * //
@ -387,12 +507,16 @@ void Foam::refinementFeatures::findNearestEdge
const pointField& samples,
const scalarField& nearestDistSqr,
labelList& nearFeature,
List<pointIndexHit>& nearInfo
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)
{
@ -418,8 +542,18 @@ void Foam::refinementFeatures::findNearestEdge
if (info.hit())
{
nearInfo[sampleI] = info;
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;
}
}
}
@ -427,18 +561,71 @@ void Foam::refinementFeatures::findNearestEdge
}
//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,
labelList& nearIndex
List<pointIndexHit>& nearInfo
) const
{
nearFeature.setSize(samples.size());
nearFeature = -1;
nearIndex.setSize(samples.size());
nearIndex = -1;
nearInfo.setSize(samples.size());
nearInfo = pointIndexHit();
forAll(pointTrees_, featI)
{
@ -453,14 +640,7 @@ void Foam::refinementFeatures::findNearestPoint
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);
distSqr = magSqr(nearInfo[sampleI].hitPoint()-sample);
}
else
{
@ -472,7 +652,12 @@ void Foam::refinementFeatures::findNearestPoint
if (info.hit())
{
nearFeature[sampleI] = featI;
nearIndex[sampleI] = info.index();
nearInfo[sampleI] = pointIndexHit
(
info.hit(),
info.hitPoint(),
tree.shapes().pointLabels()[info.index()]
);
}
}
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -35,7 +35,7 @@ SourceFiles
#ifndef refinementFeatures_H
#define refinementFeatures_H
#include "featureEdgeMesh.H"
#include "extendedFeatureEdgeMesh.H"
#include "indexedOctree.H"
#include "treeDataEdge.H"
#include "treeDataPoint.H"
@ -51,7 +51,7 @@ namespace Foam
class refinementFeatures
:
public PtrList<featureEdgeMesh>
public PtrList<extendedFeatureEdgeMesh>
{
private:
@ -86,6 +86,20 @@ private:
labelList& maxLevel
) const;
protected:
const PtrList<indexedOctree<treeDataEdge> >& edgeTrees() const
{
return edgeTrees_;
}
const PtrList<indexedOctree<treeDataPoint> >& pointTrees() const
{
return pointTrees_;
}
public:
// Constructors
@ -97,15 +111,6 @@ public:
const PtrList<dictionary>& featDicts
);
//- Construct from description and do geometric analysis to determine
// feature points
refinementFeatures
(
const objectRegistry& io,
const PtrList<dictionary>& featDicts,
const scalar minCos
);
// Member Functions
@ -123,42 +128,38 @@ public:
return distances_;
}
const PtrList<indexedOctree<treeDataEdge> >& edgeTrees() const
{
return edgeTrees_;
}
const PtrList<indexedOctree<treeDataPoint> >& pointTrees() const
{
return pointTrees_;
}
// Query
//- Highest distance of all features
scalar maxDistance() const;
//- Find nearest point on nearest feature edge
//- Find nearest point on nearest feature edge. Sets
// - nearFeature: index of feature mesh
// - nearInfo : location on feature edge and edge index
// (note: not feature edge index but index into
// edges() directly)
// - nearNormal : local feature edge normal
void findNearestEdge
(
const pointField& samples,
const scalarField& nearestDistSqr,
labelList& nearFeature,
List<pointIndexHit>& nearInfo
List<pointIndexHit>& nearInfo,
vectorField& nearNormal
) const;
//- Find nearest feature point. Is an index into feature points
// which itself is an index into the edgeMesh points.
// So the point index is
// pointTrees()[nearFeature].shapes().pointLabels()[nearIndex]
// Wip.
//- Find nearest feature point. Sets
// - nearFeature: index of feature mesh
// - nearInfo : location on feature point and point index.
// (note: not index into shapes().pointLabels() but
// index into points() directly)
void findNearestPoint
(
const pointField& samples,
const scalarField& nearestDistSqr,
labelList& nearFeature,
labelList& nearIndex
List<pointIndexHit>& nearInfo
) const;
//- Find shell level higher than ptLevel