Skip to content
GitLab
Commit ba127239 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

extract LocalFunction Geometry to curvedgeometry project

parent 23317987
Expand all Hide whitespace changes
Inline Side-by-side
dune/curvedsurfacegrid/CMakeLists.txt
View file @ ba127239
install(FILES install(FILES
backuprestore.hh backuprestore.hh
capabilities.hh capabilities.hh
concepts.hh concepts.hh
...@@ -19,9 +19,9 @@ install(FILES ...@@ -19,9 +19,9 @@ install(FILES
intersection.hh intersection.hh
intersectioniterator.hh intersectioniterator.hh
iterator.hh iterator.hh
localgeometrywrapper.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/curvedsurfacegrid) DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/curvedsurfacegrid)
add_subdirectory(gridfunctions) add_subdirectory(gridfunctions)
add_subdirectory(surfacedistance) add_subdirectory(surfacedistance)
add_subdirectory(test) add_subdirectory(test)
\ No newline at end of file
dune/curvedsurfacegrid/entity.hh
View file @ ba127239
...@@ -264,6 +264,9 @@ namespace Dune ...@@ -264,6 +264,9 @@ namespace Dune
//! type of corresponding geometry //! type of corresponding geometry
using Geometry = typename Traits::template Codim<0>::Geometry; using Geometry = typename Traits::template Codim<0>::Geometry;
//! type of corresponding local geometry
using LocalGeometry = typename Traits::template Codim<0>::LocalGeometry;
public: public:
using Super::Super; using Super::Super;
...@@ -283,7 +286,7 @@ namespace Dune ...@@ -283,7 +286,7 @@ namespace Dune
if (!geo_) { if (!geo_) {
auto localFct = localFunction(Super::gridFunction()); auto localFct = localFunction(Super::gridFunction());
localFct.bind(Super::hostEntity()); localFct.bind(Super::hostEntity());
auto fakeDefaultGeometry = Dune::CGeo::Impl::DefaultLocalGeometry{}; auto fakeDefaultGeometry = Dune::DefaultLocalGeometry<typename Super::ctype, Super::mydimension, Super::mydimension>{};
geo_ = std::make_shared<GeometryImpl>(Super::type(), localFct, fakeDefaultGeometry); geo_ = std::make_shared<GeometryImpl>(Super::type(), localFct, fakeDefaultGeometry);
} }
......
dune/curvedsurfacegrid/geometry.hh
View file @ ba127239
...@@ -17,6 +17,7 @@ ...@@ -17,6 +17,7 @@
#include <dune/common/std/type_traits.hh> #include <dune/common/std/type_traits.hh>
#include <dune/curvedgeometry/curvedgeometry.hh> #include <dune/curvedgeometry/curvedgeometry.hh>
#include <dune/curvedgeometry/localfunctiongeometry.hh>
#include <dune/geometry/affinegeometry.hh> #include <dune/geometry/affinegeometry.hh>
#include <dune/geometry/multilineargeometry.hh> #include <dune/geometry/multilineargeometry.hh>
...@@ -24,132 +25,26 @@ ...@@ -24,132 +25,26 @@
#include <dune/geometry/referenceelements.hh> #include <dune/geometry/referenceelements.hh>
#include <dune/geometry/type.hh> #include <dune/geometry/type.hh>
#include "localgeometrywrapper.hh"
namespace Dune namespace Dune
{ {
namespace CGeo namespace CGeo
{ {
namespace Impl
{
struct IdentityMatrix {};
template <class K, int n, int m, int l>
FieldMatrix<K,n,m> ABt(const FieldMatrix<K,n,l>& A, const FieldMatrix<K,m,l>& B)
{
FieldMatrix<K,n,m> ABt;
for (int i = 0; i < n; ++i)
for (int j = 0; j < m; ++j) {
ABt[i][j] = 0;
for (int k = 0; k < l; ++k)
ABt[i][j] += A[i][k] * B[j][k];
}
return ABt;
}
template <class K, int n, int m>
FieldMatrix<K,n,m> ABt(const DiagonalMatrix<K,n>& A, const FieldMatrix<K,m,n>& B)
{
FieldMatrix<K,n,m> ABt;
for (int i = 0; i < n; ++i)
for (int j = 0; j < m; ++j) {
ABt[i][j] = A[i][i] * B[j][i];
}
return ABt;
}
template <class K, int n, int m, int l>
FieldMatrix<K,n,m> AB(const FieldMatrix<K,n,l>& A, const FieldMatrix<K,l,m>& B)
{
FieldMatrix<K,n,m> ABt;
for (int i = 0; i < n; ++i)
for (int j = 0; j < m; ++j) {
ABt[i][j] = 0;
for (int k = 0; k < l; ++k)
ABt[i][j] += A[i][k] * B[k][j];
}
return ABt;
}
template <class K, int n, int m>
FieldMatrix<K,n,m> AB(const DiagonalMatrix<K,n>& A, const FieldMatrix<K,n,m>& B)
{
FieldMatrix<K,n,m> ABt;
for (int i = 0; i < n; ++i)
for (int j = 0; j < m; ++j) {
ABt[i][j] = A[i][i] * B[i][j];
}
return ABt;
}
template <class K, int n, int m>
const FieldMatrix<K,n,m>& AB(const IdentityMatrix& A, const FieldMatrix<K,n,m>& B)
{
return B;
}
// efficient implementation of a trivial local geometry for elements
struct DefaultLocalGeometry
{
template <class LocalCoordinate>
decltype(auto) global(LocalCoordinate&& local) const
{
return std::forward<LocalCoordinate>(local);
}
template <class LocalCoordinate>
IdentityMatrix jacobianTransposed(LocalCoordinate&& local) const
{
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 a differentiable localFunction * Parametrization of the geometry by a differentiable localFunction
* *
* \tparam ct coordinate type
* \tparam mydim geometry dimension
* \tparam cdim coordinate dimension
* \tparam LF localFunction parametrizing the geometry * \tparam LF localFunction parametrizing the geometry
* \tparam LG localGeometry for coordinate transform from subEntity to element, * \tparam LG localGeometry for coordinate transform from subEntity to element,
* see \ref Impl::DefaultLocalGeometry and \ref Impl::LocalGeometryInterface * see \ref Impl::DefaultLocalGeometry and \ref Impl::LocalGeometryInterface
* \tparam order Polynomial order of lagrange parametrization. If order < 0 use localFunction. * \tparam ORDER Polynomial order of lagrange parametrization. If order < 0 use localFunction.
*/ */
template <class ct, int mydim, int cdim, class LF, class LG, int order> template <class ct, int mydim, int cdim, class LF, class LG, int ORDER = -1>
class Geometry class Geometry
: public Dune::LagrangeCurvedGeometry<ct, mydim, cdim, (order > 0 ? order : 1)> : public LagrangeCurvedGeometry<ct, mydim, cdim, (ORDER > 0 ? ORDER : 1)>
{ {
using Super = Dune::LagrangeCurvedGeometry<ct, mydim, cdim, (order > 0 ? order : 1)>; using Super = LagrangeCurvedGeometry<ct, mydim, cdim, (ORDER > 0 ? ORDER : 1)>;
//! type of the localFunction representation the geometry parametrization //! type of the localFunction representation the geometry parametrization
using LocalFunction = LF; using LocalFunction = LF;
...@@ -192,57 +87,22 @@ namespace Dune ...@@ -192,57 +87,22 @@ namespace Dune
}; };
// Specialization for the case that ORDER < 0: Use LocalFunction `LF` directly as parametrization
template <class ct, int mydim, int cdim, class LF, class LG> template <class ct, int mydim, int cdim, class LF, class LG>
class Geometry<ct,mydim,cdim,LF,LG,-1> class Geometry<ct,mydim,cdim,LF,LG,-1>
: public LocalFunctionGeometry<LF,LG>
{ {
public: using Super = LocalFunctionGeometry<LF,LG>;
/// coordinate type
using ctype = ct;
/// geometry dimension
static const int mydimension = mydim;
/// coordinate dimension
static const int coorddimension = cdim;
/// type of local coordinates
using LocalCoordinate = FieldVector<ctype, mydimension>;
/// type of global coordinates
using GlobalCoordinate = FieldVector<ctype, coorddimension>;
/// type of volume
using Volume = ctype;
/// type of jacobian transposed
using JacobianTransposed = FieldMatrix<ctype, mydimension, coorddimension>;
/// type of jacobian inverse transposed
using JacobianInverseTransposed = FieldMatrix<ctype, coorddimension, mydimension>;
public:
/// type of reference element
using ReferenceElements = Dune::ReferenceElements<ctype, mydimension>;
using ReferenceElement = typename ReferenceElements::ReferenceElement;
protected:
//! helper structure containing some matrix routines. See affinegeometry.hh
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 //! type of the localFunction representation the geometry parametrization
using LocalFunction = LF; using LocalFunction = LF;
//! type of the LocalFunction representing the derivative //! type of coordinate transformation for subEntities to codim=0 entities
using DerivativeLocalFunction = std::decay_t<decltype(derivative(std::declval<LF>()))>; using LocalGeometry = LG;
//! tolerance to numerical algorithms
static ct tolerance () { return ct(16) * std::numeric_limits<ct>::epsilon(); }
//! maximal number of Newton iteration in `geometry.local(global)` /// type of reference element
static int maxIteration () { return 100; } using ReferenceElements = Dune::ReferenceElements<ct, mydim>;
using ReferenceElement = typename ReferenceElements::ReferenceElement;
public: public:
/// \brief constructor from localFunction to parametrize the geometry /// \brief constructor from localFunction to parametrize the geometry
...@@ -255,9 +115,7 @@ namespace Dune ...@@ -255,9 +115,7 @@ namespace Dune
*/ */
template <class... Args> template <class... Args>
Geometry (const ReferenceElement& refElement, const LocalFunction& localFunction, Args&&... args) Geometry (const ReferenceElement& refElement, const LocalFunction& localFunction, Args&&... args)
: refElement_(refElement) : Super(refElement, localFunction, LocalGeometry{std::forward<Args>(args)...})
, localFunction_(localFunction)
, 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
...@@ -269,280 +127,8 @@ namespace Dune ...@@ -269,280 +127,8 @@ namespace Dune
*/ */
template <class... Args> template <class... Args>
Geometry (Dune::GeometryType gt, const LocalFunction& localFunction, Args&&... args) Geometry (Dune::GeometryType gt, const LocalFunction& localFunction, Args&&... args)
: Geometry(ReferenceElements::general(gt), localFunction, std::forward<Args>(args)...) : Super(gt, localFunction, LocalGeometry{std::forward<Args>(args)...})
{} {}
/// \brief is this mapping affine? No!
bool affine () const
{
return false;
}
/// \brief obtain the element type from the reference element
Dune::GeometryType type () const
{
return refElement_.type();
}
/// \brief obtain number of corners of the corresponding reference element
int corners () const
{
return refElement_.size(mydimension);
}
/// \brief obtain coordinates of the i-th corner
GlobalCoordinate corner (int i) const
{
assert( (i >= 0) && (i < corners()) );
return global(refElement_.position(i, mydimension));
}
/// \brief obtain the centroid of the mapping's image
GlobalCoordinate center () const
{
return global(refElement_.position(0, 0));
}
/// \brief Evaluate the coordinate mapping
/**
* Evaluate the local function in local coordinates
*
* \param[in] local local coordinates
*
* \returns corresponding global coordinate
*/
GlobalCoordinate global (const LocalCoordinate& local) const
{
return localFunction()(localGeometry().global(local));
}
/// \brief Evaluate the inverse coordinate mapping
/**
* \param[in] globalCoord global coordinate to map
*
* \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
* the following function over the local coordinate space spanned by the reference element.
* \code
* (global( x ) - y).two_norm()
* \endcode
*/
LocalCoordinate local (const GlobalCoordinate& globalCoord) const
{
auto localCoord = checkedLocal(globalCoord);
if (!localCoord)
DUNE_THROW(Dune::Exception,
"local coordinate can not be recovered from given global coordinate " << globalCoord);
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
{
LocalCoordinate x = flatGeometry().local(globalCoord);
LocalCoordinate dx;
for (int i = 0; i < maxIteration(); ++i)
{
// Newton's method: DF^n dx^n = F^n, x^{n+1} -= dx^n
const GlobalCoordinate dglobal = global(x) - globalCoord;
const bool invertible = MatrixHelper::xTRightInvA(jacobianTransposed(x), dglobal, dx);
// break if jacobian is not invertible
if (!invertible)
return {};
// update x with correction
x -= dx;
// break if tolerance is reached.
if (dx.two_norm2() < tolerance())
break;
}
if (dx.two_norm2() > tolerance())
return {};
return x;
}
/// \brief Construct a global coordinate normal of the curvilinear element evaluated at
/// a given local coordinate
/**
* \note Implemented for codim=1 entities only, i.e. edges in 2D and faces in 3D
**/
GlobalCoordinate normal (const LocalCoordinate& local) const
{
if ((mydim == 1) && (cdim == 2)) { return normal1D(local, jacobianTransposed(local)); }
if ((mydim == 2) && (cdim == 3)) { return normal2D(local, jacobianTransposed(local)); }
DUNE_THROW(Dune::NotImplemented,
"ERROR: normal() method only defined for edges in 2D and faces in 3D");
return GlobalCoordinate(0);
}
/// \brief Obtain the integration element
/**
* If the Jacobian of the mapping is denoted by $J(x)$, the integration
* integration element \f$\mu(x)\f$ is given by
* \f[ \mu(x) = \sqrt{|\det (J^T(x) J(x))|}.\f]
*
* \param[in] local local coordinate to evaluate the integration element in
*
* \returns the integration element \f$\mu(x)\f$.
*/
ctype integrationElement (const LocalCoordinate& local) const
{
return MatrixHelper::sqrtDetAAT(jacobianTransposed(local));
}
/// \brief Obtain the volume of the mapping's image
/**
* Calculates the volume of the entity by numerical integration
*
* \todo Check the quadrature order of the volume integral.
*/
Volume volume () const
{
const int p = 10; // TODO: use adaptive quadrature rule
const auto& quadRule = QuadratureRules<ctype, mydimension>::rule(type(), p);
Volume vol(0);
for (auto const& qp : quadRule)
vol += integrationElement(qp.position()) * qp.weight();
return vol;
}
/// \brief Obtain the transposed of the Jacobian
/**
* \param[in] local local coordinate to evaluate Jacobian in
*
* \returns the matrix corresponding to the transposed of the Jacobian
*/
JacobianTransposed jacobianTransposed (const LocalCoordinate& local) const
{
if (!derivativeLocalFunction_) {
derivativeLocalFunction_.emplace(derivative(localFunction()));
derivativeLocalFunction_->bind(localFunction().localContext());
}
// coordinate in the localContext of the localFunction
auto&& elementLocal = localGeometry().global(local);
auto jLocal = localGeometry().jacobianTransposed(local);
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
/**
* The Jacobian's inverse is defined as a pseudo-inverse. If we denote
* the Jacobian by \f$J(x)\f$, the following condition holds:
* \f[J^{-1}(x) J(x) = I.\f]
*/
JacobianInverseTransposed jacobianInverseTransposed (const LocalCoordinate& local) const
{
JacobianInverseTransposed out;
MatrixHelper::rightInvA(jacobianTransposed(local), out);
return out;
}
/// \brief obtain the reference-element related to this geometry
friend ReferenceElement referenceElement (const Geometry& geometry)
{
return geometry.refElement();
}
protected:
// the internal stored reference element
const ReferenceElement& refElement () const
{
return refElement_;
}
// normal vector to an edge line-element
GlobalCoordinate normal1D (const LocalCoordinate& /*local*/, const JacobianTransposed& J) const
{
GlobalCoordinate res{
J[0][1],
-J[0][0]};
return res /= res.two_norm();
}
// normal vector to a triangle or quad face element
GlobalCoordinate normal2D (const LocalCoordinate& /*local*/, const JacobianTransposed& J) const
{
GlobalCoordinate res{
J[0][1] * J[1][2] - J[0][2] * J[1][1],
J[0][2] * J[1][0] - J[0][0] * J[1][2],
J[0][0] * J[1][1] - J[0][1] * J[1][0]};
return res /= res.two_norm();
}
private:
const LocalFunction& localFunction () const
{
return localFunction_;
}
const LocalGeometry& localGeometry () const
{
return localGeometry_;
}
using FlatGeometry = MultiLinearGeometry<ctype, mydim, cdim>;
const FlatGeometry& flatGeometry () const
{
if (!flatGeometry_) {
std::vector<GlobalCoordinate> corners;
corners.reserve(refElement_.size(mydimension));
for (int i = 0; i < refElement_.size(mydimension); ++i)
corners.push_back(global(refElement_.position(i, mydimension)));
flatGeometry_.emplace(refElement_, std::move(corners));
}
return *flatGeometry_;
}
using HostGeometry = std::decay_t<decltype(std::declval<LF>().localContext().geometry())>;
const HostGeometry& hostGeometry () const
{
if (!hostGeometry_)
hostGeometry_.emplace(localFunction().localContext().geometry());
return *hostGeometry_;
}
private:
//! Reference of the geometry
ReferenceElement refElement_;
//! local parametrization of the element
LocalFunction localFunction_;
//! transformation of local coordinates to element-local coordinates
LocalGeometry localGeometry_;
// some data optionally provided
mutable Std::optional<DerivativeLocalFunction> derivativeLocalFunction_;
mutable Std::optional<FlatGeometry> flatGeometry_;
mutable Std::optional<HostGeometry> hostGeometry_;
}; };
......
dune/curvedsurfacegrid/gridfamily.hh
View file @ ba127239
...@@ -16,6 +16,7 @@ ...@@ -16,6 +16,7 @@
#include <dune/curvedsurfacegrid/iterator.hh> #include <dune/curvedsurfacegrid/iterator.hh>
#include <dune/curvedsurfacegrid/idset.hh> #include <dune/curvedsurfacegrid/idset.hh>
#include <dune/curvedsurfacegrid/indexsets.hh> #include <dune/curvedsurfacegrid/indexsets.hh>
#include <dune/curvedsurfacegrid/localgeometrywrapper.hh>
#include <dune/curvedsurfacegrid/gridfunctions/gridfunction.hh> #include <dune/curvedsurfacegrid/gridfunctions/gridfunction.hh>