hierarchicshapefunctions.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-2020 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef HIERARCHICSHAPEFUNCTIONS_HH
#define HIERARCHICSHAPEFUNCTIONS_HH
 
#include "fem/barycentric.hh"
#include "fem/gridcombinatorics.hh"
#include "fem/pshapefunctions.hh"
#include "linalg/dynamicMatrix.hh"
#include "utilities/combinatorics.hh"
 
#include <dune/common/fvector.hh>
#include <dune/geometry/referenceelements.hh>
 
#include <algorithm>
#include <cmath>
#include <functional>
#include <map>
#include <memory>
#include <numeric>
#include <tuple>
 
//---------------------------------------------------------------------
namespace Kaskade
{
  // forward declarations
  template <class,int,class> class HierarchicSimplexShapeFunctionSet;
  template <class,int,class> class HierarchicExtensionSimplexShapeFunctionSet;
 
  namespace HierarchicDetail {
 
    int nSubentities(int d, int c);
 
 
    int cdimSize(int d, int p, int c);
 
 
    int size(int d, int p);
 
 
    std::pair<int,int> codim(int d, int p, int k);
 
 
    std::tuple<int,int,int> tupleIndex(int d, int p, int k);
 
 
    double shapeFunction(int d, int p, int k, double const b[]);
 
 
    void dShapeFunction(int d, int p, int k, double const b[], double* grad);
 
 
    int actualOrder(int d, int p, int k);
 
 
    int entityPermutation(int d, int c, int e, int const vGlobal[]);
 
 
    std::tuple<int,int> sfPermutation(int d, int p, int c, int e, int k, int const vGlobal[]);
 
 
    template <class ctype, int dimension, class Scalar, bool restricted>
    struct ChooseShapeFunctionSet
    {
      typedef HierarchicSimplexShapeFunctionSet<ctype,dimension,Scalar> type;
    };
 
    template <class ctype, int dimension, class Scalar>
    struct ChooseShapeFunctionSet<ctype,dimension,Scalar,true>
    {
      typedef RestrictedShapeFunctionSet<HierarchicSimplexShapeFunctionSet<ctype,dimension,Scalar> > type;
    };
 
    template <class ctype, int dimension, class Scalar, bool restricted>
    struct ChooseExtensionShapeFunctionSet
    {
      typedef HierarchicExtensionSimplexShapeFunctionSet<ctype,dimension,Scalar> type;
    };
 
    template <class ctype, int dimension, class Scalar>
    struct ChooseExtensionShapeFunctionSet<ctype,dimension,Scalar,true>
    {
      typedef RestrictedShapeFunctionSet<HierarchicExtensionSimplexShapeFunctionSet<ctype,dimension,Scalar> > type;
    };
 
  } // End of namespace HierarchicDetail
 
  //---------------------------------------------------------------------
 
  template <class ctype, int dimension, class Scalar>
  class HierarchicSimplexShapeFunction: public ShapeFunction<ctype,dimension,Scalar>
  {
  public:
    static int const dim = dimension;
    static int const comps = 1;
 
    typedef ctype CoordType;
    typedef Scalar ResultType;
 
    HierarchicSimplexShapeFunction() {}
 
    explicit HierarchicSimplexShapeFunction(int order, int k):
                    order_(order), number_(k)
    {
      assert(order>=0);
      assert(0<=k && k<HierarchicDetail::size(dim,order));
      std::tie(codim_,entity_,entityIndex_) = HierarchicDetail::tupleIndex(dim,order,k);
      actualOrder = HierarchicDetail::actualOrder(dim,order_,k);
    }
 
    HierarchicSimplexShapeFunction(HierarchicSimplexShapeFunction const& other)
    : actualOrder(other.actualOrder), order_(other.order_), number_(other.number_),
      codim_(other.codim_), entity_(other.entity_), entityIndex_(other.entityIndex_)
    {}
 
    virtual Dune::FieldVector<ResultType,comps> evaluateFunction(Dune::FieldVector<CoordType,dim> const& x) const
    {
      Dune::FieldVector<CoordType,dim+1> zeta = barycentric(x);
      return HierarchicDetail::shapeFunction(dim,order_,number_,&zeta[0]);
    }
 
    virtual Dune::FieldMatrix<Scalar,1,dim> evaluateDerivative(Dune::FieldVector<CoordType,dim> const& x) const
    {
      Dune::FieldVector<double,dim+1> zeta = barycentric(x);
      Dune::FieldVector<double,dim+1> bgrad;
      HierarchicDetail::dShapeFunction(dim,order_,number_,&zeta[0],&bgrad[0]);
      
      Dune::FieldMatrix<Scalar,1,dim> grad;
      for (int i=0; i<dim; ++i)
        grad[0][i] = bgrad[i]-bgrad[dim];
      
      
        // for (int i=0; i<dim; ++i) {
        //   double h = 1e-5;
        //   Dune::FieldVector<CoordType,dim> yl(x), yr(x);
        //   yl[i] -= h;
        //   yr[i] += h;
        //   double diff = (evaluateFunction(yr)-evaluateFunction(yl))/(2*h);
        //
        //   if (std::abs(diff-grad[0][i]) > 1e-4) {
          //     std::cout << "i=" << i << ": numdiff=" << diff << "  bdiff=" << grad[0][i] << '\n';
        //   }
        //   grad[0][i] = diff;
        // }
        
        
      return grad;
    }
 
    virtual Tensor<Scalar,1,dim,dim> evaluate2ndDerivative(Dune::FieldVector<CoordType,dim> const& x) const
    {
      static bool visited = false;
      if (!visited) 
      {
        visited = true;
        std::cerr << "Not implemented: Hierarchic shape function 2nd derivative at " << __FILE__ << ":" << __LINE__ << "\n";
      }
      return 0;
    }
 
    virtual std::tuple<int,int,int,int> location() const
    {
      return std::make_tuple(order_,codim_,entity_,entityIndex_);
    }
 
    int order() const { return actualOrder; }
 
 
  private:
    int actualOrder;
 
    int order_;
    int number_;
    int codim_;
    int entity_;
    int entityIndex_;
  };
 
 
  //---------------------------------------------------------------------
 
  template <class ctype, int dimension, class Scalar>
  class HierarchicSimplexShapeFunctionSet : public ShapeFunctionSet<ctype,dimension,Scalar>
  {
  public:
    typedef HierarchicSimplexShapeFunction<ctype,dimension,Scalar> value_type;
    typedef typename ShapeFunctionSet<ctype,dimension,Scalar>::Matrix Matrix;
 
    HierarchicSimplexShapeFunctionSet(int ord)
     : ShapeFunctionSet<ctype,dimension,Scalar>(Dune::GeometryType(Dune::GeometryType::simplex,dimension))
    {
      this->order_ = -1;
 
      // Empty dummy set.
      if (ord<0) {
        this->size_ = 0;
        return;
      }
 
      for (int o=0; o<=ord; ++o) {
        int n = HierarchicDetail::size(dimension,o);
        for (int i=0; i<n; ++i) {
          sf.push_back(value_type(o,i));
          this->order_ = std::max(this->order_,sf.back().order());
        }
      }
 
      // Use equidistant "interpolation" nodes. Note that we may have
      // more shape functions than nodes with the same order. To enable
      // a faithful reconstruction during grid transfer, we use
      // interpolation nodes for the maximal ACTUAL order of the shape
      // functions and a least squares "interpolation" approach.
      //
      // TODO: Evaluate the alternative of using interpolation nodes for
      // the maximal *NOMINAL* order of the shape functions and a least
      // coefficient l2-norm interpolation.  This does not allow a
      // faithful reconstruction (except for the lower dimensional
      // subspace of nominal order), but is faster due to the lower
      // number of interpolation points. Note that for some applications
      // the exact reproduction is necessary.
      int n = binomial(dimension+this->order_,dimension);
      this->iNodes.insert(this->iNodes.end(),n,Dune::FieldVector<ctype,dimension>());
      Dune::FieldVector<int,dimension> b(0);
      for (int i=0; i<n; ++i) {
        for (int j=0; j<dimension; ++j)
          this->iNodes[i][j] = static_cast<ctype>(b[j])/this->order_;
        next_multiindex(b.begin(),b.end(),this->order_);
      }
 
      // Compute the pseudoinverse for use in interpolation.
      computePseudoInverse();
 
      this->size_ = sf.size();
    }
 
    HierarchicSimplexShapeFunctionSet(HierarchicSimplexShapeFunctionSet const& other) = default;
 
    HierarchicSimplexShapeFunctionSet(HierarchicSimplexShapeFunctionSet&& other) = default;
 
    virtual ~HierarchicSimplexShapeFunctionSet(){}
 
    virtual value_type const& operator[](int i) const{ return sf[i]; }
 
    virtual Dune::GeometryType type() const { return Dune::GeometryType(Dune::GeometryType::simplex,dimension); }
 
    virtual void interpolate(typename ShapeFunctionSet<ctype,dimension,Scalar>::SfValueArray const& A,
                             Matrix& IA) const
    {
      // TODO: implement without temporaries
      IA = pinv * A;
    }
 
    void removeShapeFunction(size_t index)
    {
      assert(index>=0 && index < sf.size());
      sf.erase(sf.begin()+index);
      this->size_ = sf.size();
      computePseudoInverse();
    }
 
  protected:
    std::vector<value_type> sf;
    DynamicMatrix<Dune::FieldMatrix<Scalar,1,1>> pinv;
 
  private:
    void computePseudoInverse()
    {
      pinv.resize(sf.size(),this->iNodes.size());
      DynamicMatrix<Scalar> A(this->iNodes.size(),sf.size());
      for (int j=0; j<sf.size(); ++j)
        for (int i=0; i<this->iNodes.size(); ++i)
          A[i][j] = sf[j].evaluateFunction(this->iNodes[i]);
 
      pinv = pseudoinverse(A);
    }
  };
 
  //---------------------------------------------------------------------
  //---------------------------------------------------------------------
  //---------------------------------------------------------------------
 
 
  template <class ctype, int dimension, class Scalar, bool restricted = false>
  class HierarchicShapeFunctionSetContainer: public ShapeFunctionSetContainer<ctype,dimension,Scalar>
  {
  public:
    typedef ShapeFunctionSet<ctype,dimension,Scalar> value_type;
 
  private:
    typedef std::map<std::pair<Dune::GeometryType,int>,value_type const*> Container;
    typedef typename HierarchicDetail::ChooseShapeFunctionSet<ctype,dimension,Scalar,restricted>::type ShapeFunctionType;
 
  public:
 
    HierarchicShapeFunctionSetContainer(int maxOrder)
    {
      // Fill the container with shape function sets. This is not done
      // on demand in operator(), since then operator() is not easily
      // and efficiently made thread-safe.
      //
      // Currently only shape functions for simplices are defined.
      Dune::GeometryType simplex(Dune::GeometryType::simplex,dimension);
 
      value_type* sfsPrevious = 0;
 
      for (int i=-1; i<=maxOrder; ++i) {
        value_type* sfs = new ShapeFunctionType(i);
        sfs->initHierarchicalProjection(sfsPrevious);
        container.insert(typename Container::value_type(std::make_pair(simplex,i),sfs));
        sfsPrevious = sfs;
      }
    }
 
    virtual ~HierarchicShapeFunctionSetContainer()
    {
      for (typename Container::iterator i=container.begin(); i!=container.end(); ++i)
        delete i->second;
    }
 
    virtual value_type const& operator()(Dune::GeometryType type, int order) const
    {
      typename Container::const_iterator i = container.find(std::make_pair(type,order));
      assert(i!=container.end()); // should throw, not assert!
      return *i->second;
    }
 
  private:
    // A map storing Hierarchic shape function sets associated with
    // (reference element type, order). Due to runtime polymorphism,
    // pointers are stored in the map. This class owns the shape
    // function sets pointed to, and has to delete them in the
    // destructor.
    Container container;
  };
 
 
  template <class ctype, int dimension, class Scalar>
  typename HierarchicShapeFunctionSetContainer<ctype,dimension,Scalar>::value_type const&
  hierarchicShapeFunctionSet(Dune::GeometryType type, int order)
  {
    // This static variable will be destructed on program shutdown,
    // causing its contents to be deleted properly from the heap by
    // the destructor. No memory leaks should be induced.
    static HierarchicShapeFunctionSetContainer<ctype,dimension,Scalar> container(7);
 
    return container(type,order);
  }
 
  //---------------------------------------------------------------------
  //---------------------------------------------------------------------
 
 
  template <class ctype, int dimension, class Scalar>
  class HierarchicExtensionSimplexShapeFunctionSet: public ShapeFunctionSet<ctype,dimension,Scalar>
  {
  public:
    typedef HierarchicSimplexShapeFunction<ctype,dimension,Scalar> value_type;
    typedef typename ShapeFunctionSet<ctype,dimension,Scalar>::Matrix Matrix;
 
    HierarchicExtensionSimplexShapeFunctionSet(int ord)
    : ShapeFunctionSet<ctype,dimension,Scalar>(Dune::GeometryType(Dune::GeometryType::simplex,
                                                                  dimension))
    {
      this->order_ = -1;
 
      // Empty dummy set.
      if (ord<0) {
        this->size_ = 0;
        return;
      }
 
      int n = HierarchicDetail::size(dimension,ord);
      if (n==0) {
        this->size_ = 0;
        return;
      }
 
      for (int i=0; i<n; ++i) {
        sf.push_back(value_type(ord,i));
        this->order_ = std::max(this->order_,sf.back().order());
      }
 
      // Use equidistant "interpolation" nodes. To enable a faithful
      // reconstruction during grid transfer, we use interpolation nodes
      // for the maximal ACTUAL order of the shape functions and a least
      // squares "interpolation" approach.
      //
      // TODO: Since this is only an *extension*, we usually have fewer
      // shape functions. Check whether this can be exploited by using
      // fewer interpolation nodes.
      int m = binomial(dimension+this->order_,dimension);
      this->iNodes.resize(m);
      Dune::FieldVector<int,dimension> b(0);
      for (int i=0; i<m; ++i) {
        for (int j=0; j<dimension; ++j)
          this->iNodes[i][j] = static_cast<ctype>(b[j])/this->order_;
        next_multiindex(b.begin(),b.end(),this->order_);
      }
 
      // Compute the pseudoinverse for use in interpolation.
      computePseudoInverse();
 
      this->size_ = sf.size();
    }
 
    HierarchicExtensionSimplexShapeFunctionSet(HierarchicExtensionSimplexShapeFunctionSet const& other) = default;
 
    HierarchicExtensionSimplexShapeFunctionSet(HierarchicExtensionSimplexShapeFunctionSet&& other) = default;
 
    virtual value_type const& operator[](int i) const{ return sf[i]; }
 
    virtual Dune::GeometryType type() const
    {
      return Dune::GeometryType(Dune::GeometryType::simplex,dimension);
    }
 
    virtual void interpolate(typename ShapeFunctionSet<ctype,dimension,Scalar>::SfValueArray const& A,
        Matrix& IA) const
    {
      // TODO: implement without temporaries
      IA = pinv * A;
    }
 
    void removeShapeFunction(size_t index)
    {
      assert(index>=0 && index < sf.size());
      sf.erase(sf.begin()+index);
      this->size_ = sf.size();
      computePseudoInverse();
    }
 
  protected:
    std::vector<value_type> sf;
    DynamicMatrix<Dune::FieldMatrix<Scalar,1,1>> pinv;
 
  private:
    void computePseudoInverse()
    {
      DynamicMatrix<Scalar> A(this->iNodes.size(),sf.size());
      for (int j=0; j<sf.size(); ++j)
        for (int i=0; i<this->iNodes.size(); ++i)
          A[i][j] = sf[j].evaluateFunction(this->iNodes[i]);
 
      pinv = pseudoinverse(A);
    }
  };
 
  //---------------------------------------------------------------------
 
  template <class ctype, int dimension, class Scalar, bool restricted = false>
  class HierarchicExtensionShapeFunctionSetContainer: public ShapeFunctionSetContainer<ctype,dimension,Scalar>
  {
  public:
    typedef ShapeFunctionSet<ctype,dimension,Scalar> value_type;
 
  private:
    typedef std::map<std::pair<Dune::GeometryType,int>,value_type const*> Container;
    typedef typename HierarchicDetail::ChooseExtensionShapeFunctionSet<ctype,dimension,Scalar,restricted>::type ShapeFunctionType;
 
  public:
 
    HierarchicExtensionShapeFunctionSetContainer(int maxOrder)
    {
      // Fill the container with shape function sets. This is not done
      // on demand in operator(), since then operator() is not easily
      // and efficiently made thread-safe.
      //
      // Currently only shape functions for simplices are defined.
      Dune::GeometryType simplex(Dune::GeometryType::simplex,dimension);
 
      value_type* sfsPrevious = 0;
 
      for (int i=-1; i<=maxOrder; ++i) {
        value_type* sfs = new ShapeFunctionType(i);
        sfs->initHierarchicalProjection(sfsPrevious);
        container.insert(typename Container::value_type(std::make_pair(simplex,i),sfs));
        sfsPrevious = sfs;
      }
    }
 
    virtual ~HierarchicExtensionShapeFunctionSetContainer()
    {
      for (typename Container::iterator i=container.begin(); i!=container.end(); ++i)
        delete i->second;
    }
 
    virtual value_type const& operator()(Dune::GeometryType type, int order) const
    {
      typename Container::const_iterator i = container.find(std::make_pair(type,order));
      assert(i!=container.end()); // should throw, not assert!
      return *i->second;
    }
 
  private:
    // A map storing Hierarchic shape function sets associated with
    // (reference element type, order). Due to runtime polymorphism,
    // pointers are stored in the map. This class owns the shape
    // function sets pointed to, and has to delete them in the
    // destructor.
    Container container;
  };
 
 
  template <class ctype, int dimension, class Scalar>
  typename HierarchicExtensionShapeFunctionSetContainer<ctype,dimension,Scalar>::value_type const&
  hierarchicExtensionShapeFunctionSet(Dune::GeometryType type, int order)
  {
    // This static variable will be destructed on program shutdown,
    // causing its contents to be deleted properly from the heap by
    // the destructor. No memory leaks should be induced.
    static HierarchicExtensionShapeFunctionSetContainer<ctype,dimension,Scalar> container(7);
 
    return container(type,order);
  }
 
} // End of namespace Kaskade
 
#endif