mappedgeometry.hh

// -*- 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/affinegeometry.hh> // for FieldMatrixHelper
#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>
 
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