pshapefunctions.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-2024 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef PSHAPEFUNCTIONS_HH
#define PSHAPEFUNCTIONS_HH
 
#include "fem/barycentric.hh"
#include "fem/fixdune.hh"
#include "fem/mllgeometry.hh"
#include "linalg/dynamicMatrix.hh"
#include "linalg/dynamicMatrixOps.hh"
#include "linalg/tensor.hh"
 
#include <boost/multi_array.hpp>
 
#include <tuple>
#include <vector>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <memory>
 
 
 
namespace Kaskade
{
  template <class ctype, int dim, class T=double, int comp=1> 
  class ShapeFunction
  {
  public:
    virtual ~ShapeFunction() {}
 
    virtual Dune::FieldVector<T,comp> evaluateFunction(Dune::FieldVector<ctype,dim> const& xi) const = 0;
 
    virtual Dune::FieldMatrix<T,comp,dim> evaluateDerivative(Dune::FieldVector<ctype,dim> const& xi) const = 0;
 
    virtual Tensor<T,comp,dim,dim> evaluate2ndDerivative(Dune::FieldVector<ctype,dim> const& xi) const
    {
      std::cerr << "NOT IMPLEMENTED: " << __FILE__ << ":" << __LINE__ << "\n"; 
      abort();
    }
 
    virtual std::tuple<int,int,int,int> location() const = 0;
  };
 
  //---------------------------------------------------------------------
  template <class ctype, int dimension, class T, int comp=1>   
  class ShapeFunctionSet
  {
  public:
    typedef T Scalar;
    typedef ShapeFunction<ctype,dimension,T,comp> value_type;
    
    typedef DynamicMatrix<Dune::FieldMatrix<T,1,1>> Matrix;
 
    ShapeFunctionSet(Dune::GeometryType gt)
    : order_(-1), size_(-1),
      gt_(gt)
    {}
 
    ShapeFunctionSet(ShapeFunctionSet const& other) = default;
 
    static int const comps = comp;
 
 
    virtual ~ShapeFunctionSet() {}
 
 
    virtual value_type const& operator[](int i) const = 0;
 
    virtual int size() const { return size_; }
 
    int order() const { return order_; }
 
    Dune::GeometryType type() const { return gt_; }
 
    auto referenceElement() const
    {
      return Dune::referenceElement<ctype,dimension>(gt_);
    }
 
    typedef std::vector<Dune::FieldVector<ctype,dimension>> InterpolationNodes;
 
    typedef DynamicMatrix<Dune::FieldMatrix<T,comp,1>> SfValueArray;
 
    void initHierarchicalProjection(ShapeFunctionSet<ctype,dimension,T,comp> const* sfl)
    {
      // If no lower order shape function set is given, we assume an
      // empty set, i.e. the projection is just zero.
      if (!sfl || size()==0) {
        projection.setSize(size(),size());
//       projection = 0 // causes compilation error with dune-2.4.0 and clang++ on OS X (Darwin)
        projection.fill(0);
        return;
      }
 
      // Compute the hierarchic projection P = I_h E_l(x_h) I_l E_h(x_l),
      // where E_h(x_l) the evaluation of high order shape functions at
      // low order nodes, I_l the low order interpolation, E_l(x_h) the
      // evaluation of low order shape functions at high order nodes,
      // and I_h the high order shape function evaluation.
 
      // Compute the local restriction matrix A
      SfValueArray Ehxl;
      evaluate(sfl->interpolationNodes(),Ehxl);  // E_h(x_l)
      Matrix A;
      sfl->interpolate(Ehxl,A);                  // I_l  E_h(x_l)
 
      // Compute the local prolongation matrix B
      SfValueArray Elxh;
      sfl->evaluate(interpolationNodes(),Elxh);  // E_l(x_h)
      Matrix B;
      interpolate(Elxh,B);                       // I_h E_l(x_h)
 
      // Compute projection P=B*A
      projection = B*A;
    }
 
 
 
    InterpolationNodes const& interpolationNodes() const
    {
      return iNodes;
    }
 
 
    virtual void interpolate(SfValueArray const& A, Matrix& IA) const = 0;
 
    void evaluate(InterpolationNodes const& iNodes, SfValueArray& phi) const
    {
      int s = size();
      phi.setSize(iNodes.size(),s);
 
      for (int j=0; j<s; ++j) {
        value_type const& sf = (*this)[j];
        for (int i=0; i<iNodes.size(); ++i) {
          Dune::FieldVector<T,comp> v = sf.evaluateFunction(iNodes[i]);
          for (int k=0; k<comp; ++k)
            phi[i][j][k] = v[k];  // todo: improve efficiency! cache directly!
        }
      }
    }
 
    Matrix const& hierarchicProjection() const
    {
      // check that the projection has been properly initialized.
      assert(projection.N()==size() && projection.M()==size());
 
      return projection;
    }
 
    virtual void removeShapeFunction(size_t index) {}
 
  protected:
    InterpolationNodes iNodes;
    Matrix             projection;
    int                order_;
    int                size_;
 
  private:
    Dune::GeometryType gt_;
  };
 
  //---------------------------------------------------------------------
   
  template <class ctype, int dimension, class T, int comp=1>   
  class EmptyShapeFunctionSet: public ShapeFunctionSet<ctype,dimension,T,comp>
  {
    using Base = ShapeFunctionSet<ctype,dimension,T,comp>;
  public:
    EmptyShapeFunctionSet(Dune::GeometryType const& gt)
    : Base(gt)
    {
      this->size_ = 0;
    }
    virtual typename Base::value_type const& operator[](int i) const { abort(); }
    virtual void interpolate(typename Base::SfValueArray const& A, typename Base::Matrix& IA) const 
    { 
      IA.resize(0,A.M());
    }
  };
   
  //---------------------------------------------------------------------
 
  template <class ShapeFunctionSet_>
  class RestrictedShapeFunctionSet: public ShapeFunctionSet_
  {
  public:
    //typedef typename ShapeFunctionSet_::Grid Grid; causes compiler error
    typedef typename ShapeFunctionSet_::Scalar Scalar;
    typedef typename ShapeFunctionSet_::value_type value_type;
 
    explicit RestrictedShapeFunctionSet(int order) : ShapeFunctionSet_(order) {}
 
    RestrictedShapeFunctionSet(ShapeFunctionSet_ const& other) : ShapeFunctionSet_(other) {}
    RestrictedShapeFunctionSet(RestrictedShapeFunctionSet const& other) : ShapeFunctionSet_(other) {}
    RestrictedShapeFunctionSet(RestrictedShapeFunctionSet const& other, std::vector<int>* ids_) : ShapeFunctionSet_(other) {}
 
    virtual ~RestrictedShapeFunctionSet(){}
 
    void setRestriction(std::vector<int> const& ids_)
    {
      if(ids_.size()==0)
      {
        this->sf.clear();
        this->iNodes.clear();
        this->size_ = 0;
        return;
      }
 
 
      for(int i=this->sf.size()-1; i>-1; --i)
        if(std::find(ids_.begin(),ids_.end(),i)==ids_.end()) this->removeShapeFunction(i);
 
      this->size_ = ids_.size();
    }
  };
 
  //---------------------------------------------------------------------
 
  template <class ctype, int dimension, class T, int comp=1>
  class ShapeFunctionSetContainer
  {
  public:
    typedef ShapeFunctionSet<ctype,dimension,T,comp> value_type;
 
    virtual ~ShapeFunctionSetContainer() {}
 
 
    virtual value_type const& operator() (Dune::GeometryType type, int order) const = 0;
  };
 
  //---------------------------------------------------------------------
  //---------------------------------------------------------------------
 
  template <class ShapeFunctionSet>
  void computeSimplexSfPermutation(ShapeFunctionSet const& sfs,
      int const* vPermutation,
      int*       sfPermutation,
      int*       sign)
  {
    assert(vPermutation);
    assert(sfPermutation);
    assert(sign);
 
    // Ok, we rely on the evaluation & interpolation of the shape
    // function set. First we evaluate the shape functions on the
    // interpolation nodes transformed by f, then we interpolate these
    // values. This gives us a matrix which, if the shape function set
    // indeed allows simple coupling, is a signed permutation
    // matrix. Finally we extract the signs and the permutation from
    // this matrix.
 
    // Get the interpolation nodes. Make a copy because we need to
    // modify them.
    typedef typename ShapeFunctionSet::InterpolationNodes InterpolationNodes;
    InterpolationNodes in = sfs.interpolationNodes();
 
 
    // Transform the nodes. This is done by permuting their barycentric
    // coordinates. Remember that the interpolation nodes are given in
    // Cartesian, not in barycentric coordinates. This determines the
    // loop termination for j.
    for (int i=0; i<in.size(); ++i) {
      Dune::FieldVector<typename ShapeFunctionSet::Grid::ctype,ShapeFunctionSet::Grid::dimension+1> b = barycentric(in[i]);
      for (int j=0; j<ShapeFunctionSet::Grid::dimension; ++j)
        in[i][j] = b[vPermutation[j]];
    }
 
 
    // Evaluate the shape functions at the transformed nodes.
    typename ShapeFunctionSet::SfValueArray phi(in.size(),sfs.size());
    sfs.evaluate(in,phi);
 
    // Interpolate the shape function values.
    typename ShapeFunctionSet::Matrix P(sfs.size(),sfs.size());
    sfs.interpolate(phi,P);
 
 
#ifndef NDEBUG
    // Check that the permutation matrix is indeed a signed permutation
    bool ok = true;
    for (int i=0; i<sfs.size(); ++i) {
      double rowSum = 0;
      double colSum = 0;
      for (int j=0; j<sfs.size(); ++j) {
        if (std::abs(P[i][j])>1e-8 && std::abs(1-std::abs(P[i][j]))>1e-8) ok = false;
        rowSum += std::abs(P[i][j]);
        colSum += std::abs(P[j][i]);
      }
      if (std::abs(rowSum-1)>1e-8 || std::abs(colSum-1)>1e-8) ok = false;
    }
    if (!ok) {
      std::cout << "\nPermutation matrix:\n";
      for (int i=0; i<P.N(); ++i) {
        for (int j=0; j<P.M(); ++j)
          std::cout << std::setw(2) << P[i][j] << " ";
        std::cout << "\n";
      }
    }
#endif
 
 
    // Extract the permutation and the sign.
    for (int i=0; i<sfs.size(); ++i)
      for (int j=0; j<sfs.size(); ++j)
        if (std::abs(P[j][i])>0.5) {
          sfPermutation[i] = j; // XXX or the other way round?
          sign[i] = P[j][i]>0? 1: -1;
          break;
        }
 
  }
} // end of namespace Kaskade
 
#endif