5#ifndef DUNE_GEOMETRY_LOCALFINITEELEMENTGEOMETRY_HH
6#define DUNE_GEOMETRY_LOCALFINITEELEMENTGEOMETRY_HH
14#include <dune/common/fmatrix.hh>
15#include <dune/common/fvector.hh>
16#include <dune/common/math.hh>
17#include <dune/common/typetraits.hh>
18#include <dune/common/std/type_traits.hh>
38template <
class LFE,
int cdim>
41 using LocalFiniteElement = LFE;
42 using LocalBasis =
typename LFE::Traits::LocalBasisType;
43 using LocalBasisTraits =
typename LocalBasis::Traits;
47 using ctype =
typename LocalBasisTraits::DomainFieldType;
65 using Jacobian = FieldMatrix<ctype, coorddimension, mydimension>;
104 const LocalFiniteElement& localFE,
105 std::vector<GlobalCoordinate> vertices)
106 : refElement_(refElement)
108 , vertices_(
std::move(vertices))
110 assert(localFE_.size() == vertices_.size());
126 template <
class Param,
127 std::enable_if_t<std::is_invocable_r_v<GlobalCoordinate,Param,LocalCoordinate>,
int> = 0>
129 const LocalFiniteElement& localFE,
130 Param&& parametrization)
131 : refElement_(refElement)
134 localFE_.localInterpolation().interpolate(parametrization, vertices_);
143 template <
class... Args>
162 affine_.emplace(affineImpl());
169 return refElement_.type();
181 assert( (i >= 0) && (i <
corners()) );
188 return global(refElement_.position(0, 0));
204 thread_local std::vector<typename LocalBasisTraits::RangeType> shapeValues;
206 assert(shapeValues.size() == vertices_.size());
209 for (std::size_t i = 0; i < shapeValues.size(); ++i)
210 out.axpy(shapeValues[i], vertices_[i]);
235 Impl::GaussNewtonErrorCode err = Impl::gaussNewton(
241 if (err != Impl::GaussNewtonErrorCode::OK)
242 DUNE_THROW(Dune::Exception,
243 "Local coordinate can not be recovered from global coordinate, error code = " <<
int(err) <<
"\n"
244 <<
" (global(x) - y).two_norm() = " << (
global(x) - y).two_norm()
245 <<
" > tol = " << opts.absTol);
283 for (
int p = 2; p < opts.maxIt; ++p) {
285 if (abs(vol1 - vol0) < opts.absTol)
297 for (
const auto& qp : quadRule)
309 thread_local std::vector<typename LocalBasisTraits::JacobianType> shapeJacobians;
311 assert(shapeJacobians.size() == vertices_.size());
314 for (std::size_t i = 0; i < shapeJacobians.size(); ++i) {
315 for (
int j = 0; j < Jacobian::rows; ++j) {
316 shapeJacobians[i].umtv(vertices_[i][j], out[j]);
361 return geometry.refElement_;
379 return localFE_.localBasis();
384 bool affineImpl ()
const
405 auto refSimplex = referenceElement<ctype,mydimension>(GeometryTypes::simplex(
mydimension));
406 auto simplexIndices = Dune::range(refSimplex.size(
mydimension));
407 auto simplexCorners = Dune::transformedRangeView(simplexIndices,
409 AffineGeometry<ctype,mydimension,coorddimension> affineGeo(refSimplex,simplexCorners);
411 const ctype tol = sqrt(std::numeric_limits<ctype>::epsilon());
412 for (
int i = 0; i <
corners(); ++i) {
427 LocalFiniteElement localFE_{};
430 std::vector<GlobalCoordinate> vertices_{};
432 mutable std::optional<bool> affine_ = std::nullopt;
439using LocalCoordinate_t
440 = FieldVector<
typename LFE::Traits::LocalBasisType::Traits::DomainFieldType,
441 LFE::Traits::LocalBasisType::Traits::dimDomain>;
447template <
class I,
class LFE,
class GlobalCoordinate>
451template <
class I,
class LFE,
class F,
452 class Range = std::invoke_result_t<F,Impl::LocalCoordinate_t<LFE>>>
456template <
class LFE,
class GlobalCoordinate>
460template <
class LFE,
class F,
461 class Range = std::invoke_result_t<F,Impl::LocalCoordinate_t<LFE>>>
An implementation of the Geometry interface for affine geometries.
A unique label for each type of element that can occur in a grid.
Definition affinegeometry.hh:21
This class provides access to geometric and topological properties of a reference element.
Definition referenceelement.hh:52
Class providing access to the singletons of the reference elements.
Definition referenceelements.hh:128
typename Container::ReferenceElement ReferenceElement
The reference element type.
Definition referenceelements.hh:146
Geometry implementation based on local-basis function parametrization.
Definition localfiniteelementgeometry.hh:40
GeometryType type() const
Obtain the name of the reference element.
Definition localfiniteelementgeometry.hh:167
decltype(power(std::declval< ctype >(), mydimension)) Volume
type of volume
Definition localfiniteelementgeometry.hh:62
LocalCoordinate local(const GlobalCoordinate &y, Impl::GaussNewtonOptions< ctype > opts={}) const
Evaluate the inverse coordinate mapping.
Definition localfiniteelementgeometry.hh:232
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
JacobianTransposed jacobianTransposed(const LocalCoordinate &local) const
Obtain the transposed of the Jacobian.
Definition localfiniteelementgeometry.hh:327
const LocalBasis & localBasis() const
The local basis of the stored local finite-element.
Definition localfiniteelementgeometry.hh:377
typename LocalBasisTraits::DomainFieldType ctype
coordinate type
Definition localfiniteelementgeometry.hh:47
JacobianInverseTransposed jacobianInverseTransposed(const LocalCoordinate &local) const
Obtain the transposed of the Jacobian's inverse.
Definition localfiniteelementgeometry.hh:353
friend ReferenceElement referenceElement(const LocalFiniteElementGeometry &geometry)
Obtain the reference-element related to this geometry.
Definition localfiniteelementgeometry.hh:359
LocalFiniteElementGeometry()=default
Default constructed geometry results in an empty/invalid representation.
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
Jacobian jacobian(const LocalCoordinate &local) const
Obtain the Jacobian.
Definition localfiniteelementgeometry.hh:307
Volume volume(const QuadratureRule< ctype, mydimension > &quadRule) const
Obtain the volume of the mapping's image by given quadrature rules.
Definition localfiniteelementgeometry.hh:294
const std::vector< GlobalCoordinate > & coefficients() const
Obtain the coefficients of the parametrization.
Definition localfiniteelementgeometry.hh:371
LocalFiniteElementGeometry(GeometryType gt, Args &&... args)
Constructor, forwarding to the other constructors that take a reference-element.
Definition localfiniteelementgeometry.hh:144
static const int coorddimension
coordinate dimension
Definition localfiniteelementgeometry.hh:53
int corners() const
Obtain number of corners of the corresponding reference element.
Definition localfiniteelementgeometry.hh:173
GlobalCoordinate center() const
Obtain the centroid of the mapping's image.
Definition localfiniteelementgeometry.hh:186
const LocalFiniteElement & finiteElement() const
Obtain the local finite-element.
Definition localfiniteelementgeometry.hh:365
GlobalCoordinate global(const LocalCoordinate &local) const
Evaluate the coordinate mapping.
Definition localfiniteelementgeometry.hh:202
Impl::FieldMatrixHelper< ctype > MatrixHelper
Definition localfiniteelementgeometry.hh:82
typename ReferenceElements::ReferenceElement ReferenceElement
Definition localfiniteelementgeometry.hh:79
GlobalCoordinate corner(int i) const
Obtain coordinates of the i-th corner.
Definition localfiniteelementgeometry.hh:179
Volume volume(Impl::ConvergenceOptions< ctype > opts={}) const
Obtain the volume of the mapping's image.
Definition localfiniteelementgeometry.hh:276
static const int mydimension
geometry dimension
Definition localfiniteelementgeometry.hh:50
JacobianInverse jacobianInverse(const LocalCoordinate &local) const
Obtain the Jacobian's inverse.
Definition localfiniteelementgeometry.hh:339
FieldVector< ctype, coorddimension > GlobalCoordinate
type of global coordinates
Definition localfiniteelementgeometry.hh:59
FieldMatrix< ctype, coorddimension, mydimension > Jacobian
type of jacobian
Definition localfiniteelementgeometry.hh:65
FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed
type of jacobian inverse transposed
Definition localfiniteelementgeometry.hh:74
FieldMatrix< ctype, mydimension, coorddimension > JacobianInverse
type of jacobian inverse
Definition localfiniteelementgeometry.hh:71
ctype integrationElement(const LocalCoordinate &local) const
Obtain the integration element.
Definition localfiniteelementgeometry.hh:260
int order() const
Obtain the (highest) polynomial order of the parametrization.
Definition localfiniteelementgeometry.hh:149
FieldVector< ctype, mydimension > LocalCoordinate
type of local coordinates
Definition localfiniteelementgeometry.hh:56
FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed
type of jacobian transposed
Definition localfiniteelementgeometry.hh:68
LocalFiniteElementGeometry(const ReferenceElement &refElement, const LocalFiniteElement &localFE, Param &¶metrization)
Constructor from a local parametrization function, mapping local to (curved) global coordinates.
Definition localfiniteelementgeometry.hh:128
Abstract base class for quadrature rules.
Definition quadraturerules.hh:214
A container for all quadrature rules of dimension dim
Definition quadraturerules.hh:260
Unique label for each type of entities that can occur in DUNE grids.
Definition type.hh:114
constexpr bool isSimplex() const
Return true if entity is a simplex of any dimension.
Definition type.hh:319