Commit 8b343c72 authored by Janick Gerstenberger's avatar Janick Gerstenberger

[geometry] modernize `intersect` implementations

use optional to denote empty sets
parent 4e0bdcf4
set(HEADERS
intersect.hh
upwindpolygon.hh
)
......
#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 {makeLine(P)};
}
}
} // namespace Impl
} // namespace VoF
} // namespace Dune
#endif // #ifndef DUNE_VOF_GEOMETRY_1D_INTERSECT_HH
......@@ -6,23 +6,27 @@
#include <algorithm>
#include <limits>
#include <type_traits>
#include <utility>
#include <vector>
#include <dune/vof/geometry/halfspace.hh>
#include <dune/vof/geometry/1d/point.hh>
#include <dune/vof/geometry/2d/polygon.hh>
#include <dune/common/rangeutilities.hh>
#include <dune/common/std/optional.hh>
#include <dune/common/typetraits.hh>
namespace Dune {
#include <dune/vof/geometry/halfspace.hh>
#include <dune/vof/geometry/shapes/line.hh>
#include <dune/vof/geometry/shapes/polygon.hh>
namespace VoF {
/**
* \brief namespace containing specific implementations
*/
namespace __impl {
namespace Dune
{
namespace VoF
{
namespace Impl
{
/**
* \ingroup geo2d
......@@ -31,26 +35,21 @@ 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;
......@@ -63,39 +62,41 @@ 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 {makePolygon(std::move(container))};
}
......@@ -106,109 +107,50 @@ 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 {makeLine(std::move(container))};
}
} // namespace __impl
} // namespace Impl
} // namespace VoF
......
This diff is collapsed.
......@@ -3,6 +3,7 @@
#include <type_traits>
#include <dune/vof/geometry/1d/intersect.hh>
#include <dune/vof/geometry/2d/intersect.hh>
#include <dune/vof/geometry/3d/intersect.hh>
......@@ -17,16 +18,16 @@ namespace Dune {
*/
template<class A, class B>
inline static auto intersect (const A& a, const B& b)
-> decltype(__impl::intersect(a, b))
-> decltype(Impl::intersect(a, b))
{
return __impl::intersect( a, b );
return Impl::intersect( a, b );
}
template<class A, class B>
inline static auto intersect (const A& a, const B& b)
-> std::enable_if_t<!std::is_same<A, B>::value, decltype(_impl::intersect(b, a))>
-> std::enable_if_t<!std::is_same<A, B>::value, decltype(Impl::intersect(b, a))>
{
return __impl::intersect(b, a);
return Impl::intersect(b, a);
}
} // namespace VoF
......
#ifndef DUNE_VOF_GEOMETRY_POLYTOPE_HH
#define DUNE_VOF_GEOMETRY_POLYTOPE_HH
#include <functional>
#include <type_traits>
#include <utility>
#include <dune/common/std/optional.hh>
#include <dune/vof/geometry/shapes.hh>
#include <dune/vof/geometry/1d/point.hh>
#include <dune/vof/geometry/2d/polygon.hh>
#include <dune/vof/geometry/3d/polyhedron.hh>
namespace Dune
{
......@@ -16,64 +17,43 @@ namespace Dune
// makePolytope
// ------------
template < class Geometry >
static inline auto makePolytope ( const Geometry& geometry )
-> typename std::enable_if< Geometry::GlobalCoordinate::dimension == 1,
Line< typename Geometry::GlobalCoordinate > >::type
{
return makeLine ( geometry );
}
template < class Geometry >
static inline auto makePolytope ( const Geometry& geometry )
-> typename std::enable_if< Geometry::GlobalCoordinate::dimension == 2 && Geometry::mydimension == 2,
Polygon< typename Geometry::GlobalCoordinate > >::type
template<class Geometry>
static inline auto makePolytope (const Geometry& geometry)
-> Std::optional<decltype(makeLine(geometry))>
{
return makePolygon ( geometry );
return {makeLine(geometry)};
}
template < class Geometry >
static inline auto makePolytope ( const Geometry& geometry )
-> typename std::enable_if< Geometry::GlobalCoordinate::dimension == 2 && Geometry::mydimension == 1,
Line< typename Geometry::GlobalCoordinate > >::type
template<class Geometry>
static inline auto makePolytope (const Geometry& geometry)
-> Std::optional<decltype(makePolygon(geometry))>
{
return makeLine ( geometry );
return {makePolygon(geometry)};
}
template < class Geometry, class Map >
static inline auto makePolytope ( const Geometry& geometry, Map&& map )
-> typename std::enable_if< Geometry::GlobalCoordinate::dimension == 2,
Polygon< typename Geometry::GlobalCoordinate > >::type
template<class Geometry, class Map>
static inline auto makePolytope (const Geometry& geometry, Map&& map)
-> Std::optional<decltype(makePolygon(geometry, std::forward<Map>(map)))>
{
return makePolygon ( geometry, std::forward< Map >( map ) );
return {makePolygon(geometry, std::forward<Map>(map))};
}
template < class Geometry >
static inline auto makePolytope ( const Geometry& geometry )
-> typename std::enable_if< Geometry::GlobalCoordinate::dimension == 3, Polyhedron< typename Geometry::GlobalCoordinate > >::type
template<class Geometry>
static inline auto makePolytope (const Geometry& geometry)
-> Std::optional<decltype(makePolyhedron(geometry))>
{
return makePolyhedron( geometry );
return {makePolyhedron(geometry)};
}
template < class Geometry, class Map >
static inline auto makePolytope ( const Geometry& geometry, Map&& map )
-> typename std::enable_if< Geometry::GlobalCoordinate::dimension == 3, Polyhedron< typename Geometry::GlobalCoordinate > >::type
template<class Geometry, class Map>
static inline auto makePolytope (const Geometry& geometry, Map&& map)
-> Std::optional<decltype(makePolyhedron(geometry, std::forward<Map>(map)))>
{
return makePolyhedron( geometry, std::forward< Map >( map ) );
return {makePolyhedron(geometry, std::forward<Map>(map))};
}
// return the interface
// --------------------
template< class Entity, class ReconstructionSet >
auto interface( const Entity &entity, const ReconstructionSet &reconstructions )
{
auto polygon = makePolytope( entity.geometry() );
auto it = intersect( std::cref( polygon ), reconstructions[ entity ].boundary() );
auto interface = static_cast< typename decltype( it )::Result > ( it );
return interface;
}
} // namespace VoF
} // namespace Dune
......
......@@ -2,11 +2,13 @@
//- C++ includes
#include <cmath>
#include <sstream>
#include <vector>
//- dune-common includes
#include <dune/common/exceptions.hh>
#include <dune/common/parametertreeparser.hh>
#include <dune/common/std/optional.hh>
//- dune-vof includes
#include <dune/vof/geometry/halfspace.hh>
......@@ -14,74 +16,173 @@
#include <dune/vof/geometry/polytope.hh>
int main(int argc, char** argv)
try {
namespace VoF = Dune::VoF;
namespace Std = Dune::Std;
Dune::MPIHelper::instance( argc, argv );
using limits = std::numeric_limits<double>;
using GridType = Dune::GridSelector::GridType;
using Coordinate = Dune::FieldVector< double, GridType::dimension >;
using Polytope = typename std::conditional< Coordinate::dimension == 2, Dune::VoF::Polygon< Coordinate >, Dune::VoF::Polyhedron< Coordinate > >::type;
using Grid = Dune::GridSelector::GridType;
using Coordinate = Dune::FieldVector<double, Grid::dimension>;
using HalfSpace = VoF::HalfSpace<Coordinate>;
// create grid
std::stringstream gridFile;
gridFile << GridType::dimension << "dgrid.dgf";
Dune::GridPtr< GridType > gridPtr( gridFile.str() );
gridPtr->loadBalance();
GridType& grid = *gridPtr;
inline static std::string dgfUnitCube ( int dimWorld )
{
std::string dgf = "DGF\nINTERVAL\n";
for( int i = 0; i < dimWorld; ++i )
dgf += " 0";
dgf += "\n";
for( int i = 0; i < dimWorld; ++i )
dgf += " 1";
dgf += "\n";
for( int i = 0; i < dimWorld; ++i )
dgf += " 1";
dgf += "\n#\n";
return dgf;
}
int level = 2;
const int refineStepsForHalf = Dune::DGFGridInfo< GridType >::refineStepsForHalf();
template<class Polytope>
auto printPolytopeInfo(const Polytope& p, std::ostream& out = std::cout)
-> std::enable_if_t<Polytope::dimension == 2>
{
out << "- Polygon\n"
<< "edges:\t" << p.size(1) << "\n"
<< "nodes:\t" << p.size(2) << "\n";
grid.globalRefine( refineStepsForHalf * level );
for (auto i : Dune::range(p.size(2)))
out << " [" << p.vertex(i) << "]";
out << "\n";
int numRuns = 1;
out << "centroid:\t[" << centroid(p) << "]\n"
<< "volume: \t" << volume(p) << "\n";
}
for ( int i = 0; i < numRuns; ++i )
template<class Polytope>
auto printPolytopeInfo(const Polytope& p, std::ostream& out = std::cout)
-> std::enable_if_t<Polytope::dimension == 3>
{
out << "- Polyhedron\n"
<< "faces:\t" << p.size(1) << "\n";
for (auto i : Dune::range(p.size(1)))
{
const auto& entity = *grid.leafGridView().begin< 0 >();
const auto& geoEn = entity.geometry();
const auto& polytope = Dune::VoF::makePolytope( geoEn );
out << " {";
for (const auto& j : p.face(i).edges())
out << j << ", ";
out << "}";
}
out << "\n";
Coordinate normal ( 0.0 );
normal[ 0 ] = -1.0;
out << "edges:\t" << p.size(2) << "\n";
std::cout << "Checking intersection with empty part..." << std::endl;
Dune::VoF::HalfSpace< Coordinate > hs1 ( normal, geoEn.corner(0) );
Polytope interface1 = Dune::VoF::intersect( polytope, hs1 );
assert( interface1.volume() < std::numeric_limits< double >::epsilon() );
for (auto i : Dune::range(p.size(2)))
{
out << " {";
for (const auto& j : p.edge(i).nodes())
out << j << ", ";
out << "}";
}
out << "\n";
std::cout << "Checking intersection with half part..." << std::endl;
Dune::VoF::HalfSpace< Coordinate > hs2 ( normal, geoEn.center() );
Polytope interface2 = Dune::VoF::intersect( polytope, hs2 );
assert( std::abs( interface2.volume() - 0.5 * geoEn.volume() ) < std::numeric_limits< double >::epsilon() );
out << "nodes:\t" << p.size(3) << "\n";
std::cout << "Checking intersection with full part..." << std::endl;
Dune::VoF::HalfSpace< Coordinate > hs3 ( normal, geoEn.corner(1) );
Polytope interface3 = Dune::VoF::intersect( polytope, hs3 );
assert( std::abs( interface3.volume() - 1.0 * geoEn.volume() ) < std::numeric_limits< double >::epsilon() );
for (auto i : Dune::range(p.size(3)))
out << " [" << p.node(i) << "]";
out << "\n";
out << "centroid:\t[" << centroid(p) << "]\n"
<< "volume: \t" << volume(p) << "\n";
}
int main(int argc, char** argv)
try {
Dune::MPIHelper::instance(argc, argv);
// create grid
std::istringstream dgf(dgfUnitCube(Grid::dimensionworld));
// read & construct grid
auto grid = Dune::GridPtr<Grid>(dgf);
grid