diff --git a/dune/subgrid/CMakeLists.txt b/dune/subgrid/CMakeLists.txt
index 5f078aad77f62ec6cb17e1a9db05778c01c5b56a..c09b0c42fa656bcabbfd5ae6fc2d9088e1427a0a 100644
--- a/dune/subgrid/CMakeLists.txt
+++ b/dune/subgrid/CMakeLists.txt
@@ -1,3 +1,4 @@
+add_subdirectory(common)
add_subdirectory(facetools)
add_subdirectory(subgrid)
add_subdirectory(test)
diff --git a/dune/subgrid/common/CMakeLists.txt b/dune/subgrid/common/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..a31ec91c4cfb2205d700956f5a1633d72fe8b067
--- /dev/null
+++ b/dune/subgrid/common/CMakeLists.txt
@@ -0,0 +1,5 @@
+set(commonincludedir ${CMAKE_INSTALL_INCLUDEDIR}/dune/subgrid/common)
+set(commoninclude_HEADERS variant.hh)
+# include not needed for CMake
+# include $(top_srcdir)/am/global-rules
+install(FILES ${commoninclude_HEADERS} DESTINATION ${commonincludedir})
diff --git a/dune/subgrid/common/variant.hh b/dune/subgrid/common/variant.hh
new file mode 100644
index 0000000000000000000000000000000000000000..e2132ea6c66de8266669eb87f8b9a44fe6143837
--- /dev/null
+++ b/dune/subgrid/common/variant.hh
@@ -0,0 +1,387 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_COMMON_STD_VARIANT_HH
+#define DUNE_COMMON_STD_VARIANT_HH
+#include <tuple>
+#include <memory>
+#include <dune/common/hybridutilities.hh>
+#include <dune/common/exceptions.hh>
+
+namespace Dune {
+namespace Std {
+namespace Impl {
+
+ /* helper constructs to find position of a type T in a pack Ts... */
+ template <typename T, typename... Ts>
+ struct index_in_pack;
+
+ template <typename T, typename... Ts>
+ struct index_in_pack<T, T, Ts...> : std::integral_constant<std::size_t, 0> {};
+
+ template <typename T, typename U, typename... Ts>
+ struct index_in_pack<T, U, Ts...> : std::integral_constant<std::size_t, 1 + index_in_pack<T, Ts...>::value> {};
+
+ /* end helper constructs to find position of a type T in a pack Ts... */
+
+ template<typename Tp>
+ struct Buffer_ : std::aligned_storage<sizeof(Tp)> {
+ using Storage = typename std::aligned_storage_t<sizeof(Tp)>::type;
+ Storage storage_;
+
+ void* addr() {
+ return static_cast<void*>(&storage_);
+ }
+
+ const void* addr() const {
+ return static_cast<const void*>(&storage_);
+ }
+
+ Tp* ptr() {
+ return static_cast<Tp*>(addr());
+ }
+
+ const Tp* ptr() const {
+ return static_cast<const Tp*>(addr());
+ }
+ };
+
+ template<typename Tp, bool isTrivial>
+ struct TypeStorage_ { };
+
+ template<typename Tp>
+ struct TypeStorage_<Tp, true> {
+ TypeStorage_(Tp t) :
+ tp_(t) {}
+
+ template<typename... Args>
+ TypeStorage_(Args... args) :
+ tp_(args...) {}
+
+ auto& get() {
+ return tp_;
+ }
+ const auto& get() const {
+ return tp_;
+ }
+ private:
+ Tp tp_;
+ };
+
+ template<typename Tp>
+ struct TypeStorage_<Tp, false> {
+ TypeStorage_(Tp t) {
+ ::new (&tp_) Tp(t);
+ }
+
+ template<typename... Args>
+ TypeStorage_(Args... args) {
+ ::new (&tp_) Tp(std::forward<Args>(args)...);
+ }
+
+ auto& get() {
+ return *(tp_.ptr());
+ }
+ const auto& get() const {
+ return *(tp_.ptr());
+ }
+
+ private:
+ Buffer_<Tp> tp_;
+ };
+
+ template<typename... T>
+ union variant_union_ {};
+
+ template<typename Head_, typename... Tail_>
+ union variant_union_<Head_, Tail_...> {
+ constexpr variant_union_() :
+ tail_() {}
+
+ template<typename... Args>
+ constexpr variant_union_(std::integral_constant<size_t, 0>, Args&&... args) :
+ head_(std::forward<Args...>(args)...) {}
+
+ template<size_t N, typename... Args>
+ constexpr variant_union_(std::integral_constant<size_t, N>, Args&&... args) :
+ tail_(std::integral_constant<size_t, N-1>(), std::forward<Args...>(args)...) {}
+
+ auto& getByIndex(std::integral_constant<size_t, 0>) {
+ return head_.get();
+ }
+
+ const auto& getByIndex(std::integral_constant<size_t, 0>) const {
+ return head_.get();
+ }
+
+
+ template<size_t N>
+ auto& getByIndex(std::integral_constant<size_t, N>) {
+ return tail_.getByIndex(std::integral_constant<size_t, N-1>());
+ }
+
+ template<size_t N>
+ const auto& getByIndex(std::integral_constant<size_t, N>) const {
+ return tail_.getByIndex(std::integral_constant<size_t, N-1>());
+ }
+
+ template<typename Tp>
+ void set(Tp obj) {
+ Dune::Hybrid::ifElse(std::is_same<Tp, Head_>(),
+ [&](auto&& id) { head_=std::move(id(obj)); },
+ [&](auto&& id) { return id(tail_).set(std::move(obj)); }
+ );
+ }
+
+ constexpr size_t size() const {
+ return sizeof...(Tail_)+1;
+ }
+
+ private:
+ TypeStorage_<Head_, std::is_literal_type<Head_>::value> head_;
+ variant_union_<Tail_...> tail_;
+ };
+
+ template<typename...T>
+ struct variant_{
+
+ constexpr variant_() :
+ unions_() {}
+
+ template<typename Tp>
+ constexpr variant_(Tp obj) :
+ unions_(),
+ index_(index_in_pack<Tp, T...>::value)
+ {
+ unions_.set(std::move(obj));
+ }
+
+ template<typename Tp>
+ auto& get() {
+ constexpr size_t idx = index_in_pack<Tp, T...>::value;
+ if (index_ != idx)
+ DUNE_THROW(Dune::Exception, "Bad variant access.");
+
+ return get<idx>();
+ }
+
+ template<typename Tp>
+ const auto& get() const {
+ constexpr size_t idx = index_in_pack<Tp, T...>::value;
+ if (index_ != idx)
+ DUNE_THROW(Dune::Exception, "Bad variant access.");
+
+ return get<idx>();
+ }
+
+ template<typename Tp>
+ Tp* get_if() {
+ if (not holds_alternative<Tp>())
+ return (Tp*) nullptr;
+ else
+ return &(get<Tp>());
+ }
+
+ template<typename Tp>
+ const Tp* get_if() const {
+ if (not holds_alternative<Tp>())
+ return (Tp*) nullptr;
+ else
+ return &(get<Tp>());
+ }
+
+ template<size_t N>
+ auto* get_if() {
+ using Tp = std::decay_t<decltype(get<N>())>;
+ if (not holds_alternative<N>())
+ return (Tp*) nullptr;
+ else
+ return &(get<Tp>());
+ }
+
+ template<size_t N>
+ const auto* get_if() const {
+ using Tp = std::decay_t<decltype(get<N>())>;
+ if (not holds_alternative<N>())
+ return (Tp*) nullptr;
+ else
+ return &(get<Tp>());
+ }
+
+ template<size_t N>
+ auto& get() {
+ if (index_ != N)
+ DUNE_THROW(Dune::Exception, "Bad variant access.");
+ return unions_.template getByIndex(std::integral_constant<size_t, N>());
+ }
+ template<size_t N>
+ const auto& get() const {
+ if (index_ != N)
+ DUNE_THROW(Dune::Exception, "Bad variant access.");
+ return unions_.template getByIndex(std::integral_constant<size_t, N>());
+ }
+
+ template<typename Tp>
+ constexpr Tp& operator=(Tp obj) {
+ unions_.set(std::move(obj));
+ constexpr auto index = index_in_pack<Tp, T...>::value;
+ index_=index;
+ return unions_.getByIndex(std::integral_constant<size_t,index>());
+ }
+
+ constexpr std::size_t index() const noexcept {
+ return index_;
+ }
+
+ constexpr auto size() const {
+ return sizeof...(T);
+ }
+
+ /* \brief Apply visitor to the active variant.
+ *
+ * visit assumes that the result of
+ * func(T) has the same type for all types T
+ * in this variant.
+ */
+ template<typename F>
+ auto visit(F&& func) {
+ using namespace Dune::Hybrid;
+
+ using Result = decltype(func(unions_.getByIndex(std::integral_constant<size_t, 0>())));
+
+ return ifElse(std::is_same<Result, void>(), [&, this](auto id) {
+ constexpr auto tsize = size_;
+ Dune::Hybrid::forEach(Dune::Hybrid::integralRange(std::integral_constant<size_t, tsize>()), [&](auto i) {
+ if (i==this->index_)
+ func(id(unions_).getByIndex(std::integral_constant<size_t, i>()));
+ });
+ return;},
+ [&func,this](auto id) {
+ constexpr auto tsize = size_;
+
+ auto result = std::unique_ptr<Result>();
+
+ Dune::Hybrid::forEach(Dune::Hybrid::integralRange(std::integral_constant<size_t, tsize>()), [&, this](auto i) {
+ if (i==this->index_)
+ result = std::make_unique<Result>(func(id(this->unions_).getByIndex(std::integral_constant<size_t, i>())));
+ });
+ return *result;
+ });
+ }
+
+ template<typename F>
+ auto visit(F&& func) const {
+ using namespace Dune::Hybrid;
+
+ using Result = decltype(func(unions_.getByIndex(std::integral_constant<size_t, 0>())));
+
+ return ifElse(std::is_same<Result, void>(), [&, this](auto id) {
+ constexpr auto tsize = size_;
+ Dune::Hybrid::forEach(Dune::Hybrid::integralRange(std::integral_constant<size_t, tsize>()), [&](auto i) {
+ if (i==this->index_)
+ func(id(unions_).getByIndex(std::integral_constant<size_t, i>()));
+ });
+ return;},
+ [&func,this](auto id) {
+ constexpr auto tsize = size_;
+
+ auto result = std::unique_ptr<Result>();
+
+ Dune::Hybrid::forEach(Dune::Hybrid::integralRange(std::integral_constant<size_t, tsize>()), [&, this](auto i) {
+ if (i==this->index_)
+ result = std::make_unique<Result>(func(id(this->unions_).getByIndex(std::integral_constant<size_t, i>())));
+ });
+ return *result;
+ });
+ }
+
+ /** \brief Check if a given type is the one that is currently active in the variant. */
+ template<typename Tp>
+ constexpr bool holds_alternative() const {
+ // I have no idea how this could be really constexpr, but for STL-conformity,
+ // I'll leave the modifier there.
+ return (index_in_pack<Tp, T...>::value == index_);
+ }
+
+ /** \brief Check if a given type is the one that is currently active in the variant. */
+ template<size_t N>
+ constexpr bool holds_alternative() const {
+ // I have no idea how this could be really constexpr, but for STL-conformity,
+ // I'll leave the modifier there.
+ return (N == index_);
+ }
+
+ private:
+ variant_union_<T...> unions_;
+ std::size_t index_;
+ constexpr static auto size_ = sizeof...(T);
+ };
+
+} // end namespace Impl
+
+ /** \brief Incomplete re-implementation of C++17's std::variant. */
+ template<typename ...T>
+ using variant = Impl::variant_<T...>;
+
+ template<size_t N, typename... T>
+ auto& get(variant<T...>& var) {
+ return var.template get<N>();
+ }
+
+ template<size_t N, typename... T>
+ const auto& get(const variant<T...>& var) {
+ return var.template get<N>();
+ }
+
+ template<typename F, typename... T>
+ auto visit(F&& visitor, variant<T...>& var) {
+ return var.visit(std::forward<F>(visitor));
+ }
+
+ template<typename F, typename... T>
+ auto visit(F&& visitor, const variant<T...>& var) {
+ return var.visit(std::forward<F>(visitor));
+ }
+
+ template<typename Tp, typename ...T>
+ auto& get(variant<T...>& var) {
+ return var.template get<Tp>();
+ }
+
+ template<typename Tp, typename ...T>
+ const auto& get(const variant<T...>& var) {
+ return var.template get<Tp>();
+ }
+
+ template<typename Tp, typename ...T>
+ const auto* get_if(const variant<T...>& var) {
+ return var.template get_if<Tp>();
+ }
+
+ template<typename Tp, typename ...T>
+ auto* get_if(variant<T...>& var) {
+ return var.template get_if<Tp>();
+ }
+
+ template<size_t N, typename ...T>
+ const auto* get_if(const variant<T...>& var) {
+ return var.template get_if<N>();
+ }
+
+ template<size_t N, typename ...T>
+ auto* get_if(variant<T...>& var) {
+ return var.template get_if<N>();
+ }
+
+ template<typename Tp, typename ...T>
+ constexpr bool holds_alternative(const variant<T...>& var) {
+ return var.template holds_alternative<Tp>();
+ }
+
+ template <typename... T>
+ constexpr auto variant_size_v(const variant<T...>&) {
+ return std::integral_constant<std::size_t,sizeof...(T)>::value;
+ }
+
+} // end namespace Std
+} // end namespace Dune
+#endif
diff --git a/dune/subgrid/subgrid/subgridgeometry.hh b/dune/subgrid/subgrid/subgridgeometry.hh
index 26d700e3dd228a6e0c5f53ae1a19a003596345a5..578ca06328afd5d4d24e7543e45b28124259296b 100644
--- a/dune/subgrid/subgrid/subgridgeometry.hh
+++ b/dune/subgrid/subgrid/subgridgeometry.hh
@@ -5,9 +5,11 @@
* \brief The SubGridGeometry class and its specializations
*/
#include <type_traits>
+#include <dune/subgrid/common/variant.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/typetraits.hh>
+#include <dune/geometry/multilineargeometry.hh>
namespace Dune {
@@ -22,13 +24,10 @@ template<int mydim, int coorddim, class GridImp> class SubGridLocalGeometry;
because the Dune grid interface specifies the number and types of template
parameters of that class.
*/
-template<int mydim, int coorddim, class GridImp, bool isLocal>
-class SubGridGeometryBase
+template<int mydim, int coorddim, class GridImp>
+class SubGridGeometry
{
- private:
-
- typedef typename GridImp::ctype ctype;
-
+ using ctype = typename GridImp::ctype;
public:
@@ -36,12 +35,7 @@ class SubGridGeometryBase
enum {CodimInHostGrid = GridImp::HostGridType::dimension - mydim};
enum {DimensionWorld = GridImp::HostGridType::dimensionworld};
- // Select the type we are wrapping depending on the isLocal flag
- // - isLocal == true: LocalGeometry
- // - isLocal == false: Geometry
- typedef typename std::conditional<isLocal,
- typename GridImp::HostGridType::Traits::template Codim<CodimInHostGrid>::LocalGeometry,
- typename GridImp::HostGridType::Traits::template Codim<CodimInHostGrid>::Geometry>::type HostGridGeometry;
+ using HostGridGeometry = typename GridImp::HostGridType::Traits::template Codim<CodimInHostGrid>::Geometry;
//! The type of the Jacobian matrix of the mapping from the reference element to this element
typedef typename HostGridGeometry::JacobianTransposed JacobianTransposed;
@@ -50,7 +44,7 @@ class SubGridGeometryBase
/** Constructor for a given host grid geometry
*/
- SubGridGeometryBase(const HostGridGeometry& hostGeometry)
+ SubGridGeometry(const HostGridGeometry& hostGeometry)
: hostGeometry_(hostGeometry)
{
}
@@ -120,59 +114,135 @@ class SubGridGeometryBase
return hostGeometry_.center();
}
-
//! The Jacobian matrix of the mapping from the reference element to this element
const JacobianTransposed jacobianTransposed (const FieldVector<ctype, mydim>& local) const {
return hostGeometry_.jacobianTransposed(local);
}
-
//! The inverse of the Jacobian matrix of the mapping from the reference element to this element
const JacobianInverseTransposed jacobianInverseTransposed (const FieldVector<ctype, mydim>& local) const {
return hostGeometry_.jacobianInverseTransposed(local);
}
-
+ private:
const HostGridGeometry hostGeometry_;
-
};
-/** \brief Encapsulate a host grid Geometry */
+/** \brief Local SubGrid geometries. They can be either a wrapped hostgrid geometry or a multilineargeometry. */
template<int mydim, int coorddim, class GridImp>
-class SubGridGeometry :
- public SubGridGeometryBase<mydim,coorddim,GridImp,false>
+class SubGridLocalGeometry
{
+
+ using ctype = typename GridImp::ctype;
+
public:
- enum {CodimInHostGrid = GridImp::HostGridType::dimension - mydim};
+ constexpr static int hostDim = GridImp::HostGridType::dimension;
+ constexpr static int dim = hostDim;
+ constexpr static int CodimInHostGrid = hostDim - mydim;
- // select appropiate hostgrid geometry via typeswitch
- typedef typename GridImp::HostGridType::Traits::template Codim<CodimInHostGrid>::Geometry HostGridGeometry;
+ private:
+ // the local geometry is either a wrapped hostgrid local geometry or a multilineargeometry
+ // TODO does this work for arbitratry co-dimensions?
+ using HostGridLocalGeometry = typename GridImp::HostGridType::Traits::template Codim<CodimInHostGrid>::LocalGeometry;
+ using MultiLinGeometry = MultiLinearGeometry<ctype, mydim, coorddim>;
- /** Constructor for a given host grid geometry
- */
- SubGridGeometry(const HostGridGeometry& hostGeometry)
- : SubGridGeometryBase<mydim,coorddim,GridImp,false>(hostGeometry)
+ // use a variant to store either a HostGridLocalGeometry or a MultiLinearGeometry
+ using GeometryContainer = Std::variant<HostGridLocalGeometry, MultiLinGeometry>;
+
+ public:
+
+ // The codimension of this entitypointer wrt the host grid
+ //enum {CodimInHostGrid = GridImp::HostGridType::dimension - mydim};
+ enum {DimensionWorld = GridImp::HostGridType::dimensionworld};
+
+
+ //! The type of the Jacobian matrix of the mapping from the reference element to this element
+ using JacobianTransposed = typename std::common_type<FieldMatrix< ctype, mydim, coorddim>,
+ typename HostGridLocalGeometry::JacobianTransposed>::type;
+
+ //! The type of the nverse of the Jacobian matrix of the mapping from the reference element to this element
+ using JacobianInverseTransposed = typename std::common_type<Dune::FieldMatrix<ctype,coorddim, mydim>,
+ typename HostGridLocalGeometry::JacobianInverseTransposed>::type;
+
+ /** \brief Construct from a given host grid geometry*/
+ SubGridLocalGeometry(HostGridLocalGeometry hostLocalGeometry)
+ : localGeometry_(std::move(hostLocalGeometry))
{}
-};
+ /** \brief Construct from MultiLinearGeometry */
+ SubGridLocalGeometry(MultiLinGeometry multiLinGeometry)
+ : localGeometry_(std::move(multiLinGeometry))
+ {}
-/** \brief Encapsulate a host grid LocalGeometry */
-template<int mydim, int coorddim, class GridImp>
-class SubGridLocalGeometry :
- public SubGridGeometryBase<mydim,coorddim,GridImp,true>
-{
- public:
- enum {CodimInHostGrid = GridImp::HostGridType::dimension - mydim};
+ /** \brief Return the element type identifier
+ */
+ GeometryType type () const {
+ return Std::visit([&](auto&& geom) {return geom.type();}, localGeometry_);
+ }
+
+ /** \brief Return true if the geometry mapping is affine and false otherwise */
+ bool affine() const {
+ return Std::visit([&](auto&& geom) {return geom.affine();}, localGeometry_);
+ }
+
+ //! return the number of corners of this element. Corners are numbered 0...n-1
+ int corners () const {
+ return Std::visit([&](auto&& geom) {return geom.corners();}, localGeometry_);
+ }
- // select appropiate hostgrid geometry via typeswitch
- typedef typename GridImp::HostGridType::Traits::template Codim<CodimInHostGrid>::LocalGeometry HostGridLocalGeometry;
+ //! access to coordinates of corners. Index is the number of the corner
+ FieldVector<ctype, coorddim> corner (int i) const {
+ return Std::visit([&](auto&& geom) {return geom.corner(i);}, localGeometry_);
+ }
- /** Constructor for a given host grid local geometry
+ /** \brief Maps a local coordinate within reference element to
+ * global coordinate in element */
+ FieldVector<ctype, coorddim> global (const FieldVector<ctype, mydim>& local) const{
+ return Std::visit([&](auto&& geom) {return geom.global(local);}, localGeometry_);
+ }
+
+ /** \brief Maps a global coordinate within the element to a
+ * local coordinate in its reference element */
+ FieldVector<ctype, mydim> local (const FieldVector<ctype, coorddim>& global) const {
+ return Std::visit([&](auto&& geom) {return geom.local(global);}, localGeometry_);
+ }
+
+ /**
*/
- SubGridLocalGeometry(const HostGridLocalGeometry& hostGeometry)
- : SubGridGeometryBase<mydim,coorddim,GridImp,true>(hostGeometry)
- {}
+ ctype integrationElement (const FieldVector<ctype, mydim>& local) const {
+ return Std::visit([&](auto&& geom) {return geom.integrationElement(local);}, localGeometry_);
+ }
+
+ /** \brief return volume of geometry */
+ ctype volume () const {
+ return Std::visit([&](auto&& geom) {return geom.volume();}, localGeometry_);
+ }
+ /** \brief return center of geometry
+ * Note that this method is still subject to a change of name and semantics.
+ * At the moment, the center is not required to be the centroid of the
+ * geometry, or even the centroid of its corners. This makes the current
+ * default implementation acceptable, which maps the centroid of the
+ * reference element to the geometry.
+ * We may change the name (and semantic) of the method to centroid() if we
+ * find reasonably efficient ways to implement it properly.
+ */
+ FieldVector<ctype, coorddim> center () const {
+ return Std::visit([&](auto&& geom) {return geom.center();}, localGeometry_);
+ }
+
+ //! The Jacobian matrix of the mapping from the reference element to this element
+ const JacobianTransposed jacobianTransposed (const FieldVector<ctype, mydim>& local) const {
+ return Std::visit([&](auto&& geom) {return static_cast<JacobianTransposed>(geom.jacobianTransposed(local));}, localGeometry_);
+ }
+
+ //! The inverse of the Jacobian matrix of the mapping from the reference element to this element
+ const JacobianInverseTransposed jacobianInverseTransposed (const FieldVector<ctype, mydim>& local) const {
+ return Std::visit([&](auto&& geom) {return static_cast<JacobianInverseTransposed>(geom.jacobianInverseTransposed(local));}, localGeometry_);
+ }
+
+ private:
+ const GeometryContainer localGeometry_;
};
} // namespace Dune
diff --git a/dune/subgrid/subgrid/subgridintersection.hh b/dune/subgrid/subgrid/subgridintersection.hh
index f8cef36b5326e8291add18e65061a800548b65dc..8e231df7b8757bd28611d59d1b21448414978207 100644
--- a/dune/subgrid/subgrid/subgridintersection.hh
+++ b/dune/subgrid/subgrid/subgridintersection.hh
@@ -5,6 +5,7 @@
#include <dune/subgrid/subgrid/subgridgeometry.hh>
#include <dune/grid/common/exceptions.hh>
+#include <dune/geometry/multilineargeometry.hh>
namespace Dune {
@@ -220,7 +221,7 @@ public:
//! Return the global geometry of the intersection
Geometry geometry () const {
- if (outsideIntersect_.inside().level() >= insideIntersect_.inside().level())
+ if (outsideIntersect_.inside().level() > insideIntersect_.inside().level())
return Geometry(GeometryImpl(outsideIntersect_.geometry()));
else
return Geometry(GeometryImpl(insideIntersect_.geometry()));
@@ -231,21 +232,74 @@ public:
return boundary() or ((outside_.level() == inside_.level()) and insideIntersect_.conforming());
}
- //! intersection of codimension 1 of this neighbor with element where
- //! iteration started.
- //! Here returned element is in LOCAL coordinates of the element
- //! where iteration started.
+ /** \brief Return the local geometry, mapping local coordinates of the intersection to
+ * local coordinates of the inside element.
+ * If the intersection is not part of the host grid, construct one using MultiLinearGeometry.
+ */
LocalGeometry geometryInInside () const {
if (boundary())
return LocalGeometry(LocalGeometryImpl(insideIntersect_.geometryInInside()));
- else
- DUNE_THROW(NotImplemented, "SubGridLeafIntersection::geometryInInside()");
+
+ // Case 0: inside_'s the smaller element, or both elements are on the same level.
+ if (inside_.level() >= outside_.level())
+ return LocalGeometry(LocalGeometryImpl(insideIntersect_.geometryInInside()));
+
+ // Case 1: outsideIntersect_ is the intersection we want, embedded in inside_'s geometry
+
+ auto outsideIsGeo = outsideIntersect_.geometry();
+ auto insideGeo = inside_.geometry();
+
+ using Coords = std::decay_t<decltype(insideGeo.local(outsideIsGeo.corner(0)))>;
+ auto corners = std::vector<Coords>(outsideIsGeo.corners());
+
+ // Map the coordinates of the corners from the outside's local coordinates into inside's
+ // local coordinates
+ for (size_t i = 0; i < corners.size(); i++)
+ corners[i] = insideGeo.local(outsideIsGeo.corner(i));
+
+ // create a new geometry with the coordinates we just computed
+ auto geo = Dune::MultiLinearGeometry<ctype, dim-1, dim>(outsideIsGeo.type(), corners);
+ return LocalGeometry(LocalGeometryImpl(std::move(geo)));
}
- //! intersection of codimension 1 of this neighbor with element where iteration started.
- //! Here returned element is in LOCAL coordinates of neighbor
+ /** \brief Return the local geometry, mapping local coordinates of the intersection to
+ * local coordinates of the outside element.
+ */
LocalGeometry geometryInOutside () const {
- DUNE_THROW(NotImplemented, "SubGridLeafIntersection::geometryInOutside()");
+ if(boundary())
+ DUNE_THROW(Dune::Exception, "No outside geometry available for boundary intersections");
+
+ // Case 0: outside_'s the smaller element, or both elements are on the same level.
+ if (inside_.level() < outside_.level()) {
+
+ auto outsideIsGeo = outsideIntersect_.geometry();
+ auto outsideGeo = outside_.geometry();
+
+ using Coords = std::decay_t<decltype(outsideGeo.local(outsideIsGeo.corner(0)))>;
+ auto corners = std::vector<Coords>(outsideIsGeo.corners());
+ for (size_t i = 0; i < corners.size(); i++)
+ corners[i] = outsideGeo.local(outsideIsGeo.corner(i));
+
+ // create a new geometry with the coordinates we just computed
+ auto geo = Dune::MultiLinearGeometry<ctype, dim-1, dim>(outsideIsGeo.type(), corners);
+ return LocalGeometry(LocalGeometryImpl(geo));
+ }
+
+ // Case 1: outsideIntersect_ is the intersection we want, embedded in inside_'s geometry
+ auto insideIsGeo = insideIntersect_.geometry();
+ auto outsideGeo = outside_.geometry();
+
+ using Coords = std::decay_t<decltype(outsideGeo.local(insideIsGeo.corner(0)))>;
+ auto corners = std::vector<Coords>(insideIsGeo.corners());
+
+ // Map the coordinates of the corners from the insides's local coordinates into outside's
+ // local coordinates
+ for (size_t i = 0; i < corners.size(); i++)
+ corners[i] = outsideGeo.local(insideIsGeo.corner(i));
+
+ // create a new geometry with the coordinates we just computed
+ auto geo = Dune::MultiLinearGeometry<ctype, dim-1, dim>(insideIsGeo.type(), corners);
+ return LocalGeometry(LocalGeometryImpl(std::move(geo)));
}
bool equals(const SubGridLeafIntersection<GridImp> & other) const {
diff --git a/dune/subgrid/test/CMakeLists.txt b/dune/subgrid/test/CMakeLists.txt
index 1cb2c78ee5dd68fa6288e0e83b9effc9bfeb25ce..735fbd88d618686fdf69086a3f87a445c943cdae 100644
--- a/dune/subgrid/test/CMakeLists.txt
+++ b/dune/subgrid/test/CMakeLists.txt
@@ -1,5 +1,5 @@
-set(TESTPROGS test-w-onedgrid test-w-yaspgrid test-w-alugrid)
+set(TESTPROGS test-w-onedgrid test-w-yaspgrid test-w-alugrid testvariant)
if(ALBERTA_FOUND)
list(APPEND TESTPROGS test-w-albertagrid2d)
diff --git a/dune/subgrid/test/testvariant.cc b/dune/subgrid/test/testvariant.cc
new file mode 100644
index 0000000000000000000000000000000000000000..ca0393f89e1c53ab03a65ce361ddeea4fe45d7bb
--- /dev/null
+++ b/dune/subgrid/test/testvariant.cc
@@ -0,0 +1,107 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+
+#ifdef HAVE_CONFIG_H
+# include "config.h"
+#endif
+#include <iostream>
+#include <dune/common/parallel/mpihelper.hh> // An initializer of MPI
+#include <dune/common/exceptions.hh> // We use exceptions
+#include <dune/common/test/testsuite.hh>
+
+#include <dune/subgrid/common/variant.hh>
+
+// some non-default constructible type
+struct F {
+ int i;
+ F() = delete;
+ F(int j) :
+ i(j) {}
+};
+
+Dune::TestSuite testVariant() {
+ using namespace Dune;
+ TestSuite suite;
+
+ int i = 42;
+ double d = 3.14;
+ F f(13);
+ using V = std::vector<int>;
+
+ auto variant = Std::variant<int, double, F, V>();
+
+ suite.check(Std::variant_size_v(variant) == 4, "Test variant_size_v");
+
+
+ variant = d;
+ suite.check(Std::holds_alternative<double>(variant), "Test holds_alternative");
+
+ variant = f;
+ suite.check(Std::holds_alternative<F>(variant), "Test holds_alternative");
+
+ variant = i;
+ suite.check(Std::holds_alternative<int>(variant), "Test holds_alternative");
+
+ suite.check(Std::get<int>(variant) == i, "Test get<Type>");
+ suite.check(Std::get<0>(variant) == i, "Test get<Index>");
+
+ suite.check(Std::get_if<int>(variant) != nullptr, "Test get_if on right type");
+ suite.check(Std::get_if<double>(variant) == nullptr, "Test get_if on wrong type");
+
+ suite.check(Std::get_if<0>(variant) != nullptr, "Test get_if on right index");
+ suite.check(Std::get_if<1>(variant) == nullptr, "Test get_if on wrong index");
+
+ // test if get<Type> throws if one tries to get the wrong type:
+ try {
+ // currently hold type is still int, so double should throw
+ Std::get<double>(variant);
+ suite.check(false, "Test get<Type> on wrong type should have thrown");
+ }
+ catch (...) {
+ suite.check(true, "Test get<Type> on wrong type has thrown");
+ }
+
+
+ variant = V(1);
+ suite.check(Std::get<V>(variant).size() == 1, "Test with non-trivial type");
+
+ variant = f;
+
+ suite.check(variant.index() == 2, "Test index()"); // we're at type F, which has position 2
+
+ // Demonstrate visit concept and using vector as an example of a non-POD type
+ using V2 = std::vector<double>;
+ Std::variant<V, V2> variant2;
+ variant2 = V(1);
+ auto size = [](auto&& v) {return v.size();};
+ suite.check(Std::visit(size, variant2)== 1, "Test visit");
+ variant2 = V2(2);
+ suite.check(Std::visit(size, variant2)== 2, "Test visit");
+
+ // try on a const reference:
+ const auto& constv2 = variant2;
+ suite.check(Std::visit(size, constv2)== 2, "Test const visit");
+ suite.check(Std::get_if<V2>(constv2) != nullptr, "Test const get_if");
+
+
+ return suite;
+}
+
+int main(int argc, char** argv)
+{
+ try{
+ // Maybe initialize MPI
+ Dune::MPIHelper::instance(argc, argv);
+
+ Dune::TestSuite suite;
+ suite.subTest(testVariant());
+
+ return suite.exit();
+ }
+ catch (Dune::Exception &e){
+ std::cerr << "Dune reported error: " << e << std::endl;
+ }
+ catch (...){
+ std::cerr << "Unknown exception thrown!" << std::endl;
+ }
+}