Merge remote-tracking branch 'Customer-VWG/wp3-directional-refinement' into develop

This commit is contained in:
mattijs
2018-10-04 13:43:33 +01:00
45 changed files with 3200 additions and 499 deletions

View File

@ -3148,7 +3148,6 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement2
bitSet refineCell(mesh_.nCells(), newCellsToRefine);
const bitSet savedRefineCell(refineCell);
label nChanged = faceConsistentRefinement(true, cellLevel_, refineCell);
{

View File

@ -1053,7 +1053,6 @@ public:
const scalar maxLoadUnbalance
);
//- Calculate list of cells to directionally refine
labelList directionalRefineCandidates
(

View File

@ -595,6 +595,9 @@ Foam::shellSurfaces::shellSurfaces
distances_.setSize(shellI);
levels_.setSize(shellI);
dirLevels_.setSize(shellI);
smoothDirection_.setSize(shellI);
nSmoothExpansion_.setSize(shellI);
nSmoothPosition_.setSize(shellI);
extendedGapLevel_.setSize(shellI);
extendedGapMode_.setSize(shellI);
@ -671,6 +674,19 @@ Foam::shellSurfaces::shellSurfaces
}
}
// Directional smoothing
// ~~~~~~~~~~~~~~~~~~~~~
nSmoothExpansion_[shellI] = 0;
nSmoothPosition_[shellI] = 0;
smoothDirection_[shellI] =
dict.lookupOrDefault("smoothDirection", vector::zero);
if (smoothDirection_[shellI] != vector::zero)
{
dict.lookup("nSmoothExpansion") >> nSmoothExpansion_[shellI];
dict.lookup("nSmoothPosition") >> nSmoothPosition_[shellI];
}
// Gap specification
@ -797,6 +813,24 @@ Foam::labelPairList Foam::shellSurfaces::directionalSelectLevel() const
}
const Foam::labelList& Foam::shellSurfaces::nSmoothExpansion() const
{
return nSmoothExpansion_;
}
const Foam::vectorField& Foam::shellSurfaces::smoothDirection() const
{
return smoothDirection_;
}
const Foam::labelList& Foam::shellSurfaces::nSmoothPosition() const
{
return nSmoothPosition_;
}
void Foam::shellSurfaces::findHigherLevel
(
const pointField& pt,

View File

@ -86,8 +86,20 @@ private:
//- Per shell per distance the refinement level
labelListList levels_;
//- Per shell any additional directional refinement
List<Tuple2<labelPair,labelVector>> dirLevels_;
// Directional
//- Per shell any additional directional refinement
List<Tuple2<labelPair,labelVector>> dirLevels_;
//- Per shell the smoothing direction
vectorField smoothDirection_;
//- Per shell the directional smoothing iterations
labelList nSmoothExpansion_;
//- Per shell the positional smoothing iterations
labelList nSmoothPosition_;
// Gap level refinement
@ -231,6 +243,15 @@ public:
const direction dir,
labelList& shell
) const;
//- Per shell the smoothing direction
const vectorField& smoothDirection() const;
//- Per shell the directional smoothing iterations
const labelList& nSmoothExpansion() const;
//- Per shell the positional smoothing iterations
const labelList& nSmoothPosition() const;
};

View File

@ -0,0 +1,258 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2018 OpenCFD Ltd.
\\/ 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/>.
Class
Foam::PointIntegrateData
Description
Integrate along selected edges using PointEdgeWave.
SourceFiles
PointIntegrateDataI.H
\*---------------------------------------------------------------------------*/
#ifndef PointIntegrateData_H
#define PointIntegrateData_H
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
class Istream;
class Ostream;
template<class DataType>
class PointIntegrateData;
// Forward declaration of friend functions and operators
template<class DataType>
Ostream& operator<<(Ostream&, const PointIntegrateData<DataType>&);
template<class DataType>
Istream& operator>>(Istream&, PointIntegrateData<DataType>&);
/*---------------------------------------------------------------------------*\
Class PointIntegrateData declaration
\*---------------------------------------------------------------------------*/
template<class DataType>
class PointIntegrateData
{
private:
// Private data
//- Valid flag
bool valid_;
//- Integrated data
DataType data_;
public:
//- Class used to pass extra data
class trackingData
{
public:
UList<DataType>& edgeData_;
trackingData(UList<DataType>& edgeData)
:
edgeData_(edgeData)
{}
};
// Constructors
//- Construct null
inline PointIntegrateData();
//- Construct from data
inline PointIntegrateData(const DataType& data);
// Member Functions
// Access
//- Const access the data
inline const DataType& data() const
{
return data_;
};
// Needed by PointEdgeWave
//- Check whether original (invalid) value.
template<class TrackingData>
inline bool valid(TrackingData& td) const;
//- Check for identical geometrical data. Used for cyclics checking.
template<class TrackingData>
inline bool sameGeometry
(
const PointIntegrateData<DataType>&,
const scalar tol,
TrackingData& td
) const;
//- Convert origin to relative vector to leaving point
// (= point coordinate)
template<class TrackingData>
inline void leaveDomain
(
const polyPatch& patch,
const label patchPointi,
const point& pos,
TrackingData& td
);
//- Convert relative origin to absolute by adding entering point
template<class TrackingData>
inline void enterDomain
(
const polyPatch& patch,
const label patchPointi,
const point& pos,
TrackingData& td
);
//- Apply rotation matrix to the data
template<class TrackingData>
inline void transform
(
const tensor& rotTensor,
TrackingData& td
);
//- Influence of edge on point
template<class TrackingData>
inline bool updatePoint
(
const polyMesh& mesh,
const label pointI,
const label edgeI,
const PointIntegrateData<DataType>& edgeInfo,
const scalar tol,
TrackingData& td
);
//- Influence of different value on same point.
// Merge new and old info.
template<class TrackingData>
inline bool updatePoint
(
const polyMesh& mesh,
const label pointI,
const PointIntegrateData<DataType>& newPointInfo,
const scalar tol,
TrackingData& td
);
//- Influence of different value on same point.
// No information about current position whatsoever.
template<class TrackingData>
inline bool updatePoint
(
const PointIntegrateData<DataType>& newPointInfo,
const scalar tol,
TrackingData& td
);
//- Influence of point on edge.
template<class TrackingData>
inline bool updateEdge
(
const polyMesh& mesh,
const label edgeI,
const label pointI,
const PointIntegrateData<DataType>& pointInfo,
const scalar tol,
TrackingData& td
);
//- Same (like operator==)
template<class TrackingData>
inline bool equal
(
const PointIntegrateData<DataType>&,
TrackingData& td
) const;
// Member Operators
// Needed for List IO
inline bool operator==(const PointIntegrateData<DataType>&) const;
inline bool operator!=(const PointIntegrateData<DataType>&) const;
// IOstream Operators
friend Ostream& operator<< <DataType>
(
Ostream&,
const PointIntegrateData<DataType>&
);
friend Istream& operator>> <DataType>
(
Istream&,
PointIntegrateData<DataType>&
);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//- Data associated with PointIntegrateData types is contiguous
template<>
inline bool contiguous<PointIntegrateData<scalar>>()
{
return true;
}
template<>
inline bool contiguous<PointIntegrateData<vector>>()
{
return true;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "PointIntegrateDataI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,304 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2018 OpenCFD Ltd.
\\/ 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 "polyMesh.H"
#include "transform.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class DataType>
inline Foam::PointIntegrateData<DataType>::PointIntegrateData()
:
valid_(false)
{}
template<class DataType>
inline Foam::PointIntegrateData<DataType>::PointIntegrateData
(
const DataType& data
)
:
valid_(true),
data_(data)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class DataType>
template<class TrackingData>
inline bool Foam::PointIntegrateData<DataType>::valid(TrackingData& td) const
{
return valid_;
}
template<class DataType>
template<class TrackingData>
inline bool Foam::PointIntegrateData<DataType>::sameGeometry
(
const PointIntegrateData<DataType>&,
const scalar tol,
TrackingData& td
) const
{
return true;
}
template<class DataType>
template<class TrackingData>
inline void Foam::PointIntegrateData<DataType>::leaveDomain
(
const polyPatch& patch,
const label patchPointi,
const point& pos,
TrackingData& td
)
{}
template<class DataType>
template<class TrackingData>
inline void Foam::PointIntegrateData<DataType>::enterDomain
(
const polyPatch& patch,
const label patchPointi,
const point& pos,
TrackingData& td
)
{}
template<class DataType>
template<class TrackingData>
inline void Foam::PointIntegrateData<DataType>::transform
(
const tensor& rotTensor,
TrackingData& td
)
{
this->data_ = Foam::transform(rotTensor, this->data_);
}
template<class DataType>
template<class TrackingData>
inline bool Foam::PointIntegrateData<DataType>::updatePoint
(
const polyMesh& mesh,
const label pointI,
const label edgeI,
const PointIntegrateData<DataType>& edgeInfo,
const scalar tol,
TrackingData& td
)
{
// Update point from an edge
if (!valid_)
{
if (!edgeInfo.valid_)
{
FatalErrorInFunction<< "edgeInfo:" << edgeInfo << exit(FatalError);
}
this->operator=(edgeInfo);
return true;
}
else
{
return false;
}
}
template<class DataType>
template<class TrackingData>
inline bool Foam::PointIntegrateData<DataType>::updatePoint
(
const polyMesh& mesh,
const label pointI,
const PointIntegrateData<DataType>& newPointInfo,
const scalar tol,
TrackingData& td
)
{
// Update point from coupled point
if (!valid_)
{
if (!newPointInfo.valid_)
{
FatalErrorInFunction<< "newPointInfo:" << newPointInfo
<< exit(FatalError);
}
this->operator=(newPointInfo);
return true;
}
else
{
return false;
}
}
template<class DataType>
template<class TrackingData>
inline bool Foam::PointIntegrateData<DataType>::updatePoint
(
const PointIntegrateData<DataType>& newPointInfo,
const scalar tol,
TrackingData& td
)
{
if (!valid_)
{
if (!newPointInfo.valid_)
{
FatalErrorInFunction<< "newPointInfo:" << newPointInfo
<< exit(FatalError);
}
this->operator=(newPointInfo);
return true;
}
else
{
return false;
}
}
template<class DataType>
template<class TrackingData>
inline bool Foam::PointIntegrateData<DataType>::updateEdge
(
const polyMesh& mesh,
const label edgeI,
const label pointI,
const PointIntegrateData<DataType>& pointInfo,
const scalar tol,
TrackingData& td
)
{
if (!valid_)
{
if (!pointInfo.valid(td))
{
FatalErrorInFunction<< "problem: invalid point:"
<< mesh.points()[pointI]
<< " data:" << pointInfo
<< exit(FatalError);
}
this->data_ = pointInfo.data_ + td.edgeData_[edgeI];
this->valid_ = true;
return true;
}
else
{
return false;
}
}
//- Same (like operator==)
template<class DataType>
template<class TrackingData>
inline bool Foam::PointIntegrateData<DataType>::equal
(
const PointIntegrateData<DataType>& pi,
TrackingData& td
) const
{
if (!valid_)
{
return false;
}
else if (!pi.valid_)
{
FatalErrorInFunction << "pi:" << pi
<< exit(FatalError);
return false;
}
else
{
return this->data_ == pi.data_;
}
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
template<class DataType>
inline bool Foam::PointIntegrateData<DataType>::operator==
(
const Foam::PointIntegrateData<DataType>& rhs
)
const
{
return this->data_ == rhs.data_;
}
template<class DataType>
inline bool Foam::PointIntegrateData<DataType>::operator!=
(
const Foam::PointIntegrateData<DataType>& rhs
)
const
{
return !(*this == rhs);
}
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
template<class DataType>
inline Foam::Ostream& Foam::operator<<
(
Ostream& os,
const PointIntegrateData<DataType>& pd
)
{
if (os.format() == IOstream::ASCII)
{
return os << pd.valid_ << token::SPACE << pd.data();
}
else
{
return os << pd.valid_ << pd.data();
}
}
template<class DataType>
inline Foam::Istream& Foam::operator>>
(
Istream& is,
PointIntegrateData<DataType>& pd
)
{
return is >> pd.valid_ >> pd.data_;
}
// ************************************************************************* //

View File

@ -1,53 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 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 "pointData.H"
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
Foam::Ostream& Foam::operator<<(Ostream& os, const pointData& wDist)
{
if (os.format() == IOstream::ASCII)
{
return os
<< static_cast<const pointEdgePoint&>(wDist)
<< token::SPACE << wDist.s() << token::SPACE << wDist.v();
}
else
{
return os
<< static_cast<const pointEdgePoint&>(wDist)
<< wDist.s() << wDist.v();
}
}
Foam::Istream& Foam::operator>>(Istream& is, pointData& wDist)
{
return is >> static_cast<pointEdgePoint&>(wDist) >> wDist.s_ >> wDist.v_;
}
// ************************************************************************* //

View File

@ -1,193 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2015 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/>.
Class
Foam::pointData
Description
Variant of pointEdgePoint with some transported additional data.
WIP - should be templated on data like wallDistData.
Passive vector v_ is not a coordinate (so no enterDomain/leaveDomain
transformation needed)
SourceFiles
pointDataI.H
pointData.C
\*---------------------------------------------------------------------------*/
#ifndef pointData_H
#define pointData_H
#include "pointEdgePoint.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of friend functions and operators
class pointData;
Istream& operator>>(Istream&, pointData&);
Ostream& operator<<(Ostream&, const pointData&);
/*---------------------------------------------------------------------------*\
Class pointData Declaration
\*---------------------------------------------------------------------------*/
class pointData
:
public pointEdgePoint
{
// Private data
//- Additional information.
scalar s_;
//- Additional information.
vector v_;
public:
// Constructors
//- Construct null
inline pointData();
//- Construct from origin, distance
inline pointData
(
const point& origin,
const scalar distSqr,
const scalar s,
const vector& v
);
//- Construct as copy
inline pointData(const pointData&);
// Member Functions
// Access
inline scalar s() const;
inline const vector& v() const;
// Needed by meshWave
//- Apply rotation matrix to origin
template<class TrackingData>
inline void transform
(
const tensor& rotTensor,
TrackingData& td
);
//- Influence of edge on point
template<class TrackingData>
inline bool updatePoint
(
const polyMesh& mesh,
const label pointi,
const label edgeI,
const pointData& edgeInfo,
const scalar tol,
TrackingData& td
);
//- Influence of different value on same point.
// Merge new and old info.
template<class TrackingData>
inline bool updatePoint
(
const polyMesh& mesh,
const label pointi,
const pointData& newPointInfo,
const scalar tol,
TrackingData& td
);
//- Influence of different value on same point.
// No information about current position whatsoever.
template<class TrackingData>
inline bool updatePoint
(
const pointData& newPointInfo,
const scalar tol,
TrackingData& td
);
//- Influence of point on edge.
template<class TrackingData>
inline bool updateEdge
(
const polyMesh& mesh,
const label edgeI,
const label pointi,
const pointData& pointInfo,
const scalar tol,
TrackingData& td
);
// Member Operators
// Needed for List IO
inline bool operator==(const pointData&) const;
inline bool operator!=(const pointData&) const;
// IOstream Operators
friend Ostream& operator<<(Ostream&, const pointData&);
friend Istream& operator>>(Istream&, pointData&);
};
//- Data associated with pointData as contiguous as pointEdgePoint
template<>
inline bool contiguous<pointData>()
{
return contiguous<pointEdgePoint>();
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "pointDataI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -1,237 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 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 "polyMesh.H"
#include "transform.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Null constructor
inline Foam::pointData::pointData()
:
pointEdgePoint(),
s_(GREAT),
v_(point::max)
{}
// Construct from origin, distance
inline Foam::pointData::pointData
(
const point& origin,
const scalar distSqr,
const scalar s,
const vector& v
)
:
pointEdgePoint(origin, distSqr),
s_(s),
v_(v)
{}
// Construct as copy
inline Foam::pointData::pointData(const pointData& wpt)
:
pointEdgePoint(wpt),
s_(wpt.s()),
v_(wpt.v())
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
inline Foam::scalar Foam::pointData::s() const
{
return s_;
}
inline const Foam::vector& Foam::pointData::v() const
{
return v_;
}
template<class TrackingData>
inline void Foam::pointData::transform
(
const tensor& rotTensor,
TrackingData& td
)
{
pointEdgePoint::transform(rotTensor, td);
v_ = Foam::transform(rotTensor, v_);
}
// Update this with information from connected edge
template<class TrackingData>
inline bool Foam::pointData::updatePoint
(
const polyMesh& mesh,
const label pointi,
const label edgeI,
const pointData& edgeInfo,
const scalar tol,
TrackingData& td
)
{
if
(
pointEdgePoint::updatePoint
(
mesh,
pointi,
edgeI,
edgeInfo,
tol,
td
)
)
{
s_ = edgeInfo.s_;
v_ = edgeInfo.v_;
return true;
}
else
{
return false;
}
}
// Update this with new information on same point
template<class TrackingData>
inline bool Foam::pointData::updatePoint
(
const polyMesh& mesh,
const label pointi,
const pointData& newPointInfo,
const scalar tol,
TrackingData& td
)
{
if
(
pointEdgePoint::updatePoint
(
mesh,
pointi,
newPointInfo,
tol,
td
)
)
{
s_ = newPointInfo.s_;
v_ = newPointInfo.v_;
return true;
}
else
{
return false;
}
}
// Update this with new information on same point. No extra information.
template<class TrackingData>
inline bool Foam::pointData::updatePoint
(
const pointData& newPointInfo,
const scalar tol,
TrackingData& td
)
{
if (pointEdgePoint::updatePoint(newPointInfo, tol, td))
{
s_ = newPointInfo.s_;
v_ = newPointInfo.v_;
return true;
}
else
{
return false;
}
}
// Update this with information from connected point
template<class TrackingData>
inline bool Foam::pointData::updateEdge
(
const polyMesh& mesh,
const label edgeI,
const label pointi,
const pointData& pointInfo,
const scalar tol,
TrackingData& td
)
{
if
(
pointEdgePoint::updateEdge
(
mesh,
edgeI,
pointi,
pointInfo,
tol,
td
)
)
{
s_ = pointInfo.s_;
v_ = pointInfo.v_;
return true;
}
else
{
return false;
}
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
inline bool Foam::pointData::operator==(const Foam::pointData& rhs)
const
{
return
pointEdgePoint::operator==(rhs)
&& (s() == rhs.s())
&& (v() == rhs.v());
}
inline bool Foam::pointData::operator!=(const Foam::pointData& rhs)
const
{
return !(*this == rhs);
}
// ************************************************************************* //

View File

@ -41,6 +41,8 @@ License
#include "labelVector.H"
#include "profiling.H"
#include "searchableSurfaces.H"
#include "fvMeshSubset.H"
#include "interpolationTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -1642,7 +1644,6 @@ Foam::label Foam::snappyRefineDriver::directionalShellRefine
for (direction dir = 0; dir < vector::nComponents; dir++)
{
// Select the cells that need to be refined in certain direction:
//
// - cell inside/outside shell
// - original cellLevel (using mapping) mentioned in levelIncrement
// - dirCellLevel not yet up to cellLevel+levelIncrement
@ -1800,6 +1801,517 @@ Foam::label Foam::snappyRefineDriver::directionalShellRefine
}
void Foam::snappyRefineDriver::mergeAndSmoothRatio
(
const scalarList& allSeedPointDist,
const label nSmoothExpansion,
List<Tuple2<scalar, scalar>>& keyAndValue
)
{
// Merge duplicate distance from coupled locations to get unique
// distances to operate on, do on master
SortableList<scalar> unmergedDist(allSeedPointDist);
DynamicList<scalar> mergedDist;
scalar prevDist = GREAT;
forAll(unmergedDist, i)
{
scalar curDist = unmergedDist[i];
scalar difference = mag(curDist - prevDist);
if (difference > meshRefiner_.mergeDistance())
//if (difference > 0.01)
{
mergedDist.append(curDist);
prevDist = curDist;
}
}
// Sort the unique distances
SortableList<scalar> sortedDist(mergedDist);
labelList indexSet = sortedDist.indices();
// Get updated position starting from original (undistorted) mesh
scalarList seedPointsNewLocation = sortedDist;
scalar initResidual = 0.0;
scalar prevIterResidual = GREAT;
for (label iter = 0; iter < nSmoothExpansion; iter++)
{
// Position based edge averaging algorithm operated on
// all seed plane locations in normalized form.
//
// 0 1 2 3 4 5 6 (edge numbers)
// ---x---x---x---x---x---x---
// 0 1 2 3 4 5 (point numbers)
//
// Average of edge 1-3 in terms of position
// = (point3 - point0)/3
// Keeping points 0-1 frozen, new position of point 2
// = position2 + (average of edge 1-3 as above)
for(label i = 2; i<mergedDist.size()-1; i++)
{
scalar oldX00 = sortedDist[i-2];
scalar oldX1 = sortedDist[i+1];
scalar curX0 = seedPointsNewLocation[i-1];
seedPointsNewLocation[i] = curX0 + (oldX1 - oldX00)/3;
}
const scalarField residual(seedPointsNewLocation-sortedDist);
{
scalar res(sumMag(residual));
if (iter == 0)
{
initResidual = res;
}
res /= initResidual;
if (mag(prevIterResidual - res) < SMALL)
{
if (debug)
{
Pout<< "Converged with iteration " << iter
<< " initResidual: " << initResidual
<< " final residual : " << res << endl;
}
break;
}
else
{
prevIterResidual = res;
}
}
// Update the field for next iteration, avoid moving points
sortedDist = seedPointsNewLocation;
}
keyAndValue.setSize(mergedDist.size());
forAll(mergedDist, i)
{
keyAndValue[i].first() = mergedDist[i];
label index = indexSet[i];
keyAndValue[i].second() = seedPointsNewLocation[index];
}
}
Foam::label Foam::snappyRefineDriver::directionalSmooth
(
const refinementParameters& refineParams
)
{
addProfiling(split, "snappyHexMesh::refine::smooth");
Info<< nl
<< "Directional expansion ratio smoothing" << nl
<< "-------------------------------------" << nl
<< endl;
fvMesh& baseMesh = meshRefiner_.mesh();
const searchableSurfaces& geometry = meshRefiner_.surfaces().geometry();
const shellSurfaces& shells = meshRefiner_.shells();
label iter = 0;
forAll(shells.nSmoothExpansion(), shellI)
{
if
(
shells.nSmoothExpansion()[shellI] > 0
|| shells.nSmoothPosition()[shellI] > 0
)
{
label surfi = shells.shells()[shellI];
const vector& userDirection = shells.smoothDirection()[shellI];
// Extract inside points
labelList pointLabels;
{
// Get inside points
List<volumeType> volType;
geometry[surfi].getVolumeType(baseMesh.points(), volType);
label nInside = 0;
forAll(volType, pointi)
{
if (volType[pointi] == volumeType::INSIDE)
{
nInside++;
}
}
pointLabels.setSize(nInside);
nInside = 0;
forAll(volType, pointi)
{
if (volType[pointi] == volumeType::INSIDE)
{
pointLabels[nInside++] = pointi;
}
}
//bitSet isInsidePoint(baseMesh.nPoints());
//forAll(volType, pointi)
//{
// if (volType[pointi] == volumeType::INSIDE)
// {
// isInsidePoint.set(pointi);
// }
//}
//pointLabels = isInsidePoint.used();
}
// Mark all directed edges
bitSet isXEdge(baseMesh.edges().size());
{
const edgeList& edges = baseMesh.edges();
forAll(edges, edgei)
{
const edge& e = edges[edgei];
vector eVec(e.vec(baseMesh.points()));
eVec /= mag(eVec);
if (mag(eVec&userDirection) > 0.9)
{
isXEdge.set(edgei);
}
}
}
// Get the extreme of smoothing region and
// normalize all points within
const scalar totalLength =
geometry[surfi].bounds().span()
& userDirection;
const scalar startPosition =
geometry[surfi].bounds().min()
& userDirection;
scalarField normalizedPosition(pointLabels.size(), 0);
forAll(pointLabels, i)
{
label pointi = pointLabels[i];
normalizedPosition[i] =
(
((baseMesh.points()[pointi]&userDirection) - startPosition)
/ totalLength
);
}
// Sort the normalized position
labelList order;
sortedOrder(normalizedPosition, order);
DynamicList<scalar> seedPointDist;
// Select points from finest refinement (one point-per plane)
scalar prevDist = GREAT;
forAll(order, i)
{
label pointi = order[i];
scalar curDist = normalizedPosition[pointi];
if (mag(curDist - prevDist) > meshRefiner_.mergeDistance())
{
seedPointDist.append(curDist);
prevDist = curDist;
}
}
// Collect data from all processors
scalarList allSeedPointDist;
{
List<scalarList> gatheredDist(Pstream::nProcs());
gatheredDist[Pstream::myProcNo()] = seedPointDist;
Pstream::gatherList(gatheredDist);
// Combine processor lists into one big list.
allSeedPointDist =
ListListOps::combine<scalarList>
(
gatheredDist, accessOp<scalarList>()
);
}
// Pre-set the points not to smooth (after expansion)
bitSet isFrozenPoint(baseMesh.nPoints(), true);
{
scalar minSeed = min(allSeedPointDist);
Pstream::scatter(minSeed);
scalar maxSeed = max(allSeedPointDist);
Pstream::scatter(maxSeed);
forAll(normalizedPosition, posI)
{
const scalar pos = normalizedPosition[posI];
if
(
(mag(pos-minSeed) < meshRefiner_.mergeDistance())
|| (mag(pos-maxSeed) < meshRefiner_.mergeDistance())
)
{
// Boundary point: freeze
isFrozenPoint.set(pointLabels[posI]);
}
else
{
// Internal to moving region
isFrozenPoint.unset(pointLabels[posI]);
}
}
}
Info<< "Smoothing " << geometry[surfi].name() << ':' << nl
<< " Direction : " << userDirection << nl
<< " Number of points : "
<< returnReduce(pointLabels.size(), sumOp<label>())
<< " (out of " << baseMesh.globalData().nTotalPoints()
<< ")" << nl
<< " Smooth expansion iterations : "
<< shells.nSmoothExpansion()[shellI] << nl
<< " Smooth position iterations : "
<< shells.nSmoothPosition()[shellI] << nl
<< " Number of planes : "
<< allSeedPointDist.size()
<< endl;
// Make lookup from original normalized distance to new value
List<Tuple2<scalar, scalar>> keyAndValue(allSeedPointDist.size());
// Filter unique seed distances and iterate for user given steps
// or until convergence. Then get back map from old to new distance
if (Pstream::master())
{
mergeAndSmoothRatio
(
allSeedPointDist,
shells.nSmoothExpansion()[shellI],
keyAndValue
);
}
Pstream::scatter(keyAndValue);
// Construct an iterpolation table for further queries
// - although normalized values are used for query,
// it might flow out of bounds due to precision, hence clamped
const interpolationTable<scalar> table
(
keyAndValue,
bounds::repeatableBounding::CLAMP,
"undefined"
);
// Now move the points directly on the baseMesh.
pointField baseNewPoints(baseMesh.points());
forAll(pointLabels, i)
{
label pointi = pointLabels[i];
const point& curPoint = baseMesh.points()[pointi];
scalar curDist = normalizedPosition[i];
scalar newDist = table(curDist);
scalar newPosition = startPosition + newDist*totalLength;
baseNewPoints[pointi] +=
userDirection * (newPosition - (curPoint &userDirection));
}
// Moving base mesh with expansion ratio smoothing
vectorField disp(baseNewPoints-baseMesh.points());
syncTools::syncPointList
(
baseMesh,
disp,
maxMagSqrEqOp<vector>(),
point::zero
);
baseMesh.movePoints(baseMesh.points()+disp);
if (debug&meshRefinement::MESH)
{
const_cast<Time&>(baseMesh.time())++;
Pout<< "Writing directional expansion ratio smoothed"
<< " mesh to time " << meshRefiner_.timeName() << endl;
meshRefiner_.write
(
meshRefinement::debugType(debug),
meshRefinement::writeType
(
meshRefinement::writeLevel()
| meshRefinement::WRITEMESH
),
baseMesh.time().path()/meshRefiner_.timeName()
);
}
// Now we have moved the points in user specified region. Smooth
// them with neighbour points to avoid skewed cells
// Instead of moving actual mesh, operate on copy
pointField baseMeshPoints(baseMesh.points());
scalar initResidual = 0.0;
scalar prevIterResidual = GREAT;
for (iter = 0; iter < shells.nSmoothPosition()[shellI]; iter++)
{
{
const edgeList& edges = baseMesh.edges();
const labelListList& pointEdges = baseMesh.pointEdges();
pointField unsmoothedPoints(baseMeshPoints);
scalarField sumOther(baseMesh.nPoints(), 0.0);
labelList nSumOther(baseMesh.nPoints(), 0);
labelList nSumXEdges(baseMesh.nPoints(), 0);
forAll(edges, edgei)
{
const edge& e = edges[edgei];
sumOther[e[0]] +=
(unsmoothedPoints[e[1]]&userDirection);
nSumOther[e[0]]++;
sumOther[e[1]] +=
(unsmoothedPoints[e[0]]&userDirection);
nSumOther[e[1]]++;
if (isXEdge[edgei])
{
nSumXEdges[e[0]]++;
nSumXEdges[e[1]]++;
}
}
syncTools::syncPointList
(
baseMesh,
nSumXEdges,
plusEqOp<label>(),
0
);
forAll(pointLabels, i)
{
label pointi = pointLabels[i];
if (nSumXEdges[pointi] < 2)
{
// Hanging node. Remove the (single!) X edge so it
// will follow points above or below instead
const labelList& pEdges = pointEdges[pointi];
forAll(pEdges, pE)
{
label edgei = pEdges[pE];
if (isXEdge[edgei])
{
const edge& e = edges[edgei];
label otherPt = e.otherVertex(pointi);
nSumOther[pointi]--;
sumOther[pointi] -=
(
unsmoothedPoints[otherPt]
& userDirection
);
}
}
}
}
syncTools::syncPointList
(
baseMesh,
sumOther,
plusEqOp<scalar>(),
0.0
);
syncTools::syncPointList
(
baseMesh,
nSumOther,
plusEqOp<label>(),
0
);
forAll(pointLabels, i)
{
label pointi = pointLabels[i];
if ((nSumOther[pointi] >= 2) && !isFrozenPoint[pointi])
{
scalar smoothPos =
0.5
*(
(unsmoothedPoints[pointi]&userDirection)
+sumOther[pointi]/nSumOther[pointi]
);
vector& v = baseNewPoints[pointi];
v += (smoothPos-(v&userDirection))*userDirection;
}
}
const vectorField residual(baseNewPoints - baseMeshPoints);
{
scalar res(gSum(mag(residual)));
if (iter == 0)
{
initResidual = res;
}
res /= initResidual;
if (mag(prevIterResidual - res) < SMALL)
{
Info<< "Converged smoothing in iteration " << iter
<< " initResidual: " << initResidual
<< " final residual : " << res << endl;
break;
}
else
{
prevIterResidual = res;
}
}
// Just copy new location instead of moving base mesh
baseMeshPoints = baseNewPoints;
}
}
// Move mesh to new location
vectorField dispSmooth(baseMeshPoints-baseMesh.points());
syncTools::syncPointList
(
baseMesh,
dispSmooth,
maxMagSqrEqOp<vector>(),
point::zero
);
baseMesh.movePoints(baseMesh.points()+dispSmooth);
if (debug&meshRefinement::MESH)
{
const_cast<Time&>(baseMesh.time())++;
Pout<< "Writing positional smoothing iteration "
<< iter << " mesh to time " << meshRefiner_.timeName()
<< endl;
meshRefiner_.write
(
meshRefinement::debugType(debug),
meshRefinement::writeType
(
meshRefinement::writeLevel()
| meshRefinement::WRITEMESH
),
baseMesh.time().path()/meshRefiner_.timeName()
);
}
}
}
return iter;
}
void Foam::snappyRefineDriver::baffleAndSplitMesh
(
const refinementParameters& refineParams,
@ -2307,6 +2819,16 @@ void Foam::snappyRefineDriver::doRefine
100 // maxIter
);
if
(
max(meshRefiner_.shells().nSmoothExpansion()) > 0
|| max(meshRefiner_.shells().nSmoothPosition()) > 0
)
{
directionalSmooth(refineParams);
}
// Introduce baffles at surface intersections. Remove sections unreachable
// from keepPoint.
baffleAndSplitMesh

View File

@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
\\/ M anipulation | Copyright (C) 2015-2018 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -36,6 +36,8 @@ SourceFiles
#include "wordPairHashTable.H"
#include "labelList.H"
#include "scalarField.H"
#include "Tuple2.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -49,6 +51,7 @@ class snapParameters;
class meshRefinement;
class decompositionMethod;
class fvMeshDistribute;
class fvMesh;
/*---------------------------------------------------------------------------*\
Class snappyRefineDriver Declaration
@ -142,7 +145,7 @@ class snappyRefineDriver
const label maxIter
);
// Directional refinement
// Directional refinement and smoothing
//- Refine (directional) all cells inside/outside shell
label directionalShellRefine
@ -151,6 +154,18 @@ class snappyRefineDriver
const label maxIter
);
//- Calculate local edge length from cell volumes
void mergeAndSmoothRatio
(
const scalarList& allSeedPointDist,
const label nSmoothExpansion,
List<Tuple2<scalar, scalar>>& keyAndValue
);
//- Smooth the directional expansion ratio
label directionalSmooth(const refinementParameters& refineParams);
//- Add baffles and remove unreachable cells
void baffleAndSplitMesh
(
@ -229,7 +244,6 @@ public:
const refinementParameters& refineParams,
const HashTable<Pair<word>>& faceZoneToPatches
);
};

View File

@ -28,6 +28,7 @@ License
#include "meshTools.H"
#include "mapDistribute.H"
#include "flipOp.H"
#include "profiling.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -100,10 +101,9 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::projectPointsToSurface
pointField& pts
) const
{
if (debug)
{
Info<< "AMI: projecting points to surface" << endl;
}
addProfiling(ami, "AMIInterpolation::projectPointsToSurface");
DebugInfo<< "AMI: projecting points to surface" << endl;
List<pointIndexHit> nearInfo;
@ -120,8 +120,8 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::projectPointsToSurface
}
else
{
pts[i] = pts[i];
nMiss++;
// Point remains unchanged
++nMiss;
}
}
@ -148,6 +148,8 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::normaliseWeights
const scalar lowWeightTol
)
{
addProfiling(ami, "AMIInterpolation::normaliseWeights");
// Normalise the weights
wghtSum.setSize(wght.size(), 0.0);
label nLowWeight = 0;
@ -233,6 +235,8 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::agglomerate
autoPtr<mapDistribute>& tgtMap
)
{
addProfiling(ami, "AMIInterpolation::agglomerate");
label sourceCoarseSize =
(
sourceRestrictAddressing.size()
@ -329,7 +333,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::agglomerate
{
oldToNew[coarseElem] = newi;
newSubMap[newi] = coarseElem;
newi++;
++newi;
}
}
newSubMap.setSize(newi);
@ -407,7 +411,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::agglomerate
oldToNew[coarseElem] = newi;
tgtCompactMap[fineElem] = compacti;
newConstructMap[newi] = compacti++;
newi++;
++newi;
}
else
{
@ -843,6 +847,8 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
const TargetPatch& tgtPatch
)
{
addProfiling(ami, "AMIInterpolation::update");
label srcTotalSize = returnReduce(srcPatch.size(), sumOp<label>());
label tgtTotalSize = returnReduce(tgtPatch.size(), sumOp<label>());
@ -1057,6 +1063,8 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::append
const TargetPatch& tgtPatch
)
{
addProfiling(ami, "AMIInterpolation::append");
// Create a new interpolation
autoPtr<AMIInterpolation<SourcePatch, TargetPatch>> newPtr
(
@ -1283,6 +1291,8 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::interpolateToTarget
const UList<Type>& defaultValues
) const
{
addProfiling(ami, "AMIInterpolation::interpolateToTarget");
if (fld.size() != srcAddress_.size())
{
FatalErrorInFunction
@ -1368,6 +1378,8 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::interpolateToSource
const UList<Type>& defaultValues
) const
{
addProfiling(ami, "AMIInterpolation::interpolateToSource");
if (fld.size() != tgtAddress_.size())
{
FatalErrorInFunction

View File

@ -24,6 +24,7 @@ License
\*---------------------------------------------------------------------------*/
#include "faceAreaWeightAMI.H"
#include "profiling.H"
// * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * //
@ -38,6 +39,8 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::calcAddressing
label tgtFacei
)
{
addProfiling(ami, "faceAreaWeightAMI::calcAddressing");
// construct weights and addressing
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@ -123,6 +126,8 @@ bool Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::processSourceFace
List<DynamicList<scalar>>& tgtWght
)
{
addProfiling(ami, "faceAreaWeightAMI::processSourceFace");
if (tgtStartFacei == -1)
{
return false;
@ -197,6 +202,8 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::setNextFaces
bool errorOnNotFound
) const
{
addProfiling(ami, "faceAreaWeightAMI::setNextFaces");
const labelList& srcNbrFaces = this->srcPatch_.faceFaces()[srcFacei];
// initialise tgtFacei
@ -303,6 +310,8 @@ Foam::scalar Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::interArea
const label tgtFacei
) const
{
addProfiling(ami, "faceAreaWeightAMI::interArea");
scalar area = 0;
const pointField& srcPoints = this->srcPatch_.points();
@ -444,6 +453,8 @@ restartUncoveredSourceFace
List<DynamicList<scalar>>& tgtWght
)
{
addProfiling(ami, "faceAreaWeightAMI::restartUncoveredSourceFace");
// Collect all src faces with a low weight
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@ -581,6 +592,8 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::calculate
label tgtFacei
)
{
addProfiling(ami, "faceAreaWeightAMI::calculate");
bool ok =
this->initialise
(