...
 
Commits (48)
......@@ -20,7 +20,6 @@ set(HEADERS
interpolation.hh
reconstruction.hh
scheme.hh
utility.hh
velocity.hh
)
......
set(HEADERS
commoperation.hh
iterator.hh
typetraits.hh
utility.hh
)
......
#ifndef DUNE_VOF_COMMON_ITERATOR_HH
#define DUNE_VOF_COMMON_ITERATOR_HH
#include <iterator>
#include <type_traits>
#include <dune/common/iteratorrange.hh>
#include <dune/vof/common/typetraits.hh>
namespace Dune
{
namespace VoF
{
// TransformIterator
// -----------------
template<class I, class F>
class TransformIterator
{
using This = TransformIterator<I, F>;
using Traits = std::iterator_traits<I>;
public:
using iterator_category = typename Traits::iterator_category;
using difference_type = typename Traits::difference_type;
using pointer = void;
using reference = InvokeResult_t<F, typename Traits::reference>;
using value_type = std::remove_reference_t<reference>;
static_assert(!std::is_same<iterator_category, std::output_iterator_tag>::value, "Invalid 'iterator_category'.");
constexpr TransformIterator (const I& i, F f) : i_(i), f_(f) {}
constexpr auto operator* () const -> reference { return f_(*i_); }
constexpr auto operator[] (difference_type n) const -> reference { return f_(i_[n]); }
auto operator++ () -> This& { ++i_; return *this; }
auto operator++ (int) -> This { This tmp = *this; ++i_; return tmp; }
auto operator-- () -> This& { --i_; return *this; }
auto operator-- (int) -> This { This tmp = *this; --i_; return tmp; }
auto operator+= (difference_type n) -> This& { i_ += n; return *this; }
auto operator-= (difference_type n) -> This& { i_ -= n; return *this; }
friend constexpr auto operator+ (const This& it, difference_type n) -> This { return {it.i_ + n, it.f_}; }
friend constexpr auto operator+ (difference_type n, const This& it) -> This { return {it.i_ + n, it.f_}; }
friend constexpr auto operator- (const This& it, difference_type n) -> This { return {it.i_ - n, it.f_}; }
constexpr auto operator- (const This& other) const -> difference_type { return i_ - other.i_; }
constexpr bool operator== (const This& other) const { return i_ == other.i_; }
constexpr bool operator!= (const This& other) const { return i_ != other.i_; }
constexpr bool operator< (const This& other) const { return i_ < other.i_; }
constexpr bool operator<= (const This& other) const { return i_ <= other.i_; }
constexpr bool operator> (const This& other) const { return i_ > other.i_; }
constexpr bool operator>= (const This& other) const { return i_ >= other.i_; }
private:
I i_;
F f_;
};
// transformIterator
// -----------------
template<class Iterator, class F>
static constexpr auto transformIterator(const Iterator& it, F f)
-> TransformIterator<Iterator, F>
{
return {it, f};
}
// transformRange
// --------------
template<class Iterator, class F>
static constexpr auto transformRange (const Iterator& begin, const Iterator& end, F f)
-> IteratorRange<TransformIterator<Iterator, F>>
{
return {transformIterator(begin, f), transformIterator(end, f)};
}
template<class Range, class F>
static constexpr auto transformRange (const Range& r, F f)
-> IteratorRange<TransformIterator<decltype(r.begin()), F>>
{
using std::begin; using std::end;
return transformRange(begin(r), end(r), f);
}
} // namespace VoF
} // namespace Dune
#endif // #ifndef DUNE_VOF_COMMON_ITERATOR_HH
......@@ -2,6 +2,7 @@
#define DUNE_VOF_COMMON_TYPETRAITS_HH
#include <type_traits>
#include <utility>
namespace Dune
......@@ -10,6 +11,56 @@ namespace Dune
namespace VoF
{
namespace Impl
{
template<class T>
struct InvokeImpl
{
template<class F, class... Args>
static auto call (F&& f, Args&& ...args) -> decltype(std::forward<F>(f)(std::forward<Args>(args)...));
};
template<class B, class MT>
struct InvokeImpl<MT B::*>
{
template<class T, class = std::enable_if_t<std::is_base_of<B, std::decay_t<T>>::value>>
static auto get (T&& t) -> T&&;
template<class T, class = std::enable_if_t<!std::is_base_of<B, std::decay_t<T>>::value>>
static auto get (T&& t) -> decltype(*std::forward<T>(t));
template<class T, class... Args, class MT1, class = std::enable_if_t<std::is_function<MT1>::value>>
static auto call (MT1 B::*pmf, T&& t, Args&& ...args) -> decltype((get(std::forward<T>(t)).*pmf)(std::forward<Args>(args)...));
template<class T>
static auto call (MT B::*pmd, T&& t) -> decltype(get(std::forward<T>(t)).*pmd);
};
template<class F, class... Args>
auto INVOKE (F&& f, Args&& ...args) -> decltype(InvokeImpl<std::decay_t<F>>::call(std::forward<F>(f), std::forward<Args>(args)...));
template<class SFINAE, class, class...>
struct InvokeResultImpl {};
template<class F, class... Args>
struct InvokeResultImpl<decltype(void(INVOKE(std::declval<F>(), std::declval<Args>()...))), F, Args...>
{
using type = decltype(INVOKE(std::declval<F>(), std::declval<Args>()...));
};
}
// InvokeResult
// ------------
template<class F, class... Args>
struct InvokeResult : Impl::InvokeResultImpl<void, F, Args...> {};
template<class F, class... Args>
using InvokeResult_t = typename InvokeResult<F, Args...>::type;
namespace Impl
{
......
set(HEADERS
algorithm.hh
intersect.hh
upwindpolygon.hh
)
......
#ifndef DUNE_VOF_GEOMETRY_1D_ALGORITHM_HH
#define DUNE_VOF_GEOMETRY_1D_ALGORITHM_HH
#include <cassert>
#include <cmath>
#include <limits>
#include <dune/common/typetraits.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/properties.hh>
namespace Dune
{
namespace VoF
{
// locateHalfSpace
// ---------------
/**
* \ingroup Geometry
* \brief locate half space with a given normal that intersects a line so it has a given volume fraction
*
* \param line line
* \param normal normal vector
* \param fraction volume fraction
* \tparam Coord global coordinate type
*/
template<class Coordinate>
auto locateHalfSpace (const Std::optional<Line<Coordinate>>& line, const Coordinate& normal, double fraction)
-> Std::optional<HalfSpace<Coordinate>>
{
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))}};
}
} // namespace VoF
} // namespace Dune
#endif // #ifndef DUNE_VOF_GEOMETRY_1D_ALGORITHM_HH
......@@ -76,7 +76,7 @@ namespace Dune
P.axpy(-l1 / (l0 - l1), L->vertex(0));
P.axpy( l0 / (l0 - l1), L->vertex(1));
return {{P}};
return {makeLine(P)};
}
}
......
set(HEADERS
algorithm.hh
intersect.hh
upwindpolygon.hh
)
......
#ifndef DUNE_VOF_GEOMETRY_2D_ALGORITHM_HH
#define DUNE_VOF_GEOMETRY_2D_ALGORITHM_HH
#include <cassert>
#include <cmath>
#include <limits>
#include <dune/common/rangeutilities.hh>
#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>
#include <dune/vof/geometry/shapes/polygon.hh>
#include <dune/vof/geometry/shapes/properties.hh>
namespace Dune
{
namespace VoF
{
// locateHalfSpace
// ---------------
/**
* \ingroup Geometry
* \brief locate half space with a given normal that intersects a polygon so it has a given volume fraction
*
* \param polygon polygon
* \param normal normal vector
* \param fraction volume fraction
* \tparam Coord global coordinate type
*/
template<class Coordinate>
auto locateHalfSpace (const Std::optional<Polygon<Coordinate>>& polygon, const Coordinate& normal, double fraction)
-> Std::optional<HalfSpace<Coordinate>>
{
using ctype = real_t<Coordinate>;
using limits = std::numeric_limits<ctype>;
if (!polygon)
return {};
ctype pMin = 0;
ctype pMax = 0;
ctype volMin = limits::lowest();
ctype volMax = limits::max();
auto objective = [&polygon, &normal, &fraction] (ctype p) -> ctype {
return volume(intersect(polygon, Std::optional<HalfSpace<Coordinate>>{{normal, p}}))/volume(polygon) - fraction;
};
// bracketing
for (auto i : range(polygon->size()))
{
auto p = -(normal * polygon->vertex(i));
const auto vol = objective(p);
if ((vol <= volMax) && (vol >= 0))
{
pMax = p;
volMax = vol;
}
if((vol >= volMin) && (vol <= 0))
{
pMin = p;
volMin = vol;
}
}
assert( volMin <= volMax );
if (volMax == limits::max())
return {{normal, pMin}};
else if (volMin == limits::lowest())
return {{normal, pMax}};
else if (volMin == volMax)
return {{normal, pMin}};
else
return {{normal, brentsMethod(objective, pMin, pMax, 1e-12)}};
}
} // namespace VoF
} // namespace Dune
#endif // #ifndef DUNE_VOF_GEOMETRY_2D_ALGORITHM_HH
......@@ -14,7 +14,6 @@
#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/polygon.hh>
......@@ -37,7 +36,8 @@ namespace Dune
* \return the intersection polygon
*/
template<class Coordinate>
auto intersect (const Std::optional<Polygon<Coordinate>>& P, const Std::optional<Polygon<Coordinate>>& Q)
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 (!P || !Q)
......@@ -63,7 +63,8 @@ namespace Dune
* \return the intersection polygon
*/
template<class Coordinate>
auto intersect (const Std::optional<Polygon<Coordinate>>& P, const Std::optional<HalfSpace<Coordinate>>& H)
auto intersect (const Std::optional<Polygon<Coordinate>>& P,
const Std::optional<HalfSpace<Coordinate>>& H)
-> Std::optional<Polygon<Coordinate>>
{
if (!P || !H)
......@@ -95,7 +96,7 @@ namespace Dune
if (container.size() < 3)
return {};
else
return {{std::move(container)}};
return {makePolygon(std::move(container))};
}
......@@ -107,7 +108,8 @@ namespace Dune
* \return the intersection line
*/
template<class Coordinate>
auto intersect (const Std::optional<Polygon<Coordinate>>& P, const Std::optional<Hyperplane<Coordinate>>& H)
auto intersect (const Std::optional<Polygon<Coordinate>>& P,
const Std::optional<Hyperplane<Coordinate>>& H)
-> Std::optional<Line<Coordinate>>
{
using limits = std::numeric_limits<real_t<Coordinate>>;
......@@ -145,7 +147,7 @@ namespace Dune
else if (j == 1)
container[1] = container[0];
return {{std::move(container)}};
return {makeLine(std::move(container))};
}
} // namespace Impl
......
......@@ -27,9 +27,9 @@ namespace Dune
Std::optional<Polygon< typename Geometry::GlobalCoordinate>>>
{
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 }}};
return {makePolygon(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 {makePolygon(geometry.corner(1), geometry.corner(0), geometry.corner(0) - v, geometry.corner(1) - v)};
}
} // namespace VoF
......
set(HEADERS
face.hh
algorithm.hh
intersect.hh
polygonwithdirections.hh
referenceframe.hh
rotation.hh
upwindpolygon.hh
)
......
This diff is collapsed.
#ifndef DUNE_VOF_GEOMETRY_3D_FACE_HH
#define DUNE_VOF_GEOMETRY_3D_FACE_HH
/* c++ includes */
#include <cmath>
#include <vector>
namespace Dune {
namespace VoF {
template < class Coord >
struct Face
{
Face () {};
Face ( const std::vector< Coord >& nodes ) : nodes_ ( nodes )
{
assert( nodes.size() > 2 );
faceUnitNormal_ = generalizedCrossProduct( vertex( 2 ) - vertex( 1 ), vertex( 0 ) - vertex( 1 ) );
faceUnitNormal_ /= faceUnitNormal_.two_norm();
}
operator bool () const { return !nodes_.empty(); }
bool operator== ( const Face< Coord > &other) const
{
if ( other.size() != this->size() )
return false;
std::size_t i0 = 0;
for ( std::size_t i = 0; i < other.size(); ++i )
if ( other.vertex( i ) == this->vertex( 0 ) )
{
i0 = i;
break;
}
for ( std::size_t i = 0; i < other.size(); ++i )
if ( other.vertex( (i + i0)%other.size() ) != this->vertex( i ) )
return false;
return true;
}
std::size_t size () const { return nodes_.size(); }
const Coord& vertex ( const std::size_t index ) const { return nodes_[ index % size() ]; }
double triangleVolume( const std::vector< Coord > &n ) const
{
assert ( n.size() == 3 );
double sum = 0;
for( std::size_t i = 0; i < 3; ++i )
{
Coord center = n[ (i+1)%3 ] ;
center += n[ i ];
center *= 0.5;
Coord edge = n[ (i+1)%3 ];
edge -= n[ i ];
sum += center * generalizedCrossProduct( edge, faceUnitNormal_ );
}
return std::abs( sum ) / 2.0;
}
Coord centroid () const
{
Coord c ( 0.0 );
for ( std::size_t i = 0; i < size(); ++i )
c += nodes_[ i ];
c /= static_cast< double >( size() );
Coord centroid ( 0.0 );
for( std::size_t i = 0; i < size(); ++i )
{
double volTriang = triangleVolume( { vertex( i ), vertex( i+1 ), c } );
Coord cTriang = c;
cTriang += vertex( i );
cTriang += vertex( i+1 );
cTriang /= 3.0;
centroid.axpy( volTriang, cTriang );
}
centroid /= volume();
return centroid;
}
double volume () const
{
double sum = 0;
for( std::size_t i = 0; i < size(); ++i )
sum += edgeCenter( i ) * edgeOuterNormal( i );
return std::abs( sum ) / 2.0;
}
private:
const Coord edgeOuterNormal( const std::size_t i ) const
{
Coord edge = vertex( i+1 );
edge -= vertex( i );
return generalizedCrossProduct( edge, faceUnitNormal_ );
}
const Coord edgeCenter( const std::size_t i ) const
{
Coord center = vertex( i+1 );
center += vertex( i );
center *= 0.5;
return center;
}
const std::vector< Coord > nodes_;
Coord faceUnitNormal_;
};
} // namespace VoF
} // namespace Dune
#endif // DUNE_VOF_GEOMETRY_3D_FACE_HH
This diff is collapsed.
This diff is collapsed.
#ifndef DUNE_VOF_GEOMETRY_3D_REFERENCEFRAME_HH
#define DUNE_VOF_GEOMETRY_3D_REFERENCEFRAME_HH
#include <cassert>
#include <limits>
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/common/typetraits.hh>
namespace Dune
{
namespace VoF
{
template<class Coord>
struct ReferenceFrame
{
using Coordinate = Coord;
using field = field_t<Coordinate>;
using real = real_t<Coordinate>;
static constexpr int dimension = 3;
using ProjectionVector = FieldVector<field, 2>;
using RejectionVector = FieldVector<field, 1>;
using Rotation = FieldMatrix<field, 3, 3>;
static_assert(dimension == 3, "...");
static constexpr auto up () -> Coordinate { return {0, 0, 1}; }
static constexpr auto project (const Coordinate& v) -> ProjectionVector { return {v[0], v[1]}; }
static constexpr auto reject (const Coordinate& v) -> RejectionVector { return {v[2]}; }
/**
* \brief return the (proper) rotation matrix that rotates v onto up.
*
* \param v vector
*/
static constexpr auto rotation (const Coordinate& v) -> Rotation
{
using limits = std::numeric_limits<real>;
using std::abs;
assert(abs(v.two_norm2() - real{1.0}) < 3*dimension * limits::epsilon());
if ((v - up()).two_norm2() < 3*dimension * limits::epsilon())
return {{1, 0, 0},
{0, 1, 0},
{0, 0, 1}};
else if ((v + up()).two_norm2() < 3*dimension * limits::epsilon())
return {{ 1, 0, 0},
{ 0, -1, 0},
{ 0, 0, -1}};
real g = (1 - v[2]) / (v[0]*v[0] + v[1]*v[1]);
return {{v[2] + v[1]*v[1]*g, -v[0]*v[1]*g, -v[0]},
{ -v[0]*v[1]*g, v[2] + v[0]*v[0]*g, -v[1]},
{ v[0], v[1], v[2]}};
}
};
} // namespace VoF
} // namespace Dune
#endif // #ifndef DUNE_VOF_GEOMETRY_3D_REFERENCEFRAME_HH
#ifndef DUNE_VOF_GEOMETRY_3D_ROTATION_HH
#define DUNE_VOF_GEOMETRY_3D_ROTATION_HH
#include <vector>
/* dune includes */
#include <dune/common/fmatrix.hh>
/* local includes */
#include <dune/vof/geometry/shapes/polyhedron.hh>
namespace Dune
{
namespace VoF
{
template < class Coord >
const Polyhedron< Coord > rotateToReferenceFrame ( const Coord& n, const Polyhedron< Coord >& cell )
{
// If n == e_z, there is nothing to do.
if ( n == Coord ( { 0.0, 0.0, 1.0 } ) )
return cell;
Dune::FieldMatrix< double, 3, 3 > rotationMatrix;
// If n == - e_z, rotate of angle pi around x-axis.
if ( n == Coord ( { 0.0, 0.0, -1.0 } ) )
{
rotationMatrix[0][0] = 1.0;
rotationMatrix[1][1] = -1.0;
rotationMatrix[2][2] = -1.0;
}
else
{
double gamma = ( 1 - n[2] ) / ( n[0] * n[0] + n[1] * n[1] );
rotationMatrix[0][0] = n[1] * n[1] * gamma + n[2];
rotationMatrix[0][1] = - n[0] * n[1] * gamma;
rotationMatrix[1][0] = rotationMatrix[0][1];
rotationMatrix[0][2] = - n[0];
rotationMatrix[1][1] = n[0] * n[0] * gamma + n[2];
rotationMatrix[1][2] = - n[1];
rotationMatrix[2][0] = n[0];
rotationMatrix[2][1] = n[1];
rotationMatrix[2][2] = n[2];
}
std::vector< Coord > newNodes ( cell.nodes().size() );
for ( std::size_t i = 0; i < cell.nodes().size(); ++i )
rotationMatrix.mv( cell.node( i ), newNodes[ i ] );
return Polyhedron< Coord > ( cell, newNodes );
}
} // namespace VoF
} // namespace Dune
#endif // #ifndef DUNE_VOF_GEOMETRY_3D_ROTATION_HH
#ifndef DUNE_VOF_GEOMETRY_ALGORITHM_HH
#define DUNE_VOF_GEOMETRY_ALGORITHM_HH
#include <cassert>
#include <functional>
#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>
#include <dune/vof/geometry/shapes.hh>
#include <dune/vof/geometry/3d/polygonwithdirections.hh>
#include <dune/vof/geometry/3d/rotation.hh>
namespace Dune {
namespace VoF {
/**
* \ingroup Geometry
* \brief calculate fraction of the volume of the intersection of a half space with a polygon
*
* \tparam Coord global coordinate type
*/
template<template<class> class Shape, class Coordinate>
auto getVolumeFraction (const Std::optional<Shape<Coordinate>>& shape, const Std::optional<HalfSpace<Coordinate>>& halfSpace)
-> real_t<Coordinate>
{
if (!shape)
return 0;
auto intersection = intersect(shape, halfSpace);
return volume(intersection) / volume(shape);
}
// locateHalfSpace
// ---------------
/**
* \ingroup Geometry
* \brief locate half space with a given normal that intersects a polygon so it has a given volume fraction
*
* \param polygon polygon
* \param normal normal vector
* \param fraction volume fraction
* \tparam Coord global coordinate type
*/
template<class Coordinate>
auto locateHalfSpace (const Std::optional<Line<Coordinate>>& line, const Coordinate& normal, double fraction)
-> Std::optional<HalfSpace<Coordinate>>
{
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 Coordinate>
auto locateHalfSpace (const Std::optional<Polygon<Coordinate>>& polygon, const Coordinate& normal, double fraction)
-> Std::optional<HalfSpace<Coordinate>>
{
using ctype = real_t<Coordinate>;
using limits = std::numeric_limits<ctype>;
if (!polygon)
return {};
ctype pMin = 0, pMax = 0;
ctype volMin = limits::lowest(), // lower bound
volMax = limits::max(); // upper bound
auto objective = [&polygon, &normal, &fraction] (ctype p) -> ctype {
return (getVolumeFraction(polygon, {{normal, p}}) - fraction);
};
// bracketing
for (auto i : range(polygon->size()))
{
auto p = -(normal * polygon->vertex(i));
const auto vol = objective(p);
if ((vol <= volMax) && (vol >= 0))
{
pMax = p;
volMax = vol;
}
if((vol >= volMin) && (vol <= 0))
{
pMin = p;
volMin = vol;
}
}
assert( volMin <= volMax );
if (volMax == limits::max())
return {{normal, pMin}};
else if (volMin == limits::lowest())
return {{normal, pMax}};
else if (volMin == volMax)
return {{normal, pMin}};
else
return {{normal, brentsMethod(objective, pMin, pMax, 1e-12)}};
}
template< class Coord >
auto locateHalfSpace ( const Polyhedron< Coord >& cell, const Coord& innerNormal, double fraction ) -> HalfSpace< Coord >
{
using std::abs;
assert( abs( innerNormal.two_norm() - static_cast< double > ( 1.0 ) ) < 1e-12 );
Coord outerNormal ( innerNormal );
outerNormal *= -1.0;
fraction = ( fraction > 1.0 ) ? 1.0 : fraction;
fraction = ( fraction < 0.0 ) ? 0.0 : fraction;
bool inverseMode = false;
/*
if ( fraction > 0.6 )
{
fraction = 1.0 - fraction;
outerNormal *= -1.0;
inverseMode = true;
}
*/
double givenVolume = fraction * cell.volume();
Polyhedron< Coord > polyhedron ( rotateToReferenceFrame ( outerNormal, cell ) );
// Compute and sort plane constants
// --------------------------------
const std::size_t N = polyhedron.nodes().size();
std::vector< double > d ( N );
std::vector< std::size_t > sortedNodesIds ( N );
for ( std::size_t i = 0; i < N; ++i )
{
d[ i ] = polyhedron.node( i )[ 2 ];
sortedNodesIds[ i ] = i;
}
std::sort( sortedNodesIds.begin(), sortedNodesIds.end(), [ &d ]( const std::size_t i, const std::size_t j ) { return d[ i ] < d[ j ]; } );
std::sort( d.begin(), d.end() );
std::vector< double > dUnique ( d );
auto it = std::unique( dUnique.begin(), dUnique.end(), []( const double x, const double y ) { return std::abs( x - y ) < 1e-10; } );
dUnique.resize( std::distance( dUnique.begin(), it ) );
PolygonWithDirections< Polyhedron< Coord > > polygon ( polyhedron );
polygon.initialize( sortedNodesIds, d );
// Bracketing
// ----------
double Vd_k = 0, Vd_k1;
for ( std::size_t k = 0; k < dUnique.size() - 1; ++k )
{
// Compute cooefficients
const std::size_t Nn = polygon.nodes().size();
double A = 0.0;
for ( std::size_t n = 0; n < Nn; ++n )
{
assert ( !polygon.directions( n ).empty() );
const std::size_t epsn = polygon.directions( n ).size() - 1;
Coord v1 = polygon.directionVector( n, epsn );
Coord v2 = polygon.directionVector( (n+1) % Nn, 0 );
A -= ( v1[0] * v2[1] - v1[1] * v2[0] ) / ( v1[2] * v2[2] );
for ( std::size_t l = 0; l < epsn; ++l )
{
Coord w1 = polygon.directionVector( n, l );
Coord w2 = polygon.directionVector( n, l+1 );
A -= ( w1[0] * w2[1] - w1[1] * w2[0] ) / ( w1[2] * w2[2] );
}
}
assert ( A == A );
double B = 0.0;
for ( std::size_t i = 0; i < polygon.edges().size(); ++i )
{
Coord normal = polygon.correspondingFace( i ).outerNormal();
double norm = project( normal ).two_norm();
if ( norm > 0 )
B -= polygon.edge( i ).volume() * normal[2] / norm;
}
assert ( B == B );
double C = polygon.volume();
assert ( C == C );
double h = dUnique[ k+1 ] - dUnique[ k ];
auto V = [ A, B, C ] ( const double h ) { return C * h + B * h * h / 2.0 + A * h * h * h / 6.0; };
Vd_k1 = Vd_k + V( h );
if ( std::abs( Vd_k1 - givenVolume ) < 1e-14 )
{
double distance = inverseMode ? -dUnique[ k+1 ] : dUnique[ k+1 ];
return HalfSpace< Coord > ( innerNormal, distance );
}
else if ( Vd_k1 > givenVolume )
{
const auto& Vh_V = [ &V, givenVolume, Vd_k ] ( const double h ) { return V( h ) - ( givenVolume - Vd_k ); };
double distance = dUnique[ k ] + Dune::VoF::brentsMethod ( Vh_V, 0.0, h, 1e-14 );
distance = inverseMode ? -distance : distance;
return HalfSpace< Coord > ( innerNormal, distance );
}
else
polygon.evolveToNextPolygon( k, h, dUnique );
Vd_k = Vd_k1;
}
return HalfSpace< Coord > ( innerNormal, inverseMode ? -dUnique[ dUnique.size() - 1 ] : dUnique[ dUnique.size() - 1 ] );
}
} // namespace VoF
} // namespace Dune
#endif // #ifndef DUNE_VOF_GEOMETRY_ALGORITHM_HH
#include <dune/vof/geometry/1d/algorithm.hh>
#include <dune/vof/geometry/2d/algorithm.hh>
#include <dune/vof/geometry/3d/algorithm.hh>
......@@ -51,7 +51,7 @@ namespace Dune
: innerNormal_(normal), distanceToOrigin_(dist)
{
using std::abs;
assert(abs(innerNormal_.two_norm()-ctype(1.0)) < dimension * std::numeric_limits<ctype>::epsilon());
assert(abs(innerNormal_.two_norm()-ctype{1.0}) < 3*dimension * std::numeric_limits<ctype>::epsilon());
}
HalfSpace (const Coordinate& normal, const Coordinate& point) : HalfSpace(normal, -(normal * point)) {}
......
......@@ -2,4 +2,4 @@
#include <dune/vof/geometry/shapes/point.hh>
#include <dune/vof/geometry/shapes/polygon.hh>
#include <dune/vof/geometry/shapes/polyhedron.hh>
#include <dune/vof/geometry/shapes/utility.hh>
#include <dune/vof/geometry/shapes/properties.hh>
......@@ -3,7 +3,7 @@ set(HEADERS
point.hh
polygon.hh
polyhrdron.hh
utility.hh
properties.hh
)
install(FILES ${HEADERS} DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/vof/geometry/shapes)
......@@ -34,8 +34,8 @@ namespace Dune
static constexpr int dimensionworld = Coordinate::dimension;
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]}} { assert(v.size() == 2); }
explicit Line (Container vertices) : vertices_(std::move(vertices)) {}
explicit 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_; }
......@@ -47,9 +47,17 @@ namespace Dune
};
// makeLine
// --------
template<class Coordinate>
static inline auto makeLine (std::array<Coordinate, 2> vertices)
-> Line<Coordinate>
{
return Line<Coordinate>{std::move(vertices)};
}
template<class Geometry>
static inline auto makeLine(const Geometry& geometry)
-> std::enable_if_t<Geometry::mydimension == 1, Line<typename Geometry::GlobalCoordinate>>
......
......@@ -27,7 +27,7 @@ namespace Dune {
static constexpr int dimension = 0;
static constexpr int dimensionworld = Coordinate::dimension;
Point (Coordinate vertex) : vertex_(std::move(vertex)) {}
explicit Point (Coordinate vertex) : vertex_(std::move(vertex)) {}
bool operator== (const This &other) const { return vertex_ == other.vertex(); }
......@@ -46,6 +46,19 @@ namespace Dune {
Coordinate vertex_;
};
// makePoint
// ---------
template<class Coordinate>
static inline auto makePoint (Coordinate vertex)
-> Point<Coordinate>
{
return Point<Coordinate>{std::move(vertex)};
}
} // namespace VoF
} // namespace Dune
......
......@@ -11,6 +11,7 @@
#include <dune/common/exceptions.hh>
#include <dune/common/rangeutilities.hh>
#include <dune/vof/common/typetraits.hh>
#include <dune/vof/geometry/shapes/line.hh>
......@@ -39,7 +40,7 @@ namespace Dune {
static constexpr int dimensionworld = Coordinate::dimension;
static_assert(dimension <= dimensionworld, "dimension larger than dimensionworld.");
Polygon (Container vertices) : vertices_(std::move(vertices)) { assert(vertices_.size() >= 3); }
explicit Polygon (Container vertices) : vertices_(std::move(vertices)) { assert(vertices_.size() >= 3); }
int size (int codim = 2) const { assert( codim <= 2 ); return (codim == 0) ? 1 : vertices().size(); }
......@@ -64,8 +65,8 @@ namespace Dune {
/**
* \brief i-th vertex
*/
auto vertex (int i) -> Coordinate& { assert(i < size()); return vertices_[ i ]; }
auto vertex (int i) const -> const Coordinate& { assert(i < size()); return vertices()[i]; }
auto vertex (int i) -> Coordinate& { assert(i < size()); return vertices_[i]; }
auto vertex (int i) const -> const Coordinate& { return vertices()[i]; }
/**
* \brief i-th edge
......@@ -79,6 +80,7 @@ namespace Dune {
};
/**
* \ingroup geo2d
* \brief generate polygon
......@@ -87,10 +89,17 @@ namespace Dune {
* \tparam Coord global coordinate type
*/
template<class Coordinate>
static inline auto makePolygon(std::vector<Coordinate> vertices)
static inline auto makePolygon (std::vector<Coordinate> vertices)
-> Polygon<Coordinate>
{
return Polygon<Coordinate>(std::move(vertices));
return Polygon<Coordinate>{std::move(vertices)};
}
template<class... Coords>
static inline auto makePolygon (const Coords& ...coords)
-> std::enable_if_t<sizeof...(Coords) >= 3, Polygon<isHomogeneous_t<Coords...>>>
{
return makePolygon(std::vector<isHomogeneous_t<Coords...>>{coords...});
}
/**
......@@ -100,31 +109,29 @@ namespace Dune {
* \param geometry dune geometry
*/
template<class Geometry>
static inline auto makePolygon(const Geometry& geometry)
static inline auto makePolygon (const Geometry& geometry)
-> std::enable_if_t<Geometry::mydimension == 2, Polygon<typename Geometry::GlobalCoordinate>>
{
using Container = std::vector<typename Geometry::GlobalCoordinate>;
auto type = geometry.type();
if(type.isTriangle())
return makePolygon(Container{geometry.corner(0), geometry.corner(1), geometry.corner(2)});
return makePolygon(geometry.corner(0), geometry.corner(1), geometry.corner(2));
else if (type.isQuadrilateral())
return makePolygon(Container{geometry.corner(0), geometry.corner(1), geometry.corner(3), geometry.corner(2)});
return makePolygon(geometry.corner(0), geometry.corner(1), geometry.corner(3), geometry.corner(2));
else
DUNE_THROW(InvalidStateException, "Invalid GeometryType.");
}
template<class Geometry, class Map>
static inline auto makePolygon(const Geometry& geometry, Map&& map)
static inline auto makePolygon (const Geometry& geometry, Map&& map)
-> std::enable_if_t<Geometry::mydimension == 2, Polygon<typename Geometry::GlobalCoordinate>>
{
using Container = std::vector<typename Geometry::GlobalCoordinate>;
auto type = geometry.type();
if(type.isTriangle())
return makePolygon(Container{map(geometry.corner(0)), map(geometry.corner(1)), map(geometry.corner(2))});
return makePolygon(map(geometry.corner(0)), map(geometry.corner(1)), map(geometry.corner(2)));
else if (type.isQuadrilateral())
return makePolygon(Container{map(geometry.corner(0)), map(geometry.corner(1)), map(geometry.corner(3)), map(geometry.corner(2))});
return makePolygon(map(geometry.corner(0)), map(geometry.corner(1)), map(geometry.corner(3)), map(geometry.corner(2)));
else
DUNE_THROW(InvalidStateException, "Invalid GeometryType.");
}
......
This diff is collapsed.
#ifndef DUNE_VOF_GEOMETRY_SHAPES_UTILITY_HH
#define DUNE_VOF_GEOMETRY_SHAPES_UTILITY_HH
#ifndef DUNE_VOF_GEOMETRY_SHAPES_PROPERTIES_HH
#define DUNE_VOF_GEOMETRY_SHAPES_PROPERTIES_HH
#include <cassert>
#include <type_traits>
......@@ -22,13 +22,16 @@ namespace Dune
// - implement centroid of Polygon in 3D
// forward declarations
// --------------------
template<class> class Point;
template<class> class Line;
template<class> class Polygon;
template<class> struct Polyhedron;
template<class> class Polyhedron;
// volume(...)
// -----------
......@@ -40,6 +43,7 @@ namespace Dune
return shape ? volume(*shape) : 0;
}
template<class Coordinate>
auto volume (DUNE_UNUSED const Point<Coordinate>& point)
-> real_t<Coordinate>
......@@ -47,6 +51,7 @@ namespace Dune
return 1;
}
template<class Coordinate>
auto volume (const Line<Coordinate>& line)
-> real_t<Coordinate>
......@@ -54,6 +59,7 @@ namespace Dune
return (line.vertex(0)-line.vertex(1)).two_norm();
}
template<class Coordinate>
auto volume (const Polygon<Coordinate>& polygon)
-> std::enable_if_t<Coordinate::dimension == 2, real_t<Coordinate>>
......@@ -65,6 +71,7 @@ namespace Dune
return 0.5 * vol;
}
template<class Coordinate>
auto volume (const Polygon<Coordinate>& polygon)
-> std::enable_if_t<Coordinate::dimension == 3, real_t<Coordinate>>
......@@ -80,13 +87,19 @@ namespace Dune
return 0.5 * vol;
}
template<class Coordinate>
auto volume (const Polyhedron<Coordinate>& polyhedron)
-> real_t<Coordinate>
{
return polyhedron.volume();
real_t<Coordinate> vol = 0;
for (const auto& f : polyhedron.faces())
vol += volume(f) * (centroid(f) * outerNormal(f));
return vol /= 3;
}
// centroid(...)
// -------------
......@@ -130,18 +143,38 @@ namespace Dune
-> std::enable_if_t<Coordinate::dimension == 3, Coordinate>
{
const int n = polygon.size();
const auto& o = polygon.vertex(0);
auto nu = generalizedCrossProduct(polygon.vertex(1) - o, polygon.vertex(n-1) - o);
normalize(nu);
auto planeNormal = generalizedCrossProduct(polygon.vertex(1) - polygon.vertex(0), polygon.vertex(n-1) - polygon.vertex(0));
normalize(planeNormal);
return {};
Coordinate c;
real_t<Coordinate> vol = 0;
for (auto i : range(1, n-1))
{
auto voli = (generalizedCrossProduct(o, polygon.vertex(i))
+ generalizedCrossProduct(polygon.vertex(i), polygon.vertex(i+1))
+ generalizedCrossProduct(polygon.vertex(i+1), o)) * nu;
c.axpy(voli, o + polygon.vertex(i) + polygon.vertex(i+1));
vol += voli;
}
return c /= 3*vol;
}
template<class Coordinate>
inline static auto centroid (const Polyhedron<Coordinate>& polyhedron)
-> Coordinate
{
return polyhedron.center();
Coordinate c;
real_t<Coordinate> vol = 0;
for (const auto& f : polyhedron.faces())
{
auto ci = centroid(f);
auto voli = volume(f) * (ci * outerNormal(f));
c.axpy(voli, ci);
vol += voli;
}
return c *= 0.75/vol;
}
......@@ -150,4 +183,4 @@ namespace Dune
} // Dune
#endif // #ifndef DUNE_VOF_GEOMETRY_SHAPES_UTILITY_HH
#endif // #ifndef DUNE_VOF_GEOMETRY_SHAPES_PROPERTIES_HH
set(HEADERS
capabilities.hh
datahandle.hh
dataset.hh
declaration.hh
entity.hh
entityseed.hh
......
......@@ -15,6 +15,8 @@
#include <dune/vof/dataset/flag.hh>
#include <dune/vof/dataset/reconstruction.hh>
#include <dune/vof/flagging.hh>
#include <dune/vof/geometry/intersect.hh>
#include <dune/vof/geometry/polytope.hh>
#include <dune/vof/interfacegrid/geometry.hh>
#include <dune/vof/interfacegrid/mixedcellmapper.hh>
......@@ -41,7 +43,7 @@ namespace Dune
typedef typename GridView::template Codim< 0 >::Entity Element;
typedef typename ReconstructionSet::DataType::Coordinate GlobalCoordinate;
typedef typename ReconstructionSet::DataType::value_type::Coordinate GlobalCoordinate;
template< class ColorFunction, class... Args >
explicit BasicInterfaceGridDataSet ( const ColorFunction &colorFunction, Args &&... args )
......@@ -66,7 +68,11 @@ namespace Dune
reconstruction_( colorFunction, reconstructionSet_, flags_ );
}
const GlobalCoordinate &normal ( const Element &element ) const { return reconstructionSet()[ element ].innerNormal(); }
const GlobalCoordinate &normal ( const Element &element ) const
{
assert( reconstructionSet()[ element ] );
return reconstructionSet()[ element ]->innerNormal();
}
protected:
Reconstruction reconstruction_;
......@@ -92,6 +98,8 @@ namespace Dune
using Base::reconstructionSet;
using Base::normal;
typedef typename Base::ReconstructionSet ReconstructionSet;
typedef typename Base::Element Element;
typedef typename Base::GridView GridView;
typedef typename Base::GlobalCoordinate GlobalCoordinate;
......@@ -109,7 +117,7 @@ namespace Dune
explicit InterfaceGridDataSet ( const ColorFunction &colorFunction, Args &&... args )
: Base( colorFunction, std::forward< Args >( args )... ), indices_( flags() )
{
getInterfaceVertices( reconstructionSet(), flags(), vertices_, offsets_ );
getInterfaceVertices( reconstructionSet(), vertices_, offsets_ );
}
template< class ColorFunction >
......@@ -117,7 +125,7 @@ namespace Dune
{
Base::update( colorFunction );
indices_.update( flags() );
getInterfaceVertices( reconstructionSet(), flags(), vertices_, offsets_ );
getInterfaceVertices( reconstructionSet(), vertices_, offsets_ );
}
void covariantOuterNormal ( const Element &element, std::size_t i, FieldVector< ctype, 2 > &n ) const
......@@ -166,7 +174,6 @@ namespace Dune
const Offsets &offsets () const { return offsets_; }
private:
template<class ReconstructionSet>
static inline void getInterfaceVertices (const ReconstructionSet& reconstructions, Vertices& vertices, Offsets& offsets)
{