Loading...
Searching...
No Matches
localfiniteelementgeometry.hh
Go to the documentation of this file.
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
#ifndef DUNE_GEOMETRY_LOCALFINITEELEMENTGEOMETRY_HH
#define DUNE_GEOMETRY_LOCALFINITEELEMENTGEOMETRY_HH
#include <cassert>
#include <functional>
#include <limits>
#include <type_traits>
#include <vector>
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/common/math.hh>
#include <dune/common/typetraits.hh>
#include <dune/common/std/type_traits.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/geometry/referenceelements.hh>
#include <dune/geometry/type.hh>
#include <dune/geometry/utility/algorithms.hh>
#include <dune/geometry/utility/convergence.hh>
#include <dune/geometry/utility/defaultmatrixhelper.hh>
namespace Dune {
template <class LFE, int cdim>
class LocalFiniteElementGeometry
{
using LocalFiniteElement = LFE;
using LocalBasis = typename LFE::Traits::LocalBasisType;
using LocalBasisTraits = typename LocalBasis::Traits;
public:
using ctype = typename LocalBasisTraits::DomainFieldType;
static const int mydimension = LocalBasisTraits::dimDomain;
static const int coorddimension = cdim;
using LocalCoordinate = FieldVector<ctype, mydimension>;
using GlobalCoordinate = FieldVector<ctype, coorddimension>;
using Volume = decltype(power(std::declval<ctype>(),mydimension));
using Jacobian = FieldMatrix<ctype, coorddimension, mydimension>;
using JacobianTransposed = FieldMatrix<ctype, mydimension, coorddimension>;
using JacobianInverse = FieldMatrix<ctype, mydimension, coorddimension>;
using JacobianInverseTransposed = FieldMatrix<ctype, coorddimension, mydimension>;
public:
using ReferenceElements = Dune::ReferenceElements<ctype, mydimension>;
using ReferenceElement = typename ReferenceElements::ReferenceElement;
protected:
using MatrixHelper = Impl::FieldMatrixHelper<ctype>;
public:
LocalFiniteElementGeometry () = default;
LocalFiniteElementGeometry (const ReferenceElement& refElement,
const LocalFiniteElement& localFE,
std::vector<GlobalCoordinate> vertices)
: refElement_(refElement)
, localFE_(localFE)
, vertices_(std::move(vertices))
{
assert(localFE_.size() == vertices_.size());
}
template <class Param,
std::enable_if_t<std::is_invocable_r_v<GlobalCoordinate,Param,LocalCoordinate>, int> = 0>
LocalFiniteElementGeometry (const ReferenceElement& refElement,
const LocalFiniteElement& localFE,
Param&& parametrization)
: refElement_(refElement)
, localFE_(localFE)
{
localFE_.localInterpolation().interpolate(parametrization, vertices_);
}
template <class... Args>
explicit LocalFiniteElementGeometry (GeometryType gt, Args&&... args)
: LocalFiniteElementGeometry(ReferenceElements::general(gt), std::forward<Args>(args)...)
{}
int order () const
{
return localBasis().order();
}
bool affine () const
{
if (!affine_)
affine_.emplace(affineImpl());
return *affine_;
}
GeometryType type () const
{
return refElement_.type();
}
int corners () const
{
return refElement_.size(mydimension);
}
GlobalCoordinate corner (int i) const
{
assert( (i >= 0) && (i < corners()) );
return global(refElement_.position(i, mydimension));
}
GlobalCoordinate center () const
{
return global(refElement_.position(0, 0));
}
GlobalCoordinate global (const LocalCoordinate& local) const
{
thread_local std::vector<typename LocalBasisTraits::RangeType> shapeValues;
localBasis().evaluateFunction(local, shapeValues);
assert(shapeValues.size() == vertices_.size());
GlobalCoordinate out(0);
for (std::size_t i = 0; i < shapeValues.size(); ++i)
out.axpy(shapeValues[i], vertices_[i]);
return out;
}
LocalCoordinate local (const GlobalCoordinate& y, Impl::GaussNewtonOptions<ctype> opts = {}) const
{
LocalCoordinate x = refElement_.position(0,0);
Impl::GaussNewtonErrorCode err = Impl::gaussNewton(
[&](const LocalCoordinate& local) { return this->global(local); },
[&](const LocalCoordinate& local) { return this->jacobianTransposed(local); },
y, x, opts
);
if (err != Impl::GaussNewtonErrorCode::OK)
DUNE_THROW(Dune::Exception,
"Local coordinate can not be recovered from global coordinate, error code = " << int(err) << "\n"
<< " (global(x) - y).two_norm() = " << (global(x) - y).two_norm()
<< " > tol = " << opts.absTol);
return x;
}
ctype integrationElement (const LocalCoordinate& local) const
{
return MatrixHelper::sqrtDetAAT(jacobianTransposed(local));
}
Volume volume (Impl::ConvergenceOptions<ctype> opts = {}) const
{
Volume vol0 = volume(QuadratureRules<ctype, mydimension>::rule(type(), 1));
if (affine())
return vol0;
using std::abs;
for (int p = 2; p < opts.maxIt; ++p) {
Volume vol1 = volume(QuadratureRules<ctype, mydimension>::rule(type(), p));
if (abs(vol1 - vol0) < opts.absTol)
return vol1;
vol0 = vol1;
}
return vol0;
}
Volume volume (const QuadratureRule<ctype, mydimension>& quadRule) const
{
Volume vol(0);
for (const auto& qp : quadRule)
vol += integrationElement(qp.position()) * qp.weight();
return vol;
}
Jacobian jacobian (const LocalCoordinate& local) const
{
thread_local std::vector<typename LocalBasisTraits::JacobianType> shapeJacobians;
localBasis().evaluateJacobian(local, shapeJacobians);
assert(shapeJacobians.size() == vertices_.size());
Jacobian out(0);
for (std::size_t i = 0; i < shapeJacobians.size(); ++i) {
for (int j = 0; j < Jacobian::rows; ++j) {
shapeJacobians[i].umtv(vertices_[i][j], out[j]);
}
}
return out;
}
JacobianTransposed jacobianTransposed (const LocalCoordinate& local) const
{
return jacobian(local).transposed();
}
JacobianInverse jacobianInverse (const LocalCoordinate& local) const
{
JacobianInverse out;
MatrixHelper::leftInvA(jacobian(local), out);
return out;
}
JacobianInverseTransposed jacobianInverseTransposed (const LocalCoordinate& local) const
{
return jacobianInverse(local).transposed();
}
friend ReferenceElement referenceElement (const LocalFiniteElementGeometry& geometry)
{
return geometry.refElement_;
}
const LocalFiniteElement& finiteElement () const
{
return localFE_;
}
const std::vector<GlobalCoordinate>& coefficients () const
{
return vertices_;
}
const LocalBasis& localBasis () const
{
return localFE_.localBasis();
}
private:
bool affineImpl () const
{
if constexpr(mydimension == 0)
// point geometries are always affine mappings
return true;
else {
if (order() > 1)
// higher-order parametrizations are by definition not affine
return false;
if constexpr(mydimension == 1)
// linear line geometries are affine mappings
return true;
else {
if (type().isSimplex())
// linear simplex geometries are affine mappings
return true;
// multi-linear mappings on non-simplex geometries might be affine
// as well. This is tested explicitly for all vertices by constructing
// an affine mapping from dim+1 affine-independent corners and evaluating
// at the other corners.
auto refSimplex = referenceElement<ctype,mydimension>(GeometryTypes::simplex(mydimension));
auto simplexIndices = Dune::range(refSimplex.size(mydimension));
auto simplexCorners = Dune::transformedRangeView(simplexIndices,
[&](int i) { return this->global(refSimplex.position(i,mydimension)); });
AffineGeometry<ctype,mydimension,coorddimension> affineGeo(refSimplex,simplexCorners);
using std::sqrt;
const ctype tol = sqrt(std::numeric_limits<ctype>::epsilon());
for (int i = 0; i < corners(); ++i) {
const auto corner = refElement_.position(i,mydimension);
if ((affineGeo.global(corner) - global(corner)).two_norm() > tol)
return false;
}
return true;
}
}
}
private:
ReferenceElement refElement_{};
LocalFiniteElement localFE_{};
std::vector<GlobalCoordinate> vertices_{};
mutable std::optional<bool> affine_ = std::nullopt;
};
namespace Impl {
// extract the LocalCoordinate type from a LocalFiniteElement
template <class LFE>
using LocalCoordinate_t
= FieldVector<typename LFE::Traits::LocalBasisType::Traits::DomainFieldType,
LFE::Traits::LocalBasisType::Traits::dimDomain>;
} // end namespace Impl
// deduction guides
template <class I, class LFE, class GlobalCoordinate>
LocalFiniteElementGeometry (Geo::ReferenceElement<I>, const LFE&, std::vector<GlobalCoordinate>)
-> LocalFiniteElementGeometry<LFE, GlobalCoordinate::dimension>;
template <class I, class LFE, class F,
class Range = std::invoke_result_t<F,Impl::LocalCoordinate_t<LFE>>>
LocalFiniteElementGeometry (Geo::ReferenceElement<I>, const LFE&, const F&)
-> LocalFiniteElementGeometry<LFE, Range::dimension>;
template <class LFE, class GlobalCoordinate>
LocalFiniteElementGeometry (GeometryType, const LFE& localFE, std::vector<GlobalCoordinate>)
-> LocalFiniteElementGeometry<LFE, GlobalCoordinate::dimension>;
template <class LFE, class F,
class Range = std::invoke_result_t<F,Impl::LocalCoordinate_t<LFE>>>
LocalFiniteElementGeometry (GeometryType, const LFE&, const F&)
-> LocalFiniteElementGeometry<LFE, Range::dimension>;
} // namespace Dune
#endif // DUNE_GEOMETRY_LOCALFINITEELEMENTGEOMETRY_HH
referenceelements.hh
convergence.hh
defaultmatrixhelper.hh
algorithms.hh
type.hh
A unique label for each type of element that can occur in a grid.
quadraturerules.hh
Dune::GeometryTypes::simplex
constexpr GeometryType simplex(unsigned int dim)
Returns a GeometryType representing a simplex of dimension dim.
Definition type.hh:453
std
STL namespace.
Dune
Definition affinegeometry.hh:22
Dune::LocalFiniteElementGeometry
LocalFiniteElementGeometry(Geo::ReferenceElement< I >, const LFE &, std::vector< GlobalCoordinate >) -> LocalFiniteElementGeometry< LFE, GlobalCoordinate::dimension >
Dune::Geo::ReferenceElement
This class provides access to geometric and topological properties of a reference element.
Definition referenceelement.hh:52
Dune::Geo::ReferenceElements
Class providing access to the singletons of the reference elements.
Definition referenceelements.hh:128
Dune::Geo::ReferenceElements::ReferenceElement
typename Container::ReferenceElement ReferenceElement
The reference element type.
Definition referenceelements.hh:146
Dune::LocalFiniteElementGeometry
Geometry implementation based on local-basis function parametrization.
Definition localfiniteelementgeometry.hh:40
Dune::LocalFiniteElementGeometry::type
GeometryType type() const
Obtain the name of the reference element.
Definition localfiniteelementgeometry.hh:167
Dune::LocalFiniteElementGeometry::Volume
decltype(power(std::declval< ctype >(), mydimension)) Volume
type of volume
Definition localfiniteelementgeometry.hh:62
Dune::LocalFiniteElementGeometry::local
LocalCoordinate local(const GlobalCoordinate &y, Impl::GaussNewtonOptions< ctype > opts={}) const
Evaluate the inverse coordinate mapping.
Definition localfiniteelementgeometry.hh:232
Dune::LocalFiniteElementGeometry::LocalFiniteElementGeometry
LocalFiniteElementGeometry(const ReferenceElement &refElement, const LocalFiniteElement &localFE, std::vector< GlobalCoordinate > vertices)
Constructor from a vector of coefficients of the LocalBasis parametrizing the geometry.
Definition localfiniteelementgeometry.hh:103
Dune::LocalFiniteElementGeometry::jacobianTransposed
JacobianTransposed jacobianTransposed(const LocalCoordinate &local) const
Obtain the transposed of the Jacobian.
Definition localfiniteelementgeometry.hh:327
Dune::LocalFiniteElementGeometry::localBasis
const LocalBasis & localBasis() const
The local basis of the stored local finite-element.
Definition localfiniteelementgeometry.hh:377
Dune::LocalFiniteElementGeometry::ctype
typename LocalBasisTraits::DomainFieldType ctype
coordinate type
Definition localfiniteelementgeometry.hh:47
Dune::LocalFiniteElementGeometry::jacobianInverseTransposed
JacobianInverseTransposed jacobianInverseTransposed(const LocalCoordinate &local) const
Obtain the transposed of the Jacobian's inverse.
Definition localfiniteelementgeometry.hh:353
Dune::LocalFiniteElementGeometry::referenceElement
friend ReferenceElement referenceElement(const LocalFiniteElementGeometry &geometry)
Obtain the reference-element related to this geometry.
Definition localfiniteelementgeometry.hh:359
Dune::LocalFiniteElementGeometry::LocalFiniteElementGeometry
LocalFiniteElementGeometry()=default
Default constructed geometry results in an empty/invalid representation.
Dune::LocalFiniteElementGeometry::affine
bool affine() const
Is this mapping affine? Geometries of order 1 might be affine, but it needs to be checked on non-simp...
Definition localfiniteelementgeometry.hh:159
Dune::LocalFiniteElementGeometry::jacobian
Jacobian jacobian(const LocalCoordinate &local) const
Obtain the Jacobian.
Definition localfiniteelementgeometry.hh:307
Dune::LocalFiniteElementGeometry::volume
Volume volume(const QuadratureRule< ctype, mydimension > &quadRule) const
Obtain the volume of the mapping's image by given quadrature rules.
Definition localfiniteelementgeometry.hh:294
Dune::LocalFiniteElementGeometry::coefficients
const std::vector< GlobalCoordinate > & coefficients() const
Obtain the coefficients of the parametrization.
Definition localfiniteelementgeometry.hh:371
Dune::LocalFiniteElementGeometry::LocalFiniteElementGeometry
LocalFiniteElementGeometry(GeometryType gt, Args &&... args)
Constructor, forwarding to the other constructors that take a reference-element.
Definition localfiniteelementgeometry.hh:144
Dune::LocalFiniteElementGeometry::coorddimension
static const int coorddimension
coordinate dimension
Definition localfiniteelementgeometry.hh:53
Dune::LocalFiniteElementGeometry::corners
int corners() const
Obtain number of corners of the corresponding reference element.
Definition localfiniteelementgeometry.hh:173
Dune::LocalFiniteElementGeometry::center
GlobalCoordinate center() const
Obtain the centroid of the mapping's image.
Definition localfiniteelementgeometry.hh:186
Dune::LocalFiniteElementGeometry::finiteElement
const LocalFiniteElement & finiteElement() const
Obtain the local finite-element.
Definition localfiniteelementgeometry.hh:365
Dune::LocalFiniteElementGeometry::global
GlobalCoordinate global(const LocalCoordinate &local) const
Evaluate the coordinate mapping.
Definition localfiniteelementgeometry.hh:202
Dune::LocalFiniteElementGeometry::MatrixHelper
Impl::FieldMatrixHelper< ctype > MatrixHelper
Definition localfiniteelementgeometry.hh:82
Dune::LocalFiniteElementGeometry::ReferenceElement
typename ReferenceElements::ReferenceElement ReferenceElement
Definition localfiniteelementgeometry.hh:79
Dune::LocalFiniteElementGeometry::corner
GlobalCoordinate corner(int i) const
Obtain coordinates of the i-th corner.
Definition localfiniteelementgeometry.hh:179
Dune::LocalFiniteElementGeometry::volume
Volume volume(Impl::ConvergenceOptions< ctype > opts={}) const
Obtain the volume of the mapping's image.
Definition localfiniteelementgeometry.hh:276
Dune::LocalFiniteElementGeometry::mydimension
static const int mydimension
geometry dimension
Definition localfiniteelementgeometry.hh:50
Dune::LocalFiniteElementGeometry::jacobianInverse
JacobianInverse jacobianInverse(const LocalCoordinate &local) const
Obtain the Jacobian's inverse.
Definition localfiniteelementgeometry.hh:339
Dune::LocalFiniteElementGeometry::GlobalCoordinate
FieldVector< ctype, coorddimension > GlobalCoordinate
type of global coordinates
Definition localfiniteelementgeometry.hh:59
Dune::LocalFiniteElementGeometry::Jacobian
FieldMatrix< ctype, coorddimension, mydimension > Jacobian
type of jacobian
Definition localfiniteelementgeometry.hh:65
Dune::LocalFiniteElementGeometry::JacobianInverseTransposed
FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed
type of jacobian inverse transposed
Definition localfiniteelementgeometry.hh:74
Dune::LocalFiniteElementGeometry::ReferenceElements
Dune::ReferenceElements< ctype, mydimension > ReferenceElements
type of reference element
Definition localfiniteelementgeometry.hh:78
Dune::LocalFiniteElementGeometry::JacobianInverse
FieldMatrix< ctype, mydimension, coorddimension > JacobianInverse
type of jacobian inverse
Definition localfiniteelementgeometry.hh:71
Dune::LocalFiniteElementGeometry::integrationElement
ctype integrationElement(const LocalCoordinate &local) const
Obtain the integration element.
Definition localfiniteelementgeometry.hh:260
Dune::LocalFiniteElementGeometry::order
int order() const
Obtain the (highest) polynomial order of the parametrization.
Definition localfiniteelementgeometry.hh:149
Dune::LocalFiniteElementGeometry::LocalCoordinate
FieldVector< ctype, mydimension > LocalCoordinate
type of local coordinates
Definition localfiniteelementgeometry.hh:56
Dune::LocalFiniteElementGeometry::JacobianTransposed
FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed
type of jacobian transposed
Definition localfiniteelementgeometry.hh:68
Dune::LocalFiniteElementGeometry::LocalFiniteElementGeometry
LocalFiniteElementGeometry(const ReferenceElement &refElement, const LocalFiniteElement &localFE, Param &&parametrization)
Constructor from a local parametrization function, mapping local to (curved) global coordinates.
Definition localfiniteelementgeometry.hh:128
Dune::QuadratureRule
Abstract base class for quadrature rules.
Definition quadraturerules.hh:214
Dune::QuadratureRules::rule
static const QuadratureRule & rule(const GeometryType &t, int p, QuadratureType::Enum qt=QuadratureType::GaussLegendre)
select the appropriate QuadratureRule for GeometryType t and order p
Definition quadraturerules.hh:326
Dune::GeometryType
Unique label for each type of entities that can occur in DUNE grids.
Definition type.hh:114