LocalBasisCache.hpp

#pragma once

#include <vector>

#include <dune/common/hash.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/geometry/type.hh>

#include <amdis/utility/ConcurrentCache.hpp>

namespace AMDiS
{
  template <class LocalBasisType>
  class LocalBasisCache
  {
  public:
    using Traits = typename LocalBasisType::Traits;

    using DomainType = typename Traits::DomainType;
    using RangeType = typename Traits::RangeType;
    using RangeFieldType = typename Traits::RangeFieldType;
    using JacobianType = typename Traits::JacobianType;

    using QuadratureRule = Dune::QuadratureRule<RangeFieldType,Traits::dimDomain>;
    using QuadratureRules = Dune::QuadratureRules<RangeFieldType,Traits::dimDomain>;

  private:

    /// Pair of GeometryType and quadrature order and size
    template <class Tag>
    struct QuadKey
    {
      unsigned int id;    // topologyId
      int p;              // quadrature order
      std::size_t size;   // number of quadrature points

      struct hasher
      {
        std::size_t operator()(QuadKey const& t) const
        {
          std::size_t seed = 0;
          Dune::hash_combine(seed, t.id);
          Dune::hash_combine(seed, t.p);
          Dune::hash_combine(seed, t.size);
          return seed;
        }
      };

      friend bool operator==(QuadKey const& lhs, QuadKey const& rhs)
      {
        return std::tie(lhs.id,lhs.p,lhs.size) == std::tie(rhs.id,rhs.p,rhs.size);
      }
    };

    /// Pair of GeometryType and local coordinates
    template <class Tag>
    struct CoordsKey
    {
      unsigned int id;    // topologyId
      DomainType local;   // local coordinate

      struct hasher
      {
        std::size_t operator()(CoordsKey const& t) const
        {
          std::size_t seed = 0;
          Dune::hash_combine(seed, t.id);
          Dune::hash_range(seed, t.local.begin(), t.local.end());
          return seed;
        }
      };

      friend bool operator==(CoordsKey const& lhs, CoordsKey const& rhs)
      {
        return std::tie(lhs.id,lhs.local) == std::tie(rhs.id,rhs.local);
      }
    };

    using LocalValuesKey = CoordsKey<struct ValuesTag>;
    using LocalGradientsKey = CoordsKey<struct GradientsTag>;

    using ValuesKey = QuadKey<struct ValuesTag>;
    using GradientsKey = QuadKey<struct GradientsTag>;

  public:
    using ShapeValues = std::vector<RangeType>;
    using ShapeGradients = std::vector<JacobianType>;

    using ShapeValuesContainer = std::vector<ShapeValues>;
    using ShapeGradientsContainer = std::vector<ShapeGradients>;


    LocalBasisCache(LocalBasisType const& localBasis)
      : localBasis_(&localBasis)
    {}

    // evaluate local basis-function at all quadrature points
    template <class LocalContext>
    ShapeValuesContainer const& evaluateFunctionAtQp(LocalContext const& context, QuadratureRule const& quad) const
    {
      ValuesKey key{context.type().id(),quad.order(),quad.size()};
      return ShapeValuesContainerCache::get(key, [&](ValuesKey const&)
      {
        ShapeValuesContainer data(quad.size());
        for (std::size_t iq = 0; iq < quad.size(); ++iq)
          localBasis_->evaluateFunction(context.local(quad[iq].position()), data[iq]);
        return data;
      });
    }

    // evaluate local basis-gradients at all quadrature points
    template <class LocalContext>
    ShapeGradientsContainer const& evaluateJacobianAtQp(LocalContext const& context, QuadratureRule const& quad) const
    {
      GradientsKey key{context.type().id(),quad.order(),quad.size()};
      return ShapeGradientsContainerCache::get(key, [&](GradientsKey const&)
      {
        ShapeGradientsContainer data(quad.size());
        for (std::size_t iq = 0; iq < quad.size(); ++iq)
          localBasis_->evaluateJacobian(context.local(quad[iq].position()), data[iq]);
        return data;
      });
    }


    // evaluate local basis-function at point
    template <class Element>
    ShapeValues const& evaluateFunction(Element const& element, DomainType const& local) const
    {
      LocalValuesKey key{element.type().id(),local};
      return ShapeValuesCache::get(key, [&](LocalValuesKey const&)
      {
        ShapeValues data;
        localBasis_->evaluateFunction(local, data);
        return data;
      });
    }

    // evaluate local basis-gradients at point
    template <class Element>
    ShapeGradients const& evaluateJacobian(Element const& element, DomainType const& local) const
    {
      LocalGradientsKey key{element.type().id(),local};
      return ShapeGradientsCache::get(key, [&](LocalGradientsKey const&)
      {
        ShapeGradients data;
        localBasis_->evaluateJacobian(local, data);
        return data;
      });
    }

  private:
    using ShapeValuesCache = ConcurrentCache<LocalValuesKey, ShapeValues, ThreadLocalPolicy,
      std::unordered_map<LocalValuesKey, ShapeValues, typename LocalValuesKey::hasher>>;
    using ShapeGradientsCache = ConcurrentCache<LocalGradientsKey, ShapeGradients, ThreadLocalPolicy,
      std::unordered_map<LocalGradientsKey, ShapeGradients, typename LocalGradientsKey::hasher>>;

    using ShapeValuesContainerCache = ConcurrentCache<ValuesKey, ShapeValuesContainer, ThreadLocalPolicy,
      std::unordered_map<ValuesKey, ShapeValuesContainer, typename ValuesKey::hasher>>;
    using ShapeGradientsContainerCache = ConcurrentCache<GradientsKey, ShapeGradientsContainer, ThreadLocalPolicy,
      std::unordered_map<GradientsKey, ShapeGradientsContainer, typename GradientsKey::hasher>>;

    LocalBasisType const* localBasis_;
  };

} // end namespace AMDiS