shapefunctioncache.hh

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*  This file is part of the library KASKADE 7                               */
/*  https://www.zib.de/research/projects/kaskade7-finite-element-toolbox     */
/*                                                                           */
/*  Copyright (C) 2002-2019 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef SHAPEFUNCTIONCACHE_HH
#define SHAPEFUNCTIONCACHE_HH
#include <map>
#include <tuple>
 
#include "boost/multi_array.hpp"
#include "boost/signals2.hpp"
#include "boost/thread.hpp"
 
#include <boost/fusion/algorithm.hpp>
#include <boost/fusion/sequence.hpp>
#include <boost/mpl/int.hpp>
#include <boost/mpl/range_c.hpp>
 
#include "dune/geometry/quadraturerules.hh"
#include "dune/istl/bvector.hh"
 
#include "fem/fetransfer.hh"
#include "fem/pshapefunctions.hh"
#include "linalg/tensor.hh"
#include "utilities/threading.hh"
 
 
//---------------------------------------------------------------------
namespace Kaskade
{
  template <typename Scalar, int dim, int components=1>
  struct VariationalArg
  {
    VariationalArg() {}
    
    using ValueType = Dune::FieldVector<Scalar,components>;
    
    explicit VariationalArg(ValueType const& value_)
    : value(value_), derivative(0), hessian(0) {}
    
    
    VariationalArg(ValueType const& value_,
                   Dune::FieldMatrix<Scalar,components,dim> const& derivative_)
    : value(value_), derivative(derivative_),  hessian(0) {}
    
    VariationalArg(ValueType const& value_,
                   Dune::FieldMatrix<Scalar,components,dim> const& derivative_,
                   Tensor<Scalar,components,dim,dim> const& hessian_)
    : value(value_), derivative(derivative_), hessian(hessian_) {}
    
    
    ValueType value;
 
    Dune::FieldMatrix<Scalar,components,dim> derivative;
    Tensor<Scalar,components,dim,dim>        hessian;
  };
 
  //------------------------------------------------------------------------------------------------
  //------------------------------------------------------------------------------------------------
 
  template <class G, class T, int ComponentsEnd=4>
  class ShapeFunctionCache
  {
 
  private:
    typedef typename boost::mpl::range_c<int,1,ComponentsEnd>::type ComponentRange;
 
    // This is a functor that, given a number of
    // components (a boost::mpl integer), defines an associative container type
    // with a tuple (shape function set address, quadrature rule address) as key
    // and a twodimensional boost::multi_array of pairs of vector and
    // matrix with given components as value type.
    struct MakeMapType
    {
      template <class Int>
      auto operator()(Int i) 
      {
        static int const nComp = Int::value;
 
        typedef VariationalArg<T,G::dimension,nComp> EntryType;
        typedef boost::multi_array<EntryType,2> DataType;
        typedef ShapeFunctionSet<typename G::ctype,G::dimension,T,Int::value> Sfs;
        typedef std::tuple<Sfs const*,void const*> KeyType;
        return std::map<KeyType,DataType>();
      }
    };
 
    typedef typename boost::fusion::result_of::as_vector<
                typename boost::fusion::result_of::transform<ComponentRange,MakeMapType>::type
        >::type MapsType;
 
  public:
 
    template <int nComp>
    struct DataType 
    {
      typedef typename std::result_of<MakeMapType(boost::mpl::int_<nComp>)>::type::mapped_type type;
    };
    
    template <int nComp>
    struct LocalDataType 
    {
      typedef typename DataType<nComp>::type::template const_array_view<1>::type type;
    };
 
    
    template <int nComp, int subDim>
    Dune::FieldVector<T,nComp> evaluateFunction(ShapeFunctionSet<typename G::ctype,G::dimension,T,nComp> const& sfs, int isf,
                                                Dune::QuadratureRule<typename G::ctype,subDim> const& qr, int ip, int subId)
    {
      assert(0<=isf && isf<sfs.size());
      assert(0<=ip && ip<qr.size());
      return get(sfs,qr,subId)[ip][isf].value;
    }
 
    template <int nComp, int subDim>
    Dune::FieldMatrix<T,nComp,G::dimension> const& evaluateDerivative(ShapeFunctionSet<typename G::ctype,G::dimension,T,nComp> const& sfs, int isf,
                                                                      Dune::QuadratureRule<typename G::ctype,subDim> const& qr, int ip,
                                                                      int subId)
    {
      assert(0<=isf && isf<sfs.size());
      assert(0<=ip && ip<qr.size());
      return get(sfs,qr,subId)[ip][isf].derivative;
    }
 
    template <int nComp, int subDim>
    typename LocalDataType<nComp>::type
    evaluate(ShapeFunctionSet<typename G::ctype,G::dimension,T,nComp> const& sfs,
             Dune::QuadratureRule<typename G::ctype,subDim> const& qr, int ip, int subId)
    {
      typedef typename DataType<nComp>::type DType;
      typedef typename DType::index_range range;
      typename DType::index_gen indices;
      assert(sfs.size()>=0);
      return get(sfs,qr,subId)[indices[ip][range(0,sfs.size())]];
    }
 
  
   template <int nComp, int subDim>
   typename DataType<nComp>::type const& evaluate(ShapeFunctionSet<typename G::ctype,G::dimension,T,nComp> const& sfs,
                                                  Dune::QuadratureRule<typename G::ctype,subDim> const& qr, int subId) 
   {
     return get(sfs,qr,subId);
   }
   
   
 
  int maximumSubentityIndex() const 
  {
    return maxSubId;
  }
 
  private:
    // TODO: use boost::container::flat_map from boost 1.49 on
    static int const maxSubId = 11; // 12 edges in a cube - that's the upper limit (up to 3d, of course...)
    
    // For each subentity number, we maintain a different cache for direct access
    MapsType caches[maxSubId+1];
    
    template <int nComp, int subDim>
    typename DataType<nComp>::type const&
    get(ShapeFunctionSet<typename G::ctype,G::dimension,T,nComp> const& sfs, Dune::QuadratureRule<typename G::ctype,subDim> const& qr, int subId)
    {
      assert(0<=subId && subId<=maxSubId);
      using MapType = typename std::result_of<MakeMapType(boost::mpl::int_<nComp>)>::type;
      typedef typename MapType::key_type KeyType;
      typedef typename DataType<nComp>::type DType;
      MapType& cache = boost::fusion::at_c<nComp-1>(caches[subId]);
 
      typename MapType::const_iterator i = cache.find(KeyType(&sfs,&qr));
      
      // check whether the shape functions' values are already available
      if (i==cache.end())
      {
        // No, they are not. Evaluate shape functions and store the result.
        Dune::FieldVector<T,G::dimension> pos;
        auto const& refElem = sfs.referenceElement();
        assert(sfs.size()>=0);
        DType data(boost::extents[qr.size()][sfs.size()]);
 
        if (G::dimension != subDim)
          for (int isf=0; isf<sfs.size(); ++isf)
            for (int ip=0; ip<qr.size(); ++ip)
            {
              {
                boost::mutex::scoped_lock lock(refElementMutex);
                pos = refElem.template geometry<G::dimension-subDim>(subId).global(qr[ip].position());
              }
              data[ip][isf].value = sfs[isf].evaluateFunction(pos);
              data[ip][isf].derivative = sfs[isf].evaluateDerivative(pos);
              data[ip][isf].hessian = sfs[isf].evaluate2ndDerivative(pos);
            }
        else
          for (int isf=0; isf<sfs.size(); ++isf)
            for (int ip=0; ip<qr.size(); ++ip)
            {
              pos = refElem.template geometry<G::dimension-subDim>(subId).global(qr[ip].position());
              data[ip][isf].value = sfs[isf].evaluateFunction(pos);
              data[ip][isf].derivative = sfs[isf].evaluateDerivative(pos);
              data[ip][isf].hessian = sfs[isf].evaluate2ndDerivative(pos);
            }
 
        // store the values, obtaining an iterator pointing to the newly inserted values
        i = cache.insert(cache.begin(),typename MapType::value_type(KeyType(&sfs,&qr),data));
      }
      
      // From here on, the values are in the map (at the position pointed to by i).
      return i->second;
    }
  };
} // end of namespace Kaskade
//---------------------------------------------------------------------
//---------------------------------------------------------------------
 
 
#endif