From ceaf95c12677d041d1ba392abeb708d20bccae9f Mon Sep 17 00:00:00 2001
From: Lasse Hinrichsen <lh1887@mi.fu-berlin.de>
Date: Wed, 18 Oct 2017 16:29:14 +0200
Subject: [PATCH] Add implementation for geometryInInside and geometryInOutside
---
dune/subgrid/subgrid/subgridgeometry.hh | 152 ++++++++++++++------
dune/subgrid/subgrid/subgridintersection.hh | 74 ++++++++--
2 files changed, 175 insertions(+), 51 deletions(-)
diff --git a/dune/subgrid/subgrid/subgridgeometry.hh b/dune/subgrid/subgrid/subgridgeometry.hh
index 26d700e..578ca06 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 f8cef36..8e231df 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 {
--
GitLab