...
 
Commits (14)
add_subdirectory(common)
add_subdirectory(curvature)
add_subdirectory(dataset)
add_subdirectory(evolution)
add_subdirectory(geometry)
add_subdirectory(interfacegrid)
......@@ -9,19 +10,16 @@ add_subdirectory(stencil)
add_subdirectory(test)
set(HEADERS
algorithm.hh
brents.hh
colorfunction.hh
curvatureset.hh
curvature.hh
dataset.hh
eoc.hh
evolution.hh
flagging.hh
flagset.hh
mixedcellmapper.hh
interpolation.hh
reconstruction.hh
reconstructionset.hh
utility.hh
scheme.hh
velocity.hh
)
......
#ifndef COLORFUNCTION_HH
#define COLORFUNCTION_HH
#ifndef DUNE_VOF_COLORFUNCTION_HH
#define DUNE_VOF_COLORFUNCTION_HH
//- dune-grid includes
#include <dune/grid/common/mcmgmapper.hh>
#include <cassert>
#include <dune/vof/dataset.hh>
#include <dune/vof/dataset/dataset.hh>
#include <dune/vof/dataset/localview.hh>
// ColorFunction
......@@ -16,44 +16,54 @@ namespace Dune
namespace VoF
{
template< class GV >
struct ColorFunction
: public Dune::VoF::DataSet < GV, double >
template<class GV>
class ColorFunction
: public DataSet<GV, double>
{
using GridView = GV;
using ThisType = ColorFunction< GridView >;
using BaseType = typename Dune::VoF::DataSet< GridView, double >;
using ctype = double;
using BaseType = DataSet<GV, double>;
using ThisType = ColorFunction<GV>;
public:
ColorFunction ( const GridView &gridView ) : BaseType( gridView ) {}
using ctype = double;
using BaseType::BaseType;
using BaseType::gridView;
using BaseType::size;
void axpy ( const ctype a, ThisType &x )
auto axpy (const ctype a, const ThisType& x) -> ThisType&
{
assert( x.size() == this->size() );
assert(x.size() == size());
for ( const auto& entity : elements( this->gridView() ) )
this->operator[]( entity ) += a * x[ entity ];
for (const auto& entity : elements(gridView()))
this->operator[](entity) += a * x[entity];
return *this;
}
template< class BinaryInStream >
void write ( BinaryInStream &in ) const
template<class BinaryInStream>
void write (BinaryInStream& in) const
{
for ( const auto &entity : elements( this->gridView() ) )
in.writeDouble( this->operator[]( entity ) );
for (const auto& value : *this)
in.writeDouble(value);
}
template< class BinaryOutStream >
void read ( BinaryOutStream &out )
template<class BinaryOutStream>
void read (BinaryOutStream& out)
{
for ( const auto &entity : elements( this->gridView() ) )
out.readDouble( this->operator[]( entity ) );
for (auto& value : *this)
out.readDouble(value);
}
};
template<class GridView>
auto localView (const ColorFunction<GridView>& color) -> DataSetLocalView<GridView, double>
{
return DataSetLocalView<GridView, double>(color);
}
} // namespace VoF
} // namespace Dune
#endif // #ifndef COLORFUNCTION_HH
#endif // #ifndef DUNE_VOF_COLORFUNCTION_HH
......@@ -2,6 +2,7 @@ set(HEADERS
commoperation.hh
iterator.hh
typetraits.hh
utility.hh
)
install(FILES ${HEADERS} DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/vof/common)
#ifndef DUNE_VOF_COMMON_COMMOPERATION_HH
#define DUNE_VOF_COMMON_COMMOPERATION_HH
#include <dune/common/unused.hh>
namespace Dune
{
namespace VoF
......@@ -17,7 +20,7 @@ namespace Dune
struct Copy
{
template< class T >
T operator()( T a, T b ) const { return a; }
T operator()( T a, DUNE_UNUSED T b ) const { return a; }
};
}
......
#ifndef DUNE_VOF_COMMON_UTILITY_HH
#define DUNE_VOF_COMMON_UTILITY_HH
#include <cmath>
#include <algorithm>
#include <limits>
#include <tuple>
#include <dune/common/fvector.hh>
#include <dune/common/rangeutilities.hh>
#include <dune/common/typetraits.hh>
#include <dune/vof/stencil/heightfunction.hh>
namespace Dune
{
namespace VoF
{
// clamp
// -----
template<class K>
auto clamp(K v, K lo, K hi) -> K
{
using std::min; using std::max;
return max(min(v, hi), lo);
}
// getHeightValues
// ---------------
template<class ColorFunction>
static inline auto getHeightValues (const ColorFunction& color,
const HeightFunctionStencil<typename ColorFunction::GridView>& stencil)
-> FieldVector<typename ColorFunction::ctype, HeightFunctionStencil<typename ColorFunction::GridView>::noc>
{
using field = typename ColorFunction::ctype;
using Values = FieldVector<field, HeightFunctionStencil<typename ColorFunction::GridView>::noc>;
const field TOL = 1e-10;
Values values;
for (auto c : range(stencil.columns()))
{
if (!stencil.valid(c, 0))
continue;
field u0 = clamp(color[stencil(c, 0)], 0.0, 1.0);
values[c] += u0;
// upwards sweep
field lastU = u0;
for (auto t : range(1, 1+stencil.tup()))
{
if (!stencil.valid(c, t))
break;
field u = clamp(color[stencil(c, t)], 0.0, 1.0);
if (u > lastU - TOL && !(u > 1.0 - TOL))
break;
values[c] += u;
lastU = u;
}
// downwards sweep
lastU = u0;
for (int t : range(1, 1-stencil.tdown()))
{
if (!stencil.valid(c, -t))
break;
field u = clamp(color[stencil(c, -t)], 0.0, 1.0);
if (u < lastU + TOL && !(u < TOL))
u = 1.0;
values[c] += u;
lastU = u;
}
}
return values;
}
// getOrientation
// --------------
template<class Coordinate>
static inline auto getOrientation(const Coordinate& normal)
-> std::tuple<int, int>
{
using field = field_t<Coordinate>;
using std::abs; using std::sqrt;
int dir = 0;
field max = std::numeric_limits<field>::min();
for (auto i : range(Coordinate::dimension))
if(abs(normal[i]) > max)
{
dir = i;
max = abs(normal[i]);
}
if (abs(max - 1.0/sqrt(2.0)) < 1e-3)
dir = Coordinate::dimension-1;
return {dir, (normal[dir] > 0) ? -1.0 : 1.0};
}
} // namespace VoF
} // namespace Dune
#endif // #ifndef DUNE_VOF_COMMON_UTILITY_HH
#ifndef DUNE_VOF_CURVATURE_HH
#define DUNE_VOF_CURVATURE_HH
#include <dune/vof/curvature/cartesianheightfunctioncurvature.hh>
#include <dune/vof/curvature/swartzcurvature.hh>
#include <dune/vof/curvature/heightfunction.hh>
#include <dune/vof/curvature/swartz.hh>
namespace Dune
{
......@@ -10,10 +10,11 @@ namespace Dune
namespace VoF
{
template< class Stencils >
static inline auto curvature ( Stencils& stencils ) -> SwartzCurvature< typename Stencils::GridView, Stencils >
template<class Stencils>
static inline auto curvature (const Stencils& stencils)
-> HeightFunctionCurvature<typename Stencils::GridView, Stencils>
{
return CartesianHeightFunctionCurvature< typename Stencils::GridView, Stencils >( stencils );
return HeightFunctionCurvature<typename Stencils::GridView, Stencils>(stencils);
}
} // namespace VoF
......
set(HEADERS
cartesianheightfunctioncurvature.hh
swartzcurvature.hh
heightfunction.hh
swartz.hh
)
install(FILES ${HEADERS} DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/vof/curvature)
#ifndef DUNE_VOF_CARTESIANHEIGHTFUNCTIONCURVATURE_HH
#define DUNE_VOF_CARTESIANHEIGHTFUNCTIONCURVATURE_HH
#include <algorithm>
#include <cmath>
#include <utility>
//- dune-vof includes
#include <dune/vof/dataset.hh>
#include <dune/vof/geometry/intersect.hh>
#include <dune/vof/geometry/utility.hh>
#include <dune/vof/stencil/heightfunctionstencil.hh>
#include <dune/vof/reconstruction/modifiedyoungs.hh>
#include <dune/vof/reconstruction/heightfunction.hh>
#include <dune/vof/utility.hh>
namespace Dune
{
namespace VoF
{
// Calculate the curvature of the inferface in each mixed cell by the height function method.
/**
* \ingroup Method
* \brief set of curvatures
*
* \tparam GV grid view
* \tparam VNST vertex neighbor stencils
*/
template< class GV, class VNST >
struct CartesianHeightFunctionCurvature
: public HeightFunctionReconstruction< GV, VNST, ModifiedYoungsReconstruction< GV, VNST > >
{
private:
using BaseType = HeightFunctionReconstruction< GV, VNST, ModifiedYoungsReconstruction< GV, VNST > >;
public:
using GridView = typename BaseType::GridView;
using StencilSet = typename BaseType::StencilSet;
using Stencil = Dune::VoF::HeightFunctionStencil< GridView >;
using Entity = typename decltype(std::declval< GridView >().template begin< 0 >())::Entity;
using Coordinate = typename Entity::Geometry::GlobalCoordinate;
static constexpr int dim = GridView::dimension;
using Heights = Dune::FieldVector< double, Stencil::noc >;
using Orientation = std::tuple< int, int >;
public:
explicit CartesianHeightFunctionCurvature ( const StencilSet &vertexNeighborStencils )
: BaseType( vertexNeighborStencils ) {}
template< class ColorFunction, class ReconstructionSet, class Flags, class CurvatureSet >
void operator() ( const ColorFunction &color, const ReconstructionSet &reconstructions, const Flags &flags,
CurvatureSet &curvature, bool communicate = true ) const
{
for ( const auto& entity : elements( color.gridView(), Partitions::interiorBorder ) )
{
curvature[ entity ] = 0.0;
satisfiesConstraint_[ entity ] = 0;
if ( !flags.isMixed( entity ) )
continue;
applyLocal( entity, color, reconstructions, curvature );
}
curvature.communicate();
satisfiesConstraint_.communicate();
for ( int i = 0; i < 1; ++i )
{
CurvatureSet newCurvature ( color.gridView() );
for ( const auto& entity : elements( color.gridView(), Partitions::interiorBorder ) )
{
if ( !flags.isMixed( entity ) )
continue;
averageCurvature( entity, curvature, flags, newCurvature, (i == 0) );
}
curvature = newCurvature;
}
if ( communicate )
curvature.communicate();
}
private:
using BaseType::satisfiesConstraint_;
using BaseType::vertexStencil;
template< class ColorFunction, class ReconstructionSet, class CurvatureSet >
void applyLocal ( const Entity &entity, const ColorFunction &color, const ReconstructionSet &reconstructions,
CurvatureSet &curvature ) const
{
const Coordinate &normal = reconstructions[ entity ].innerNormal();
const Orientation orientation = this->getOrientation( normal );
const auto entityInfo = GridView::Grid::getRealImplementation( entity ).entityInfo();
const Stencil stencil ( color.gridView(), entityInfo, orientation, 1 );
Heights heights = this->getHeightValues( color, stencil );
// Constraint
double uMid = heights[ ( heights.size() - 1 ) / 2 ];
int effTdown = stencil.effectiveTdown();
if ( uMid < effTdown || uMid > effTdown + 1 )
return;
const Dune::FieldMatrix< double, dim, dim > &jac = entity.geometry().jacobianTransposed( Coordinate( 0.0 ) );
const int o = std::get< 0 >( orientation );
const double deltaX = jac[ dim - o - 1 ][ dim - o - 1 ];
// TODO: use different deltaXs for the two directions in 3D
for( std::size_t i = 0; i < decltype( stencil )::noc; ++i )
if ( heights[ i ] == 0.0 )
return;
satisfiesConstraint_[ entity ] = 1;
curvature[ entity ] = kappa< dim >( heights, deltaX );
}
template < int dimension >
auto kappa ( const Heights &heights, const double dx ) const
-> typename std::enable_if< dimension == 1, double >::type
{
return 0;
}
template < int dimension >
auto kappa ( const Heights &heights, const double dx ) const
-> typename std::enable_if< dimension == 2, double >::type
{
double Hx = ( heights[ 2 ] - heights[ 0 ] ) / 2.0;
double Hxx = ( heights[ 2 ] - 2 * heights[ 1 ] + heights[ 0 ] ) / dx;
return - Hxx / std::pow( 1.0 + Hx * Hx, 3.0 / 2.0 );
}
template < int dimension >
auto kappa ( const Heights &heights, const double dx ) const
-> typename std::enable_if< dimension == 3, double >::type
{
double Hx = ( heights[ 5 ] - heights[ 3 ] ) / 2.0;
double Hy = ( heights[ 7 ] - heights[ 1 ] ) / 2.0;
double Hxx = ( heights[ 5 ] - 2.0 * heights[ 4 ] + heights[ 3 ] ) / dx;
double Hyy = ( heights[ 7 ] - 2.0 * heights[ 4 ] + heights[ 1 ] ) / dx;
double Hxy = ( heights[ 8 ] - heights[ 2 ] - heights[ 6 ] + heights[ 0 ] ) / ( 4.0 * dx );
return - ( Hxx + Hyy + Hxx * Hy * Hy + Hyy * Hx * Hx - 2.0 * Hxy * Hx * Hy ) / ( std::pow( 1.0 + Hx * Hx + Hy * Hy, 3.0 / 2.0 ) );
}
template< class CurvatureSet, class Flags >
void averageCurvature( const Entity &entity, const CurvatureSet &curvature, const Flags &flags, CurvatureSet &newCurvature, bool firstTurn ) const
{
newCurvature[ entity ] = 0;
int n = 0;
if ( firstTurn && satisfiesConstraint_[ entity ] )
{
newCurvature[ entity ] = curvature[ entity ];
return;
}
if ( satisfiesConstraint_[ entity ] )
{
newCurvature[ entity ] += curvature[ entity ];
n++;
}
for( const auto& neighbor : vertexStencil( entity ) )
{
if ( !flags.isMixed( neighbor ) )
continue;
if ( satisfiesConstraint_[ neighbor ] )
{
newCurvature[ entity ] += curvature[ neighbor ];
n++;
}
}
if ( n > 0 )
newCurvature[ entity ] /= static_cast< double >( n );
}
};
} // namespace VoF
} // namespace Dune
#endif // #ifndef DUNE_VOF_CARTESIANHEIGHTFUNCTIONCURVATURE_HH
#ifndef DUNE_VOF_CURVATURE_HEIGHTFUNCTION_HH
#define DUNE_VOF_CURVATURE_HEIGHTFUNCTION_HH
#include <algorithm>
#include <cmath>
#include <tuple>
#include <utility>
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/common/typetraits.hh>
//- dune-vof includes
#include <dune/vof/common/utility.hh>
#include <dune/vof/dataset/curvature.hh>
#include <dune/vof/dataset/dataset.hh>
#include <dune/vof/dataset/flag.hh>
#include <dune/vof/dataset/reconstruction.hh>
#include <dune/vof/stencil/heightfunction.hh>
namespace Dune
{
namespace VoF
{
// Calculate the curvature of the inferface in each mixed cell by the height function method.
/**
* \ingroup Method
* \brief set of curvatures
*
* \tparam GV grid view
* \tparam VNST vertex neighbor stencils
*/
template< class GV, class StS>
struct HeightFunctionCurvature
{
using GridView = GV;
using Stencils = StS;
using Reconstructions = ReconstructionSet<GridView>;
using Curvatures = CurvatureSet<GridView>;
using Flags = FlagSet<GridView>;
using Curvature = typename Curvatures::DataType;
private:
using Stencil = typename Stencils::Stencil;
using Support = HeightFunctionStencil<GridView>;
using Entity = typename GridView::template Codim<0>::Entity;
using Coordinate = typename Entity::Geometry::GlobalCoordinate;
using field = field_t<Coordinate>;
using ctype = typename Curvature::value_type;
static constexpr int dim = Coordinate::dimension;
using Jacobian = FieldMatrix<field, dim, dim>;
using Heights = FieldVector<field, Support::noc>;
using Orientation = std::tuple<int, int>;
using Reconstruction = typename Reconstructions::DataType;
public:
explicit HeightFunctionCurvature (const Stencils& stencils)
: stencils_(stencils)
{}
template<class ColorFunction>
void operator() (const ColorFunction& color, const Reconstructions& reconstructions, const Flags& flags,
Curvatures &curvatures, bool communicate = true) const
{
curvatures.clear();
for (const auto& entity : elements(color.gridView(), Partitions::interiorBorder))
curvatures[entity] = applyLocal(entity, color, reconstructions[entity]);
curvatures.communicate();
Curvatures avgCurvatures (curvatures);
for (const auto& entity : elements(color.gridView(), Partitions::interiorBorder))
if (reconstructions[entity] && !curvatures[entity])
avgCurvatures[entity] = average(stencil(entity), reconstructions, curvatures);
curvatures = avgCurvatures;
if (communicate)
curvatures.communicate();
}
private:
template<class ColorFunction>
static inline auto applyLocal (const Entity& entity, const ColorFunction& color, const Reconstruction& reconstruction)
-> Curvature
{
if (!reconstruction)
return {};
const auto orientation = getOrientation(reconstruction->innerNormal());
const auto entityInfo = GridView::Grid::getRealImplementation(entity).entityInfo();
const auto support = Support(color.gridView(), entityInfo, orientation, 1);
auto values = getHeightValues(color, support);
// Constraint
auto uMid = values[(Support::noc - 1) / 2];
int effTdown = support.effectiveTdown();
if (uMid < effTdown || uMid > effTdown + 1)
return {};
const Jacobian& jac = entity.geometry().jacobianTransposed(Coordinate{});
const int o = std::get<0>(orientation);
const auto dx = jac[dim-o-1][dim-o-1];
// TODO: use different deltaXs for the two directions in 3D
if (std::any_of(values.begin(), values.end(), [](const auto& v){ return (v == 0.0); }))
return {};
return kappa<dim>(values, dx);
}
static inline auto average (const Stencil& stencil, const Reconstructions& reconstructions, const Curvatures& curvatures)
-> Curvature
{
int n = 0;
ctype kappa = 0.0;
for (const auto& neighbor : stencil)
{
if (!reconstructions[neighbor] || !curvatures[neighbor])
continue;
kappa += *curvatures[neighbor];
n++;
}
if (n > 0)
return {kappa / n};
else
return {};
}
template<int dimension>
static inline auto kappa (const Heights& heights, const field dx)
-> std::enable_if_t<dimension == 1, ctype>
{
return 0;
}
template<int dimension>
static inline auto kappa (const Heights& heights, const field dx)
-> std::enable_if_t<dimension == 2, ctype>
{
auto Hx = (heights[2] - heights[ 0 ]) / 2.0;
auto Hxx = (heights[2] - 2*heights[1] + heights[0]) / dx;
using std::pow;
return -Hxx / pow(1.0 + Hx * Hx, 3.0 / 2.0);
}
template <int dimension>
static inline auto kappa (const Heights& heights, const field dx)
-> std::enable_if_t<dimension == 3, ctype>
{
double Hx = (heights[5] - heights[3]) / 2.0;
double Hy = (heights[7] - heights[1]) / 2.0;
double Hxx = (heights[5] - 2*heights[4] + heights[3]) / dx;
double Hyy = (heights[7] - 2*heights[4] + heights[1]) / dx;
double Hxy = (heights[8] - heights[2] - heights[6] + heights[0]) / (4*dx);
using std::pow;
return -(Hxx + Hyy + Hxx*Hy*Hy + Hyy*Hx*Hx - 2*Hxy*Hx*Hy) / pow(1.0 + Hx*Hx + Hy*Hy, 3.0 / 2.0);
}
auto stencil (const Entity& entity) const -> const Stencil& { return stencils_[entity]; }
const Stencils& stencils_;
};
} // namespace VoF
} // namespace Dune
#endif // #ifndef DUNE_VOF_CURVATURE_HEIGHTFUNCTION_HH
#ifndef DUNE_VOF_CURVATURE_SWARTZ_HH
#define DUNE_VOF_CURVATURE_SWARTZ_HH
#include <cmath>
#include <utility>
#include <dune/common/unused.hh>
#include <dune/common/typetraits.hh>
#include <dune/grid/common/rangegenerators.hh>
#include <dune/vof/dataset/curvature.hh>
#include <dune/vof/dataset/flag.hh>
#include <dune/vof/dataset/reconstruction.hh>
#include <dune/vof/geometry/intersect.hh>
#include <dune/vof/geometry/polytope.hh>
#include <dune/vof/geometry/utility.hh>
namespace Dune
{
namespace VoF
{
/**
* \ingroup Method
* \brief curvature
*
* \tparam GV grid view
* \tparam ST stencils
* \tparam RS reconstruction set
* \tparam FL flags
*/
template<class GV, class StS>
struct SwartzCurvature
{
using GridView = GV;
using Stencils = StS;
using Reconstructions = ReconstructionSet<GridView>;
using Curvatures = CurvatureSet<GridView>;
using Flags = FlagSet<GridView>;
using Curvature = typename Curvatures::DataType;
private:
using Stencil = typename Stencils::Stencil;
using Entity = typename GridView::template Codim<0>::Entity;
using Coordinate = typename Entity::Geometry::GlobalCoordinate;
using real = real_t<Coordinate>;
static_assert(Coordinate::dimension == 2, "...");
public:
explicit SwartzCurvature (const Stencils& stencils)
: stencils_( stencils )
{}
template<class ColorFunction>
void operator() (DUNE_UNUSED const ColorFunction& color, const Reconstructions& reconstructions, DUNE_UNUSED const Flags& flags,
Curvatures& curvatures, bool communicate = true) const
{
curvatures.clear();
for (const auto& entity : elements(reconstructions.gridView(), Partitions::interiorBorder))
curvatures[entity] = applyLocal(entity, reconstructions);
curvatures.communicate();
Curvatures avgCurvatures(curvatures);
for ( const auto& entity : elements(reconstructions.gridView(), Partitions::interiorBorder))
if (reconstructions[entity] && !curvatures[entity])
avgCurvatures[entity] = average(stencil(entity), reconstructions, curvatures);
curvatures = avgCurvatures;
if (communicate)
curvatures.communicate();
}
private:
auto applyLocal (const Entity& entity, const Reconstructions& reconstructions) const
-> Curvature
{
auto halfSpace = reconstructions[entity];
auto interfaceEn = intersect(makePolytope(entity.geometry()), boundary(reconstructions[entity]));
auto centroidEn = centroid(interfaceEn);
auto volumeEn = volume(interfaceEn);
if (!centroidEn)
return {};
//second moment: mu = |I|^2/12
auto muEn = volumeEn * volumeEn / 12.0;
real weights = 0.0;
real kappa = 0.0;
for(const auto& neighbor1 : stencil(entity))
{
auto interfaceNb1 = intersect(makePolytope(neighbor1.geometry()), boundary(reconstructions[neighbor1]));
auto centroidNb1 = centroid(interfaceNb1);
auto volumeNb1 = volume(interfaceNb1);
if (!centroidNb1)
continue;
real muNb1 = volumeNb1 * volumeNb1 / 12.0;
Coordinate midwayNb1 = *centroidNb1 + *centroidEn;
midwayNb1 *= 0.5;
Coordinate tanNb1 = *centroidEn - *centroidNb1;
auto triangleNb1 = tanNb1.two_norm();
tanNb1 /= triangleNb1;
Coordinate nu12 = generalizedCrossProduct(tanNb1);
if (nu12 * halfSpace->innerNormal() < 0)
continue;
for(const auto& neighbor2 : stencil(entity))
{
if (neighbor1 == neighbor2)
continue;
auto interfaceNb2 = intersect(makePolytope(neighbor2.geometry()), boundary(reconstructions[neighbor2]));
auto centroidNb2 = centroid(interfaceNb2);
auto volumeNb2 = volume(interfaceNb2);
if (!centroidNb2)
continue;
real muNb2 = volumeNb2 * volumeNb2 / 12.0;
Coordinate midwayNb2 = *centroidNb2 + *centroidEn;
midwayNb2 *= 0.5;
Coordinate tanNb2 = *centroidNb2 - *centroidEn;
auto triangleNb2 = tanNb2.two_norm();
tanNb2 /= triangleNb2;
Coordinate nu23 = generalizedCrossProduct(tanNb2);
if (nu23 * halfSpace->innerNormal() < 0)
continue;
Coordinate midwayDiff = midwayNb2 - midwayNb1;
auto deltaS = midwayDiff.two_norm();
if (midwayDiff * (tanNb1 + tanNb2) < 0)
continue;
// taking all together
real M = ((muNb2 - muEn) / triangleNb2 - (muEn - muNb1) / triangleNb1) / (2.0 * deltaS);
kappa += determinant(nu12, nu23) / ((1.0 + M) * deltaS);
weights += 1.0;
}
}
if (weights > 0)
return {kappa / weights};
else
return {};
}
static inline auto average (const Stencil& stencil, const Reconstructions& reconstructions, const Curvatures& curvatures)
-> Curvature
{
int n = 0;
real kappa = 0.0;
for (const auto& neighbor : stencil)
{
if (!reconstructions[neighbor])
continue;
kappa += curvatures[neighbor].value_or(0.0);
n++;
}
if (n > 0)
return {kappa / n};
else
return {};
}
static inline auto determinant (const Coordinate& a, const Coordinate& b)
-> real_t<Coordinate>
{
return a[0] * b[1] - a[1] * b[0];
}
auto stencil (const Entity& entity) const -> const Stencil& { return stencils_[entity]; }
const Stencils& stencils_;
};
} // namespace VoF
} // namespace Dune
#endif // #ifndef DUNE_VOF_CURVATURE_SWARTZ_HH
#ifndef DUNE_VOF_SWARTZCURVATURE_HH
#define DUNE_VOF_SWARTZCURVATURE_HH
#include <algorithm>
#include <cmath>
#include <utility>
#include <dune/grid/common/rangegenerators.hh>
#include <dune/vof/geometry/intersect.hh>
#include <dune/vof/geometry/utility.hh>
namespace Dune
{
namespace VoF
{
/**
* \ingroup Method
* \brief curvature
*
* \tparam GV grid view
* \tparam ST stencils
* \tparam RS reconstruction set
* \tparam FL flags
*/
template< class GV, class ST >
struct SwartzCurvature
{
using GridView = GV;
using StencilSet = ST;
using Entity = typename decltype(std::declval< GridView >().template begin< 0 >())::Entity;
using Coordinate = typename Entity::Geometry::GlobalCoordinate;
public:
explicit SwartzCurvature ( const StencilSet &stencils )
: gridView_( stencils.gridView() ), stencils_( stencils )
{}
template< class DF, class ReconstructionSet, class Flags, class CurvatureSet >
void operator() ( const DF &color, const ReconstructionSet &reconstructions, const Flags &flags,
CurvatureSet &curvatureSet, bool communicate = true ) const
{
for ( const auto& entity : elements( gridView(), Partitions::interiorBorder ) )
{
curvatureSet[ entity ] = 0.0;
if ( !flags.isMixed( entity ) )
continue;
applyLocal( entity, reconstructions, flags, curvatureSet );
}
curvatureSet.communicate();
CurvatureSet newCurvature ( curvatureSet );
for ( const auto& entity : elements( gridView(), Partitions::interiorBorder ) )
{
if ( !flags.isMixed( entity ) )
continue;
if ( curvatureSet[ entity ] == 0.0 )
averageCurvature( entity, curvatureSet, flags, newCurvature );
}
curvatureSet = newCurvature;
if ( communicate )
curvatureSet.communicate();
}
private:
template< class ReconstructionSet, class Flags, class CurvatureSet >
void applyLocal ( const Entity &entity, const ReconstructionSet &reconstructions, const Flags &flags, CurvatureSet &curvatureSet ) const
{
double weights = 0.0;
auto interfaceEn = interface( entity, reconstructions );
Coordinate centroidEn = interfaceEn.centroid();
//second moment: mu = |I|^2/12
double muEn = interfaceEn.volume();
muEn *= muEn / 12.0;
for( const auto& neighbor1 : stencil( entity ) )
{
if ( !flags.isMixed( neighbor1 ) )
continue;
auto interfaceNb1 = interface( neighbor1, reconstructions );
Coordinate centroidNb1 = interfaceNb1.centroid();
Coordinate midwayNb1 = centroidNb1 + centroidEn;
midwayNb1 *= 0.5;
double triangleNb1 = ( centroidNb1 - centroidEn ).two_norm();
double muNb1 = interfaceNb1.volume();
muNb1 *= muNb1 / 12.0;
Coordinate tanNb1 = centroidEn - centroidNb1;
tanNb1 /= tanNb1.two_norm();
Coordinate nu12 = generalizedCrossProduct( tanNb1 );
if ( nu12 * reconstructions[ entity ].innerNormal() < 0 )
continue;
for( const auto& neighbor2 : stencil( entity ) )
{
if ( !flags.isMixed( neighbor2 ) )
continue;
if ( neighbor1 == neighbor2 )
continue;
auto interfaceNb2 = interface( neighbor2, reconstructions );
Coordinate centroidNb2 = interfaceNb2.centroid();
Coordinate midwayNb2 = centroidNb2 + centroidEn;
midwayNb2 *= 0.5;
double triangleNb2 = ( centroidNb2 - centroidEn ).two_norm();
double muNb2 = interfaceNb2.volume();
muNb2 *= muNb2 / 12.0;
Coordinate midwayDiff = midwayNb2 - midwayNb1;
double deltaS = midwayDiff.two_norm();
Coordinate tanNb2 = centroidNb2 - centroidEn;
tanNb2 /= tanNb2.two_norm();
Coordinate nu23 = generalizedCrossProduct( tanNb2 );
if ( nu23 * reconstructions[ entity ].innerNormal() < 0 )
continue;
if ( midwayDiff * ( tanNb1 + tanNb2 ) < 0 )
continue;
static_assert( GridView::dimension == 2, "Only implemented in 2D yet!" );
auto times = []( const Coordinate &a, const Coordinate &b ) -> double { return a[0] * b[1] - a[1] * b[0]; };
// taking all together
double M = ( ( muNb2 - muEn ) / triangleNb2 - ( muEn - muNb1 ) / triangleNb1 ) / ( 2.0 * deltaS );
double weight = 1.0; //interfaceNb1.volume() * interfaceNb2.volume();
curvatureSet[ entity ] += weight * times( nu12, nu23 ) / ( ( 1.0 + M ) * deltaS );
weights += weight;
}
}
if ( weights > 0 )
curvatureSet[ entity ] /= weights;
}
template< class CurvatureSet, class Flags >
void averageCurvature( const Entity &entity, const CurvatureSet &curvatureSet, const Flags &flags,
CurvatureSet &newCurvature ) const
{
int n = 0;
for ( const auto& neighbor : stencil( entity ) )
{
if ( !flags.isMixed( neighbor ) )
continue;
newCurvature[ entity ] += curvatureSet[ neighbor ];
n++;
}
if ( n > 0 )
newCurvature[ entity ] /= static_cast< double >( n );
}
const GridView &gridView () const { return gridView_; }
const auto &stencil ( const Entity &entity ) const { return stencils_[ entity ]; }
GridView gridView_;
const StencilSet &stencils_;
};
} // namespace VoF
} // namespace Dune
#endif // #ifndef DUNE_VOF_SWARTZCURVATURE_HH
#ifndef DUNE_VOF_DATASET__HH
#define DUNE_VOF_DATASET__HH
#include <algorithm>
#include <vector>
#include <dune/grid/common/gridenums.hh>
#include <dune/grid/common/datahandleif.hh>
namespace Dune
{
namespace VoF
{
// DataSet
// -------
/**
* \ingroup Other
* \brief storage of grid data
*
* \tparam GV grid view
* \tparam T data type
*/
template< class GV, class T >
struct DataSet
{
using GridView = GV;
using DataType = T;
using Entity = typename GridView::template Codim< 0 >::Entity;
private:
using IndexSet = typename GridView::IndexSet;
public:
using iterator = typename std::vector< DataType >::iterator;
using const_iterator = typename std::vector< DataType >::const_iterator;
using Index = typename IndexSet::IndexType;
/**
* \brief MPI communication handler
*/
template < class Reduce >
struct Exchange;
explicit DataSet ( GridView gridView )
: gridView_( gridView ),
dataSet_( indexSet().size( 0 ) )
{}
const DataType& operator[] ( const Entity &entity ) const { return dataSet_[ indexSet().index( entity ) ]; }
DataType& operator[] ( const Entity &entity ) { return dataSet_[ indexSet().index( entity ) ]; }
const DataType& operator[] ( const Index &index ) const { return dataSet_[ index ]; }
DataType& operator[] ( const Index &index ) { return dataSet_[ index ]; }
iterator begin () { return dataSet_.begin(); }
iterator end () { return dataSet_.end(); }
const_iterator begin () const { return dataSet_.begin(); }
const_iterator end () const { return dataSet_.end(); }
void clear() { std::fill( dataSet_.begin(), dataSet_.end(), DataType() ); }
std::size_t size() const { return dataSet_.size(); }
const GridView &gridView () const { return gridView_; }
template< class Reduce >
void communicate ( Dune::InterfaceType interface, Reduce reduce )
{
auto exchange = Exchange< Reduce > ( *this, std::move( reduce ) );
gridView_.communicate( exchange, interface, Dune::ForwardCommunication );
}
void communicate ()
{
auto reduce = [] ( DataType a, DataType b ) { return a; };
auto exchange = Exchange< decltype( reduce ) > ( *this, std::move( reduce ) );
gridView_.communicate( exchange, Dune::InteriorBorder_All_Interface, Dune::ForwardCommunication );
}
private:
const IndexSet &indexSet () const { return gridView_.indexSet(); }
GridView gridView_;
std::vector< DataType > dataSet_;
};
// Exchange class for MPI
template< class GV, class T >
template < class Reduce >
struct DataSet< GV, T >::Exchange
: public Dune::CommDataHandleIF < Exchange< Reduce >, T >
{
Exchange ( DataSet &dataSet, Reduce reduce ) : dataSet_ ( dataSet ), reduce_( std::move( reduce ) ) {}
const bool contains ( const int dim, const int codim ) const { return ( codim == 0 ); }
const bool fixedsize ( const int dim, const int codim ) const { return true; }
template < class Entity >
const size_t size ( const Entity &e ) const { return 1; }
template < class MessageBuffer, class Entity, typename std::enable_if< Entity::codimension == 0, int >::type = 0 >
void gather ( MessageBuffer &buff, const Entity &e ) const
{
buff.write( dataSet_[ e ] );
}
template < class MessageBuffer, class Entity, typename std::enable_if< Entity::codimension != 0, int >::type = 0 >
void gather ( MessageBuffer &buff, const Entity &e ) const
{}
template < class MessageBuffer, class Entity, typename std::enable_if< Entity::codimension == 0, int >::type = 0 >
void scatter ( MessageBuffer &buff, const Entity &e, std::size_t n )
{
T x ;
buff.read( x );
T &y = dataSet_[ e ];
y = reduce_( x, y );
}
template < class MessageBuffer, class Entity, typename std::enable_if< Entity::codimension != 0, int >::type = 0 >
void scatter ( MessageBuffer &buff, const Entity &e, std::size_t n )
{}
private:
DataSet &dataSet_;
Reduce reduce_;
};
} // namespace VoF
} // namespace Dune
#endif // #ifndef DUNE_VOF_DATASET__HH
#include <dune/vof/dataset/curvature.hh>
#include <dune/vof/dataset/dataset.hh>
#include <dune/vof/dataset/flag.hh>
#include <dune/vof/dataset/localview.hh>
#include <dune/vof/dataset/reconstruction.hh>
set(HEADERS
curvature.hh
dataset.hh
flag.hh
localview.hh
reconstruction.hh
)
install(FILES ${HEADERS} DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/vof/dataset)
#ifndef DUNE_VOF_CURVATURESET_HH
#define DUNE_VOF_CURVATURESET_HH
#ifndef DUNE_VOF_DATASET_CURVATURE_HH
#define DUNE_VOF_DATASET_CURVATURE_HH
#include <limits>
#include <dune/common/std/optional.hh>
#include <dune/common/typetraits.hh>
#include <dune/vof/dataset.hh>
#include <dune/vof/dataset/dataset.hh>
#include <dune/vof/dataset/localview.hh>
namespace Dune
......@@ -21,12 +25,12 @@ namespace Dune
*
* \tparam GridView grid view
*/
template< class GridView >
using CurvatureSet = DataSet< GridView, real_t< typename GridView::Grid::ctype > >;
template<class GridView>
using CurvatureSet = DataSet<GridView, Std::optional<real_t<typename GridView::template Codim<0>::Geometry::GlobalCoordinate>>>;
} // namespace VoF
} // namespace Dune
#endif // #ifndef DUNE_VOF_CURVATURESET_HH
#endif // #ifndef DUNE_VOF_DATASET_CURVATURE_HH
#ifndef DUNE_VOF_DATASET_DATASET_HH
#define DUNE_VOF_DATASET_DATASET_HH
#include <algorithm>
#include <vector>
#include <dune/common/unused.hh>
#include <dune/grid/common/gridenums.hh>
#include <dune/grid/common/datahandleif.hh>
namespace Dune
{
namespace VoF
{
// DataSet
// -------
/**
* \ingroup Other
* \brief storage of grid data
*
* \tparam GV grid view
* \tparam T data type
*/
template< class GV, class T >
struct DataSet
{
using GridView = GV;
using DataType = T;
using Entity = typename GridView::template Codim< 0 >::Entity;
private:
using IndexSet = typename GridView::IndexSet;
public:
using iterator = typename std::vector< DataType >::iterator;
using const_iterator = typename std::vector< DataType >::const_iterator;
using Index = typename IndexSet::IndexType;
/**
* \brief MPI communication handler
*/
template < class Reduce >
struct Exchange;
explicit DataSet ( GridView gridView )
: gridView_( gridView ),
dataSet_( indexSet().size( 0 ) )
{}
const DataType& operator[] ( const Entity &entity ) const { return dataSet_[ indexSet().index( entity ) ]; }
DataType& operator[] ( const Entity &entity ) { return dataSet_[ indexSet().index( entity ) ]; }
const DataType& operator[] ( const Index &index ) const { return dataSet_[ index ]; }
DataType& operator[] ( const Index &index ) { return dataSet_[ index ]; }
iterator begin () { return dataSet_.begin(); }
iterator end () { return dataSet_.end(); }
const_iterator begin () const { return dataSet_.begin(); }
const_iterator end () const { return dataSet_.end(); }
void clear() { std::fill( dataSet_.begin(), dataSet_.end(), DataType() ); }
std::size_t size() const { return dataSet_.size(); }
const GridView &gridView () const { return gridView_; }
template< class Reduce >
void communicate ( Dune::InterfaceType interface, Reduce reduce )
{
auto exchange = Exchange< Reduce > ( *this, std::move( reduce ) );
gridView_.communicate( exchange, interface, Dune::ForwardCommunication );
}
void communicate ()
{
auto reduce = [] ( DataType a, DUNE_UNUSED DataType b ) { return a; };
auto exchange = Exchange< decltype( reduce ) > ( *this, std::move( reduce ) );
gridView_.communicate( exchange, Dune::InteriorBorder_All_Interface, Dune::ForwardCommunication );
}
private:
const IndexSet &indexSet () const { return gridView_.indexSet(); }
GridView gridView_;
std::vector< DataType > dataSet_;
};
// Exchange class for MPI
template< class GV, class T >
template < class Reduce >
struct DataSet< GV, T >::Exchange
: public Dune::CommDataHandleIF < Exchange< Reduce >, T >
{
Exchange ( DataSet &dataSet, Reduce reduce ) : dataSet_ ( dataSet ), reduce_( std::move( reduce ) ) {}
bool contains ( DUNE_UNUSED const int dim, const int codim ) const { return ( codim == 0 ); }
bool fixedsize ( DUNE_UNUSED const int dim, DUNE_UNUSED const int codim ) const { return true; }
template < class Entity >
size_t size ( DUNE_UNUSED const Entity &e ) const { return 1; }
template < class MessageBuffer, class Entity, typename std::enable_if< Entity::codimension == 0, int >::type = 0 >
void gather ( MessageBuffer &buff, const Entity &e ) const
{
buff.write( dataSet_[ e ] );
}
template < class MessageBuffer, class Entity, typename std::enable_if< Entity::codimension != 0, int >::type = 0 >
void gather ( DUNE_UNUSED MessageBuffer &buff, DUNE_UNUSED const Entity &e ) const
{}
template < class MessageBuffer, class Entity, typename std::enable_if< Entity::codimension == 0, int >::type = 0 >
void scatter ( MessageBuffer &buff, const Entity &e, DUNE_UNUSED std::size_t n )
{
T x ;
buff.read( x );
T &y = dataSet_[ e ];
y = reduce_( x, y );
}
template < class MessageBuffer, class Entity, typename std::enable_if< Entity::codimension != 0, int >::type = 0 >
void scatter ( DUNE_UNUSED MessageBuffer &buff, DUNE_UNUSED const Entity &e, DUNE_UNUSED std::size_t n )
{}
private:
DataSet &dataSet_;
Reduce reduce_;
};
} // namespace VoF
} // namespace Dune
#endif // #ifndef DUNE_VOF_DATASET_DATASET_HH
......@@ -5,7 +5,8 @@
#include <cmath>
#include <utility>
#include <dune/vof/dataset.hh>
#include <dune/vof/dataset/dataset.hh>
#include <dune/vof/dataset/localview.hh>
namespace Dune
{
......@@ -34,58 +35,56 @@ namespace Dune
* \tparam GridView grid view
*/
template< class GridView >
class FlagSet : public DataSet< GridView, Flag >
class FlagSet : public DataSet<GridView, Flag>
{
private:
using This = FlagSet< GridView >;
using Base = DataSet< GridView, Flag >;
using This = FlagSet<GridView>;
using Base = DataSet<GridView, Flag>;
using Entity = typename Base::Entity;