...
 
Commits (26)
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
interpolation.hh
mixedcellmapper.hh
reconstruction.hh
reconstructionset.hh
scheme.hh
utility.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
set(HEADERS
commoperation.hh
typetraits.hh
utility.hh
)
install(FILES ${HEADERS} DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/vof/common)
#ifndef DUNE_VOF_COMMON_UTILITY_HH
#define DUNE_VOF_COMMON_UTILITY_HH
#include <algorithm>
namespace Dune
{
namespace VoF
{
template<class K>
auto clamp(K v, K lo, K hi) -> K
{
using std::min; using std::max;
return max(min(v, hi), lo);
}
} // namespace VoF
} // namespace Dune
#endif // #ifndef DUNE_VOF_COMMON_UTILITY_HH
......@@ -6,13 +6,13 @@
#include <utility>
//- dune-vof includes
#include <dune/vof/dataset.hh>
#include <dune/vof/dataset/dataset.hh>
#include <dune/vof/geometry/intersect.hh>
#include <dune/vof/geometry/utility.hh>
#include <dune/vof/stencil/heightfunctionstencil.hh>
#include <dune/vof/stencil/heightfunction.hh>
#include <dune/vof/reconstruction/modifiedyoungs.hh>
#include <dune/vof/reconstruction/heightfunction.hh>
#include <dune/vof/utility.hh>
#include <dune/vof/common/utility.hh>
namespace Dune
......
#ifndef DUNE_VOF_DATASET__HH
#define DUNE_VOF_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__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
......@@ -22,12 +25,12 @@ namespace Dune
*
* \tparam GridView grid view
*/
template< class GridView >
using CurvatureSet = DataSet< GridView, Std::optional<real_t<typename GridView::template Codim<0>::Geometry::GlobalCoordinate>>>;
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;
template< std::size_t Lower, std::size_t Upper >
template<std::size_t Lower, std::size_t Upper>
struct Range
{
constexpr Range () = default;
static constexpr bool contains ( std::size_t val ) { return ( val - Lower <= Upper - Lower ); }
static constexpr bool contains (std::size_t val) { return (val - Lower <= Upper - Lower); }
template< class T, class = std::enable_if_t< std::is_enum< T >{} >,
class = std::enable_if_t< std::is_integral< std::underlying_type_t< T > >{} > >
static constexpr bool contains ( T val )
template<class T, class = std::enable_if_t<std::is_enum<T>::value>,
class = std::enable_if_t<std::is_integral<std::underlying_type_t<T>>::value>>
static constexpr bool contains (T val)
{
return contains( static_cast< std::underlying_type_t< T > >( val ) );
return contains(static_cast<std::underlying_type_t<T>>(val));
}
};
public:
FlagSet ( GridView gridView )
: Base( gridView )
{}
using Base::Base;
using Base::operator[];
double operator[] ( const typename Base::Index &index ) const { return static_cast< double >( this->Base::operator[]( index ) ); }
using Empty = Range< static_cast< std::size_t >( Flag::empty ),
static_cast< std::size_t >( Flag::empty ) >;
using Mixed = Range< static_cast< std::size_t >( Flag::mixed ),
static_cast< std::size_t >( Flag::mixedfull ) >;
double operator[] (const typename Base::Index& index) const { return static_cast<double>(this->Base::operator[](index)); }
using Full = Range< static_cast< std::size_t >( Flag::full ),
static_cast< std::size_t >( Flag::full ) >;
using Empty = Range<static_cast<std::size_t>(Flag::empty), static_cast<std::size_t>(Flag::empty)>;
using Mixed = Range<static_cast<std::size_t>(Flag::mixed), static_cast<std::size_t>(Flag::mixedfull)>;
using Full = Range<static_cast<std::size_t>(Flag::full), static_cast<std::size_t>(Flag::full)>;
bool isEmpty ( const Entity& entity ) const { return this->Base::operator[]( entity ) == Flag::empty; }
bool isMixed ( const Entity& entity ) const { return inRange( entity, Mixed{} ); }
bool isFull ( const Entity& entity ) const { return this->Base::operator[]( entity ) == Flag::full; }
bool isEmpty (const Entity& entity) const { return this->operator[](entity) == Flag::empty; }
bool isMixed (const Entity& entity) const { return inRange(entity, Mixed{}); }
bool isFull (const Entity& entity) const { return this->operator[](entity) == Flag::full; }
private:
template< class _Range >
bool inRange ( const Entity& en, _Range = {} ) const
template<class _Range>
bool inRange (const Entity& en, _Range = {}) const
{
return _Range::contains( this->operator[]( en ) );
return _Range::contains(this->operator[](en));
}
};
template<class GridView>
auto localView (const FlagSet<GridView>& flags) -> DataSetLocalView<GridView, Flag>
{
return DataSetLocalView<GridView, Flag>(flags);
}
} // namespace VoF
} // namespace Dune
......
#ifndef DUNE_VOF_DATASET_LOCALVIEW_HH
#define DUNE_VOF_DATASET_LOCALVIEW_HH
#include <dune/common/std/optional.hh>
#include <dune/common/unused.hh>
#include <dune/vof/dataset/dataset.hh>
namespace Dune
{
namespace VoF
{
template<class GV, class T>
struct DataSetLocalView
{
using GridView = GV;
using DataType = T;
using Entity = typename GridView::template Codim<0>::Entity;
using Coordinate = typename Entity::Geometry::LocalCoordinate;
explicit DataSetLocalView (const DataSet<GridView, DataType>& data) : data_(data) {}
void bind (const Entity& entity) { value_ = data_[entity]; }
void unbind () { value_.reset(); }
auto operator() (DUNE_UNUSED const Coordinate& x) const -> DataType { assert(value_); return *value_; }
private:
const DataSet<GridView, T>& data_;
Std::optional<T> value_ = {};
};
template<class GV, class T>
struct DataSetLocalView<GV, Std::optional<T>>
{
using GridView = GV;
using DataType = T;
using Entity = typename GridView::template Codim<0>::Entity;
using Coordinate = typename Entity::Geometry::LocalCoordinate;
explicit DataSetLocalView (const DataSet<GridView, Std::optional<DataType>>& data) : data_(data) {}
DataSetLocalView (const DataSet<GridView, Std::optional<DataType>>& data, DataType fallback)
: data_(data), fallback_(fallback)
{}
void bind (const Entity& entity) { value_ = data_[entity].value_or(fallback_); }
void unbind () { value_.reset(); }
auto operator() (DUNE_UNUSED const Coordinate& x) const -> DataType { assert(value_); return *value_; }
private:
const DataSet<GridView, Std::optional<DataType>>& data_;
DataType fallback_ = std::numeric_limits<DataType>::quiet_NaN();
Std::optional<DataType> value_ = {};
};
template<class GridView, class T>
auto localView (const DataSet<GridView, T>& data)
-> DataSetLocalView<GridView, T>
{
return DataSetLocalView<GridView, T>(data);
}
template<class GridView, class T>
auto localView (const DataSet<GridView, Std::optional<T>>& data, T fallback)
-> DataSetLocalView<GridView, Std::optional<T>>
{
return DataSetLocalView<GridView, Std::optional<T>>(data, fallback);
}
} // namespace VoF
} // namespace Dune
#endif // #ifndef DUNE_VOF_DATASET_LOCALVIEW_HH
#ifndef DUNE_VOF_RECONSTRUCTIONSET_HH
#define DUNE_VOF_RECONSTRUCTIONSET_HH
#ifndef DUNE_VOF_DATASET_RECONSTRUCTION_HH
#define DUNE_VOF_DATASET_RECONSTRUCTION_HH
#include <dune/common/std/optional.hh>
#include <dune/vof/dataset.hh>
#include <dune/vof/dataset/dataset.hh>
#include <dune/vof/dataset/localview.hh>
#include <dune/vof/geometry/halfspace.hh>
......@@ -29,4 +30,4 @@ namespace Dune
} // namespace Dune
#endif // #ifndef DUNE_VOF_RECONSTRUCTIONSET_HH
#endif // #ifndef DUNE_VOF_DATASET_RECONSTRUCTION_HH
......@@ -7,7 +7,7 @@
#include <dune/grid/common/partitionset.hh>
#include <dune/vof/flagset.hh>
#include <dune/vof/dataset/flag.hh>
namespace Dune
{
......
......@@ -13,6 +13,7 @@ set(HEADERS
intersection.hh
intersectioniterator.hh
iterator.hh
mixedcellmapper.hh
persistentcontainer.hh
)
......
......@@ -7,15 +7,16 @@
#include <utility>
#include <vector>
#include <dune/common/rangeutilities.hh>
#include <dune/geometry/dimension.hh>
#include <dune/vof/colorfunction.hh>
#include <dune/vof/flagset.hh>
#include <dune/vof/dataset/flag.hh>
#include <dune/vof/dataset/reconstruction.hh>
#include <dune/vof/flagging.hh>
#include <dune/vof/interfacegrid/geometry.hh>
#include <dune/vof/mixedcellmapper.hh>
#include <dune/vof/reconstructionset.hh>
#include <dune/vof/utility.hh>
#include <dune/vof/interfacegrid/mixedcellmapper.hh>
namespace Dune
{
......@@ -165,6 +166,31 @@ namespace Dune
const Offsets &offsets () const { return offsets_; }
private:
template<class ReconstructionSet>
static inline void getInterfaceVertices (const ReconstructionSet& reconstructions, Vertices& vertices, Offsets& offsets)
{
vertices.clear();
offsets.clear();
std::size_t offset = 0;
offsets.push_back(offset);
for (const auto &entity : elements(reconstructions.gridView()))
{
auto segment = interface(entity, reconstructions);
if (!segment)
continue;
auto size = segment->size();
offset += size;
offsets.push_back(offset);
for (auto i : range(size))
vertices.push_back(segment->vertex(i));
}
}
Indices indices_;
Vertices vertices_;
Offsets offsets_;
......
#ifndef DUNE_VOF_MIXEDCELLMAPPER_HH
#define DUNE_VOF_MIXEDCELLMAPPER_HH
#ifndef DUNE_VOF_INTERFACEGRID_MIXEDCELLMAPPER_HH
#define DUNE_VOF_INTERFACEGRID_MIXEDCELLMAPPER_HH
#include <cassert>
#include <cstddef>
......@@ -10,7 +10,7 @@
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/grid/common/partitionset.hh>
#include <dune/vof/flagset.hh>
#include <dune/vof/dataset/flag.hh>
namespace Dune
......@@ -97,4 +97,4 @@ namespace Dune
} // namespace Dune
#endif // #ifndef DUNE_VOF_MIXEDCELLMAPPER_HH
#endif // #ifndef DUNE_VOF_INTERFACEGRID_MIXEDCELLMAPPER_HH
......@@ -4,8 +4,7 @@
#include <dune/vof/reconstruction/heightfunction.hh>
#include <dune/vof/reconstruction/modifiedswartz.hh>
#include <dune/vof/reconstruction/modifiedyoungs.hh>
#include <dune/vof/reconstruction/swartz.hh>
#include <dune/vof/reconstructionset.hh>
#include <dune/vof/dataset/reconstruction.hh>
namespace Dune
{
......
set(HEADERS
heightfunction.hh
generalheightfunction.hh
swartz.hh
modifiedswartz.hh
modifiedyoungs.hh
modifiedyoungssecondorder.hh
......
#ifndef DUNE_VOF_RECONSTRUCTION_GENERALHEIGHTFUNCTION_HH
#define DUNE_VOF_RECONSTRUCTION_GENERALHEIGHTFUNCTION_HH
#include <cmath>
#include <algorithm>
#include <limits>
#include <vector>
#include <tuple>
#include <dune/common/fmatrix.hh>
#include <dune/grid/common/rangegenerators.hh>
#include <dune/vof/geometry/algorithm.hh>
#include <dune/vof/geometry/polytope.hh>
#include <dune/vof/geometry/2d/intersect.hh>
#include <dune/vof/geometry/utility.hh>
namespace Dune
{
namespace VoF
{
/**
* \ingroup Reconstruction
* \brief general height function reconstruction operator
*
* \tparam DF discrete function type
* \tparam RS reconstruction set type
* \tparam IR initial reconstruction type
*/
template< class GV, class VSt, class IR >
struct GeneralHeightFunctionReconstruction
{
using GridView = GV;
using StencilSet = VSt;
using InitialReconstruction = IR;
private:
using Entity = typename decltype(std::declval< GridView >().template begin< 0 >())::Entity;
using Coordinate = typename Entity::Geometry::GlobalCoordinate;
using Reconstruction = HalfSpace< Coordinate >;
static constexpr int dim = Coordinate::dimension;
static constexpr int noc = static_cast< int > ( std::pow( 3, dim-1 ) );
static constexpr double length = 3;
static constexpr double width = 1;
using Line = Dune::VoF::Line< Coordinate >;
using Lines = std::array< std::array< Line, 2 >, noc >;
using Segments = std::array< std::array< std::vector< std::array< double, 3 > >, 2 >, noc >;
using Heights = Dune::FieldVector< double, noc >;
public:
explicit GeneralHeightFunctionReconstruction ( const StencilSet &vertexStencilSet ) : initializer_( vertexStencilSet ) {}
/**
* \brief (global) operator application
*
* \tparam ColorFunction
* \tparam ReconstructionSet
* \tparam Flags
* \param color color function
* \param reconstructions set of interface
* \param flags set of flags
*/
template< class ColorFunction, class ReconstructionSet, class Flags >
void operator() ( const ColorFunction &color, ReconstructionSet &reconstructions, const Flags &flags, bool communicate = true )
{
initializer_( color, reconstructions, flags );
for ( const auto &entity : elements( color.gridView(), Partitions::interiorBorder ) )
{
if ( !flags.isMixed( entity ) )
continue;
applyLocal( entity, color, flags, reconstructions );
}
if ( communicate )
reconstructions.communicate();
}
/**
* \brief (local) operator application
*
* \tparam ColorFunction
* \tparam Flags
* \tparam ReconstructionSet
* \param entity current element
* \param color color functions
* \param flags set of flags
* \param reconstructions reconstruction set
*/
template< class ColorFunction, class Flags, class ReconstructionSet >
void applyLocal ( const Entity &entity, const ColorFunction &color, const Flags &flags, ReconstructionSet &reconstructions )
{
const Coordinate &normal = reconstructions[ entity ].innerNormal();
const Coordinate centroid = interface( entity, reconstructions ).centroid();
const double h = std::pow( entity.geometry().volume(), 1.0 / dim );
Coordinate shift = generalizedCrossProduct( normal );
Lines lines;
for( int i = 0; i < 3; ++i )
{
Coordinate mid = centroid;
mid.axpy( ( i-1 ) * width * h, shift );
Coordinate up = mid;
Coordinate down = mid;
up.axpy( -length * h, normal );
down.axpy( length * h, normal );
lines[ i ][ 0 ] = Line( mid, up );
lines[ i ][ 1 ] = Line( mid, down );
}
// Get every height function value as tuple of double: length, colorvalue, levelset of centroid
Segments segments = computeSegments( color, lines, centroid, h, reconstructions[ entity ] );
Heights heights ( 0.0 );
double TOL = std::numeric_limits< double >::epsilon();
for( std::size_t i = 0; i < noc; ++i )
for( std::size_t r = 0; r < 2; ++r )
{
std::vector< std::array< double, 3 > > &segs = segments[ i ][ r ];
double lastU = segs[ 0 ][ 1 ];
heights[ i ] += lastU * segs[ 0 ][ 0 ];
for ( std::size_t s = 1; s < segs.size(); ++s )
{
double u = segs[ s ][ 1 ];
if ( r == 0 && u > lastU - TOL )
u = 0.0;
if ( r == 1 && u < lastU + TOL )
u = 1.0;
heights[ i ] += u * segs[ s ][ 0 ];
lastU = u;
}
}
Coordinate newNormal = computeNormal( heights, h, normal );
reconstructions[ entity ] = locateHalfSpace( makePolytope( entity.geometry() ), newNormal, color[ entity ] );
}
private:
template< class ColorFunction >
Segments computeSegments( const ColorFunction &color, const Lines &lines, const Coordinate& center, const double h, const Reconstruction &reconstruction ) const
{
Segments segments;
for ( const auto &entity : elements( color.gridView() ) ) // TODO: Do not run over the whole grid!
{
if ( ( entity.geometry().center() - center ).two_norm() > std::sqrt( ( length * h ) * ( length * h ) + h * h ) + h )
continue;
for ( int i = 0; i < noc; ++i )
for ( int r = 0; r < 2; ++r )
{
Line segment = getSegment( lines[ i ][ r ], entity, color.gridView() );
double segmentVolume = volume( segment );
if ( segmentVolume < std::numeric_limits< double >::epsilon() )
continue;
segments[ i ][ r ].push_back( {{ segmentVolume, color[ entity ], reconstruction.levelSet( segment.centroid() ) }} );
}
}
auto compDown = []( const std::array< double, 3 > &a, const std::array< double, 3 > &b ) -> bool { return a[ 2 ] < b[ 2 ]; };
for ( int i = 0; i < noc; ++i )
std::sort( segments[ i ][ 1 ].begin(), segments[ i ][ 1 ].end(), compDown );
auto compUp = []( const std::array< double, 3 > &a, const std::array< double, 3 > &b ) -> bool { return a[ 2 ] > b[ 2 ]; };
for ( int i = 0; i < noc; ++i )
std::sort( segments[ i ][ 0 ].begin(), segments[ i ][ 0 ].end(), compUp );
return segments;
};
#if GRIDDIM == 2
Coordinate computeNormal ( const Heights &heights, const double h, const Coordinate &normal ) const
{
double Hx;
if ( heights[ 0 ] == 0.0 )
Hx = ( heights[ 2 ] - heights[ 1 ] ) / ( width * h );
else if ( heights[ 2 ] == 0.0 )
Hx = ( heights[ 1 ] - heights[ 0 ] ) / ( width * h );
else
Hx = ( heights[ 2 ] - heights[ 0 ] ) / ( width * 2.0 * h );
Coordinate n ( { Hx, -1.0 } );
normalize( n );
// Rotate to reference frame
Coordinate newNormal ( 0 );
newNormal[ 0 ] = - n[ 0 ] * normal[ 1 ] - n[ 1 ] * normal[ 0 ];
newNormal[ 1 ] = n[ 0 ] * normal[ 0 ] - n[ 1 ] * normal[ 1 ];
normalize( newNormal );
return newNormal;
}
#elif GRIDDIM == 3
Coordinate computeNormal ( const Heights &heights, const Coordinate &normal ) const
{
double Hx;
if ( heights[ 5 ] == 0.0 )
Hx = ( heights[ 4 ] - heights[ 3 ] );
else if ( heights[ 3 ] == 0.0 )
Hx = ( heights[ 5 ] - heights[ 4 ] );
else
Hx = ( heights[ 5 ] - heights[ 3 ] ) / 2.0;
double Hy;
if ( heights[ 7 ] == 0.0 )
Hy = ( heights[ 4 ] - heights[ 1 ] );
else if ( heights[ 1 ] == 0.0 )
Hy = ( heights[ 7 ] - heights[ 4 ] );
else
Hy = ( heights[ 7 ] - heights[ 1 ] ) / 2.0;
int i = std::get< 0 >( orientation );
int j = std::get< 1 >( orientation );
Coordinate newNormal;
newNormal[ i ] = -j;
newNormal[ (i+1)%3 ] = Hx * -j;
newNormal[ (i+2)%3 ] = Hy;
normalize( newNormal );
return newNormal;
}
#endif
template< class Entity, class GridView >
Line getSegment ( const Line& line, const Entity& entity, const GridView &gridView ) const
{
Line segment( line );
for ( const auto &intersection : intersections( gridView, entity ) )
{
auto normal = intersection.centerUnitOuterNormal();
normal *= -1.0;
const auto center = intersection.geometry().center();
const HalfSpace< Coordinate > halfspace( normal, center );
segment = intersect( segment, halfspace );
if ( segment == Line() )
return Line();
}
return segment;
}
InitialReconstruction initializer_;
};
} // end of namespace VoF
} // end of namespace Dune
#endif
This diff is collapsed.
This diff is collapsed.
......@@ -11,12 +11,12 @@
#include <dune/grid/common/rangegenerators.hh>
#include <dune/vof/flagset.hh>
#include <dune/vof/dataset/flag.hh>
#include <dune/vof/dataset/reconstruction.hh>
#include <dune/vof/geometry/algorithm.hh>
#include <dune/vof/geometry/polytope.hh>
#include <dune/vof/geometry/utility.hh>
#include <dune/vof/reconstructionset.hh>
#include <dune/vof/utility.hh>
#include <dune/vof/common/utility.hh>
namespace Dune
......@@ -37,13 +37,13 @@ namespace Dune
{
using GridView = GV;
using Stencils = StS;
using Resonstructions = ReconstructionSet<GridView>;
using Reconstructions = ReconstructionSet<GridView>;
using Flags = FlagSet<GridView>;
using Reconstruction = typename Reconstructions::DataType;
private:
using Stencil = typename StencilSet::Stencil;
using Stencil = typename Stencils::Stencil;
using Entity = typename GridView::template Codim<0>::Entity;
using Coordinate = typename Entity::Geometry::GlobalCoordinate;
......@@ -128,7 +128,7 @@ namespace Dune
const real weight = 1.0 / d.two_norm2();
AtA.axpy(0.25 * weight * weight, outerProduct( d, d ));
Atb.axpy(weight * (0.5 - colorEn), d);
Atb.axpy(0.25 * weight * weight * (0.5 - colorEn), d);
}
}
......
......@@ -11,13 +11,13 @@
#include <dune/grid/common/rangegenerators.hh>
#include <dune/vof/flagset.hh>
#include <dune/vof/dataset/flag.hh>
#include <dune/vof/dataset/reconstruction.hh>
#include <dune/vof/geometry/algorithm.hh>
#include <dune/vof/geometry/polytope.hh>
#include <dune/vof/geometry/utility.hh>
#include <dune/vof/reconstructionset.hh>
#include <dune/vof/reconstruction/modifiedyoungs.hh>
#include <dune/vof/utility.hh>
#include <dune/vof/common/utility.hh>
namespace Dune
......@@ -39,7 +39,7 @@ namespace Dune
{
using GridView = GV;
using Stencils = StS;
using Resonstructions = ReconstructionSet<GridView>;
using Reconstructions = ReconstructionSet<GridView>;
using Flags = FlagSet<GridView>;
using Reconstruction = typename Reconstructions::DataType;
......@@ -49,8 +49,6 @@ namespace Dune
using Entity = typename GridView::template Codim<0>::Entity;
using Coordinate = typename Entity::Geometry::GlobalCoordinate;
static constexpr int dim = Coordinate::dimension;
using real = real_t<Coordinate>;
using field = field_t<Coordinate>;
......@@ -64,7 +62,7 @@ namespace Dune
public:
explicit ModifiedYoungsSecondOrderReconstruction (const Stencils& stencils)
: stencils_(stencils), fallbackReconstruction_(stencils)
: stencils_(stencils), fallback_(stencils)
{}
/**
......@@ -107,7 +105,7 @@ namespace Dune
{
const auto& stencil = stencils_[entity];
if (stencil.size() < std::pow(3, dim)-1)
return fallbackReconstruction.applyLocal(entity, color);
return fallback_.applyLocal(entity, color);
const auto& geometry = entity.geometry();
const auto dimRange = range(Coordinate::dimension);
......@@ -149,11 +147,11 @@ namespace Dune
return {};
normalize(normal);
return reconstruction = locateHalfSpace(makePolytope(geometry), normal, colorEn);
return locateHalfSpace(makePolytope(geometry), normal, colorEn);
}
const Stencils& stencils_;
const FirstOrderReconstruction fallbackReconstruction;
const FallbackReconstruction fallback_;
};
} // namespace VoF
......
#ifndef DUNE_VOF_RECONSTRUCTION_SWARTZ_HH
#define DUNE_VOF_RECONSTRUCTION_SWARTZ_HH
#include <cmath>
#include <numeric>
#include <vector>
#include <dune/common/exceptions.hh>
#include <dune/grid/common/partitionset.hh>
#include <dune/vof/geometry/algorithm.hh>
#include <dune/vof/geometry/intersect.hh>
#include <dune/vof/geometry/utility.hh>
namespace Dune
{
namespace VoF
{
// SwartzReconstruction
// --------------------
/**
* \ingroup Reconstruction
* \brief modified Swartz reconstruction operator
* \details Rider, W.J., Kothe, D.B., Reconstructing Volume Tracking, p. 15ff
*
* \tparam GV grid view
* \tparam StS stencils type
*/
template< class GV, class StS >
struct SwartzReconstruction
{
using GridView = GV;
using StencilSet = StS;
private:
using Stencil = typename StencilSet::Stencil;
using Entity = typename decltype(std::declval< GridView >().template begin< 0 >())::Entity;
using Coordinate = typename Entity::Geometry::GlobalCoordinate;
static constexpr int dim = Coordinate::dimension;
using ctype = typename Coordinate::value_type;
public:
explicit SwartzReconstruction ( const StencilSet &stencils, const std::size_t maxIterations = 50 )
: stencils_( stencils ), maxIterations_( maxIterations )
{}
/**
* \brief (global) operator application
*
* \tparam ColorFunction
* \tparam ReconstructionSet
* \tparam Flags set
* \param color color function
* \param reconstructions set of interface
* \param flags set of flags
*/
template< class ColorFunction, class ReconstructionSet, class Flags >
void operator() ( const ColorFunction &color, ReconstructionSet &reconstructions, const Flags &flags ) const
{
initializer_( color, reconstructions, flags );
for ( const auto &entity : elements( color.gridView(), Partitions::interiorBorder ) )
{
if ( !flags.isMixed( entity ) )
continue;
applyLocal( entity, color, flags, reconstructions );
}
reconstructions.communicate();
}
private:
/**
* \brief (local) operator application
*
* \tparam ColorFunction
* \tparam Flags
* \tparam ReconstructionSet
* \param entity current element
* \param flags set of flags
* \param color color functions
* \param reconstructions set of reconstruction
*/
template< class ColorFunction, class Flags, class ReconstructionSet >
void applyLocal ( const Entity &entity, const ColorFunction &color, const Flags &flags, ReconstructionSet &reconstructions ) const
{
auto &reconstruction = reconstructions[ entity ];
Coordinate normal = reconstruction.innerNormal();
Coordinate normalOld = normal;
const auto geometry = entity.geometry();
Coordinate center = geometry.center();
const double colorEn = clamp( color[ entity ], 0.0, 1.0 );
auto fInverse = [ &normalOld ] ( const auto &geometry, const auto color )
{
const auto polytope = makePolytope( geometry );
const auto hs = locateHalfSpace( polytope, normalOld, color );
return hs.levelSet( geometry.center() );
};
auto computeNormal = [ normalOld ] ( const Coordinate &x1, const Coordinate &x2, const double a1, const double a2, Coordinate &normal ) -> int
{
using Matrix = FieldMatrix< ctype, dim, dim >;
using Vector = FieldVector< ctype, dim >;
const Matrix X { x1 - Coordinate( 0.5 ), x2 - Coordinate( 0.5 ) };
Vector a { a1, a2 };
Vector one ( -1.0 );
Vector b, q;
X.solve( b, a );
X.solve( q, one );
Vector r1, r2;
X.mv( b, r1 );
X.mv( q, r2 );
double a_ = q * q;
double b_ = 2.0 * ( q * b );
double c_ = b * b - 1.0;
double diskr = b_ * b_ - 4.0 * a_ * c_;
if ( diskr < 1e-14 )
return 0;
else if ( std::abs( diskr ) < 1e-14 )
{
double s = - b_ / ( 2.0 * a_ );
normal += b;
normal.axpy( s, q );
return 1;
}
else
{
double s1 = ( - b_ + std::sqrt( diskr ) ) / ( 2.0 * a_ );
double s2 = ( - b_ - std::sqrt( diskr ) ) / ( 2.0 * a_ );
Coordinate guess = b;
guess.axpy( s1, q );
if ( normalOld * guess > 0 )
{
normal += guess;
return 1;
}
else
{
guess.axpy( s2 - s1, q );
normal += guess;
return 1;
}
}
};
HalfSpace< Coordinate > hsEn;
do
{
int n = 0;
normalOld = normal;
normal = Coordinate ( 0.0 );
double sigmaEn = fInverse( geometry, colorEn );
for ( const auto &neighbor : stencil( entity ) )
{
if ( !flags.isMixed( neighbor ) )
continue;
double colorNb = clamp( color[ neighbor ], 0.0, 1.0 );
const auto geometryNb = neighbor.geometry();
auto centerNb = geometryNb.center();
double sigmaNb = fInverse( geometryNb, colorNb );
n += computeNormal( center, centerNb, sigmaEn, sigmaNb, normal );
}
if( n == 0 )
return;
normalize( normal );
break; // makes second order too!
}
while ( std::acos( normal * normalOld ) > 1e-6 );
reconstruction = locateHalfSpace( makePolytope( geometry ), normal, colorEn );
}
const Stencil &stencil ( const Entity &entity ) const { return stencils_[ entity ]; }
const StencilSet &stencils_;
const std::size_t maxIterations_;
};
} // namespace VoF
} // namespace Dune
#endif // #ifndef DUNE_VOF_RECONSTRUCTION_SWARTZ_HH
#ifndef DUNE_VOF_ALGORITHM_HH
#define DUNE_VOF_ALGORITHM_HH
#ifndef DUNE_VOF_SCHEME_HH
#define DUNE_VOF_SCHEME_HH
// C++ includes
#include <utility>
......@@ -8,21 +8,24 @@
#include <dune/common/timer.hh>
// dune-vof includes
#include "test/errors.hh"
#include <dune/vof/dataset/flag.hh>