crouzeixraviart.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_LOCALFUNCTIONS_CROUZEIXRAVIART_HH
#define DUNE_LOCALFUNCTIONS_CROUZEIXRAVIART_HH
 
#include <array>
#include <numeric>
 
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
 
#include <dune/geometry/type.hh>
#include <dune/geometry/referenceelements.hh>
 
#include <dune/localfunctions/common/localbasis.hh>
#include <dune/localfunctions/common/localfiniteelementtraits.hh>
#include <dune/localfunctions/common/localinterpolation.hh>
#include <dune/localfunctions/common/localkey.hh>
 
namespace Dune { namespace Impl
{
  template<class D, class R, unsigned int dim>
  class CrouzeixRaviartLocalBasis
  {
  public:
    using Traits = LocalBasisTraits<D,dim,FieldVector<D,dim>,R,1,FieldVector<R,1>,FieldMatrix<R,1,dim> >;
 
    static constexpr unsigned int size ()
    {
      return dim+1;
    }
 
    void evaluateFunction(const typename Traits::DomainType& x,
                          std::vector<typename Traits::RangeType>& out) const
    {
      out.resize(size());
 
      std::fill(out.begin(), out.end()-1, 1.0);
      out.back() = 1.0-dim;
 
      for (unsigned int i=0; i<dim; i++)
      {
        out[i] -= dim * x[dim-i-1];
        out.back() += dim*x[i];
      }
    }
 
    void evaluateJacobian(const typename Traits::DomainType& x,
                          std::vector<typename Traits::JacobianType>& out) const
    {
      out.resize(size());
 
      for (unsigned i=0; i<dim; i++)
        for (unsigned j=0; j<dim; j++)
          out[i][0][j] = (i==(dim-1-j)) ? -(double)dim : 0;
 
      std::fill(out.back()[0].begin(), out.back()[0].end(), dim);
    }
 
    void partial(const std::array<unsigned int,dim>& order,
                 const typename Traits::DomainType& in,
                 std::vector<typename Traits::RangeType>& out) const
    {
      auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
 
      out.resize(size());
 
      if (totalOrder == 0) {
        evaluateFunction(in, out);
        return;
      }
 
      if (totalOrder==1)
      {
        auto direction = std::find(order.begin(), order.end(), 1)-order.begin();
 
        for (unsigned int i=0; i<dim; i++)
          out[i] = (i==(dim-1-direction)) ? -(double)dim : 0.0;
 
        out.back()[0] = dim;
      }
      else  // all higher order derivatives are zero
        std::fill(out.begin(), out.end(), 0);
    }
 
    static constexpr unsigned int order ()
    {
      return 1;
    }
  };
 
  template<unsigned int dim>
  class CrouzeixRaviartLocalCoefficients
  {
  public:
    CrouzeixRaviartLocalCoefficients ()
    : localKeys_(size())
    {
      for (std::size_t i=0; i<size(); i++)
        localKeys_[i] = LocalKey(i,dim-1,0);
    }
 
    static constexpr std::size_t size ()
    {
      return dim+1;
    }
 
    const LocalKey& localKey (std::size_t i) const
    {
      return localKeys_[i];
    }
 
  private:
    std::vector<LocalKey> localKeys_;
  };
 
  template<class LocalBasis>
  class CrouzeixRaviartLocalInterpolation
  {
  public:
 
    template<typename F, typename C>
    void interpolate (const F& ff, std::vector<C>& out) const
    {
      constexpr auto dim = LocalBasis::Traits::dimDomain;
 
      auto&& f = Impl::makeFunctionWithCallOperator<typename LocalBasis::Traits::DomainType>(ff);
 
      out.resize(LocalBasis::size());
 
      // Evaluate at the facet midpoints
      auto refSimplex = ReferenceElements<typename LocalBasis::Traits::DomainFieldType,dim>::simplex();
      for (int i=0; i<refSimplex.size(1); i++)
        out[i] = f(refSimplex.template geometry<1>(i).center());
    }
  };
 
} }   // namespace Dune::Impl
 
namespace Dune
{
  template<class D, class R, int dim>
  class CrouzeixRaviartLocalFiniteElement
  {
  public:
    using Traits = LocalFiniteElementTraits<Impl::CrouzeixRaviartLocalBasis<D,R,dim>,
                                            Impl::CrouzeixRaviartLocalCoefficients<dim>,
                                            Impl::CrouzeixRaviartLocalInterpolation<Impl::CrouzeixRaviartLocalBasis<D,R,dim> > >;
 
    const typename Traits::LocalBasisType& localBasis() const
    {
      return basis_;
    }
 
    const typename Traits::LocalCoefficientsType& localCoefficients() const
    {
      return coefficients_;
    }
 
    const typename Traits::LocalInterpolationType& localInterpolation() const
    {
      return interpolation_;
    }
 
    static constexpr std::size_t size()
    {
      return dim+1;
    }
 
    static constexpr GeometryType type()
    {
      return GeometryTypes::simplex(dim);
    }
 
  private:
    Impl::CrouzeixRaviartLocalBasis<D,R,dim> basis_;
    Impl::CrouzeixRaviartLocalCoefficients<dim> coefficients_;
    Impl::CrouzeixRaviartLocalInterpolation<Impl::CrouzeixRaviartLocalBasis<D,R,dim> > interpolation_;
  };
 
}        // namespace Dune
 
#endif   // DUNE_LOCALFUNCTIONS_CROUZEIXRAVIART_HH