Newer
Older
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\/ 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 "sampledSet.H"
#include "polyMesh.H"
#include "primitiveMesh.H"
#include "meshSearch.H"
#include "writer.H"
#include "particle.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(sampledSet, 0);
defineRunTimeSelectionTable(sampledSet, word);
}
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
void Foam::sampledSet::checkDimensions() const
{
if
(
(cells_.size() != size())
|| (faces_.size() != size())
|| (segments_.size() != size())
|| (curveDist_.size() != size())
)
{
FatalErrorInFunction
<< "sizes not equal : "
<< " points:" << size()
<< " cells:" << cells_.size()
<< " faces:" << faces_.size()
<< " segments:" << segments_.size()
<< " curveDist:" << curveDist_.size()
<< abort(FatalError);
}
}
Foam::label Foam::sampledSet::getBoundaryCell(const label facei) const
return mesh().faceOwner()[facei];
Foam::label Foam::sampledSet::getNeighbourCell(const label facei) const
Henry Weller
committed
{
if (facei >= mesh().nInternalFaces())
Henry Weller
committed
{
return mesh().faceOwner()[facei];
Henry Weller
committed
}
else
{
return mesh().faceNeighbour()[facei];
Henry Weller
committed
}
}
Henry Weller
committed
Foam::label Foam::sampledSet::pointInCell
Henry Weller
committed
const point& p,
const label samplei
) const
{
Henry Weller
committed
// Collect the face owner and neighbour cells of the sample into an array
// for convenience
const label cells[4] =
Henry Weller
committed
{
mesh().faceOwner()[faces_[samplei]],
getNeighbourCell(faces_[samplei]),
mesh().faceOwner()[faces_[samplei+1]],
getNeighbourCell(faces_[samplei+1])
};
// Find the sampled cell by checking the owners and neighbours of the
// sampled faces
Henry Weller
committed
label cellm =
Henry Weller
committed
(cells[0] == cells[2] || cells[0] == cells[3]) ? cells[0]
: (cells[1] == cells[2] || cells[1] == cells[3]) ? cells[1]
Henry Weller
committed
: -1;
Henry Weller
committed
if (cellm != -1)
{
Henry Weller
committed
// If found the sampled cell check the point is in the cell
// otherwise ignore
Henry Weller
committed
if (!mesh().pointInCell(p, cellm, searchEngine_.decompMode()))
cellm = -1;
Henry Weller
committed
if (debug)
Henry Weller
committed
WarningInFunction
<< "Could not find mid-point " << p
<< " cell " << cellm << endl;
}
}
}
Henry Weller
committed
else
Henry Weller
committed
{
Henry Weller
committed
// If the sample does not pass through a single cell check if the point
// is in any of the owners or neighbours otherwise ignore
for (label i=0; i<4; ++i)
Henry Weller
committed
{
if (mesh().pointInCell(p, cells[i], searchEngine_.decompMode()))
{
return cells[i];
}
}
if (debug)
{
WarningInFunction
<< "Could not find cell for mid-point" << nl
<< " samplei: " << samplei
<< " pts[samplei]: " << operator[](samplei)
<< " face[samplei]: " << faces_[samplei]
<< " pts[samplei+1]: " << operator[](samplei+1)
<< " face[samplei+1]: " << faces_[samplei+1]
<< " cellio: " << cells[0]
<< " cellin: " << cells[1]
<< " celljo: " << cells[2]
<< " celljn: " << cells[3]
<< endl;
}
Henry Weller
committed
}
return cellm;
}
Foam::scalar Foam::sampledSet::calcSign
(
const point& sample
) const
{
vector vec = sample - mesh().faceCentres()[facei];
scalar magVec = mag(vec);
if (magVec < VSMALL)
{
// Sample on face centre. Regard as inside
return -1;
}
vec /= magVec;
const vector n = normalised(mesh().faceAreas()[facei]);
return n & vec;
}
Foam::label Foam::sampledSet::findNearFace
(
const point& sample,
const scalar smallDist
) const
{
const cell& myFaces = mesh().cells()[celli];
forAll(myFaces, myFacei)
const face& f = mesh().faces()[myFaces[myFacei]];
pointHit inter = f.nearestPoint(sample, mesh().points());
scalar dist;
if (inter.hit())
{
dist = mag(inter.hitPoint() - sample);
}
else
{
dist = mag(inter.missPoint() - sample);
}
if (dist < smallDist)
{
return myFaces[myFacei];
}
}
return -1;
}
Foam::point Foam::sampledSet::pushIn
(
const point& facePt,
) const
{
label celli = mesh().faceOwner()[facei];
const point& cC = mesh().cellCentres()[celli];
point newPosition = facePt;
// Taken from particle::initCellFacePt()
label tetFacei;
label tetPtI;
mesh().findTetFacePt(celli, facePt, tetFacei, tetPtI);
// This is the tolerance that was defined as a static constant of the
// particle class. It is no longer used by particle, following the switch to
// barycentric tracking. The only place in which the tolerance is now used
// is here. I'm not sure what the purpose of this code is, but it probably
// wants removing. It is doing tet-searches for a particle position. This
// should almost certainly be left to the particle class.
const scalar trackingCorrectionTol = 1e-5;
if (tetFacei == -1 || tetPtI == -1)
{
newPosition = facePt;
label trap(1.0/trackingCorrectionTol + 1);
label iterNo = 0;
do
{
newPosition += trackingCorrectionTol*(cC - facePt);
mesh().findTetFacePt
(
newPosition,
tetPtI
);
++iterNo;
} while (tetFacei < 0 && iterNo <= trap);
}
if (tetFacei == -1)
FatalErrorInFunction
<< "After pushing " << facePt << " to " << newPosition
<< " it is still outside face " << facei
<< " at " << mesh().faceCentres()[facei]
<< " of cell " << celli
<< " at " << cC << endl
<< "Please change your starting point"
<< abort(FatalError);
}
return newPosition;
}
bool Foam::sampledSet::getTrackingPoint
(
const point& samplePt,
const point& bPoint,
const label bFacei,
Henry Weller
committed
const scalar smallDist,
point& trackPt,
label& trackCelli,
label& trackFacei
) const
{
bool isGoodSample = false;
if (bFacei == -1)
{
// No boundary intersection. Try and find cell samplePt is in
trackCelli = mesh().findCell(samplePt, searchEngine_.decompMode());
if
(
(trackCelli == -1)
|| !mesh().pointInCell
(
samplePt,
trackCelli,
searchEngine_.decompMode()
)
)
{
// Line samplePt - end_ does not intersect domain at all.
// (or is along edge)
trackCelli = -1;
trackFacei = -1;
isGoodSample = false;
}
else
{
// Start is inside. Use it as tracking point
trackPt = samplePt;
trackFacei = -1;
isGoodSample = true;
}
}
else if (mag(samplePt - bPoint) < smallDist)
{
// samplePt close to bPoint. Snap to it
trackPt = pushIn(bPoint, bFacei);
trackFacei = bFacei;
trackCelli = getBoundaryCell(trackFacei);
isGoodSample = true;
}
else
{
scalar sign = calcSign(bFacei, samplePt);
if (sign < 0)
{
// samplePt inside or marginally outside.
trackPt = samplePt;
trackFacei = -1;
trackCelli = mesh().findCell(trackPt, searchEngine_.decompMode());
isGoodSample = true;
}
else
{
// samplePt outside. use bPoint
trackPt = pushIn(bPoint, bFacei);
trackFacei = bFacei;
trackCelli = getBoundaryCell(trackFacei);
isGoodSample = false;
}
}
DebugInFunction
<< " samplePt:" << samplePt
<< " bPoint:" << bPoint
<< " bFacei:" << bFacei << nl
<< " Calculated first tracking point :"
<< " trackPt:" << trackPt
<< " trackCelli:" << trackCelli
<< " trackFacei:" << trackFacei
<< " isGoodSample:" << isGoodSample
<< endl;
return isGoodSample;
}
void Foam::sampledSet::setSamples
(
const List<point>& samplingPts,
const labelList& samplingCells,
const labelList& samplingFaces,
const labelList& samplingSegments,
const scalarList& samplingCurveDist
)
{
setPoints(samplingPts);
curveDist_ = samplingCurveDist;
segments_ = samplingSegments;
cells_ = samplingCells;
faces_ = samplingFaces;
checkDimensions();
}
void Foam::sampledSet::setSamples
(
List<point>&& samplingPts,
labelList&& samplingCells,
labelList&& samplingFaces,
labelList&& samplingSegments,
scalarList&& samplingCurveDist
)
{
setPoints(std::move(samplingPts));
curveDist_ = std::move(samplingCurveDist);
segments_ = std::move(samplingSegments);
cells_ = std::move(samplingCells);
faces_ = std::move(samplingFaces);
checkDimensions();
Foam::autoPtr<Foam::coordSet> Foam::sampledSet::gather
(
labelList& indexSet,
labelList& allSegments
) const
{
// Combine sampleSet from processors. Sort by curveDist. Return
// ordering in indexSet.
// Note: only master results are valid
List<point> allPts;
globalIndex::gatherOp(*this, allPts);
globalIndex::gatherOp(segments(), allSegments);
scalarList allCurveDist;
globalIndex::gatherOp(curveDist(), allCurveDist);
if (Pstream::master() && allCurveDist.empty())
{
WarningInFunction
<< "Sample set " << name()
<< " has zero points." << endl;
}
// Sort curveDist and use to fill masterSamplePts
Foam::sortedOrder(allCurveDist, indexSet); // uses stable sort
scalarList sortedDist(allCurveDist, indexSet); // with indices for mapping
allSegments = UIndirectList<label>(allSegments, indexSet)();
return autoPtr<coordSet>::New
name(),
axis(),
List<point>(UIndirectList<point>(allPts, indexSet)),
sortedDist
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::sampledSet::sampledSet
(
const word& name,
const polyMesh& mesh,
const meshSearch& searchEngine,
const coordSet::coordFormat axisType
coordSet(name, axisType),
mesh_(mesh),
searchEngine_(searchEngine),
segments_(),
cells_(),
faces_()
{}
Foam::sampledSet::sampledSet
(
const word& name,
const polyMesh& mesh,
const meshSearch& searchEngine,
const word& axis
coordSet(name, axis),
mesh_(mesh),
searchEngine_(searchEngine),
segments_(),
cells_(),
faces_()
Foam::sampledSet::sampledSet
(
const word& name,
const polyMesh& mesh,
const meshSearch& searchEngine,
const dictionary& dict
)
:
coordSet(name, dict.get<word>("axis")),
mesh_(mesh),
searchEngine_(searchEngine),
segments_(),
cells_(),
faces_()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::autoPtr<Foam::sampledSet> Foam::sampledSet::New
(
const word& name,
const polyMesh& mesh,
const meshSearch& searchEngine,
const dictionary& dict
)
{
const word sampleType(dict.get<word>("type"));
auto cstrIter = wordConstructorTablePtr_->cfind(sampleType);
FatalIOErrorInLookup
"sample",
sampleType,
*wordConstructorTablePtr_
) << exit(FatalIOError);
}
return autoPtr<sampledSet>
(
cstrIter()
(
name,
mesh,
searchEngine,
dict.optionalSubDict(sampleType + "Coeffs")
)
);
}
Foam::Ostream& Foam::sampledSet::write(Ostream& os) const
{
coordSet::write(os);
os << nl << "\t(celli)\t(facei)" << nl;
forAll(*this, samplei)
os << '\t' << cells_[samplei]
<< '\t' << faces_[samplei]
<< nl;