nedelecshapefunctions.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-2016 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef NEDELECSHAPEFUNCTIONS_HH
#define NEDELECSHAPEFUNCTIONS_HH
 
#include <map>
 
#include <dune/common/fvector.hh>
#include <dune/geometry/referenceelements.hh>
#include <dune/geometry/type.hh>
 
#include "fem/barycentric.hh"
#include "fem/lagrangeshapefunctions.hh"
#include "fem/pshapefunctions.hh"
 
#include "linalg/dynamicMatrix.hh"
 
namespace Kaskade
{
  template <class ctype, int dimension, class T>
  class NedelecSimplexShapeFunction: public ShapeFunction<ctype,dimension,T,dimension>
  {
  public:
    static int const dim = dimension;
    static int const comps = dim;
    static int const order = 1;
 
    typedef ctype CoordType;
    typedef T ResultType;
 
    NedelecSimplexShapeFunction()
    {
      *this = NedelecSimplexShapeFunction<ctype,dimension,T>(0);
    }
 
    explicit NedelecSimplexShapeFunction(int e)
    {
      // compute index local in element
      entity_ = e;
 
      // Compute barycentric dimensions of edge end points. Note that the origin, which 
      // has barycentric index dim, has the subentity number 0. Hence we need to shift and
      // wrap the subentity indices to obtain the barycentric indices.
      Dune::ReferenceElement<CoordType,dim> const &simplex = Dune::ReferenceElements<CoordType,dim>::simplex();
      int p0Id = simplex.subEntity(e,dim-1,0,dim);
      p0_ = p0Id==0? dim: p0Id-1;
      int p1Id = simplex.subEntity(e,dim-1,1,dim);
      p1_ = p1Id==0? dim: p1Id-1;
 
      // compute gradients of barycentric coordinate functions
      if (p0_==dim) dp0_ = -1;
      else {        dp0_ = 0; dp0_[p0_] = 1; }
      if (p1_==dim) dp1_ = -1;
      else {        dp1_ = 0; dp1_[p1_] = 1; }
 
      // compute position and direction of edge
      pos_ = simplex.position(e,dim-1);
      edge_ = simplex.position(p1Id,dim) - simplex.position(p0Id,dim);
    }
    
    virtual Dune::FieldVector<T,dim> evaluateFunction(Dune::FieldVector<ctype,dim> const& x) const
    {
      Dune::FieldVector<CoordType,dim+1> bc = barycentric(x);
      return bc[p0_]*dp1_ - bc[p1_]*dp0_;
    }
 
    virtual Dune::FieldMatrix<ResultType,dim,dim> evaluateDerivative(Dune::FieldVector<CoordType,dim> const& x) const
    {
      // derivative of bc[p0]*dp1 - bc[p1]*dp0 is dp1 * bc[p0]' - dp0 * bc[p1]', where the products are *outer products*
      // resulting in the proper matrices. Start with derivatives of the barycentric coordinates \lambda_i and \lambda_j
      Dune::FieldVector<T,dim> dbc0(0.0), dbc1(0.0);
      for (int dir=0; dir<dim; ++dir) 
      {
        if (p0_ == dir) dbc0 = -1.0;
        else            dbc0[p0_] = 1.0;
        if (p1_ == dir) dbc1 = -1.0;
        else            dbc1[p1_] = 1.0;
      }
      
      return outerProduct(dp1_,dbc0) - outerProduct(dp0_,dbc1);
      
      // old code, should be equivlaent...
//       // TODO: faster implementation needed
//       Dune::FieldMatrix<ResultType,dim,dim> d;
//       for (int c=0; c<dim; ++c)
//         for (int dir=0; dir<dim; ++dir) {
//           ResultType dbcp0_ddir=0, dbcp1_ddir=0;
//           if (p0_==dir) dbcp0_ddir = 1;
//           if (p0_==dim) dbcp0_ddir = -1;
//           if (p1_==dir) dbcp1_ddir = 1;
//           if (p1_==dim) dbcp1_ddir = -1;
// 
//           d[c][dir] = dbcp0_ddir*dp1_[c] - dbcp1_ddir*dp0_[c];
//         }
// 
//       return d;
    }
 
    virtual Tensor<ResultType,1,dim,dim> evaluate2ndDerivative(Dune::FieldVector<CoordType,dim> const& x) const
    {
      // Nedelec shape functions are linear -- the second derivative vanishes.
      return Tensor<ResultType,1,dim,dim>(0.0);
    }
 
    virtual std::tuple<int,int,int,int> location() const
    {
      return std::make_tuple(1,dim-1,entity_,0);
    }
 
    Dune::FieldVector<CoordType,dim> position() const { return pos_; }
 
    Dune::FieldVector<CoordType,dim> edge() const { return edge_; }
 
  private:
    int entity_;
    // p0_ and p1_ contain the barycentric coordinate into which the
    // "from" and "to" points of the edge lie, respectively.
    int p0_, p1_;
    Dune::FieldVector<ResultType,dim> dp0_, dp1_;
    Dune::FieldVector<CoordType,dim> pos_, edge_; 
  };
 
  //--------------------------------------------------
 
  template <class ctype, int dimension, class T>
  class NedelecSimplexShapeFunctionSet: public ShapeFunctionSet<ctype,dimension,T,dimension>
  {
    static int const dim = dimension;
    static int const noEdges = dim*(dim+1)/2; // number of edges
 
 
  public:
    typedef NedelecSimplexShapeFunction<ctype,dim,T> value_type;
 
    NedelecSimplexShapeFunctionSet():
      ShapeFunctionSet<ctype,dim,T,dim>(Dune::GeometryType(Dune::GeometryType::simplex,dim))
      {
     auto const& simplex = Dune::ReferenceElements<ctype,dim>::simplex();
      assert(noEdges == simplex.size(dim-1));
 
      sf.reserve(noEdges);
      for (int i=0; i<noEdges; ++i)
        sf.push_back(value_type(i));
 
      // The interpolation uses the fact that for any edge vector e_i of
      // the reference simplex and any shape function phi_j evaluated on
      // this edge it holds that <e_i,phi_j>=0 if i!=j. Thus we first
      // store the tangential component <e_i,phi_i> for all shape
      // functions.
      this->iNodes.reserve(noEdges);
      tangentialComponents.reserve(noEdges);
      for (size_t i=0; i<noEdges; ++i) {
        this->iNodes.push_back(sf[i].position());
        tangentialComponents.push_back(sf[i].edge()*sf[i].evaluateFunction(sf[i].position()));
      }
 
      this->order_ = 1;
      this->size_ = sf.size();
      }
 
    virtual value_type const& operator[](int i) const{ return sf[i]; }
    virtual Dune::GeometryType type() const { return Dune::GeometryType(Dune::GeometryType::simplex,dim); }
 
    virtual void interpolate(typename ShapeFunctionSet<ctype,dim,T,dim>::SfValueArray const& A,
                             DynamicMatrix<Dune::FieldMatrix<T,1,1> >& IA) const
    {
      assert(A.N()==noEdges);
 
      IA.setSize(noEdges,A.M());
      for (int i=0; i<noEdges; ++i)
        for (int j=0; j<A.M(); ++j)
          IA[i][j] = sf[i].edge()*asVector(A[i][j])/tangentialComponents[i];
    }
 
 
  private:
    std::vector<value_type> sf;
    std::vector<T>          tangentialComponents;
  };
 
  //---------------------------------------------------------------------
  //---------------------------------------------------------------------
 
  template <class ctype, int dimension, class T>
  class NedelecShapeFunctionSetContainer: public ShapeFunctionSetContainer<ctype,dimension,T,dimension>
  {
  public:
    typedef ShapeFunctionSet<ctype,dimension,T,dimension> value_type;
 
  private:
    typedef std::map<std::pair<Dune::GeometryType,int>,value_type const*> Container;
 
  public:
 
    NedelecShapeFunctionSetContainer()
    {
      // 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);
      container.insert(typename Container::value_type(std::make_pair(simplex,1),new NedelecSimplexShapeFunctionSet<ctype,dimension,T>));
    }
 
 
    virtual ~NedelecShapeFunctionSetContainer()
    {
      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 Nedelec 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 T>
  typename NedelecShapeFunctionSetContainer<ctype,dimension,T>::value_type const& nedelecShapeFunctionSet(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 NedelecShapeFunctionSetContainer<ctype,dimension,T> container;
    return container(type,order);
  }
  
  //-------------------------------------------------------------------------------
  //-------------------------------------------------------------------------------
  
  template <class ctype, int dimension, class T>
  class HdivSimplexShapeFunction: public ShapeFunction<ctype,dimension,T,dimension>
  {
  public:
    static int const dim = dimension;
    static int const comps = dim;
    static int const order = 1;
 
    typedef ctype CoordType;
 
    HdivSimplexShapeFunction() = default;
 
    template <class Coefficients>
    explicit HdivSimplexShapeFunction(Coefficients const& coeff_, std::tuple<int,int,int,int> loc_)
    : coeff(coeff_), loc(loc_), lsfs(lagrangeShapeFunctionSet<ctype,dimension>(Dune::GeometryType(Dune::GeometryType::simplex,dim),std::get<0>(loc)))
    {
      assert(coeff.size() == lsfs.size());
    }
    
    virtual std::unique_ptr<ShapeFunction<ctype,dim,T,dim>> clone() const
    {
      return std::unique_ptr<ShapeFunction<ctype,dim,T,dim>>(new HdivSimplexShapeFunction(*this));
    }
 
    virtual Dune::FieldVector<T,dim> evaluateFunction(Dune::FieldVector<ctype,dim> const& xi) const
    {
      Dune::FieldVector<T,dim> phi(0);
      for (int i=0; i<coeff.size(); ++i)
          phi += lsfs[i].evaluateFunction(xi)[0] * coeff[i];
      return phi;
    }
 
    virtual Dune::FieldMatrix<T,dim,dim> evaluateDerivative(Dune::FieldVector<CoordType,dim> const& xi) const
    {
      Dune::FieldMatrix<T,dim,dim> dphi(0);
      for (int i=0; i<coeff.size(); ++i)
      {
        auto lsf = lsfs[i].evaluateDerivative(xi)[0]; // has only one component
        for (int j=0; j<dim; ++j)
          dphi[j] +=  coeff[i][j] * lsf;
      }
      return dphi;
    }
 
    virtual Tensor<T,dim,dim,dim> evaluate2ndDerivative(Dune::FieldVector<CoordType,dim> const& xi) const
    {
      Tensor<T,dim,dim,dim> ddphi(0);
      for (int i=0; i<coeff.size(); ++i)
      {
        auto lsf = lsfs[i].evaluate2ndDerivative(xi)[0];// has only one component
        for (int j=0; j<dim; ++j)
          ddphi[j] +=  coeff[i][j] * lsf;
      }
      return ddphi;
    }
 
    virtual std::tuple<int,int,int,int> location() const
    {
      return loc;
    }
 
  private:
    std::vector<Dune::FieldVector<T,comps>>   coeff;  // coefficients of Lagrangian shape functions.
    std::tuple<int,int,int,int> loc;                  // location information
    ShapeFunctionSet<ctype,dim,T> const&  lsfs;       // Lagrangian shape function set
  };
 
  //--------------------------------------------------
 
  template <class ctype, int dimension, class T=double>
  class HdivSimplexShapeFunctionSet: public ShapeFunctionSet<ctype,dimension,T,dimension>
  {
    static int const dim = dimension; 
 
 
  public:
    typedef HdivSimplexShapeFunction<ctype,dim,T> value_type;
 
    HdivSimplexShapeFunctionSet(int order);
    
    virtual value_type const& operator[](int i) const{ return sf[i]; }
    
    virtual Dune::GeometryType type() const { return Dune::GeometryType(Dune::GeometryType::simplex,dim); }
 
    virtual void interpolate(typename ShapeFunctionSet<ctype,dim,T,dim>::SfValueArray const& A,
                             DynamicMatrix<Dune::FieldMatrix<T,1,1> >& IA) const
    {
      std::cerr << "not implemented\n";
      abort(); // to be implemented
    }
 
 
  private:
    std::vector<value_type> sf;
  };
 
  //---------------------------------------------------------------------
  //---------------------------------------------------------------------
 
  template <class ctype, int dimension, class T>
  class HdivShapeFunctionSetContainer: public ShapeFunctionSetContainer<ctype,dimension,T,dimension>
  {
  public:
    typedef ShapeFunctionSet<ctype,dimension,T,dimension> value_type;
 
  private:
    typedef std::map<std::pair<Dune::GeometryType,int>,value_type const*> Container;
 
  public:
 
    HdivShapeFunctionSetContainer()
    {
      // 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);
      for (int order=1; order<6; ++order)
        container.insert(typename Container::value_type(std::make_pair(simplex,order),new HdivSimplexShapeFunctionSet<ctype,dimension,T>(order)));
    }
 
 
    virtual ~HdivShapeFunctionSetContainer()
    {
      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 Hdiv 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 T=double>
  typename HdivShapeFunctionSetContainer<ctype,dimension,T>::value_type const& hdivShapeFunctionSet(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 HdivShapeFunctionSetContainer<ctype,dimension,T> container;
    return container(type,order);
  }
  
} // end of namespace Kaskade
 
 
#endif