...
 
Commits (42)
......@@ -18,6 +18,7 @@ set(HEADERS
evolution.hh
flagging.hh
flagset.hh
interpolation.hh
mixedcellmapper.hh
reconstruction.hh
reconstructionset.hh
......
set(HEADERS
commoperation.hh
typetraits.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_TYPETRAITS_HH
#define DUNE_VOF_COMMON_TYPETRAITS_HH
#include <type_traits>
namespace Dune
{
namespace VoF
{
namespace Impl
{
template<class T, class... Us>
struct _isHomogeneous;
template<class T>
struct _isHomogeneous<T> : std::true_type
{
using type = T;
};
template<class T, class ... Us>
struct _isHomogeneous<T, T, Us...> : _isHomogeneous<T, Us...> {};
template<class T, class U, class... Vs>
struct _isHomogeneous<T, U, Vs...> : std::false_type {};
}
// isHomogeneous
// -------------
/**
* \brief check if all types in a parameter pack are the same
*
* \tparam ...Ts pack of types
*/
template<class... Ts>
struct isHomogeneous : Impl::_isHomogeneous<Ts...> {};
template<>
struct isHomogeneous<> : std::false_type {};
template<class... Ts>
using isHomogeneous_t = typename isHomogeneous<Ts...>::type;
} // namespace VoF
} // namespace Dune
#endif // #ifndef DUNE_VOF_COMMON_TYPETRAITS_HH
#ifndef DUNE_VOF_CURVATURESET_HH
#define DUNE_VOF_CURVATURESET_HH
#include <dune/common/std/optional.hh>
#include <dune/common/typetraits.hh>
#include <dune/vof/dataset.hh>
......@@ -22,7 +23,7 @@ namespace Dune
* \tparam GridView grid view
*/
template< class GridView >
using CurvatureSet = DataSet< GridView, real_t< typename GridView::Grid::ctype > >;
using CurvatureSet = DataSet< GridView, Std::optional<real_t<typename GridView::template Codim<0>::Geometry::GlobalCoordinate>>>;
} // namespace VoF
......
......@@ -4,9 +4,12 @@
#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
......@@ -75,7 +78,7 @@ namespace Dune
void communicate ()
{
auto reduce = [] ( DataType a, DataType b ) { return a; };
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 );
}
......@@ -96,12 +99,12 @@ namespace Dune
{
Exchange ( DataSet &dataSet, Reduce reduce ) : dataSet_ ( dataSet ), reduce_( std::move( reduce ) ) {}
const bool contains ( const int dim, const int codim ) const { return ( codim == 0 ); }
bool contains ( DUNE_UNUSED const int dim, const int codim ) const { return ( codim == 0 ); }
const bool fixedsize ( const int dim, const int codim ) const { return true; }
bool fixedsize ( DUNE_UNUSED const int dim, DUNE_UNUSED const int codim ) const { return true; }
template < class Entity >
const size_t size ( const Entity &e ) const { return 1; }
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
......@@ -110,11 +113,11 @@ namespace Dune
}
template < class MessageBuffer, class Entity, typename std::enable_if< Entity::codimension != 0, int >::type = 0 >
void gather ( MessageBuffer &buff, const Entity &e ) const
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, std::size_t n )
void scatter ( MessageBuffer &buff, const Entity &e, DUNE_UNUSED std::size_t n )
{
T x ;
buff.read( x );
......@@ -123,7 +126,7 @@ namespace Dune
}
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 )
void scatter ( DUNE_UNUSED MessageBuffer &buff, DUNE_UNUSED const Entity &e, DUNE_UNUSED std::size_t n )
{}
private:
......
......@@ -6,6 +6,9 @@
//- dune-common includes
#include <dune/common/fvector.hh>
#include <dune/common/unused.hh>
//- dune-geometry includes
#include <dune/geometry/referenceelements.hh>
//- dune-grid includes
......@@ -93,7 +96,7 @@ namespace Dune
template< class ReconstructionSet, class Flags, class Velocity, class ColorFunction >
double applyLocal ( const Entity &entity,
const ReconstructionSet &reconstructions,
const Flags &flags,
DUNE_UNUSED const Flags &flags,
Velocity& velocity,
double time,
double deltaT,
......@@ -101,24 +104,27 @@ namespace Dune
{
auto polytope = makePolytope( entity.geometry() );
if ( flags.isMixed( entity ) )
polytope = intersect( polytope, reconstructions[ entity ] );
auto halfSpace = reconstructions[ entity ];
if ( halfSpace )
polytope = intersect( polytope, halfSpace );
for ( int i = 0; i < polytope.size(); ++i )
getNewPosition( polytope.vertex( i ), velocity, time, deltaT );
if ( polytope )
{
for ( int i = 0; i < polytope->size(); ++i )
getNewPosition( polytope->vertex( i ), velocity, time, deltaT );
double controlVolume = polytope.volume();
double controlVolume = volume( polytope );
for ( const auto &elem : elements( gridView(), Partitions::all ) )
{
const auto &geo = elem.geometry();
Polygon< Coordinate > cut = intersect( polytope, makePolytope( geo ) );
double vol = cut.volume();
color[ elem ] += vol / geo.volume();
controlVolume -= vol;
}
for ( const auto &elem : elements( gridView(), Partitions::all ) )
{
const auto &geo = elem.geometry();
double vol = volume( intersect( polytope, makePolytope( geo ) ) );
color[ elem ] += vol / geo.volume();
controlVolume -= vol;
}
assert( std::abs( controlVolume ) < 1e-14 );
assert( std::abs( controlVolume ) < 1e-14 );
}
return deltaT;
}
......
......@@ -113,15 +113,12 @@ namespace Dune
sumFluxes += intersection.geometry().volume() * abs( outerNormal * v );
v *= deltaT;
ctype flux = 0.0;
const auto &neighbor = intersection.outside();
// outflow
if ( v * outerNormal > 0 )
{
geometricFlux( intersection.inside(), intersection.geometry(), reconstructions, flags, v, flux );
auto flux = geometricFlux( intersection.inside(), intersection.geometry(), reconstructions, flags, v);
update[ entity ] -= flux / volume;
......@@ -157,20 +154,21 @@ namespace Dune
* \param flux
*/
template< class Geometry, class ReconstructionSet, class Flags >
void geometricFlux ( const Entity& entity,
auto geometricFlux ( const Entity& entity,
const Geometry& geometry,
const ReconstructionSet& reconstructions,
const Flags& flags,
const Coordinate& v,
ctype& flux ) const
const Coordinate& v) const
-> ctype
{
flux = 0.0;
auto upwind = upwindPolygon( geometry, v );
auto polytope = upwindPolygon( geometry, v );
if ( flags.isMixed( entity ) )
flux = truncVolume( upwind, reconstructions[ entity ] );
return volume( intersect( polytope, reconstructions[ entity ] ) );
else if ( flags.isFull( entity ) )
flux = upwind.volume();
return volume( polytope );
else
return 0;
}
const GridView& gridView() const { return gridView_; }
......
#ifndef DUNE_VOF_GEOMETRY_1D_INTERSECT_HH
#define DUNE_VOF_GEOMETRY_1D_INTERSECT_HH
#include <dune/common/std/optional.hh>
#include <dune/vof/geometry/halfspace.hh>
#include <dune/vof/geometry/shapes/line.hh>
#include <dune/vof/geometry/shapes/point.hh>
namespace Dune
{
namespace VoF
{
namespace Impl
{
/**
* \ingroup geo2d
* \brief implementation for an intersection between an edge and a half space
*
* \tparam Coord type of the global coordinate
* \return the intersection edge segment
*/
template<class Coordinate>
auto intersect (const Std::optional<Line<Coordinate>>& L, const Std::optional<HalfSpace<Coordinate>>& H)
-> Std::optional<Line<Coordinate>>
{
if (!L || !H)
return {};
auto l0 = H->levelSet(L->vertex(0));
auto l1 = H->levelSet(L->vertex(1));
if (l0 >= 0.0 && l1 >= 0.0)
return L;
else if (l0 < 0.0 && l1 < 0.0)
return {};
else
{
Coordinate P;
P.axpy(-l1 / (l0 - l1), L->vertex(0));
P.axpy( l0 / (l0 - l1), L->vertex(1));
if (l0 >= 0.0)
return {{L->vertex(0), P}};
else
return {{P, L->vertex(1)}};
}
}
/**
* \ingroup geo2d
* \brief implementation for an intersection between an edge and a hyperplane
*
* \tparam Coord type of the global coordinate
* \return the intersection edge segment
*/
template<class Coordinate>
auto intersect (const Std::optional<Line<Coordinate>>& L, const Std::optional<Hyperplane<Coordinate>>& H)
-> Std::optional<Point<Coordinate>>
{
if (!L || !H)
return {};
auto l0 = H->levelSet(L->vertex(0));
auto l1 = H->levelSet(L->vertex(1));
if ((l0 >= 0.0 && l1 >= 0.0) || (l0 < 0.0 && l1 < 0.0))
return {};
else
{
Coordinate P;
P.axpy(-l1 / (l0 - l1), L->vertex(0));
P.axpy( l0 / (l0 - l1), L->vertex(1));
return {{P}};
}
}
} // namespace Impl
} // namespace VoF
} // namespace Dune
#endif // #ifndef DUNE_VOF_GEOMETRY_1D_INTERSECT_HH
......@@ -3,6 +3,8 @@
#include <type_traits>
#include <dune/common/std/optional.hh>
//- local includes
#include <dune/vof/geometry/shapes/line.hh>
......@@ -18,11 +20,11 @@ namespace Dune
* \brief generate upwind polygon for 1d
*/
template< class Geomtry >
static inline auto upwindPolygon1d ( const Geomtry& geometry, const typename Geometry::GlobalCoordinate& v )
-> std::enable_if_t< Geometry::mydimension == 0, Line< typename Geometry::GlobalCoordinate > >
template<class Geometry>
static inline auto upwindPolygon1d (const Geometry& geometry, const typename Geometry::GlobalCoordinate& v)
-> std::enable_if_t<Geometry::mydimension == 0, Std::optional<Line<typename Geometry::GlobalCoordinate>>>
{
return { geometry.corner( 0 ), geometry.corner( 0 ) - v };
return {{geometry.corner(0), geometry.corner(0) - v}};
}
} // namespace VoF
......
......@@ -6,24 +6,28 @@
#include <algorithm>
#include <limits>
#include <type_traits>
#include <utility>
#include <vector>
#include <dune/common/rangeutilities.hh>
#include <dune/common/std/optional.hh>
#include <dune/common/typetraits.hh>
#include <dune/vof/geometry/1d/intersect.hh>
#include <dune/vof/geometry/halfspace.hh>
#include <dune/vof/geometry/shapes/line.hh>
#include <dune/vof/geometry/shapes/point.hh>
#include <dune/vof/geometry/shapes/polygon.hh>
namespace Dune {
namespace VoF {
/**
* \brief namespace containing specific implementations
*/
namespace __impl {
namespace Dune
{
namespace VoF
{
namespace Impl
{
/**
* \ingroup geo2d
......@@ -32,26 +36,20 @@ namespace Dune {
* \tparam Coord type of the global coordinate
* \return the intersection polygon
*/
template< class Coord >
auto intersect ( const Polygon< Coord >& first, const Polygon< Coord >& second ) -> Polygon< Coord >
template<class Coordinate>
auto intersect (const Std::optional<Polygon<Coordinate>>& P, const Std::optional<Polygon<Coordinate>>& Q)
-> std::enable_if_t<Coordinate::dimension == 2, Std::optional<Polygon<Coordinate>>>
{
if ( first.size() < 3 || second.size() < 3 )
return Polygon< Coord >();
Polygon< Coord > result = second;
if (!P || !Q)
return {};
for( int i = 0; i < first.size( 1 ); ++i )
Std::optional<Polygon<Coordinate>> result(Q);
for (int i : range(P->size(1)))
{
Coord normal = generalizedCrossProduct( first.vertex( (i+1)%first.size() ) - first.vertex( i ) );
normalize( normal );
const Coord point = first.vertex( i );
const HalfSpace< Coord > halfspace( normal, point );
result = intersect(result, makeHalfSpace(P->vertex(i), P->vertex((i+1)%P->size())));
result = intersect( result, halfspace );
if ( result == Polygon< Coord >() )
return Polygon< Coord >();
if (!result)
break;
}
return result;
......@@ -64,39 +62,40 @@ namespace Dune {
* \tparam Coord type of the global coordinate
* \return the intersection polygon
*/
template< class Coord >
auto intersect ( const Polygon< Coord >& polygon, const HalfSpace< Coord >& halfSpace ) -> Polygon< Coord >
template<class Coordinate>
auto intersect (const Std::optional<Polygon<Coordinate>>& P, const Std::optional<HalfSpace<Coordinate>>& H)
-> Std::optional<Polygon<Coordinate>>
{
if ( !halfSpace )
return Polygon< Coord >();
if (!P || !H)
return {};
auto container = typename Polygon< Coord >::Container();
container.reserve( polygon.size() );
auto container = typename Polygon<Coordinate>::Container{};
container.reserve(P->size());
for( int i = 0; i < polygon.size( 1 ); ++i )
for (auto i : range(P->size(1)))
{
auto edge = polygon.edge( i );
auto E = P->edge(i);
auto l0 = halfSpace.levelSet( edge.vertex( 0 ) );
auto l1 = halfSpace.levelSet( edge.vertex( 1 ) );
auto l0 = H->levelSet(E.vertex(0));
auto l1 = H->levelSet(E.vertex(1));
if( l0 > 0.0 )
container.push_back( edge.vertex( 0 ) );
if (l0 > 0.0)
container.push_back(E.vertex(0));
if ( ( l0 > 0.0 ) ^ ( l1 > 0.0 ) )
if ((l0 > 0.0) ^ (l1 > 0.0))
{
Coord point;
point.axpy( -l1 / ( l0 - l1 ), edge.vertex( 0 ) );
point.axpy( l0 / ( l0 - l1 ), edge.vertex( 1 ) );
Coordinate point;
point.axpy(-l1 / (l0 - l1), E.vertex(0));
point.axpy( l0 / (l0 - l1), E.vertex(1));
container.push_back( point );
container.push_back(point);
}
}
if ( container.empty() )
return Polygon< Coord >();
if (container.size() < 3)
return {};
else
return Polygon< Coord >( std::move( container ) );
return {{std::move(container)}};
}
......@@ -107,109 +106,49 @@ namespace Dune {
* \tparam Coord type of the global coordinate
* \return the intersection line
*/
template< class Coord >
auto intersect ( const Polygon< Coord >& polygon, const HyperPlane< Coord >& plane ) -> Line< Coord >
template<class Coordinate>
auto intersect (const Std::optional<Polygon<Coordinate>>& P, const Std::optional<Hyperplane<Coordinate>>& H)
-> Std::optional<Line<Coordinate>>
{
using limits = std::numeric_limits< typename Coord::value_type >;
using limits = std::numeric_limits<real_t<Coordinate>>;
using std::abs;
if ( !plane )
return Line< Coord >();
if (!P || !H)
return {};
auto container = typename Line< Coord >::Container();
auto container = typename Line<Coordinate>::Container();
std::size_t j = 0u;
for( int i = 0; i < polygon.size( 1 ); ++i )
for(auto i : range(P->size(1)))
{
if( j == 2 ) break;
if (j == 2) break;
auto edge = polygon.edge( i );
auto E = P->edge(i);
auto l0 = plane.levelSet( edge.vertex( 0 ) );
auto l1 = plane.levelSet( edge.vertex( 1 ) );
auto l0 = H->levelSet(E.vertex(0));
auto l1 = H->levelSet(E.vertex(1));
if ( ( l0 > 0.0 && l1 < 0.0 ) || ( l0 < 0.0 && l1 > 0.0 ) )
if ((l0 > 0.0 && l1 < 0.0) || (l0 < 0.0 && l1 > 0.0))
{
Coord point;
point.axpy( -l1 / ( l0 - l1 ), edge.vertex( 0 ) );
point.axpy( l0 / ( l0 - l1 ), edge.vertex( 1 ) );
Coordinate point;
point.axpy(-l1 / (l0 - l1), E.vertex(0));
point.axpy( l0 / (l0 - l1), E.vertex(1));
container[ j++ ] = point ;
container[j++] = point;
}
else if ( abs( l0 ) < limits::epsilon() )
container[ j++ ] = edge.vertex( 0 );
else if (abs(l0) < limits::epsilon())
container[j++] = E.vertex(0);
}
if ( j == 0 )
return Line< Coord >();
else if ( j == 1 )
container[ 1 ] = container[ 0 ];
return Line< Coord >( std::move( container ) );
}
/**
* \ingroup geo2d
* \brief implementation for an intersection between an edge and a half space
*
* \tparam Coord type of the global coordinate
* \return the intersection edge segment
*/
template< class Coord >
auto intersect ( const Line< Coord >& edge, const HalfSpace< Coord >& halfSpace ) -> Line< Coord >
{
if ( !halfSpace )
return Line< Coord >();
auto l0 = halfSpace.levelSet( edge.vertex( 0 ) );
auto l1 = halfSpace.levelSet( edge.vertex( 1 ) );
if ( l0 >= 0.0 && l1 >= 0.0 )
return edge;
else if ( l0 < 0.0 && l1 < 0.0 )
if (j == 0)
return {};
else
{
Coord point;
point.axpy( -l1 / ( l0 - l1 ), edge.vertex( 0 ) );
point.axpy( l0 / ( l0 - l1 ), edge.vertex( 1 ) );
if ( l0 >= 0.0 )
return { edge.vertex( 0 ), point };
else
return { point, edge.vertex( 1 ) };
}
}
/**
* \ingroup geo2d
* \brief implementation for an intersection between an edge and a hyperplane
*
* \tparam Coord type of the global coordinate
* \return the intersection edge segment
*/
template< class Coord >
Point< Coord > intersect ( const Line< Coord >& edge, const HyperPlane< Coord >& plane )
{
if ( !plane )
return Point< Coord >();
auto l0 = plane.levelSet( edge.vertex( 0 ) );
auto l1 = plane.levelSet( edge.vertex( 1 ) );
else if (j == 1)
container[1] = container[0];
if (( l0 >= 0.0 && l1 >= 0.0 ) || ( l0 < 0.0 && l1 < 0.0 ))
return Point< Coord >();
else
{
Coord point;
point.axpy( -l1 / ( l0 - l1 ), edge.vertex( 0 ) );
point.axpy( l0 / ( l0 - l1 ), edge.vertex( 1 ) );
return Point< Coord >( point );
}
return {{std::move(container)}};
}
} // namespace __impl
} // namespace Impl
} // namespace VoF
......
......@@ -3,6 +3,8 @@
#include <type_traits>
#include <dune/common/std/optional.hh>
//- local includes
#include <dune/vof/geometry/shapes/polygon.hh>
#include <dune/vof/geometry/utility.hh>
......@@ -19,15 +21,15 @@ namespace Dune
* \brief generate upwind polygon for 2d
*/
template< class Geometry >
static inline auto upwindPolygon2d ( const Geometry& geometry, const typename Geometry::GlobalCoordinate& v )
-> std::enable_if_t< Geometry::mydimension == 1 && Geometry::Globalcoordinate::dimension == 2,
Polygon< typename Geometry::GlobalCoordinate > >
template<class Geometry>
static inline auto upwindPolygon2d (const Geometry& geometry, const typename Geometry::GlobalCoordinate& v)
-> std::enable_if_t<Geometry::mydimension == 1 && Geometry::GlobalCoordinate::dimension == 2,
Std::optional<Polygon< typename Geometry::GlobalCoordinate>>>
{
if ( ( generalizedCrossProduct( geometry.corner( 1 ) - geometry.corner( 0 ) ) * v ) < 0.0 )
return {{ geometry.corner( 0 ), geometry.corner( 1 ), geometry.corner( 1 ) - v, geometry.corner( 0 ) - v }};
if ((v * generalizedCrossProduct(geometry.corner(1) - geometry.corner(0))) < 0.0)
return {{{ geometry.corner(0), geometry.corner(1), geometry.corner(1) - v, geometry.corner(0) - v }}};
else
return {{ geometry.corner( 1 ), geometry.corner( 0 ), geometry.corner( 0 ) - v, geometry.corner( 1 ) - v }};
return {{{ geometry.corner(1), geometry.corner(0), geometry.corner(0) - v, geometry.corner(1) - v }}};
}
} // namespace VoF
......
......@@ -24,7 +24,7 @@ namespace Dune {
/**
* \brief namespace containing specific implementations
*/
namespace __impl {
namespace Impl {
template< class Coord >
......@@ -272,7 +272,7 @@ namespace Dune {
* \return the intersection line
*/
template< class Coord >
auto intersect ( const Polyhedron< Coord >& polyhedron, const HyperPlane< Coord >& plane ) -> Dune::VoF::Face< Coord >
auto intersect ( const Polyhedron< Coord >& polyhedron, const Hyperplane< Coord >& plane ) -> Dune::VoF::Face< Coord >
{
if ( !plane )
return Dune::VoF::Face< Coord >();
......@@ -355,7 +355,7 @@ namespace Dune {
return typename Dune::VoF::Face< Coordinate >( nodes );
}
} // namespace __impl
} // namespace Impl
} // namespace VoF
......
......@@ -5,6 +5,8 @@
#include <type_traits>
#include <vector>
#include <dune/common/std/optional.hh>
//- local includes
#include <dune/vof/geometry/shapes/polyhedron.hh>
#include <dune/vof/geometry/utility.hh>
......@@ -21,28 +23,29 @@ namespace Dune
* \brief generate upwind polygon for 3d
*/
template< class Geometry >
inline auto upwindPolygon3d ( const Geometry& geometry, const typename Geometry::GlobalCoordinate& v )
-> std::enable_if_t< Geometry::mydimension == 2 && Geometry::GlobalCoordinate::dimension == 3,
Polyhedron< typename Geometry::GlobalCoordinate > >
template<class Geometry>
inline auto upwindPolygon3d (const Geometry& geometry, const typename Geometry::GlobalCoordinate& v)
-> std::enable_if_t<Geometry::mydimension == 2 && Geometry::GlobalCoordinate::dimension == 3,
Std::optional<Polyhedron< typename Geometry::GlobalCoordinate>>>
{
std::vector< typename Geometry::GlobalCoordinate > nodes;
std::vector<typename Geometry::GlobalCoordinate> nodes;
if ( generalizedCrossProduct( geometry.corner( 1 ) - geometry.corner( 0 ), geometry.corner( 2 ) - geometry.corner( 0 ) ) * v < 0.0 )
auto _range = range(geometry.corners());
if (v * generalizedCrossProduct( geometry.corner(1) - geometry.corner(0), geometry.corner(2) - geometry.corner(0)) < 0.0)
{
for ( int i = 0; i < geometry.corners(); ++i )
nodes.push_back( geometry.corner( i ) - v );
for (auto i : _range)
nodes.push_back(geometry.corner(i) - v);
for ( int i = 0; i < geometry.corners(); ++i )
nodes.push_back( geometry.corner( i ) );
for (auto i : _range)
nodes.push_back(geometry.corner(i));
}
else
{
for ( int i = 0; i < geometry.corners(); ++i )
nodes.push_back( geometry.corner( i ) );
for (auto i : _range)
nodes.push_back(geometry.corner(i));
for ( int i = 0; i < geometry.corners(); ++i )
nodes.push_back( geometry.corner( i ) - v );
for (auto i : _range)
nodes.push_back(geometry.corner(i) - v);
}
......@@ -57,7 +60,7 @@ namespace Dune
const std::vector< std::vector< std::size_t > > faces {
{{ 9, 3, 1, 15 }}, {{ 0, 7, 11, 13 }}, {{ 5, 2, 17, 10 }}, {{ 4, 14, 12 }}, {{ 6, 8, 16 }}
};
return { faces, edges, nodes };
return {{ faces, edges, nodes }};
}
else if( type.isQuadrilateral() )
{
......@@ -70,10 +73,10 @@ namespace Dune
const std::vector< std::vector< std::size_t > > faces {
{{ 0, 8, 14, 16 }}, {{ 5, 3, 21, 13 }}, {{ 6, 1, 22, 12 }}, {{ 19, 2, 11, 15 }}, {{ 4, 7, 17, 18 }}, {{ 10, 9, 23, 20 }}
};
return { faces, edges, nodes };
return {{ faces, edges, nodes }};
}
else
DUNE_THROW( InvalidStateException, "Invalid GeometryType." );
DUNE_THROW(InvalidStateException, "Invalid GeometryType.");
}
......
......@@ -7,6 +7,9 @@
#include <limits>
#include <vector>
#include <dune/common/typetraits.hh>
#include <dune/common/std/optional.hh>
#include <dune/vof/brents.hh>
#include <dune/vof/geometry/halfspace.hh>
#include <dune/vof/geometry/intersect.hh>
......@@ -26,20 +29,15 @@ namespace Dune {
*
* \tparam Coord global coordinate type
*/
template< class Coord >
double getVolumeFraction ( const Polygon< Coord > &polygon, const HalfSpace< Coord > &halfSpace )
{
Polygon< Coord > intersection = intersect( std::cref( polygon ), std::cref( halfSpace ) );
return intersection.volume() / polygon.volume();
}
template< class Coord >
double getVolumeFraction ( const Polyhedron< Coord > &polyhedron, const HalfSpace< Coord > &halfSpace )
template<template<class> class Shape, class Coordinate>
auto getVolumeFraction (const Std::optional<Shape<Coordinate>>& shape, const Std::optional<HalfSpace<Coordinate>>& halfSpace)
-> real_t<Coordinate>
{
Polyhedron< Coord > intersection = intersect( std::cref( polyhedron ), std::cref( halfSpace ) );
if (!shape)
return 0;
return intersection.volume() / polyhedron.volume();
auto intersection = intersect(shape, halfSpace);
return volume(intersection) / volume(shape);
}
......@@ -56,57 +54,66 @@ namespace Dune {
* \tparam Coord global coordinate type
*/
template< class Coord >
auto locateHalfSpace ( const Line< Coord >& line, const Coord& normal, double fraction ) -> HalfSpace< Coord >
template<class Coordinate>
auto locateHalfSpace (const Std::optional<Line<Coordinate>>& line, const Coordinate& normal, double fraction)
-> Std::optional<HalfSpace<Coordinate>>
{
return HalfSpace< Coord >( normal, line.centroid() + line.volume() * normal * ( 0.5 - fraction ) );
if (!line)
return {};
using std::abs;
assert(abs(normal.two_norm2() - real_t<Coordinate>{1.0}) < Coordinate::dimension * std::numeric_limits<real_t<Coordinate>>::epsilon());
return {{normal, volume(line)*(fraction-0.5) - (normal * centroid(line))}};
}
template< class Coord >
auto locateHalfSpace ( const Polygon< Coord >& polygon, const Coord& normal, double fraction ) -> HalfSpace< Coord >
template<class Coordinate>
auto locateHalfSpace (const Std::optional<Polygon<Coordinate>>& polygon, const Coordinate& normal, double fraction)
-> Std::optional<HalfSpace<Coordinate>>
{
using ctype = typename HalfSpace< Coord >::ctype;
using limits = std::numeric_limits< ctype >;
using ctype = real_t<Coordinate>;
using limits = std::numeric_limits<ctype>;
if (!polygon)
return {};
ctype pMin = 0, pMax = 0;
ctype volume,
volMin = limits::lowest(), // lower bound
ctype volMin = limits::lowest(), // lower bound
volMax = limits::max(); // upper bound
// upper/lower bound for
for( int i = 0; i < polygon.size(); ++i )
{
ctype dist = -( normal * polygon.vertex( i ) );
HalfSpace< Coord > hs( normal, dist );
auto objective = [&polygon, &normal, &fraction] (ctype p) -> ctype {
return (getVolumeFraction(polygon, {{normal, p}}) - fraction);
};
volume = getVolumeFraction( polygon, hs );
// bracketing
for (auto i : range(polygon->size()))
{
auto p = -(normal * polygon->vertex(i));
const auto vol = objective(p);
if( ( volume <= volMax ) && ( volume >= fraction ) )
if ((vol <= volMax) && (vol >= 0))
{
pMax = dist;
volMax = volume;
pMax = p;
volMax = vol;
}
if( ( volume >= volMin ) && ( volume <= fraction ) )
if((vol >= volMin) && (vol <= 0))
{
pMin = dist;
volMin = volume;
pMin = p;
volMin = vol;
}
}
assert( volMin <= volMax );
if ( volMax == limits::max() )
return HalfSpace< Coord >( normal, pMin );
else if ( volMin == limits::lowest() )
return HalfSpace< Coord >( normal, pMax );
else if ( volMin == volMax )
return HalfSpace< Coord >( normal, pMin );
if (volMax == limits::max())
return {{normal, pMin}};
else if (volMin == limits::lowest())
return {{normal, pMax}};
else if (volMin == volMax)
return {{normal, pMin}};
else
return HalfSpace< Coord >( normal, brentsMethod( [ &polygon, &normal, &fraction ] ( ctype p ) -> ctype {
HalfSpace< Coord > hs( normal, p );
return ( getVolumeFraction( polygon, hs ) - fraction );
}, pMin, pMax, 1e-12 ) );
return {{normal, brentsMethod(objective, pMin, pMax, 1e-12)}};
}
......
#ifndef DUNE_VOF_GEOMETRY_HALFSPACE_HH
#ifndef DUNE_VOF_GEOMETRY_HALFSPACE_HH
#define DUNE_VOF_GEOMETRY_HALFSPACE_HH
#include <cassert>
#include <cstddef>
#include <dune/common/deprecated.hh>
#include <limits>
#include <type_traits>
#include <vector>
#include <dune/common/std/optional.hh>
#include <dune/common/typetraits.hh>
#include <dune/vof/common/typetraits.hh>
#include <dune/vof/geometry/utility.hh>
namespace Dune
{
namespace VoF
......@@ -15,8 +23,8 @@ namespace Dune
// forward declarations
// --------------------
template< class > struct HalfSpace;
template< class > struct HyperPlane;
template<class> struct HalfSpace;
template<class> struct Hyperplane;
// HalfSpace
......@@ -28,91 +36,110 @@ namespace Dune
*
* \tparam Coord global coordinate type
*/
template < class Coord >
template<class Coord>
struct HalfSpace
{
using Coordinate = Coord;
using ctype = typename Coordinate::value_type;
using Boundary = HyperPlane< Coordinate >;
static constexpr int dimension = Coordinate::dimension;
static constexpr int dimensionworld = Coordinate::dimension;
HalfSpace ()
: innerNormal_( 0 ), distanceToOrigin_( 0.0 )
{}
using Coordinate = Coord;
using ctype = real_t<Coordinate>;
HalfSpace ( const Coordinate& normal, const ctype& dist )
: innerNormal_( normal ), distanceToOrigin_( dist )
{}
using Boundary = Hyperplane<Coordinate>;
HalfSpace ( const Boundary& plane )
: HalfSpace( plane.normal_, plane.distanceToOrigin_ )
{}
HalfSpace ( const Coordinate& normal, const Coordinate& point )
: innerNormal_( normal ), distanceToOrigin_( -1.0*( normal * point ) )
{}
static constexpr int dimension = Coordinate::dimension;
static constexpr int dimensionworld = Coordinate::dimension;
HalfSpace ( const std::vector< Coordinate > &points )
HalfSpace (const Coordinate& normal, const ctype& dist)
: innerNormal_(normal), distanceToOrigin_(dist)
{
assert( points.size() == dimension );
std::array< Coordinate, dimension-1 > vectors;
for ( std::size_t i = 1; i < dimension; ++i )
vectors[ i-1 ] = points[ i ] - points[ 0 ];
innerNormal_ = generalizedCrossProduct( vectors );
innerNormal_ /= innerNormal_.two_norm();
distanceToOrigin_ = -1.0*( innerNormal_ * points[ 0 ] );
using std::abs;
assert(abs(innerNormal_.two_norm()-ctype(1.0)) < dimension * std::numeric_limits<ctype>::epsilon());
}
HalfSpace (const Coordinate& normal, const Coordinate& point) : HalfSpace(normal, -(normal * point)) {}
HalfSpace (const Boundary& plane) : HalfSpace(plane.normal_, plane.distanceToOrigin_) {}
/**
* \brief inner normal (used in the representation)
*/
const Coordinate& innerNormal () const { return innerNormal_; }
auto innerNormal () const -> const Coordinate& { return innerNormal_; }
/**
* \brief bounding hyperplane
* \brief outer normal (returns a copy)
*/
Boundary boundary () const { return Boundary( innerNormal(), distanceToOrigin_ ); }
auto outerNormal () const -> Coordinate
{
Coordinate outer(innerNormal_);
return outer *= -1.0;
}
/**
* \brief level set function of the bounding plane
* \brief bounding hyperplane
*/
ctype levelSet ( const Coordinate& point ) const { return innerNormal() * point + distanceToOrigin_; }
auto boundary () const -> Boundary { return Boundary{innerNormal(), distanceToOrigin_}; }
/**
* \brief outer normal (returns a copy)
* \brief level set function of the bounding plane
*/
Coordinate outerNormal () const
{
Coordinate outer( innerNormal() );
outer *= -1.0;
return outer;
}
ctype levelSet ( const Coordinate& point ) const { return (innerNormal() * point) + distanceToOrigin_; }
/**
* \brief flip halfspace (returns a copy)
*/
HalfSpace< Coordinate > flip () const
auto flip () const -> HalfSpace<Coordinate> { return {outerNormal(), -distanceToOrigin_}; }
private:
Coordinate innerNormal_;
ctype distanceToOrigin_;
};
namespace Impl
{
template<class K, int N, class... Cs>
inline static auto _makeHalfSpace (const FieldVector<K, N>& v, const Cs& ...vs)
-> std::enable_if_t<sizeof...(Cs) == N-1 && isHomogeneous<FieldVector<K, N>, Cs...>::value,
Std::optional<HalfSpace<FieldVector<K, N>>>>
{
return HalfSpace< Coordinate > ( outerNormal(), - distanceToOrigin_ );
auto normal = generalizedCrossProduct((vs-v)...);
auto len = normal.two_norm2();
if (len < N * std::numeric_limits<real_t<K>>::epsilon())
return {};
return {{normal /= sqrt(len), v}};
}
explicit operator bool () const
template<class K, int N, std::size_t... Is>
inline static auto _makeHalfSpace (const std::vector<FieldVector<K, N>>& coords, std::index_sequence<Is...>)
-> std::enable_if_t<sizeof...(Is) == N, Std::optional<HalfSpace<FieldVector<K, N>>>>
{
return innerNormal_.two_norm2() > 0.5;
if (coords.size() == N)
return _makeHalfSpace(coords[Is]...);
else
return {};
}
private:
Coordinate innerNormal_;
ctype distanceToOrigin_;
};
}
// makeHalfSpace
// -------------
template<class... Args>
inline static auto makeHalfSpace (const Args& ...args)
-> decltype(Impl::_makeHalfSpace(args...))
{
return Impl::_makeHalfSpace(args...);
}
template<class K, int N>
inline static auto makeHalfSpace (const std::vector<FieldVector<K, N>>& coordinates)
-> decltype(Impl::_makeHalfSpace(coordinates, std::make_index_sequence<N>{}))
{
return Impl::_makeHalfSpace(coordinates, std::make_index_sequence<N>{});
}
// HyperPlane
// Hyperplane
// ----------
/**
......@@ -121,27 +148,21 @@ namespace Dune
*
* \tparam Coord global coordinate type
*/
template< class Coord >
struct HyperPlane
template<class Coord>
struct Hyperplane
{
using Coordinate = Coord;
using ctype = typename Coordinate::value_type;
using ctype = real_t<Coordinate>;
static constexpr int dimension = Coordinate::dimension - 1;
static constexpr int dimensionworld = Coordinate::dimension;
HyperPlane ()
: normal_( 0 ), distanceToOrigin_( 0.0 )
{}
HyperPlane ( const Coordinate& normal, const ctype& dist )
: normal_( normal ), distanceToOrigin_( dist )
{}
Hyperplane (const Coordinate& normal, const ctype& dist) : normal_(normal), distanceToOrigin_(dist) {}
/**
* \brief normal
*/
const Coordinate& normal () const { return normal_; }
auto normal () const -> const Coordinate& { return normal_; }
/**
* \brief distance to origin
......@@ -151,26 +172,69 @@ namespace Dune
/**
* \brief level set
*/
ctype levelSet ( const Coordinate& point ) const { return normal_ * point + distanceToOrigin_; }
ctype levelSet (const Coordinate& point) const { return (normal_ * point) + distanceToOrigin_; }
/**
* \brief cast into half space bounded by this plane
*/
operator HalfSpace< Coordinate > () const { return HalfSpace< Coordinate >( *this ); }
explicit operator bool () const
{
return normal_.two_norm2() > 0.5;
}
operator HalfSpace<Coordinate> () const { return HalfSpace<Coordinate>(*this); }
private:
friend HalfSpace< Coordinate >;
friend HalfSpace<Coordinate>;
Coordinate normal_;
ctype distanceToOrigin_;
};
// utility functions
// -----------------
template<class Coordinate>
auto boundary (const Std::optional<HalfSpace<Coordinate>>& H)
-> Std::optional<Hyperplane<Coordinate>>
{
if (!H)
return {};
return {H->boundary()};
}
template<class Coordinate>
auto complement (const Std::optional<HalfSpace<Coordinate>>& H)
-> Std::optional<HalfSpace<Coordinate>>
{
if (!H)
return {};
return {H->flip()};
}
template<class Coordinate>
auto levelSet (const Std::optional<HalfSpace<Coordinate>>& H, const Coordinate& x)
-> real_t<Coordinate>
{
if (!H)
return std::numeric_limits<real_t<Coordinate>>::infinity();
return H->levelSet(x);
}
template<class Coordinate>
auto levelSet (const Std::optional<Hyperplane<Coordinate>>& H, const Coordinate& x)
-> real_t<Coordinate>
{
if (!H)
return std::numeric_limits<real_t<Coordinate>>::infinity();
return H->levelSet(x);
}
} // namespace VoF
} // namespace Dune
......
#ifndef DUNE_VOF_GEOMETRY_INTERSECT_HH
#define DUNE_VOF_GEOMETRY_INTERSECT_HH
#include <cassert>
#include <cstddef>
#include <algorithm>
#include <functional>
#include <limits>
#include <type_traits>
#include <utility>
#include <vector>
#include <dune/common/exceptions.hh>
#include <dune/vof/geometry/halfspace.hh>
#include <dune/vof/geometry/1d/intersect.hh>
#include <dune/vof/geometry/2d/intersect.hh>
#include <dune/vof/geometry/3d/intersect.hh>
......@@ -22,67 +12,22 @@ namespace Dune {
namespace VoF {
/**
* \ingroup Geometry
* \brief expression template for intersection algorithms
*/