Loading...
Searching...
No Matches
mappedgeometry.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_MAPPEDGEOMETRY_HH
#define DUNE_GEOMETRY_MAPPEDGEOMETRY_HH
#include <cassert>
#include <limits>
#include <optional>
#include <stdexcept>
#include <type_traits>
#include <dune/common/copyableoptional.hh>
#include <dune/common/exceptions.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/common/math.hh>
#include <dune/common/transpose.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 Map, class Geo>
class MappedGeometry
{
public:
using LocalCoordinate = typename Geo::LocalCoordinate;
using GlobalCoordinate = std::remove_reference_t<decltype(std::declval<Map>()(std::declval<typename Geo::GlobalCoordinate>()))>;
using ctype = typename Geo::ctype;
static constexpr int mydimension = LocalCoordinate::size();
static constexpr int coorddimension = GlobalCoordinate::size();
using Volume = std::remove_reference_t<decltype(Dune::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>;
private:
using ReferenceElements = Dune::ReferenceElements<ctype, mydimension>;
using ReferenceElement = typename ReferenceElements::ReferenceElement;
protected:
using MatrixHelper = Impl::FieldMatrixHelper<ctype>;
// type of the mapping representation the geometry parametrization
using Mapping = Map;
// type of the geometry that is wrapped
using Geometry = Geo;
// type of a mapping representing the derivative of `Map` w.r.t. `GlobalCoordinate`
using DerivativeMapping = std::remove_reference_t<decltype(derivative(std::declval<Map>()))>;
public:
template <class Geo_, class Map_,
std::enable_if_t<Dune::IsInteroperable<Map, Map_>::value, int> = 0,
std::enable_if_t<Dune::IsInteroperable<Geo, Geo_>::value, int> = 0>
MappedGeometry (Map_&& mapping, Geo_&& geometry, bool affine = false)
: mapping_(std::forward<Map_>(mapping))
, dMapping_(derivative(*mapping_))
, geometry_(std::forward<Geo_>(geometry))
, affine_(affine)
{}
bool affine () const
{
return affine_;
}
GeometryType type () const
{
return geometry_.type();
}
int corners () const
{
return geometry_.corners();
}
GlobalCoordinate corner (int i) const
{
assert( (i >= 0) && (i < corners()) );
return mapping()(geometry_.corner(i));
}
GlobalCoordinate center () const
{
return mapping()(geometry_.center());
}
GlobalCoordinate global (const LocalCoordinate& local) const
{
return mapping()(geometry_.global(local));
}
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
{
auto&& jLocal = geometry_.jacobian(local);
auto&& jMapping = (*dMapping_)(geometry_.global(local));
return jMapping * jLocal;
}
JacobianTransposed jacobianTransposed (const LocalCoordinate& local) const
{
return transpose(jacobian(local));
}
JacobianInverse jacobianInverse (const LocalCoordinate& local) const
{
JacobianInverse out;
MatrixHelper::leftInvA(jacobian(local), out);
return out;
}
JacobianInverseTransposed jacobianInverseTransposed (const LocalCoordinate& local) const
{
return transpose(jacobianInverse(local));
}
friend ReferenceElement referenceElement (const MappedGeometry& geometry)
{
return geometry.refElement();
}
protected:
// the internal stored reference element
ReferenceElement refElement () const
{
return referenceElement(geometry_);
}
private:
// internal reference to the stored mapping
const Mapping& mapping () const
{
return *mapping_;
}
// internal reference to the wrapped geometry
const Geometry& geometry () const
{
return geometry_;
}
private:
CopyableOptional<Mapping> mapping_;
CopyableOptional<DerivativeMapping> dMapping_;
Geometry geometry_;
bool affine_;
};
// deduction guides
template <class Map, class Geo>
MappedGeometry (const Map&, const Geo&)
-> MappedGeometry<Map,Geo>;
template <class Map, class Geo>
MappedGeometry (const Map&, const Geo&, bool)
-> MappedGeometry<Map,Geo>;
} // end namespace Dune
#endif // DUNE_GEOMETRY_MAPPEDGEOMETRY_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
std
STL namespace.
Dune
Definition affinegeometry.hh:22
Dune::MappedGeometry
MappedGeometry(const Map &, const Geo &) -> MappedGeometry< Map, Geo >
Dune::Geo
Definition affinegeometry.hh:28
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::MappedGeometry
Geometry parametrized by a LocalFunction and a LocalGeometry.
Definition mappedgeometry.hh:65
Dune::MappedGeometry::GlobalCoordinate
std::remove_reference_t< decltype(std::declval< Map >()(std::declval< typename Geo::GlobalCoordinate >()))> GlobalCoordinate
type of global coordinates
Definition mappedgeometry.hh:71
Dune::MappedGeometry::center
GlobalCoordinate center() const
Map the center of the wrapped geometry.
Definition mappedgeometry.hh:162
Dune::MappedGeometry::jacobianInverseTransposed
JacobianInverseTransposed jacobianInverseTransposed(const LocalCoordinate &local) const
Obtain the transposed of the Jacobian's inverse.
Definition mappedgeometry.hh:310
Dune::MappedGeometry::jacobianInverse
JacobianInverse jacobianInverse(const LocalCoordinate &local) const
Obtain the Jacobian's inverse.
Definition mappedgeometry.hh:296
Dune::MappedGeometry::jacobianTransposed
JacobianTransposed jacobianTransposed(const LocalCoordinate &local) const
Obtain the transposed of the Jacobian.
Definition mappedgeometry.hh:284
Dune::MappedGeometry::MappedGeometry
MappedGeometry(Map_ &&mapping, Geo_ &&geometry, bool affine=false)
Constructor from mapping to parametrize the geometry.
Definition mappedgeometry.hh:125
Dune::MappedGeometry::JacobianTransposed
FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed
type of jacobian transposed
Definition mappedgeometry.hh:89
Dune::MappedGeometry::mydimension
static constexpr int mydimension
geometry dimension
Definition mappedgeometry.hh:77
Dune::MappedGeometry::type
GeometryType type() const
Obtain the geometry type from the reference element.
Definition mappedgeometry.hh:143
Dune::MappedGeometry::LocalCoordinate
typename Geo::LocalCoordinate LocalCoordinate
type of local coordinates
Definition mappedgeometry.hh:68
Dune::MappedGeometry::ctype
typename Geo::ctype ctype
coordinate type
Definition mappedgeometry.hh:74
Dune::MappedGeometry::local
LocalCoordinate local(const GlobalCoordinate &y, Impl::GaussNewtonOptions< ctype > opts={}) const
Evaluate the inverse coordinate mapping.
Definition mappedgeometry.hh:197
Dune::MappedGeometry::Geometry
Geo Geometry
Definition mappedgeometry.hh:108
Dune::MappedGeometry::corner
GlobalCoordinate corner(int i) const
Obtain coordinates of the i-th corner.
Definition mappedgeometry.hh:155
Dune::MappedGeometry::coorddimension
static constexpr int coorddimension
coordinate dimension
Definition mappedgeometry.hh:80
Dune::MappedGeometry::JacobianInverseTransposed
FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed
type of jacobian inverse transposed
Definition mappedgeometry.hh:95
Dune::MappedGeometry::Mapping
Map Mapping
Definition mappedgeometry.hh:105
Dune::MappedGeometry::global
GlobalCoordinate global(const LocalCoordinate &local) const
Evaluate the coordinate mapping.
Definition mappedgeometry.hh:176
Dune::MappedGeometry::refElement
ReferenceElement refElement() const
Definition mappedgeometry.hh:323
Dune::MappedGeometry::volume
Volume volume(const QuadratureRule< ctype, mydimension > &quadRule) const
Obtain the volume of the mapping's image by given quadrature rules.
Definition mappedgeometry.hh:259
Dune::MappedGeometry::corners
int corners() const
Obtain number of corners of the corresponding reference element.
Definition mappedgeometry.hh:149
Dune::MappedGeometry::JacobianInverse
FieldMatrix< ctype, mydimension, coorddimension > JacobianInverse
type of jacobian inverse
Definition mappedgeometry.hh:92
Dune::MappedGeometry::referenceElement
friend ReferenceElement referenceElement(const MappedGeometry &geometry)
Obtain the reference-element related to this geometry.
Definition mappedgeometry.hh:316
Dune::MappedGeometry::volume
Volume volume(Impl::ConvergenceOptions< ctype > opts={}) const
Obtain the volume of the mapping's image.
Definition mappedgeometry.hh:241
Dune::MappedGeometry::Volume
std::remove_reference_t< decltype(Dune::power(std::declval< ctype >(), mydimension))> Volume
type of volume
Definition mappedgeometry.hh:83
Dune::MappedGeometry::affine
bool affine() const
Is this mapping affine? Not in general, since we don't know anything about the mapping....
Definition mappedgeometry.hh:137
Dune::MappedGeometry::DerivativeMapping
std::remove_reference_t< decltype(derivative(std::declval< Map >()))> DerivativeMapping
Definition mappedgeometry.hh:111
Dune::MappedGeometry::jacobian
Jacobian jacobian(const LocalCoordinate &local) const
Obtain the Jacobian.
Definition mappedgeometry.hh:272
Dune::MappedGeometry::MatrixHelper
Impl::FieldMatrixHelper< ctype > MatrixHelper
Definition mappedgeometry.hh:102
Dune::MappedGeometry::Jacobian
FieldMatrix< ctype, coorddimension, mydimension > Jacobian
type of jacobian
Definition mappedgeometry.hh:86
Dune::MappedGeometry::integrationElement
ctype integrationElement(const LocalCoordinate &local) const
Obtain the integration element.
Definition mappedgeometry.hh:225
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