localfiniteelementgeometry.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_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/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 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