LocalToGlobalAdapter.hpp

#pragma once
 
#include <cstddef>
#include <vector>
 
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/common/typetraits.hh>
#include <dune/geometry/type.hh>
 
#include <amdis/common/DerivativeTraits.hpp>
#include <amdis/common/FieldMatVec.hpp>
 
namespace AMDiS
{
 
  template <class LocalBasisTraits, std::size_t dimGlobal>
  struct LocalToGlobalBasisAdapterTraits
  {
    static const std::size_t dimDomainLocal = LocalBasisTraits::dimDomain;
    static const std::size_t dimDomainGlobal = dimGlobal;
    static const std::size_t dimRange = LocalBasisTraits::dimRange;
 
    using DomainField = typename LocalBasisTraits::DomainFieldType;
    using DomainLocal = typename LocalBasisTraits::DomainType;
    using DomainGlobal = FieldVector<DomainField, dimDomainGlobal>;
 
    using RangeField = typename LocalBasisTraits::RangeFieldType;
    using Range = typename LocalBasisTraits::RangeType;
 
    using GradientRange = typename DerivativeTraits<Range(DomainGlobal), tag::gradient>::Range;
    using PartialRange = typename DerivativeTraits<Range(DomainGlobal), tag::partial>::Range;
  };
 
 
  template <class BasisCache, class Geometry>
  class LocalToGlobalBasisAdapter
  {
    using LocalBasis = typename BasisCache::LocalBasis;
 
    static_assert(std::is_same_v<typename LocalBasis::Traits::DomainFieldType, typename Geometry::ctype>,
      "LocalToGlobalBasisAdapter: LocalBasis must use the same ctype as Geometry");
 
    static_assert(std::size_t(LocalBasis::Traits::dimDomain) == std::size_t(Geometry::mydimension),
      "LocalToGlobalBasisAdapter: LocalBasis domain dimension must match local dimension of Geometry");
 
  public:
    using Cache = BasisCache;
    using Traits = LocalToGlobalBasisAdapterTraits<typename LocalBasis::Traits, Geometry::coorddimension>;
 
 
    LocalToGlobalBasisAdapter(Cache const& cache, Geometry const& geometry)
      : localBasis_(cache.finiteElement().localBasis())
      , cache_(cache)
      , geometry_(geometry)
      , size_(localBasis_.size())
    {}
 
    std::size_t size() const { return size_; }
 
 
    std::size_t order() const
    {
      if (geometry_.affine())
        // affine linear
        return localBasis_.order();
      else
        // assume at most order dim
        return localBasis_.order() + Traits::dimDomainGlobal - 1;
    }
 
    void evaluateFunction(typename Traits::DomainLocal const& x,
                          std::vector<typename Traits::Range>& out) const
    {
      out = cache_.localBasisValuesAt(x);
    }
 
    auto const& valuesAt(typename Traits::DomainLocal const& x) const
    {
      return cache_.localBasisValuesAt(x);
    }
 
    void evaluateGradient(typename Traits::DomainLocal const& x,
                          std::vector<typename Traits::GradientRange>& out) const
    {
      auto const& localJacobian = cache_.localBasisJacobiansAt(x);
      auto&& geoJacobian = geometry_.jacobianInverseTransposed(x);
 
      out.resize(size_);
      for (std::size_t i = 0; i < size_; ++i)
        geoJacobian.mv(Dune::MatVec::as_vector(localJacobian[i]),
                       Dune::MatVec::as_vector(out[i]));
    }
 
    auto const& gradientsAt(typename Traits::DomainLocal const& x) const
    {
      evaluateGradient(x, outGrad_);
      return outGrad_;
    }
 
    void evaluatePartial(typename Traits::DomainLocal const& x,
                         std::size_t comp,
                         std::vector<typename Traits::PartialRange>& out) const
    {
      auto const& localJacobian = cache_.localBasisJacobiansAt(x);
      // NOTE: geoJacobian might be a Dune::DiagonalMatrix with limited interface!
      auto&& geoJacobian = geometry_.jacobianInverseTransposed(x);
 
      out.resize(size_);
      typename Traits::GradientRange grad;
      auto&& grad_ = Dune::MatVec::as_vector(grad);
      for (std::size_t i = 0; i < size_; ++i) {
        geoJacobian.mv(Dune::MatVec::as_vector(localJacobian[i]), grad_);
        out[i] = grad_[comp];
      }
    }
 
    auto const& partialsAt(typename Traits::DomainLocal const& x, std::size_t comp) const
    {
      evaluatePartial(x, comp, outPartials_);
      return outPartials_;
    }
 
  private:
    LocalBasis const& localBasis_;
 
    Cache const& cache_;
 
    Geometry const& geometry_;
 
    std::size_t size_;
 
    // data return as reference
    mutable std::vector<typename Traits::GradientRange> outGrad_;
    mutable std::vector<typename Traits::PartialRange> outPartials_;
  };
 
} // namespace AMDiS