...
 
Commits (16)
#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
......@@ -9,6 +19,8 @@ namespace Dune
namespace VoF
{
// clamp
// -----
template<class K>
auto clamp(K v, K lo, K hi) -> K
......@@ -18,6 +30,89 @@ namespace Dune
}
// 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
......
#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
......
#ifndef DUNE_VOF_CARTESIANHEIGHTFUNCTIONCURVATURE_HH
#define DUNE_VOF_CARTESIANHEIGHTFUNCTIONCURVATURE_HH
#include <algorithm>
#include <cmath>
#include <utility>
//- dune-vof includes
#include <dune/vof/dataset/dataset.hh>
#include <dune/vof/geometry/intersect.hh>
#include <dune/vof/geometry/utility.hh>
#include <dune/vof/stencil/heightfunction.hh>
#include <dune/vof/reconstruction/modifiedyoungs.hh>
#include <dune/vof/reconstruction/heightfunction.hh>
#include <dune/vof/common/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
......@@ -6,7 +6,6 @@
#include <dune/common/std/optional.hh>
#include <dune/vof/geometry/intersect.hh>
#include <dune/vof/geometry/shapes.hh>
......@@ -55,15 +54,6 @@ namespace Dune
}
// return the interface
// --------------------
template<class Entity, class ReconstructionSet>
static inline auto interface (const Entity &entity, const ReconstructionSet &reconstructions)
-> decltype(auto)
{
return intersect(makePolytope(entity.geometry()), boundary(reconstructions[entity]));
}
} // namespace VoF
} // namespace Dune
......
......@@ -35,7 +35,7 @@ namespace Dune
static_assert(dimension <= dimensionworld, "dimension larger than dimensionworld.");
Line (Container vertices) : vertices_(std::move(vertices)) {}
Line (const std::vector<Coordinate>& v) : vertices_{{v[0], v[1]}} {}
Line (const std::vector<Coordinate>& v) : vertices_{{v[0], v[1]}} { assert(v.size() == 2); }
Line (const Coordinate& v0, const Coordinate& v1) : vertices_{{v0, v1}} {}
bool operator== (const This &other) const { return vertices_ == other.vertices_; }
......
......@@ -111,7 +111,7 @@ namespace Dune
inline static auto centroid (const Line<Coordinate>& line)
-> Coordinate
{
return (line.vertex(0)-line.vertex(1)) *= 0.5;
return (line.vertex(0) + line.vertex(1)) *= 0.5;
}
template<class Coordinate>
......
......@@ -9,8 +9,6 @@
#include <utility>
#include <dune/common/fvector.hh>
#include <dune/common/rangeutilities.hh>
#include <dune/common/unused.hh>
#include <dune/vof/dataset/flag.hh>
#include <dune/vof/dataset/reconstruction.hh>
......@@ -53,11 +51,10 @@ namespace Dune
using Entity = typename GridView::template Codim<0>::Entity;
using Coordinate = typename Entity::Geometry::GlobalCoordinate;
using Heights = Dune::FieldVector<double, Support::noc>;
using Orientation = std::tuple<int, int>;
static constexpr int dim = Coordinate::dimension;
static_assert(std::is_same<Stencils, typename Initializer::Stencils>::value, "...");
public:
......@@ -87,7 +84,6 @@ namespace Dune
for (const auto& entity : elements(color.gridView(), Partitions::interiorBorder))
reconstructions[entity] = applyLocal(entity, color, reconstructions[entity]);
if ( communicate )
reconstructions.communicate();
}
......@@ -114,85 +110,10 @@ namespace Dune
const auto entityInfo = GridView::Grid::getRealImplementation(entity).entityInfo();
const auto support = Support(color.gridView(), entityInfo, orientation);
Coordinate normal = computeNormal<dim>(getHeightValues(color, support), orientation);
Coordinate normal = computeNormal<Coordinate::dimension>(getHeightValues(color, support), orientation);
return locateHalfSpace(makePolytope(entity.geometry()), normal, color[entity]);
}
protected:
static inline auto getOrientation(const Coordinate &normal)
-> Orientation
{
using std::abs; using std::sqrt;
std::size_t dir = 0;
double max = std::numeric_limits< double >::min();
for (std::size_t i : range(Coordinate::dimension))
if (abs(normal[i]) > max)
{
dir = i;
max = std::abs(normal[i]);
}
// For symmetric considerations use vertical stencil in case of about 45 degrees
if (abs( max - 1.0 / sqrt(2.0)) < 1e-3)
dir = dim-1;
int sign = (normal[dir] > 0) ? -1.0 : 1.0;
return std::make_tuple(dir, sign);
}
template<class ColorFunction>
static inline auto getHeightValues (const ColorFunction& color, const Support& support)
-> Heights
{
const double TOL = 1e-10;
Heights heights;
for (auto c : range(support.columns()))
{
if (!support.valid(c, 0))
continue;
double u0 = clamp(color[support(c, 0)], 0.0, 1.0);
heights[c] += u0;
double lastU = u0;
// upwards sweep
for (auto t : range(1, support.tup()))
{
if (!support.valid(c, t))
break;
double u = clamp(color[support(c, t)], 0.0, 1.0);
if (u > lastU - TOL && !(u > 1.0 - TOL))
break;
heights[c] += u;
lastU = u;
}
lastU = u0;
// downwards sweep
for (int t : range(1, -support.tdown()))
{
if (!support.valid(c, -t))
break;
double u = clamp(color[support(c, -t)], 0.0, 1.0);
if (u < lastU + TOL && !(u < TOL))
u = 1.0;
heights[c] += u;
lastU = u;
}
}
return heights;
}
private:
template <int dimension>
static inline auto computeNormal (const Heights& heights, const Orientation& orientation)
......
......@@ -21,8 +21,8 @@
//- dune-vof includes
#include <dune/vof/colorfunction.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>
#include <dune/vof/dataset/curvature.hh>
#include <dune/vof/dataset/flag.hh>
#include <dune/vof/dataset/reconstruction.hh>
......@@ -58,11 +58,11 @@ double algorithm (const GridView& gridView, const Dune::ParameterTree& parameter
using Stencils = VoF::VertexNeighborsStencil<GridView>;
using FlagOperator = VoF::FlagOperator<GridView>;
#if GRIDDIM == 2
using CurvatureOperator = Dune::VoF::SwartzCurvature< GridView, Stencils >;
#else
using CurvatureOperator = Dune::VoF::CartesianHeightFunctionCurvature< GridView, Stencils >;
#endif
// #if GRIDDIM == 2
// using CurvatureOperator = Dune::VoF::SwartzCurvature< GridView, Stencils >;
// #else
using CurvatureOperator = Dune::VoF::HeightFunctionCurvature< GridView, Stencils >;
// #endif
using DataWriter = Dune::VTKWriter<GridView>;
using FieldInfo = VTK::FieldInfo;
......@@ -143,7 +143,7 @@ double algorithm (const GridView& gridView, const Dune::ParameterTree& parameter
normalY[entity] = normal[1];
//normalZ[entity] = normal[2];
}
vtkwriter.pwrite("vof", path.str(), "");
vtkwriter.pwrite("vof", path.str(), "", VTK::appendedraw);
}
return error;
......@@ -175,7 +175,7 @@ try {
grid.globalRefine(refineStepsForHalf * level0);
int maxLevel = 1;
int maxLevel = 3;
std::ofstream errorsFile;
errorsFile.open("errors");
......