functionspace.hh

#ifndef DUNE_FEM_FUNCTIONSPACE_HH
#define DUNE_FEM_FUNCTIONSPACE_HH
 
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
 
#include <dune/geometry/dimension.hh>
 
#include <dune/fem/space/common/functionspaceinterface.hh>
 
namespace Dune
{
 
  namespace Fem
  {
 
    // Forward declaration of
    // base class for vector valued function spaces.
    template< class DomainField, class RangeField, int dimD, int dimR>
    class FunctionSpace;
 
    template< class DomainField, class RangeField, int dimD, int dimR>
    struct VectorSpaceTraits
    {
      typedef DomainField DomainFieldType;
      typedef RangeField RangeFieldType;
 
      enum { dimDomain = dimD };
      enum { dimRange = dimR };
 
      typedef FieldVector< DomainFieldType, dimDomain > DomainType;
 
      typedef FieldVector< RangeFieldType, dimRange> RangeType;
 
      typedef FieldMatrix< RangeFieldType, dimRange, dimDomain > LinearMappingType;
 
      typedef FunctionSpace< DomainFieldType, RangeFieldType, dimDomain, 1>
        ScalarFunctionSpaceType;
    };
 
    template< class DomainField, class RangeField, int dimD, int dimR>
    class FunctionSpace
    : public FunctionSpaceInterface
        < VectorSpaceTraits< DomainField, RangeField, dimD, dimR> >
    {
      typedef FunctionSpace< DomainField, RangeField, dimD, dimR> ThisType;
 
    public:
      typedef ThisType FunctionSpaceType;
    };
 
    /* Forward declaration of
      base class for matrix valued function spaces.
    */
    template <typename DomainFieldImp,typename RangeFieldImp,int n,int m1,int m2>
    class MatrixFunctionSpace;
 
    template <typename K,int n,int m>
    class RangeMatrix : public FieldMatrix<K,n,m> {
     public:
      typedef FieldMatrix<K,n,m> BaseType;
      typedef typename BaseType::row_type RowType;
      enum {
        rows = BaseType::rows,
        cols = BaseType::cols,
        dimension = BaseType::rows*BaseType::cols
      };
      //===== constructors
      RangeMatrix () : BaseType() {}
 
      RangeMatrix (const K& k) : BaseType(k) {}
 
      K& operator()(int r,int c) {
        return static_cast<BaseType&>(*this)[r][c];
      }
      const K operator()(int r,int c) const {
        return static_cast<const BaseType&>(*this)[r][c];
      }
      const RowType& row(int r) const{
        return static_cast<BaseType&>(*this)[r];
      }
      RowType& row(int r){
        return static_cast<BaseType&>(*this)[r];
      }
 
      K& operator[](int i) {
        int r = i/cols;
        int c = i%cols;
        return (*this)(r,c);
      }
      const K operator[](int i) const {
        int r = i/cols;
        int c = i%cols;
        return (*this)(r,c);
      }
 
      K operator* (const BaseType& y)
      {
        K ret(0);
        for (int i=0; i<n; i++)
          ret += static_cast<BaseType&>(*this)[i] * y[i];
        return ret;
      }
 
      RangeMatrix& axpy (const K& a, const BaseType& y)
      {
         for (int i=0; i<n; i++)
            static_cast<BaseType&>(*this)[i].axpy(a,y[i]);
         return *this;
      }
    };
    template <typename DomainFieldImp,typename RangeFieldImp,int n,int m1,int m2>
    class MatrixMapping :
        public FieldMatrix<RangeFieldImp,m1*m2,n> {
      // Implement L : VD -> VR where VR is the space of FieldMatrix<RFI,m1,m2>
      // VD is space of FieldVector<DFT,n>
     public:
      typedef DomainFieldImp DomainFieldType;
      typedef RangeFieldImp RangeFieldType;
      typedef RangeMatrix<RangeFieldImp, m1,m2> RangeType;
      typedef FieldMatrix<RangeFieldImp,m1*m2,n> BaseType;
      //===== constructors
      MatrixMapping () : BaseType() {}
      MatrixMapping (const RangeFieldImp& k) : BaseType(k) {}
 
      FieldVector<DomainFieldImp,n>& operator[](int i) {
        return static_cast<BaseType&>(*this)[i];
      }
      const FieldVector<DomainFieldImp,n>& operator[](int i) const {
        return static_cast<const BaseType&>(*this)[i];
      }
    };
 
    template <typename DomainFieldImp,typename RangeFieldImp,int n,int m1,int m2>
    struct MatrixSpaceTraits {
      typedef DomainFieldImp DomainFieldType;
      typedef RangeFieldImp RangeFieldType;
      typedef FieldVector<DomainFieldImp, n> DomainType;
      typedef RangeMatrix<RangeFieldImp, m1,m2> RangeType;
      typedef MatrixMapping<DomainFieldImp,RangeFieldImp, n, m1,m2> LinearMappingType;
      typedef MatrixFunctionSpace<DomainFieldImp,RangeFieldImp,n,1,1> ScalarFunctionSpaceType;
      enum { dimRange = m1 * m2};
      enum { dimDomain = n };
    };
 
    template <typename DomainFieldImp,typename RangeFieldImp,int n,int m1,int m2>
    class MatrixFunctionSpace : public
          FunctionSpaceInterface<MatrixSpaceTraits<DomainFieldImp,RangeFieldImp,n,m1,m2> > {};
 
    // function space converter objects
    // --------------------------------
 
    template < class FunctionSpaceImp, int newDimDomain >
    struct ToNewDimDomainFunctionSpace;
 
    template < class FunctionSpaceImp, int newDimRange >
    struct ToNewDimRangeFunctionSpace;
 
    template< class DomainFieldImp, class RangeFieldImp, int dimDomain, int dimRange, int newDimDomain >
    struct ToNewDimDomainFunctionSpace< FunctionSpace< DomainFieldImp, RangeFieldImp, dimDomain, dimRange >, newDimDomain >
    {
      typedef FunctionSpace< DomainFieldImp, RangeFieldImp, newDimDomain, dimRange> Type;
    };
 
    template< class DomainFieldImp, class RangeFieldImp, int n, int m1, int m2, int newDimDomain >
    struct ToNewDimDomainFunctionSpace< MatrixFunctionSpace< DomainFieldImp, RangeFieldImp, n, m1, m2 >, newDimDomain >
    {
      typedef MatrixFunctionSpace< DomainFieldImp, RangeFieldImp, newDimDomain, m1, m2 > Type;
    };
 
    template< class DomainFieldImp, class RangeFieldImp, int dimDomain, int dimRange, int newDimRange >
    struct ToNewDimRangeFunctionSpace< FunctionSpace< DomainFieldImp, RangeFieldImp, dimDomain, dimRange >, newDimRange >
    {
      typedef FunctionSpace< DomainFieldImp, RangeFieldImp, dimDomain, newDimRange > Type;
    };
 
 
 
    // GridFunctionSpace
    // -----------------
 
    namespace Impl
    {
 
      template< class GridPart, class >
      struct GridFunctionSpace;
 
      template< class GridPart, class K, int dimRange >
      struct GridFunctionSpace< GridPart, FunctionSpace< typename GridPart::ctype, K, GridPart::dimensionworld, dimRange > >
      {
        typedef FunctionSpace< typename GridPart::ctype, K, GridPart::dimensionworld, dimRange > Type;
      };
 
      template< class GridPart, class K, int rows, int cols >
      struct GridFunctionSpace< GridPart, MatrixFunctionSpace< typename GridPart::ctype, K, GridPart::dimensionworld, rows, cols > >
      {
        typedef MatrixFunctionSpace< typename GridPart::ctype, K, GridPart::dimensionworld, rows, cols > Type;
      };
 
      template< class GridPart, class K, int dimRange >
      struct GridFunctionSpace< GridPart, Dune::FieldVector< K, dimRange > >
      {
        typedef FunctionSpace< typename GridPart::ctype, K, GridPart::dimensionworld, dimRange > Type;
      };
 
      template< class GridPart, class K, int rows, int cols >
      struct GridFunctionSpace< GridPart, Dune::FieldMatrix< K, rows, cols > >
      {
        typedef MatrixFunctionSpace< typename GridPart::ctype, K, GridPart::dimensionworld, rows, cols > Type;
      };
 
      template< class GridPart, int dimRange >
      struct GridFunctionSpace< GridPart, Dune::Dim< dimRange > >
      {
        typedef typename GridFunctionSpace< GridPart, Dune::FieldVector< typename GridPart::ctype, dimRange > >::Type Type;
      };
 
    } // namespace Impl
 
    template< class GridPart, class T >
    using GridFunctionSpace = typename Impl::GridFunctionSpace< GridPart, T >::Type;
 
  } // namespace Fem
 
} // namespace Dune
 
#endif // #ifndef DUNE_FEM_FUNCTIONSPACE_HH