Commit 7f8a2779 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

corrected the exact derivative of sphere function and wrapped CurvedGeometry...

corrected the exact derivative of sphere function and wrapped CurvedGeometry to support localtransformation of local coordinates
parent db972060
...@@ -205,31 +205,49 @@ namespace Dune ...@@ -205,31 +205,49 @@ namespace Dune
using Super = EntityBase<codim, G>; using Super = EntityBase<codim, G>;
using Traits = typename std::remove_const_t<G>::Traits; using Traits = typename std::remove_const_t<G>::Traits;
//! type of host elements, i.e., of host entities of codimension 0
using HostElement = typename Traits::HostGrid::template Codim<0>::Entity;
using GridFunction = typename Super::GridFunction;
public: public:
//! type of corresponding geometry //! type of corresponding geometry
using Geometry = typename Traits::template Codim<codim>::Geometry; using Geometry = typename Traits::template Codim<codim>::Geometry;
using ctype = typename Geometry::ctype;
public: public:
using Super::Super; using Super::Super;
// construct the entity from a subentity of a host-entity
EntityBaseGeometry (const GridFunction& gridFunction, const HostElement& hostElement, int i)
: Super(gridFunction, hostElement, i)
, hostElement_(hostElement)
, subEntity_(i)
{}
//! obtain the geometry of this entity //! obtain the geometry of this entity
Geometry geometry () const Geometry geometry () const
{ {
if (!geo_) { if (!geo_) {
auto localFct = localFunction(Super::gridFunction()); auto localFct = localFunction(Super::gridFunction());
// localFct.bind(Super::hostEntity()); if (hostElement_) {
localFct.bind(*hostElement_);
// construct a fake local geometry for transformations.
// auto refElem = referenceElement<typename Geometry::ctype, Super::dimension>(GeometryTypes::simplex(Super::dimension)); auto refElem = referenceElement<ctype, Super::dimension>(hostElement_->type());
geo_ = std::make_shared<GeometryImpl>(Super::type(), localFct); auto localGeometry = refElem.template geometry<codim>(subEntity_);
geo_ = std::make_shared<GeometryImpl>(Super::type(), localFct, localGeometry);
}
else {
DUNE_THROW(Dune::NotImplemented, "Geometry of entities of codim!=0 not implemented");
}
} }
return Geometry(*geo_); return Geometry(*geo_);
//return Super::hostEntity().geometry();
} }
private: private:
Std::optional<HostElement> hostElement_;
int subEntity_ = -1;
using GeometryImpl = typename Traits::template Codim<codim>::GeometryImpl; using GeometryImpl = typename Traits::template Codim<codim>::GeometryImpl;
mutable std::shared_ptr<GeometryImpl> geo_; mutable std::shared_ptr<GeometryImpl> geo_;
}; };
...@@ -265,7 +283,8 @@ namespace Dune ...@@ -265,7 +283,8 @@ namespace Dune
if (!geo_) { if (!geo_) {
auto localFct = localFunction(Super::gridFunction()); auto localFct = localFunction(Super::gridFunction());
localFct.bind(Super::hostEntity()); localFct.bind(Super::hostEntity());
geo_ = std::make_shared<GeometryImpl>(Super::type(), localFct, Dune::CGeo::Impl::DefaultLocalGeometry{}); auto fakeDefaultGeometry = Dune::CGeo::Impl::DefaultLocalGeometry{};
geo_ = std::make_shared<GeometryImpl>(Super::type(), localFct, fakeDefaultGeometry);
} }
return Geometry(*geo_); return Geometry(*geo_);
......
...@@ -16,6 +16,8 @@ ...@@ -16,6 +16,8 @@
#include <dune/common/std/optional.hh> #include <dune/common/std/optional.hh>
#include <dune/common/std/type_traits.hh> #include <dune/common/std/type_traits.hh>
#include <dune/curvedgeometry/curvedgeometry.hh>
#include <dune/geometry/affinegeometry.hh> #include <dune/geometry/affinegeometry.hh>
#include <dune/geometry/multilineargeometry.hh> #include <dune/geometry/multilineargeometry.hh>
#include <dune/geometry/quadraturerules.hh> #include <dune/geometry/quadraturerules.hh>
...@@ -90,6 +92,7 @@ namespace Dune ...@@ -90,6 +92,7 @@ namespace Dune
} }
// efficient implementation of a trivial local geometry for elements
struct DefaultLocalGeometry struct DefaultLocalGeometry
{ {
template <class LocalCoordinate> template <class LocalCoordinate>
...@@ -104,23 +107,93 @@ namespace Dune ...@@ -104,23 +107,93 @@ namespace Dune
return {}; return {};
} }
}; };
// type-erased wrapper for local geometries
template <class LG>
struct LocalGeometryInterface
{
using LocalCoordinate = typename LG::LocalCoordinate;
using GlobalCoordinate = typename LG::GlobalCoordinate;
using JacobianTransposed = typename LG::JacobianTransposed;
template <class LocalGeometry>
LocalGeometryInterface(const LocalGeometry& lg)
: global([lg](auto const& local) { return lg.global(local); })
, jacobianTransposed([lg](auto const& local) { return lg.jacobianTransposed(local); })
{}
template <class GlobalFct, class JacobianTransposedFct>
LocalGeometryInterface(GlobalFct&& g, JacobianTransposedFct&& jt)
: global(std::forward<GlobalFct>(g))
, jacobianTransposed(std::forward<JacobianTransposedFct>(jt))
{}
std::function<GlobalCoordinate(LocalCoordinate)> global;
std::function<JacobianTransposed(LocalCoordinate)> jacobianTransposed;
};
} }
/// \brief Curved geometry implementation based on local basis function parametrization /// \brief Curved geometry implementation based on local basis function parametrization
/** /**
* Parametrization of the geometry by any localfunction given by a finite-element cache * Parametrization of the geometry by a differentiable localFunction
* defined in the Traits class. See \ref LagrangeGeometry for an example instantiation
* of this class with local lagrange basis functions.
* *
* \tparam ct coordinate type * \tparam ct coordinate type
* \tparam mydim geometry dimension * \tparam mydim geometry dimension
* \tparam cdim coordinate dimension * \tparam cdim coordinate dimension
* * \tparam LF localFunction parametrizing the geometry
* The requirements on the traits are documented along with their default, * \tparam LG localGeometry for coordinate transform from subEntity to element,
* GeometryTraits. * see \ref Impl::DefaultLocalGeometry and \ref Impl::LocalGeometryInterface
* \tparam order Polynomial order of lagrange parametrization. If order < 0 use localFunction.
*/ */
template <class ct, int mydim, int cdim, class LocalFunction, class LocalGeometry = Impl::DefaultLocalGeometry> template <class ct, int mydim, int cdim, class LF, class LG, int order>
class Geometry class Geometry
: public Dune::LagrangeCurvedGeometry<ct, mydim, cdim, (order > 0 ? order : 1)>
{
using Super = Dune::LagrangeCurvedGeometry<ct, mydim, cdim, (order > 0 ? order : 1)>;
//! type of the localFunction representation the geometry parametrization
using LocalFunction = LF;
//! type of coordinate transformation for subEntities to codim=0 entities
using LocalGeometry = LG;
/// type of reference element
using ReferenceElements = Dune::ReferenceElements<ct, mydim>;
using ReferenceElement = typename ReferenceElements::ReferenceElement;
public:
/// \brief constructor from localFunction to parametrize the geometry
/**
* \param[in] refElement reference element for the geometry
* \param[in] localFunction localFunction for the parametrization of the geometry
* (stored by value)
* \param[in] args... argument to construct the local geometry
*
*/
template <class... Args>
Geometry (const ReferenceElement& refElement, const LocalFunction& localFunction, Args&&... args)
: Super(refElement, [localFunction, localGeometry=LocalGeometry{std::forward<Args>(args)...}](auto const& local)
{
return localFunction(localGeometry.global(local));
})
{}
/// \brief constructor, forwarding to the other constructor that take a reference-element
/**
* \param[in] gt geometry type
* \param[in] localFunction localFunction for the parametrization of the geometry
* (stored by value)
* \param[in] args... argument to construct the local geometry
*/
template <class... Args>
Geometry (Dune::GeometryType gt, const LocalFunction& localFunction, Args&&... args)
: Geometry(ReferenceElements::general(gt), localFunction, std::forward<Args>(args)...)
{}
};
template <class ct, int mydim, int cdim, class LF, class LG>
class Geometry<ct,mydim,cdim,LF,LG,-1>
{ {
public: public:
/// coordinate type /// coordinate type
...@@ -156,9 +229,14 @@ namespace Dune ...@@ -156,9 +229,14 @@ namespace Dune
//! helper structure containing some matrix routines. See affinegeometry.hh //! helper structure containing some matrix routines. See affinegeometry.hh
using MatrixHelper = Dune::Impl::FieldMatrixHelper<ct>; using MatrixHelper = Dune::Impl::FieldMatrixHelper<ct>;
//! type of coordinate transformation for subEntities to codim=0 entities
using LocalGeometry = LG;
//! type of the localFunction representation the geometry parametrization
using LocalFunction = LF;
//! type of the LocalFunction representing the derivative //! type of the LocalFunction representing the derivative
using DerivativeLocalFunction using DerivativeLocalFunction = std::decay_t<decltype(derivative(std::declval<LF>()))>;
= std::decay_t<decltype(derivative(std::declval<LocalFunction>()))>;
//! tolerance to numerical algorithms //! tolerance to numerical algorithms
static ct tolerance () { return ct(16) * std::numeric_limits<ct>::epsilon(); } static ct tolerance () { return ct(16) * std::numeric_limits<ct>::epsilon(); }
...@@ -167,31 +245,31 @@ namespace Dune ...@@ -167,31 +245,31 @@ namespace Dune
static int maxIteration () { return 100; } static int maxIteration () { return 100; }
public: public:
/// \brief constructor from a vector of coefficients of the LocalBasis parametrizing /// \brief constructor from localFunction to parametrize the geometry
/// the geometry
/** /**
* \param[in] refElement reference element for the geometry * \param[in] refElement reference element for the geometry
* \param[in] vertices vertices to store internally * \param[in] localFunction localFunction for the parametrization of the geometry
* (stored by value)
* \param[in] args... argument to construct the local geometry
* *
* \note The vertices are stored internally, so if possible move an external vertex
* storedge to this constructor
*/ */
Geometry (const ReferenceElement& refElement, const LocalFunction& localFunction, template <class... Args>
Std::optional<LocalGeometry> localGeometry = {}) Geometry (const ReferenceElement& refElement, const LocalFunction& localFunction, Args&&... args)
: refElement_(refElement) : refElement_(refElement)
, localFunction_(localFunction) , localFunction_(localFunction)
, localGeometry_(localGeometry) , localGeometry_(std::forward<Args>(args)...)
{} {}
/// \brief constructor, forwarding to the other constructor that take a reference-element /// \brief constructor, forwarding to the other constructor that take a reference-element
/** /**
* \param[in] gt geometry type * \param[in] gt geometry type
* \param[in] localFunction Function that can be evaluated in local coordinates and * \param[in] localFunction localFunction for the parametrization of the geometry
* can be differentiated * (stored by value)
* \param[in] args... argument to construct the local geometry
*/ */
Geometry (Dune::GeometryType gt, const LocalFunction& localFunction, template <class... Args>
Std::optional<LocalGeometry> localGeometry = {}) Geometry (Dune::GeometryType gt, const LocalFunction& localFunction, Args&&... args)
: Geometry(ReferenceElements::general(gt), localFunction, localGeometry) : Geometry(ReferenceElements::general(gt), localFunction, std::forward<Args>(args)...)
{} {}
/// \brief is this mapping affine? No! /// \brief is this mapping affine? No!
...@@ -200,7 +278,7 @@ namespace Dune ...@@ -200,7 +278,7 @@ namespace Dune
return false; return false;
} }
/// \brief obtain the name of the reference element /// \brief obtain the element type from the reference element
Dune::GeometryType type () const Dune::GeometryType type () const
{ {
return refElement_.type(); return refElement_.type();
...@@ -213,10 +291,6 @@ namespace Dune ...@@ -213,10 +291,6 @@ namespace Dune
} }
/// \brief obtain coordinates of the i-th corner /// \brief obtain coordinates of the i-th corner
/**
* \todo In some local-bases, the corners are part of the stored vertices, and sometimes
* even stored first before higher order vertices. Maybe this can be utilized somehow.
*/
GlobalCoordinate corner (int i) const GlobalCoordinate corner (int i) const
{ {
assert( (i >= 0) && (i < corners()) ); assert( (i >= 0) && (i < corners()) );
...@@ -247,6 +321,9 @@ namespace Dune ...@@ -247,6 +321,9 @@ namespace Dune
* \param[in] globalCoord global coordinate to map * \param[in] globalCoord global coordinate to map
* *
* \return corresponding local coordinate * \return corresponding local coordinate
* \throws in case of an error indicating that the local coordinate can not be obtained,
* an exception is thrown. See \ref checkedLocal for a variant that returns
* an optional instead.
* *
* \note For given global coordinate `y` the returned local coordinate `x` that minimizes * \note For given global coordinate `y` the returned local coordinate `x` that minimizes
* the following function over the local coordinate space spanned by the reference element. * the following function over the local coordinate space spanned by the reference element.
...@@ -258,11 +335,19 @@ namespace Dune ...@@ -258,11 +335,19 @@ namespace Dune
{ {
auto localCoord = checkedLocal(globalCoord); auto localCoord = checkedLocal(globalCoord);
if (!localCoord) if (!localCoord)
DUNE_THROW(Exception, "local coordinate can not be recovered from given global coordinate " << globalCoord); DUNE_THROW(Dune::Exception,
"local coordinate can not be recovered from given global coordinate " << globalCoord);
return *localCoord; return *localCoord;
} }
/// \brief Evaluate the inverse coordinate mapping
/**
* \param[in] globalCoord global coordinate to map
*
* \return corresponding local coordinate or nothing, in case the local coordinate
* could not be calculated. The return value is wrapped in an optional.
**/
Std::optional<LocalCoordinate> checkedLocal (const GlobalCoordinate& globalCoord) const Std::optional<LocalCoordinate> checkedLocal (const GlobalCoordinate& globalCoord) const
{ {
LocalCoordinate x = flatGeometry().local(globalCoord); LocalCoordinate x = flatGeometry().local(globalCoord);
...@@ -317,9 +402,6 @@ namespace Dune ...@@ -317,9 +402,6 @@ namespace Dune
* \param[in] local local coordinate to evaluate the integration element in * \param[in] local local coordinate to evaluate the integration element in
* *
* \returns the integration element \f$\mu(x)\f$. * \returns the integration element \f$\mu(x)\f$.
*
* \todo Implement caching of the jacobianTransposed matrix, since it is used for the
* integrationElement and the jacobianInverseTransposed
*/ */
ctype integrationElement (const LocalCoordinate& local) const ctype integrationElement (const LocalCoordinate& local) const
{ {
...@@ -356,11 +438,13 @@ namespace Dune ...@@ -356,11 +438,13 @@ namespace Dune
derivativeLocalFunction_->bind(localFunction().localContext()); derivativeLocalFunction_->bind(localFunction().localContext());
} }
// TODO: maybe multiply with localGeometry().jacobianTransposed() // coordinate in the localContext of the localFunction
auto j0 = localGeometry().jacobianTransposed(local); auto&& elementLocal = localGeometry().global(local);
auto jt = hostGeometry().jacobianTransposed(localGeometry().global(local));
auto jacobian = (*derivativeLocalFunction_)(localGeometry().global(local)); auto jLocal = localGeometry().jacobianTransposed(local);
return Impl::AB(j0,Impl::ABt(jt,jacobian)); auto jGlobal = hostGeometry().jacobianTransposed(elementLocal);
auto jacobian = (*derivativeLocalFunction_)(elementLocal);
return Impl::AB(jLocal,Impl::ABt(jGlobal,jacobian));
} }
/// \brief obtain the transposed of the Jacobian's inverse /// \brief obtain the transposed of the Jacobian's inverse
...@@ -410,7 +494,7 @@ namespace Dune ...@@ -410,7 +494,7 @@ namespace Dune
return res /= res.two_norm(); return res /= res.two_norm();
} }
public: private:
const LocalFunction& localFunction () const const LocalFunction& localFunction () const
{ {
return localFunction_; return localFunction_;
...@@ -418,7 +502,7 @@ namespace Dune ...@@ -418,7 +502,7 @@ namespace Dune
const LocalGeometry& localGeometry () const const LocalGeometry& localGeometry () const
{ {
return localGeometry_.value(); return localGeometry_;
} }
using FlatGeometry = MultiLinearGeometry<ctype, mydim, cdim>; using FlatGeometry = MultiLinearGeometry<ctype, mydim, cdim>;
...@@ -436,7 +520,7 @@ namespace Dune ...@@ -436,7 +520,7 @@ namespace Dune
return *flatGeometry_; return *flatGeometry_;
} }
using HostGeometry = std::decay_t<decltype(std::declval<LocalFunction>().localContext().geometry())>; using HostGeometry = std::decay_t<decltype(std::declval<LF>().localContext().geometry())>;
const HostGeometry& hostGeometry () const const HostGeometry& hostGeometry () const
{ {
if (!hostGeometry_) if (!hostGeometry_)
...@@ -453,7 +537,7 @@ namespace Dune ...@@ -453,7 +537,7 @@ namespace Dune
LocalFunction localFunction_; LocalFunction localFunction_;
//! transformation of local coordinates to element-local coordinates //! transformation of local coordinates to element-local coordinates
Std::optional<LocalGeometry> localGeometry_; LocalGeometry localGeometry_;
// some data optionally provided // some data optionally provided
mutable Std::optional<DerivativeLocalFunction> derivativeLocalFunction_; mutable Std::optional<DerivativeLocalFunction> derivativeLocalFunction_;
...@@ -465,15 +549,20 @@ namespace Dune ...@@ -465,15 +549,20 @@ namespace Dune
#ifndef DOXYGEN #ifndef DOXYGEN
// Specialization for vertex geometries // Specialization for vertex geometries
template <class ct, int cdim, class Localfunction> template <class ct, int cdim, class LF, class LG>
class Geometry<ct,0,cdim,Localfunction> class VertexGeometry
: public AffineGeometry<ct,0,cdim> : public AffineGeometry<ct,0,cdim>
{ {
using Super = AffineGeometry<ct,0,cdim>; using Super = AffineGeometry<ct,0,cdim>;
using LocalFunction = LF;
using LocalGeometry = LG;
//! tolerance to numerical algorithms //! tolerance to numerical algorithms
static ct tolerance () { return ct(16) * std::numeric_limits<ct>::epsilon(); } static ct tolerance () { return ct(16) * std::numeric_limits<ct>::epsilon(); }
struct Construct {};
public: public:
/// type of reference element /// type of reference element
using ReferenceElements = Dune::ReferenceElements<ct, 0>; using ReferenceElements = Dune::ReferenceElements<ct, 0>;
...@@ -486,17 +575,21 @@ namespace Dune ...@@ -486,17 +575,21 @@ namespace Dune
using GlobalCoordinate = FieldVector<ct, cdim>; using GlobalCoordinate = FieldVector<ct, cdim>;
protected: protected:
Geometry (const ReferenceElement& refElement, std::vector<GlobalCoordinate> vertices) VertexGeometry (Construct, const ReferenceElement& refElement, const LocalFunction& localFunction,
: Super(refElement, std::move(vertices)) const LocalGeometry& lg)
: Super(refElement, std::vector<GlobalCoordinate>{localFunction(lg.global(refElement.position(0,0)))})
{} {}
public: public:
Geometry (const ReferenceElement& refElement, const Localfunction& localFunction) template <class... Args>
: Geometry(refElement, std::vector<GlobalCoordinate>{localFunction(refElement.position(0,0))}) VertexGeometry (const ReferenceElement& refElement, const LocalFunction& lf, Args&&... args)
: VertexGeometry(Construct{}, refElement, lf, LocalGeometry(std::forward<Args>(args)...))
{} {}
Geometry (Dune::GeometryType gt, const Localfunction& localFunction) template <class... Args>
: Geometry(ReferenceElements::general(gt), localFunction) VertexGeometry (Dune::GeometryType gt, const LocalFunction& lf, Args&&... args)
: VertexGeometry(Construct{}, ReferenceElements::general(gt), lf,
LocalGeometry(std::forward<Args>(args)...))
{} {}
Std::optional<LocalCoordinate> checkedLocal (const GlobalCoordinate& globalCoord) const Std::optional<LocalCoordinate> checkedLocal (const GlobalCoordinate& globalCoord) const
...@@ -515,8 +608,28 @@ namespace Dune ...@@ -515,8 +608,28 @@ namespace Dune
return GlobalCoordinate(0); return GlobalCoordinate(0);
} }
private:
VertexGeometry const& flatGeometry () const { return *this; }
};
// Specialization for vertex geometries
template <class ct, int cdim, class LF, class LG, int order>
class Geometry<ct,0,cdim,LF,LG,order>
: public VertexGeometry<ct,cdim,LF,LG>
{
using Super = VertexGeometry<ct,cdim,LF,LG>;
public:
using Super::Super;
};
template <class ct, int cdim, class LF, class LG>
class Geometry<ct,0,cdim,LF,LG,-1>
: public VertexGeometry<ct,cdim,LF,LG>
{
using Super = VertexGeometry<ct,cdim,LF,LG>;
public: public:
Geometry const& flatGeometry () const { return *this; } using Super::Super;
}; };
#endif // DOXYGEN #endif // DOXYGEN
......
...@@ -4,7 +4,6 @@ ...@@ -4,7 +4,6 @@
#define DUNE_CURVED_SURFACE_GRID_GRIDFAMILY_HH #define DUNE_CURVED_SURFACE_GRID_GRIDFAMILY_HH
#include <dune/grid/common/grid.hh> #include <dune/grid/common/grid.hh>
#include <dune/curvedgeometry/curvedgeometry.hh>
#include <dune/curvedsurfacegrid/capabilities.hh> #include <dune/curvedsurfacegrid/capabilities.hh>
#include <dune/curvedsurfacegrid/entity.hh> #include <dune/curvedsurfacegrid/entity.hh>
#include <dune/curvedsurfacegrid/entityseed.hh> #include <dune/curvedsurfacegrid/entityseed.hh>
...@@ -90,26 +89,17 @@ namespace Dune ...@@ -90,26 +89,17 @@ namespace Dune
{ {
using LocalGeometry = typename HostGrid::template Codim<codim>::LocalGeometry; using LocalGeometry = typename HostGrid::template Codim<codim>::LocalGeometry;
using LocalTransformation = std::conditional_t<(codim == 0),
CGeo::Impl::DefaultLocalGeometry,
CGeo::Impl::LocalGeometryInterface<LocalGeometry>>;
template < int mydim, int cdim, class GridImpl > template < int mydim, int cdim, class GridImpl >
using GeometryImplTemplate = std::conditional_t<(differentiableLocalFunction || (order > 0)), using GeometryImplTemplate
CGeo::Geometry<ctype, mydim, cdim, LocalFunction, = CGeo::Geometry<ctype, mydim, cdim, LocalFunction, LocalTransformation, order>;
std::conditional_t<(codim == 0), CGeo::Impl::DefaultLocalGeometry, LocalGeometry>>,
Dune::LagrangeCurvedGeometry<ctype, mydim, cdim, (order > 0 ? order : 1)>>;
// geometry types