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 const& 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
    {
      thread_local std::vector<typename Traits::GradientRange> grad;
      evaluateGradient(x, grad);
      return grad;
    }

    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 const& 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
    {
      thread_local std::vector<typename Traits::PartialRange> d_comp;
      evaluatePartial(x, comp, d_comp);
      return d_comp;
    }

  private:
    LocalBasis const& localBasis_;

    Cache const& cache_;

    Geometry const& geometry_;

    std::size_t size_;
  };

} // namespace AMDiS