ENH: adjust displacementLayered framework

- add patch-point tracking to structured walk

- provision for cylindrical interpolation scheme

STYLE: more efficient use of bitSet
This commit is contained in:
Mark Olesen
2019-07-18 15:54:55 +02:00
committed by Andrew Heather
parent 18984f6c27
commit 55b1c57275
5 changed files with 224 additions and 191 deletions

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2015-2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2016 OpenFOAM Foundation
@ -59,6 +59,26 @@ namespace Foam
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
const Foam::coordSystem::cylindrical&
Foam::displacementLayeredMotionMotionSolver::getCylindrical
(
const label cellZoneI,
const dictionary& zoneDict
)
{
auto iter = cylSystems_.cfind(cellZoneI);
if (iter.found())
{
return *(iter.val());
}
cylSystems_.set(cellZoneI, new coordSystem::cylindrical(zoneDict));
return *cylSystems_[cellZoneI];
}
void Foam::displacementLayeredMotionMotionSolver::calcZoneMask
(
const label cellZoneI,
@ -66,68 +86,49 @@ void Foam::displacementLayeredMotionMotionSolver::calcZoneMask
bitSet& isZoneEdge
) const
{
isZonePoint.resize(mesh().nPoints());
isZoneEdge.resize(mesh().nEdges());
if (cellZoneI == -1)
{
isZonePoint.setSize(mesh().nPoints());
isZonePoint = true;
isZoneEdge.setSize(mesh().nEdges());
isZoneEdge = true;
return;
}
else
isZonePoint.reset();
isZoneEdge.reset();
const cellZone& cz = mesh().cellZones()[cellZoneI];
// Mark points, edges inside cellZone
for (const label celli : cz)
{
const cellZone& cz = mesh().cellZones()[cellZoneI];
label nPoints = 0;
forAll(cz, i)
{
const labelList& cPoints = mesh().cellPoints(cz[i]);
forAll(cPoints, cPointi)
{
if (isZonePoint.set(cPoints[cPointi]))
{
++nPoints;
}
}
}
syncTools::syncPointList
(
mesh(),
isZonePoint,
orEqOp<unsigned int>(),
0
);
// Mark edge inside cellZone
label nEdges = 0;
forAll(cz, i)
{
const labelList& cEdges = mesh().cellEdges(cz[i]);
forAll(cEdges, cEdgeI)
{
if (isZoneEdge.set(cEdges[cEdgeI]))
{
++nEdges;
}
}
}
syncTools::syncEdgeList
(
mesh(),
isZoneEdge,
orEqOp<unsigned int>(),
0
);
if (debug)
{
Info<< "On cellZone " << cz.name()
<< " marked " << returnReduce(nPoints, sumOp<label>())
<< " points and " << returnReduce(nEdges, sumOp<label>())
<< " edges." << endl;
}
isZonePoint.set(mesh().cellPoints(celli));
isZoneEdge.set(mesh().cellEdges(celli));
}
syncTools::syncPointList
(
mesh(),
isZonePoint,
orEqOp<unsigned int>(),
0
);
syncTools::syncEdgeList
(
mesh(),
isZoneEdge,
orEqOp<unsigned int>(),
0
);
DebugInfo
<< "On cellZone " << cz.name() << " marked "
<< returnReduce(isZonePoint.count(), sumOp<label>()) << " points and "
<< returnReduce(isZoneEdge.count(), sumOp<label>()) << " edges" << nl;
}
@ -140,19 +141,23 @@ void Foam::displacementLayeredMotionMotionSolver::walkStructured
const labelList& seedPoints,
const vectorField& seedData,
scalarField& distance,
vectorField& data
vectorField& data,
labelField& patchPoints
) const
{
List<pointEdgeStructuredWalk> seedInfo(seedPoints.size());
forAll(seedPoints, i)
{
const label pointi = seedPoints[i];
seedInfo[i] = pointEdgeStructuredWalk
(
points0()[seedPoints[i]], // location of data
points0()[seedPoints[i]], // previous location
points0()[pointi], // location of data
points0()[pointi], // previous location
0.0,
seedData[i]
seedData[i],
pointi // pass thru seed point id
);
}
@ -162,36 +167,31 @@ void Foam::displacementLayeredMotionMotionSolver::walkStructured
// Mark points inside cellZone.
// Note that we use points0, not mesh.points()
// so as not to accumulate errors.
forAll(isZonePoint, pointi)
for (const label pointi : isZonePoint)
{
if (isZonePoint[pointi])
{
allPointInfo[pointi] = pointEdgeStructuredWalk
(
points0()[pointi], // location of data
vector::max, // not valid
0.0,
Zero // passive data
);
}
allPointInfo[pointi] = pointEdgeStructuredWalk
(
points0()[pointi], // location of data
point::max, // not valid
0.0,
Zero // passive data
);
}
// Current info on edges
List<pointEdgeStructuredWalk> allEdgeInfo(mesh().nEdges());
// Mark edges inside cellZone
forAll(isZoneEdge, edgeI)
for (const label edgei : isZoneEdge)
{
if (isZoneEdge[edgeI])
{
allEdgeInfo[edgeI] = pointEdgeStructuredWalk
(
mesh().edges()[edgeI].centre(points0()), // location of data
vector::max, // not valid
0.0,
Zero
);
}
allEdgeInfo[edgei] = pointEdgeStructuredWalk
(
mesh().edges()[edgei].centre(points0()), // location of data
point::max, // not valid
0.0,
Zero
);
}
// Walk
@ -207,12 +207,17 @@ void Foam::displacementLayeredMotionMotionSolver::walkStructured
);
// Extract distance and passive data
forAll(allPointInfo, pointi)
for (const label pointi : isZonePoint)
{
if (isZonePoint[pointi])
const auto& pointInfo = allPointInfo[pointi];
distance[pointi] = pointInfo.dist();
data[pointi] = pointInfo.data();
// Optional information
if (patchPoints.size())
{
distance[pointi] = allPointInfo[pointi].dist();
data[pointi] = allPointInfo[pointi].data();
patchPoints[pointi] = pointInfo.index();
}
}
}
@ -229,8 +234,8 @@ Foam::displacementLayeredMotionMotionSolver::faceZoneEvaluate
const label patchi
) const
{
tmp<vectorField> tfld(new vectorField(meshPoints.size()));
vectorField& fld = tfld.ref();
auto tfld = tmp<vectorField>::New(meshPoints.size());
auto& fld = tfld.ref();
const word type(dict.get<word>("type"));
@ -265,7 +270,7 @@ Foam::displacementLayeredMotionMotionSolver::faceZoneEvaluate
// Reads name of name of patch. Then get average point displacement on
// patch. That becomes the value of fld.
const word patchName(dict.get<word>("patch"));
label patchID = mesh().boundaryMesh().findPatchID(patchName);
const label patchID = mesh().boundaryMesh().findPatchID(patchName);
pointField pdf
(
pointDisplacement_.boundaryField()[patchID].patchInternalField()
@ -275,9 +280,11 @@ Foam::displacementLayeredMotionMotionSolver::faceZoneEvaluate
else
{
FatalIOErrorInFunction(*this)
<< "Unknown faceZonePatch type " << type << " for faceZone "
<< fz.name() << exit(FatalIOError);
<< "Unknown faceZonePatch type " << type
<< " for faceZone " << fz.name() << nl
<< exit(FatalIOError);
}
return tfld;
}
@ -288,8 +295,7 @@ void Foam::displacementLayeredMotionMotionSolver::cellZoneSolve
const dictionary& zoneDict
)
{
bitSet isZonePoint(mesh().nPoints());
bitSet isZoneEdge(mesh().nEdges());
bitSet isZonePoint, isZoneEdge;
calcZoneMask(cellZoneI, isZonePoint, isZoneEdge);
const dictionary& patchesDict = zoneDict.subDict("boundaryField");
@ -303,6 +309,7 @@ void Foam::displacementLayeredMotionMotionSolver::cellZoneSolve
<< exit(FatalIOError);
}
PtrList<labelField> patchPoints(patchesDict.size());
PtrList<scalarField> patchDist(patchesDict.size());
PtrList<pointVectorField> patchDisp(patchesDict.size());
@ -311,7 +318,7 @@ void Foam::displacementLayeredMotionMotionSolver::cellZoneSolve
for (const entry& dEntry : patchesDict)
{
const word& faceZoneName = dEntry.keyword();
label zoneI = mesh().faceZones().findZoneID(faceZoneName);
const label zoneI = mesh().faceZones().findZoneID(faceZoneName);
if (zoneI == -1)
{
FatalIOErrorInFunction(*this)
@ -323,6 +330,7 @@ void Foam::displacementLayeredMotionMotionSolver::cellZoneSolve
// Determine the points of the faceZone within the cellZone
const faceZone& fz = mesh().faceZones()[zoneI];
patchPoints.set(patchi, new labelField(mesh().nPoints(), label(-1)));
patchDist.set(patchi, new scalarField(mesh().nPoints()));
patchDisp.set
(
@ -342,7 +350,7 @@ void Foam::displacementLayeredMotionMotionSolver::cellZoneSolve
)
);
patchi++;
++patchi;
}
@ -382,17 +390,16 @@ void Foam::displacementLayeredMotionMotionSolver::cellZoneSolve
patchi
);
if (debug)
{
Info<< "For cellZone:" << cellZoneI
<< " for faceZone:" << fz.name()
<< " nPoints:" << tseed().size()
<< " have patchField:"
<< " max:" << gMax(tseed())
<< " min:" << gMin(tseed())
<< " avg:" << gAverage(tseed())
<< endl;
}
DebugInfo
<< "For cellZone:" << cellZoneI
<< " for faceZone:" << fz.name()
<< " nPoints:" << tseed().size()
<< " have patchField:"
<< " max:" << gMax(tseed())
<< " min:" << gMin(tseed())
<< " avg:" << gAverage(tseed())
<< endl;
// Set distance and transported value
walkStructured
@ -404,13 +411,14 @@ void Foam::displacementLayeredMotionMotionSolver::cellZoneSolve
meshPoints,
tseed,
patchDist[patchi],
patchDisp[patchi]
patchDisp[patchi],
patchPoints[patchi]
);
// Implement real bc.
patchDisp[patchi].correctBoundaryConditions();
patchi++;
++patchi;
}
@ -435,14 +443,13 @@ void Foam::displacementLayeredMotionMotionSolver::cellZoneSolve
dimensionedScalar(dimLength, Zero)
);
forAll(distance, pointi)
for (const label pointi : isZonePoint)
{
scalar d1 = patchDist[0][pointi];
scalar d2 = patchDist[1][pointi];
const scalar d1 = patchDist[0][pointi];
const scalar d2 = patchDist[1][pointi];
if (d1 + d2 > SMALL)
{
scalar s = d1/(d1 + d2);
distance[pointi] = s;
distance[pointi] = d1/(d1 + d2);
}
}
@ -457,36 +464,43 @@ void Foam::displacementLayeredMotionMotionSolver::cellZoneSolve
if (interpolationScheme == "oneSided")
{
forAll(pointDisplacement_, pointi)
for (const label pointi : isZonePoint)
{
if (isZonePoint[pointi])
{
pointDisplacement_[pointi] = patchDisp[0][pointi];
}
pointDisplacement_[pointi] = patchDisp[0][pointi];
}
}
else if (interpolationScheme == "linear")
{
forAll(pointDisplacement_, pointi)
for (const label pointi : isZonePoint)
{
if (isZonePoint[pointi])
{
scalar d1 = patchDist[0][pointi];
scalar d2 = patchDist[1][pointi];
scalar s = d1/(d1 + d2 + VSMALL);
const scalar d1 = patchDist[0][pointi];
const scalar d2 = patchDist[1][pointi];
const scalar s = d1/(d1 + d2 + VSMALL);
const vector& pd1 = patchDisp[0][pointi];
const vector& pd2 = patchDisp[1][pointi];
const vector& pd1 = patchDisp[0][pointi];
const vector& pd2 = patchDisp[1][pointi];
pointDisplacement_[pointi] = (1 - s)*pd1 + s*pd2;
}
pointDisplacement_[pointi] = (1 - s)*pd1 + s*pd2;
}
}
else if (interpolationScheme == "cylindrical")
{
const coordSystem::cylindrical& cs =
this->getCylindrical(cellZoneI, zoneDict);
// May wish to implement alternative distance calculation here
FatalErrorInFunction
<< "cylindrical interpolation not yet available" << nl
<< "coordinate system " << cs << nl
<< exit(FatalError);
}
else
{
FatalErrorInFunction
<< "Invalid interpolationScheme: " << interpolationScheme
<< ". Valid schemes are 'oneSided' and 'linear'"
<< ". Valid schemes: (oneSided linear cylindrical)" << nl
<< exit(FatalError);
}
}
@ -518,13 +532,6 @@ displacementLayeredMotionMotionSolver
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::displacementLayeredMotionMotionSolver::
~displacementLayeredMotionMotionSolver()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::tmp<Foam::pointField>
@ -553,7 +560,7 @@ void Foam::displacementLayeredMotionMotionSolver::solve()
const word& cellZoneName = dEntry.keyword();
const dictionary& regionDict = dEntry.dict();
label zoneI = mesh().cellZones().findZoneID(cellZoneName);
const label zoneI = mesh().cellZones().findZoneID(cellZoneName);
Info<< "solving for zone: " << cellZoneName << endl;
@ -586,11 +593,11 @@ void Foam::displacementLayeredMotionMotionSolver::updateMesh
forAll(points0_, pointi)
{
label oldPointi = mpm.pointMap()[pointi];
const label oldPointi = mpm.pointMap()[pointi];
if (oldPointi >= 0)
{
label masterPointi = mpm.reversePointMap()[oldPointi];
const label masterPointi = mpm.reversePointMap()[oldPointi];
if ((masterPointi != pointi))
{

View File

@ -63,14 +63,17 @@ SourceFiles
#define displacementLayeredMotionMotionSolver_H
#include "displacementMotionSolver.H"
#include "bitSet.H"
#include "cylindricalCS.H"
#include "PtrMap.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward class declarations
// Forward Declarations
class bitSet;
/*---------------------------------------------------------------------------*\
Class displacementLayeredMotionMotionSolver Declaration
@ -80,8 +83,21 @@ class displacementLayeredMotionMotionSolver
:
public displacementMotionSolver
{
// Private Data
//- Cylindrical coordinate systems (per cellZone)
PtrMap<coordSystem::cylindrical> cylSystems_;
// Private Member Functions
//- Get existing or create new cylindrical system
const coordSystem::cylindrical& getCylindrical
(
const label cellZoneI,
const dictionary& zoneDict
);
void calcZoneMask
(
const label cellZoneI,
@ -97,7 +113,8 @@ class displacementLayeredMotionMotionSolver
const labelList& seedPoints,
const vectorField& seedData,
scalarField& distance,
vectorField& data
vectorField& data,
labelField& patchPoints
) const;
tmp<vectorField> faceZoneEvaluate
@ -137,23 +154,23 @@ public:
//- Construct from polyMesh and IOdictionary
displacementLayeredMotionMotionSolver
(
const polyMesh&,
const IOdictionary&
const polyMesh& mesh,
const IOdictionary& dict
);
//- Construct from polyMesh, IOdictionary and displacement and
// zero-time fields
//- Construct from polyMesh, IOdictionary, displacement and
//- zero-time fields
displacementLayeredMotionMotionSolver
(
const polyMesh&,
const IOdictionary&,
const pointVectorField&,
const pointIOField&
const polyMesh& mesh,
const IOdictionary& dict,
const pointVectorField& pointDisplacement,
const pointIOField& points0
);
//- Destructor
~displacementLayeredMotionMotionSolver();
~displacementLayeredMotionMotionSolver() = default;
// Member Functions
@ -165,7 +182,7 @@ public:
virtual void solve();
//- Update topology
virtual void updateMesh(const mapPolyMesh&);
virtual void updateMesh(const mapPolyMesh& mpm);
};

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2016 OpenFOAM Foundation
@ -31,24 +31,25 @@ License
Foam::Ostream& Foam::operator<<
(
Foam::Ostream& os,
const Foam::pointEdgeStructuredWalk& wDist
Ostream& os,
const pointEdgeStructuredWalk& wDist
)
{
return os
<< wDist.point0_ << wDist.previousPoint_
<< wDist.dist_ << wDist.data_;
<< wDist.dist_ << wDist.data_ << wDist.index_;
}
Foam::Istream& Foam::operator>>
(
Foam::Istream& is,
Foam::pointEdgeStructuredWalk& wDist
Istream& is,
pointEdgeStructuredWalk& wDist
)
{
return is
>> wDist.point0_ >> wDist.previousPoint_
>> wDist.dist_ >> wDist.data_;
>> wDist.dist_ >> wDist.data_ >> wDist.index_;
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2016 OpenFOAM Foundation
@ -46,13 +46,9 @@ SourceFiles
namespace Foam
{
// Forward declaration of classes
// Forward declarations
class polyPatch;
class polyMesh;
// Forward declaration of friend functions and operators
class pointEdgeStructuredWalk;
Istream& operator>>(Istream&, pointEdgeStructuredWalk&);
@ -60,12 +56,12 @@ Ostream& operator<<(Ostream&, const pointEdgeStructuredWalk&);
/*---------------------------------------------------------------------------*\
Class pointEdgeStructuredWalk Declaration
Class pointEdgeStructuredWalk Declaration
\*---------------------------------------------------------------------------*/
class pointEdgeStructuredWalk
{
// Private data
// Private Data
//- Starting location
point point0_;
@ -79,6 +75,10 @@ class pointEdgeStructuredWalk
//- Passive data
vector data_;
//- Index of passive data (optional)
label index_;
// Private Member Functions
//- Evaluate distance to point.
@ -100,22 +100,29 @@ public:
//- Construct from components
inline pointEdgeStructuredWalk
(
const point&,
const point&,
const scalar,
const vector&
const point& point0,
const point& previousPoint,
const scalar dist,
const vector& data,
const label index = -1
);
// Member Functions
// Access
// Access
inline bool inZone() const;
//- True if starting point is valid (ie, not point::max)
inline bool inZone() const;
inline scalar dist() const;
//- The distance information
inline scalar dist() const;
inline const vector& data() const;
//- Tracking data
inline const vector& data() const;
//- Index (if any) associated with data
inline label index() const;
// Needed by meshWave

View File

@ -30,7 +30,6 @@ License
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// Update this with w2.
template<class TrackingData>
inline bool Foam::pointEdgeStructuredWalk::update
(
@ -45,6 +44,7 @@ inline bool Foam::pointEdgeStructuredWalk::update
dist_ = w2.dist_ + mag(point0_-w2.previousPoint_);
previousPoint_ = point0_;
data_ = w2.data_;
index_ = w2.index_;
return true;
}
@ -55,29 +55,30 @@ inline bool Foam::pointEdgeStructuredWalk::update
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Null constructor
inline Foam::pointEdgeStructuredWalk::pointEdgeStructuredWalk()
:
point0_(vector::max),
previousPoint_(vector::max),
point0_(point::max),
previousPoint_(point::max),
dist_(0),
data_(Zero)
data_(Zero),
index_(-1)
{}
// Construct from origin, distance
inline Foam::pointEdgeStructuredWalk::pointEdgeStructuredWalk
(
const point& point0,
const point& previousPoint,
const scalar dist,
const vector& data
const vector& data,
const label index
)
:
point0_(point0),
previousPoint_(previousPoint),
dist_(dist),
data_(data)
data_(data),
index_(index)
{}
@ -85,16 +86,10 @@ inline Foam::pointEdgeStructuredWalk::pointEdgeStructuredWalk
inline bool Foam::pointEdgeStructuredWalk::inZone() const
{
return point0_ != vector::max;
return point0_ != point::max;
}
//inline const Foam::point& Foam::pointEdgeStructuredWalk::previousPoint() const
//{
// return previousPoint_;
//}
inline Foam::scalar Foam::pointEdgeStructuredWalk::dist() const
{
return dist_;
@ -107,10 +102,16 @@ inline const Foam::vector& Foam::pointEdgeStructuredWalk::data() const
}
inline Foam::label Foam::pointEdgeStructuredWalk::index() const
{
return index_;
}
template<class TrackingData>
inline bool Foam::pointEdgeStructuredWalk::valid(TrackingData& td) const
{
return previousPoint_ != vector::max;
return previousPoint_ != point::max;
}