Add implementation for geometryInInside and geometryInOutside

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

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 {