vectorexpression.hh

#ifndef __DUNE_ACFEM_VECTOREXPRESSION_HH__
#define __DUNE_ACFEM_VECTOREXPRESSION_HH__
 
#include "expressionoperations.hh"
 
namespace Dune
{
 
  namespace ACFem
  {
 
    template<class ExpressionImp>
    class VectorExpression
      : public ExpressionTemplate<ExpressionImp>
    {
    };
 
    template<class BinOP, class Vector1, class Vector2>
    class BinaryVectorExpression;
 
    template<class OP, class Vector>
    class UnaryVectorExpression;
 
    template<class BinOp, class Vector1, class Vector2>
    class BinaryVectorExpression
      : public VectorExpression<BinaryVectorExpression<BinOp, Vector1, Vector2> >
    {
      typedef BinOp BinOpType;
      typedef
      ExpressionResult<BinOpType,
                       decltype(std::declval<Vector1>()[0]),
                       decltype(std::declval<Vector2>()[0])>
      ResultType;
     public:
      static_assert(std::is_same<typename Vector1::FieldType, typename Vector2::FieldType>::value,
                    "Cannot combine vector expression over different fields.");
 
      typedef typename Vector1::FieldType FieldType;
      typedef typename ResultType::Type ElementType;
 
      BinaryVectorExpression(const VectorExpression<Vector1>& v1, const VectorExpression<Vector2>& v2)
        : vector1_(v1), vector2_(v2)
      {};
 
      ElementType operator[](size_t index) const
      {
        return ResultType::apply(vector1_()[index], vector2_()[index]);
      }
 
      size_t size() const {
        assert(vector1_().size() == vector2_().size());
        return vector1_().size();
      }
 
     protected:
      ExpressionStorage<Vector1> vector1_;
      ExpressionStorage<Vector2> vector2_;
    };
 
    template<class Scalar, class Vector>
    class BinaryVectorExpression<SMultiplyOperation, Scalar, Vector>
      : public VectorExpression<BinaryVectorExpression<SMultiplyOperation, Scalar, Vector> >
    {
      typedef SMultiplyOperation BinOpType;
      typedef
      ExpressionResult<BinOpType,
                       decltype(parameterValue(std::declval<Scalar>())),
                       decltype(std::declval<Vector>()[0])>
      ResultType;
     public:
      typedef typename Vector::FieldType FieldType;
      typedef typename ResultType::Type ElementType;
 
      BinaryVectorExpression(const Scalar& s, const VectorExpression<Vector>& v)
        : scalar_(s), vector_(v)
      {};
 
      ElementType operator[](size_t index) const
      {
        return ResultType::apply(parameterValue(scalar_()), vector_()[index]);
      }
 
      size_t size() const {
        return vector_().size();
      }
 
     protected:
      ExpressionStorage<Scalar> scalar_;
      ExpressionStorage<Vector> vector_;
    };
 
    template<class Op, class Vector>
    class UnaryVectorExpression
      : public VectorExpression<UnaryVectorExpression<Op, Vector> >
    {
      typedef Op OpType;
      typedef
      typename std::remove_reference<decltype(std::declval<Vector>()[0])>::type
      ComponentType;
      typedef ExpressionResult<OpType, ComponentType> ResultType;
     public:
      typedef typename Vector::FieldType FieldType;
      typedef typename ResultType::Type ElementType;
 
      UnaryVectorExpression(const Vector& v)
        : vector_(v)
      {};
 
      ElementType operator[](size_t index) const
      {
        return ResultType::apply(vector_()[index]);
      }
 
      size_t size() const {
        return vector_().size();
      }
 
     protected:
      ExpressionStorage<Vector> vector_;
    };
 
    template<class Vector>
    class UnaryVectorExpression<IdentityOperation, Vector>
      : public VectorExpression<UnaryVectorExpression<IdentityOperation, Vector> >
    {
      typedef IdentityOperation OpType;
      typedef
      typename std::remove_reference<decltype(std::declval<Vector>()[0])>::type
      ComponentType;
     public:
      typedef typename FieldTraits<ComponentType>::field_type FieldType;
     private:
      enum {
        hasFieldComponents = std::is_same<ComponentType, FieldType>::value
      };
     public:
      typedef
      typename std::conditional<hasFieldComponents,
                                FieldType,
                                UnaryVectorExpression<IdentityOperation,
                                                      ComponentType>
                                >::type
      ElementType;
 
      UnaryVectorExpression(const Vector& v)
        : vector_(v)
      {};
 
      ElementType operator[](size_t index) const
      {
        return ElementType(vector_()[index]);
      }
 
      size_t size() const {
        return vector_().size();
      }
 
     protected:
      ExpressionStorage<Vector> vector_;
    };
 
    template<class Field, int rows, int cols>
    class UnaryVectorExpression<IdentityOperation, FieldMatrix<Field, rows, cols> >
      : public
    VectorExpression<
      UnaryVectorExpression<IdentityOperation, FieldMatrix<Field, rows, cols> > >
    {
      typedef IdentityOperation OpType;
      typedef FieldMatrix<Field, rows, cols> MatrixType;
      typedef typename MatrixType::row_type RowType;
     public:
      typedef typename FieldTraits<Field>::field_type FieldType;
      typedef UnaryVectorExpression<IdentityOperation, RowType> ElementType;
 
      UnaryVectorExpression(const MatrixType& m)
        : matrix_(m)
      {};
 
      ElementType operator[](size_t index) const
      {
        return ElementType(matrix_()[index]);
      }
 
      size_t size() const {
        return matrix_().size();
      }
 
     protected:
      ExpressionStorage<MatrixType> matrix_;
    };
 
    template<class Expr1, class Expr2>
    static inline
    BinaryVectorExpression<PlusOperation, Expr1, Expr2>
    operator+(const VectorExpression<Expr1>& e1, const VectorExpression<Expr2>& e2)
    {
      typedef
        BinaryVectorExpression<PlusOperation, Expr1, Expr2>
        ExpressionType;
 
      return ExpressionType(e1, e2);
    }
 
    template<class Expr1, class Expr2>
    static inline
    BinaryVectorExpression<MinusOperation, Expr1, Expr2>
    operator-(const VectorExpression<Expr1>& e1, const VectorExpression<Expr2>& e2)
    {
      typedef
        BinaryVectorExpression<MinusOperation, Expr1, Expr2>
        ExpressionType;
 
      return ExpressionType(e1, e2);
    }
 
    namespace {
      // Helper for scalar products. To find an end in this recursion
      // operator*() evaluates to a Field element if the components of
      // the expressions evaluate to a field element.
      template<class Expr1, class Expr2>
      struct InnerProductHelper
      {
        template<bool haveFieldVectors, class dummy = void>
        struct ProductFactory;
 
        template<class dummy>
        struct ProductFactory<true, dummy>
        {
          typedef typename Expr1::FieldType ResultType;
 
          static ResultType apply(const Expr1& l, const Expr2& r)
          {
            typename Expr1::FieldType res(0);
 
            for (unsigned i = 0; i < l.size(); ++i) {
              res += l[i] * r[i];
            }
 
            return res;
          };
        };
 
        template<class dummy>
        struct ProductFactory<false, dummy>
        {
          typedef
          BinaryVectorExpression<MultiplyOperation, Expr1, Expr2>
          ResultType;
 
          static ResultType apply(const Expr1& l, const Expr2& r)
          {
            return ResultType(l, r);
          };
        };
 
        enum {
          sameElements = std::is_same<typename Expr1::ElementType, typename Expr2::ElementType>::value,
          fieldVectors = (sameElements &&
                          std::is_same<typename Expr1::FieldType, typename Expr1::ElementType>::value)
        };
 
        typedef typename ProductFactory<fieldVectors>::ResultType ResultType;
        static ResultType apply(const Expr1& l, const Expr2& r)
        {
          return ProductFactory<fieldVectors>::apply(l, r);
        }
      };
 
    }
 
    template<class Expr1, class Expr2>
    static inline
    typename InnerProductHelper<Expr1, Expr2>::ResultType
    operator*(const VectorExpression<Expr1>& e1, const VectorExpression<Expr2>& e2)
    {
      return InnerProductHelper<Expr1, Expr2>::apply(e1.expression(), e2.expression());
    }
 
    template<class Expr>
    static inline
    BinaryVectorExpression<SMultiplyOperation, typename Expr::FieldType, Expr>
    operator*(const typename Expr::FieldType& s, const VectorExpression<Expr>& e)
    {
      typedef typename Expr::FieldType FieldType;
      typedef BinaryVectorExpression<SMultiplyOperation, FieldType, Expr> ExpressionType;
 
      return ExpressionType(s, e);
    }
 
    template<class Expr>
    static inline
    auto
    operator*(const VectorExpression<Expr>& e, const typename Expr::FieldType& s)
      -> decltype(s * e)
    {
      return s * e;
    }
 
    template<class Expr>
    static inline
    BinaryVectorExpression<SMultiplyOperation, typename Expr::FieldType, Expr>
    operator*(const FieldVector<typename Expr::FieldType, 1>& s, const VectorExpression<Expr>& e)
    {
      typedef typename Expr::FieldType FieldType;
      typedef BinaryVectorExpression<SMultiplyOperation, FieldType, Expr> ExpressionType;
 
      return ExpressionType(s[0], e);
    }
 
    template<class Expr>
    static inline
    BinaryVectorExpression<SMultiplyOperation, typename Expr::FieldType, Expr>
    operator*(const VectorExpression<Expr>& e, const FieldVector<typename Expr::FieldType, 1>& s)
    {
      typedef typename Expr::FieldType FieldType;
      typedef BinaryVectorExpression<SMultiplyOperation, FieldType, Expr> ExpressionType;
 
      return ExpressionType(s[0], e);
    }
 
    template<class Expr>
    static inline
    UnaryVectorExpression<MinusOperation, Expr>
    operator-(const VectorExpression<Expr>& arg)
    {
      typedef UnaryVectorExpression<MinusOperation, Expr> ExpressionType;
 
      return ExpressionType(arg);
    }
 
    template<class Vector>
    static inline
    UnaryVectorExpression<IdentityOperation, Vector>
    asExpr(const Vector& arg)
    {
      return UnaryVectorExpression<IdentityOperation, Vector>(arg);
    }
 
    template<class Field, int rows, int cols>
    static inline
    UnaryVectorExpression<IdentityOperation, FieldMatrix<Field, rows, cols> >
    asExpr(const FieldMatrix<Field, rows, cols>& arg)
    {
      return UnaryVectorExpression<IdentityOperation, FieldMatrix<Field, rows, cols> >(arg);
    }
 
    template<class Op, class Expr>
    static inline
    UnaryVectorExpression<Op, Expr>
    asExpr(const UnaryVectorExpression<Op, Expr>& arg)
    {
      return arg.expression();
    }
 
    template<class Op, class Expr1, class Expr2>
    static inline
    BinaryVectorExpression<Op, Expr1, Expr2>
    asExpr(const BinaryVectorExpression<Op, Expr1, Expr2>& arg)
    {
      return arg.expression();
    }
 
    template<class Vector, class Op, class Expr1, class Expr2>
    static inline
    Vector& operator+=(Vector& v,
                       const BinaryVectorExpression<Op, Expr1, Expr2>& e)
    {
      for (unsigned i = 0; i < v.size(); ++i) {
        v[i] += e[i];
      }
 
      return v;
    }
 
    template<class Vector, class Op, class Expr>
    static inline
    Vector& operator+=(Vector& v,
                       const UnaryVectorExpression<Op, Expr>& e)
    {
      for (unsigned i = 0; i < v.size(); ++i) {
        v[i] += e[i];
      }
 
      return v;
    }
 
    template<class Vector, class Op, class Expr1, class Expr2>
    static inline
    Vector& operator-=(Vector& v,
                       const BinaryVectorExpression<Op, Expr1, Expr2>& e)
    {
      for (unsigned i = 0; i < v.size(); ++i) {
        v[i] -= e[i];
      }
 
      return v;
    }
 
    template<class Vector, class Op, class Expr>
    static inline
    Vector& operator-=(Vector& v,
                       const UnaryVectorExpression<Op, Expr>& e)
    {
      for (unsigned i = 0; i < v.size(); ++i) {
        v[i] -= e[i];
      }
 
      return v;
    }
 
    namespace {
      // assignment helper
 
      template<bool hasFieldElements, class Dummy = void>
      struct VectorAssignHelper;
 
      template<>
      struct VectorAssignHelper<true>
      {
        template<class To, class From>
        static void doAssign(To& to, const From& from)
        {
          to = from;
        }
      };
 
      template<>
      struct VectorAssignHelper<false>
      {
        template<class To, class From>
        static void doAssign(To& to, const From& from)
        {
          assign(to, from); // recurse
        }
      };
    }
 
    template<class Vector, class Expr>
    static inline
    Vector& assign(Vector& v, const VectorExpression<Expr>& e_)
    {
      const Expr& e(e_.expression());
      enum {
        hasFieldElements =
        std::is_same<typename Expr::ElementType, typename Expr::FieldType>::value
      };
 
      for (unsigned i = 0; i < v.size(); ++i) {
        VectorAssignHelper<hasFieldElements>::doAssign(v[i], e[i]);
      }
 
      return v;
    }
 
 
 
  } // namespace ACFem
 
} // namespace Dune
 
#endif // __DUNE_ACFEM_VECTOREXPRESSION_HH__