localtoglobaladaptors.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
 
#ifndef DUNE_LOCALFUNCTIONS_COMMON_LOCALTOGLOBALADAPTORS_HH
#define DUNE_LOCALFUNCTIONS_COMMON_LOCALTOGLOBALADAPTORS_HH
 
#include <cstddef>
#include <vector>
 
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/common/typetraits.hh>
 
#include <dune/geometry/type.hh>
 
namespace Dune {
 
 
  template<class LocalBasisTraits, std::size_t dimDomainGlobal_>
  struct LocalToGlobalBasisAdaptorTraits {
    typedef typename LocalBasisTraits::DomainFieldType DomainField;
    static const std::size_t dimDomainLocal = LocalBasisTraits::dimDomain;
    static const std::size_t dimDomainGlobal = dimDomainGlobal_;
    typedef typename LocalBasisTraits::DomainType DomainLocal;
    typedef FieldVector<DomainField, dimDomainGlobal> DomainGlobal;
 
    typedef typename LocalBasisTraits::RangeFieldType RangeField;
    static const std::size_t dimRange = LocalBasisTraits::dimRange;
    typedef typename LocalBasisTraits::RangeType Range;
 
    typedef FieldMatrix<RangeField, dimRange, dimDomainGlobal> Jacobian;
  };
 
 
  template<class LocalBasis, class Geometry>
  class ScalarLocalToGlobalBasisAdaptor {
    static_assert(LocalBasis::Traits::dimRange == 1,
                  "ScalarLocalToGlobalBasisAdaptor can only wrap a "
                  "scalar local basis.");
    static_assert((std::is_same<typename LocalBasis::Traits::DomainFieldType,
                           typename Geometry::ctype>::value),
                   "ScalarLocalToGlobalBasisAdaptor: LocalBasis must use "
                   "the same ctype as Geometry");
    static_assert
      ( static_cast<std::size_t>(LocalBasis::Traits::dimDomain) ==
      static_cast<std::size_t>(Geometry::mydimension),
      "ScalarLocalToGlobalBasisAdaptor: LocalBasis domain dimension must "
      "match local dimension of Geometry");
 
    const LocalBasis& localBasis;
    Geometry geometry;
 
  public:
    typedef LocalToGlobalBasisAdaptorTraits<typename LocalBasis::Traits,
        Geometry::coorddimension> Traits;
 
 
    ScalarLocalToGlobalBasisAdaptor(const LocalBasis& localBasis_,
                                    const Geometry& geometry_) :
      localBasis(localBasis_), geometry(geometry_)
    { }
 
    std::size_t size() const { return localBasis.size(); }
 
    std::size_t order() const {
      if(geometry.affine())
        // affine linear
        return localBasis.order();
      else
        // assume at most order dim
        return localBasis.order() + Traits::dimDomainGlobal - 1;
    }
 
    void evaluateFunction(const typename Traits::DomainLocal& in,
                          std::vector<typename Traits::Range>& out) const
    {
      localBasis.evaluateFunction(in, out);
    }
 
    void evaluateJacobian(const typename Traits::DomainLocal& in,
                          std::vector<typename Traits::Jacobian>& out) const
    {
      std::vector<typename LocalBasis::Traits::JacobianType>
      localJacobian(size());
      localBasis.evaluateJacobian(in, localJacobian);
 
      const typename Geometry::JacobianInverseTransposed &geoJacobian =
        geometry.jacobianInverseTransposed(in);
 
      out.resize(size());
      for(std::size_t i = 0; i < size(); ++i)
        geoJacobian.mv(localJacobian[i][0], out[i][0]);
    }
  };
 
 
  template<class LocalInterpolation, class Traits_>
  class LocalToGlobalInterpolationAdaptor {
    const LocalInterpolation& localInterpolation;
 
  public:
    typedef Traits_ Traits;
 
 
    LocalToGlobalInterpolationAdaptor
      ( const LocalInterpolation& localInterpolation_) :
      localInterpolation(localInterpolation_)
    { }
 
    template<class Function, class Coeff>
    void interpolate(const Function& function, std::vector<Coeff>& out) const
    { localInterpolation.interpolate(function, out); }
  };
 
 
  template<class LocalFiniteElement, class Geometry>
  struct ScalarLocalToGlobalFiniteElementAdaptor {
    struct Traits {
      typedef ScalarLocalToGlobalBasisAdaptor<typename LocalFiniteElement::
          Traits::LocalBasisType, Geometry> Basis;
      typedef LocalToGlobalInterpolationAdaptor<typename LocalFiniteElement::
          Traits::LocalInterpolationType, typename Basis::Traits>
      Interpolation;
      typedef typename LocalFiniteElement::Traits::LocalCoefficientsType
      Coefficients;
    };
 
  private:
    const LocalFiniteElement &localFE;
    typename Traits::Basis basis_;
    typename Traits::Interpolation interpolation_;
 
  public:
 
    ScalarLocalToGlobalFiniteElementAdaptor
      ( const LocalFiniteElement& localFE_, const Geometry &geometry) :
      localFE(localFE_),
      basis_(localFE.localBasis(), geometry),
      interpolation_(localFE.localInterpolation())
    { }
 
    const typename Traits::Basis& basis() const { return basis_; }
    const typename Traits::Interpolation& interpolation() const
    { return interpolation_; }
    const typename Traits::Coefficients& coefficients() const
    { return localFE.localCoefficients(); }
    GeometryType type() const { return localFE.type(); }
  };
 
 
  template<class LocalFiniteElement, class Geometry>
  class ScalarLocalToGlobalFiniteElementAdaptorFactory {
    const LocalFiniteElement& localFE;
 
  public:
    typedef ScalarLocalToGlobalFiniteElementAdaptor<LocalFiniteElement,
        Geometry> FiniteElement;
 
 
    ScalarLocalToGlobalFiniteElementAdaptorFactory
      (const LocalFiniteElement &localFE_) : localFE(localFE_) {}
 
 
    const FiniteElement make(const Geometry& geometry) {
      return FiniteElement(localFE, geometry);
    }
  };
 
} // namespace Dune
 
#endif // DUNE_LOCALFUNCTIONS_COMMON_LOCALTOGLOBALADAPTORS_HH