librheolef-dev 6.7-6 / usr / include / rheolef / field_expr_v2_linear.h

#ifndef _RHEOLEF_FIELD_EXPR_V2_LINEAR_H
#define _RHEOLEF_FIELD_EXPR_V2_LINEAR_H
///
/// This file is part of Rheolef.
///
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
///
/// Rheolef is free software; you can redistribute it and/or modify
/// it under the terms of the GNU General Public License as published by
/// the Free Software Foundation; either version 2 of the License, or
/// (at your option) any later version.
///
/// Rheolef is distributed in the hope that it will be useful,
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
/// GNU General Public License for more details.
///
/// You should have received a copy of the GNU General Public License
/// along with Rheolef; if not, write to the Free Software
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
/// 
/// =========================================================================
//
// field-valued affine expressions:
//
//	field wh = 2*uh - vh + 1:
//
// author: Pierre.Saramito@imag.fr
//
// date: 13 september 2015
//
// Notes; use template expressions and SFINAE techniques
//
// SUMMARY:
// 1. concept
// 2. unary operations
//    2.1. unary wrapper
//    2.2. unary calls
// 3. binary operations
//    3.1. binary wrapper
//    3.2. binary calls
// 4. field assignments
//    4.1. field assignment members
//    4.2. computed assignment
// 5. misc
//    5.1. duality product
//    5.2. form d = diag (expr)
//    5.3. output a linear expession

#include "rheolef/operators.h"
#include "rheolef/field.h"
#include "rheolef/field_component.h"
#include "rheolef/form.h"

namespace rheolef {

// -------------------------------------------------------------------
// 1. concept
// -------------------------------------------------------------------
// field_expr_v2_linear, a type concept for expression types

namespace details {

// Define a trait type for detecting field expression valid arguments
template<class T>          struct is_field_expr_v2_linear_arg                                : std::false_type {};
template<class T, class M> struct is_field_expr_v2_linear_arg <field_basic<T, M> >           : std::true_type {};
template<class T, class M> struct is_field_expr_v2_linear_arg <field_indirect       <T, M> > : std::true_type {}; 
template<class T, class M> struct is_field_expr_v2_linear_arg <field_indirect_const <T, M> > : std::true_type {};
template<class T, class M> struct is_field_expr_v2_linear_arg <field_component      <T, M> > : std::true_type {};
template<class T, class M> struct is_field_expr_v2_linear_arg <field_component_const<T, M> > : std::true_type {};

// predicate for filtering the field_basic class
template <class T>          struct is_field                     : std::false_type {};
template <class T, class M> struct is_field <field_basic<T,M> > : std::true_type {};

} // namespace details
// -------------------------------------------
// 2. unary operations
// -------------------------------------------
// 2.1. unary wrapper
// -------------------------------------------
namespace details {

template <class Op, class Expr>
struct field_expr_v2_unary {

// typedefs:

  typedef typename Expr::size_type      size_type;
  typedef typename Expr::value_type     value_type;
  typedef typename Expr::memory_type    memory_type;
  typedef typename Expr::const_iterator expr_const_iterator;
  typedef typename scalar_traits<value_type>::type    scalar_type;
  typedef typename float_traits <scalar_type>::type   float_type;

// allocatos:

  field_expr_v2_unary (const Op& op, const Expr& expr)
   : _op(op), _expr_iter(expr.begin_dof()), _space(expr.get_space()) {}

  template <class BinaryOp, class Constant>
  field_expr_v2_unary (const BinaryOp& binop, const Constant& c, const Expr& expr)
   : _op(binop,c), _expr_iter(expr.begin_dof()), _space(expr.get_space()) {}

  template <class BinaryOp, class Constant>
  field_expr_v2_unary (const BinaryOp& binop, const Expr& expr, const Constant& c)
   : _op(binop,c), _expr_iter(expr.begin_dof()), _space(expr.get_space()) {}

// accessors:

  const space_basic<scalar_type,memory_type>& get_space() const { return _space; }
  std::string stamp() const { return _space.stamp(); }

// minimal forward iterator interface:

  struct const_iterator {
    typedef std::forward_iterator_tag         iterator_category;
    typedef typename Expr::value_type         value_type;
    typedef value_type&                       reference;
    typedef value_type*                       pointer;
    typedef std::ptrdiff_t                    difference_type;
    const_iterator (Op op, expr_const_iterator expr_iter)
     : _op(op), _expr_iter (expr_iter) {}
    const_iterator& operator++ () { ++_expr_iter; return *this; }
    value_type operator* () const {  return _op (*_expr_iter); }
  protected:
    const Op                                   _op;
    expr_const_iterator                        _expr_iter;
  };
  const_iterator begin_dof() const { return const_iterator (_op, _expr_iter); }
protected:
  const Op                                   _op;
  const expr_const_iterator                  _expr_iter;
  const space_basic<scalar_type,memory_type> _space;
};
template<class Op, class Expr> struct is_field_expr_v2_linear_arg     <field_expr_v2_unary<Op,Expr> > : std::true_type {};

} // namespace details

// -------------------------------------------
// 2.2. unary calls
// -------------------------------------------
#define _RHEOLEF_field_expr_v2_unary_operator(OP, FUNCTOR) 			\
template <class Expr>								\
inline										\
typename 									\
std::enable_if<									\
  details::is_field_expr_v2_linear_arg<Expr>::value,				\
  details::field_expr_v2_unary<							\
    FUNCTOR,									\
    Expr									\
  >										\
>::type										\
operator OP (const Expr& expr)							\
{										\
  typedef details::field_expr_v2_unary <FUNCTOR, Expr>   expr_t;		\
  return expr_t (FUNCTOR(), expr); 						\
}

_RHEOLEF_field_expr_v2_unary_operator (+, details::unary_plus)
_RHEOLEF_field_expr_v2_unary_operator (-, details::negate)
#undef _RHEOLEF_field_expr_v2_unary_operator

// -------------------------------------------
// 3. binary operations
// -------------------------------------------
// 3.1. binary wrapper
// -------------------------------------------
namespace details {

template <class Op, class Expr1, class Expr2>
struct field_expr_v2_binary {

  typedef typename Expr1::size_type      size_type;
  typedef typename promote<
    typename Expr1::value_type,
    typename Expr2::value_type>::type   value_type;
  typedef typename Expr1::memory_type   memory_type; // TODO: check Expr2::memory_type
  typedef typename Expr1::const_iterator expr1_const_iterator;
  typedef typename Expr2::const_iterator expr2_const_iterator;
  typedef typename scalar_traits<value_type>::type    scalar_type;
  typedef typename float_traits <scalar_type>::type   float_type;

// allocators:

  field_expr_v2_binary (const Op& op, const Expr1& expr1, const Expr2& expr2)
   : _op (op),
     _iter1 (expr1.begin_dof()),
     _iter2 (expr2.begin_dof()),
     _space (expr1.get_space())
  {
    check_macro (expr1.stamp() == expr2.stamp(), "linear binary field expression: incompatible spaces "
                << expr1.stamp() << " and " << expr2.stamp());
  }

// accessors:

  const space_basic<scalar_type,memory_type>& get_space() const { return _space; }
  std::string stamp() const { return _space.stamp(); }

// minimal forward iterator interface:

  struct const_iterator {
    typedef std::forward_iterator_tag         iterator_category;
    typedef typename promote<
      typename Expr1::value_type,
      typename Expr2::value_type>::type       value_type;
    typedef value_type&                       reference;
    typedef value_type*                       pointer;
    typedef std::ptrdiff_t                    difference_type;
    const_iterator (Op op, expr1_const_iterator iter1, expr2_const_iterator iter2)
     : _op(op), _iter1 (iter1), _iter2 (iter2) {}
    const_iterator& operator++ () { ++_iter1; ++_iter2; return *this; }
    value_type operator* () const {  return _op (*_iter1, *_iter2); }
  protected:
    const Op                                  _op;
    expr1_const_iterator                      _iter1;
    expr2_const_iterator                      _iter2;
  };
  const_iterator begin_dof() const { return const_iterator (_op, _iter1, _iter2); }
protected:
  const Op                                    _op;
  expr1_const_iterator                        _iter1;
  expr2_const_iterator                        _iter2;
  const space_basic<scalar_type,memory_type>  _space;
};
template<class Op, class Expr1, class Expr2> struct is_field_expr_v2_linear_arg     <field_expr_v2_binary<Op,Expr1,Expr2> > : std::true_type {};

template <class Op, class Expr1, class Expr2, class Sfinae = void>
struct field_expr_v2_binary_traits { /* catch-all case */ };

// field_expr +- field_expr
template <class Op, class Expr1, class Expr2>
struct field_expr_v2_binary_traits <Op, Expr1, Expr2,
  typename std::enable_if<
    details::is_field_expr_v2_linear_arg<Expr1>::value &&
    details::is_field_expr_v2_linear_arg<Expr2>::value>::type>
{
  typedef field_expr_v2_binary <Op,Expr1,Expr2> type;
};
// constant +-* field_expr
template <class Op, class Expr1, class Expr2>
struct field_expr_v2_binary_traits <Op, Expr1, Expr2,
  typename std::enable_if<
       is_field_expr_v2_constant<Expr1>::value 
    && details::is_field_expr_v2_linear_arg<Expr2>::value
  >::type>
{
  typedef field_expr_v2_unary <details::binder_first<Op,Expr1>,Expr2> type;
};
// field_expr +-* constant
template <class Op, class Expr1, class Expr2>
struct field_expr_v2_binary_traits <Op, Expr1, Expr2,
  typename std::enable_if<
       details::is_field_expr_v2_linear_arg<Expr1>::value
    && is_field_expr_v2_constant<Expr2>::value>
  ::type>
{
  typedef field_expr_v2_unary <details::binder_second<Op,Expr2>,Expr1> type;
};

} // namespace details

// -------------------------------------------
// 3.2. binary calls
// -------------------------------------------
// uh+vh; uh+c ; c+uh
#define _RHEOLEF_field_expr_v2_binary_operator(OP, FUNCTOR) 			\
template <class Expr1, class Expr2>						\
inline										\
typename 									\
std::enable_if<									\
    (details::is_field_expr_v2_linear_arg<Expr1>::value &&			\
     details::is_field_expr_v2_linear_arg<Expr2>::value) ||			\
    (details::is_field_expr_v2_constant  <Expr1>::value &&			\
     details::is_field_expr_v2_linear_arg<Expr2>::value) ||			\
    (details::is_field_expr_v2_linear_arg<Expr1>::value &&			\
     details::is_field_expr_v2_constant  <Expr2>::value)			\
 ,typename 									\
  details::field_expr_v2_binary_traits<						\
    FUNCTOR,									\
    Expr1, Expr2								\
  >::type									\
>::type										\
operator OP (const Expr1& expr1, const Expr2& expr2)				\
{										\
  typedef typename details::field_expr_v2_binary_traits <FUNCTOR, Expr1, Expr2>::type  expr_t; \
  return expr_t (FUNCTOR(), expr1, expr2);					\
}
_RHEOLEF_field_expr_v2_binary_operator (+, details::plus)
_RHEOLEF_field_expr_v2_binary_operator (-, details::minus)
#undef _RHEOLEF_field_expr_v2_binary_operator

// c*uh ; uh*c
template <class Expr1, class Expr2>
inline
typename
std::enable_if<
    (details::is_field_expr_v2_constant  <Expr1>::value &&
     details::is_field_expr_v2_linear_arg<Expr2>::value) ||
    (details::is_field_expr_v2_linear_arg<Expr1>::value &&
     details::is_field_expr_v2_constant  <Expr2>::value)
 ,typename
  details::field_expr_v2_binary_traits<
    details::multiplies,
    Expr1, Expr2
  >::type
>::type
operator* (const Expr1& expr1, const Expr2& expr2)
{
  typedef typename details::field_expr_v2_binary_traits <details::multiplies, Expr1, Expr2>::type  expr_t;
  return expr_t (details::multiplies(), expr1, expr2);
}

// uh/c
template <class Expr1, class Expr2>
inline
typename
std::enable_if<
    (details::is_field_expr_v2_linear_arg<Expr1>::value &&
     details::is_field_expr_v2_constant  <Expr2>::value)
 ,typename
  details::field_expr_v2_binary_traits<
    details::divides,
    Expr1, Expr2
  >::type
>::type
operator/ (const Expr1& expr1, const Expr2& expr2)
{
  typedef typename details::field_expr_v2_binary_traits <details::divides, Expr1, Expr2>::type  expr_t;
  return expr_t (details::divides(), expr1, expr2);
}

// -------------------------------------------
// 4. field assignments
// -------------------------------------------
// 4.1. field assignment members
// -------------------------------------------
// 4.1.1 field
// -------------------------------------------
// uh = expr;
template<class T, class M>
template<class Expr, class Sfinae>
inline
field_basic<T,M>& 
field_basic<T,M>::operator= (const Expr& expr) {
  if (stamp() == "") {
    resize (expr.get_space());
  } else {
    check_macro (stamp() == expr.stamp(), "field = field_expression : incompatible spaces "
                << stamp() << " and " << expr.stamp());
  }
  details::assign_with_operator (begin_dof(), end_dof(), expr.begin_dof(), details::assign_op());
  return *this;
}
// field uh = expr;
template<class T, class M>
template<class Expr, class Sfinae>
inline
field_basic<T,M>::field_basic (const Expr& expr)
 : _V (),
   _u (),
   _b (),
   _dis_dof_indexes_requires_update(true),
   _dis_dof_assembly_requires_update(false)
{
    operator= (expr);
}

#ifdef TO_CLEAN
template<class T, class M>
inline
field_basic<T,M>&
field_basic<T,M>::operator=  (const field_indirect<T,M>& expr)
{
  if (stamp() == "") {
    resize (expr.get_space());
  } else {
    check_macro (stamp() == expr.stamp(), "incompatible spaces "
	<< stamp() << " and " << expr.stamp()
	<< " in field = field[domain]");
  }
  dis_dof_indexes_requires_update();
  algo::copy_n (expr.begin_dof(), expr.ndof(), begin_dof());
  return *this;
}
template<class T, class M>
inline
field_basic<T,M>&
field_basic<T,M>::operator=  (const field_indirect_const<T,M>& expr)
{
  if (stamp() == "") {
    resize (expr.get_space());
  } else {
    check_macro (stamp() == expr.stamp(), "incompatible spaces "
	<< stamp() << " and " << expr.stamp()
	<< " in field = field[domain]");
  }
  dis_dof_indexes_requires_update();
  algo::copy_n (expr.begin_dof(), expr.ndof(), begin_dof());
  return *this;
}
template<class T, class M>
inline
field_basic<T,M>&
field_basic<T,M>::operator=  (const field_component_const<T,M>& uh_comp)
{
  if (stamp() == "") {
    resize (uh_comp.get_space());
  } else {
    check_macro (stamp() == uh_comp.stamp(), "incompatible spaces "
        << stamp() << " and " << uh_comp.stamp()
        << " in field = field[i_comp]");
  }
  dis_dof_indexes_requires_update();
  std::copy (uh_comp.begin_dof(), uh_comp.end_dof(), begin_dof());
  return *this;
}
template<class T, class M>
inline
field_basic<T,M>::field_basic (const field_indirect<T,M>& expr)
 : _V(),
   _u(),
   _b(),
   _dis_dof_indexes_requires_update(true),
   _dis_dof_assembly_requires_update(false)
{
    operator= (expr);
}
template<class T, class M>
inline
field_basic<T,M>::field_basic (const field_indirect_const<T,M>& expr)
 : _V(),
   _u(),
   _b(),
   _dis_dof_indexes_requires_update(true),
   _dis_dof_assembly_requires_update(false)
{
    operator= (expr);
}
template<class T, class M>
inline
field_basic<T,M>&
field_basic<T,M>::operator=  (const field_component<T,M>& uh_comp)
{
  return operator= (field_component_const<T,M>(uh_comp));
}
template<class T, class M>
inline
field_basic<T,M>::field_basic (const field_component<T,M>& uh_comp)
 : _V(),
   _u(),
   _b(),
   _dis_dof_indexes_requires_update(true),
   _dis_dof_assembly_requires_update(false)
{
  operator= (uh_comp);
}
template<class T, class M>
inline
field_basic<T,M>::field_basic (const field_component_const<T,M>& uh_comp)
 : _V(),
   _u(),
   _b(),
   _dis_dof_indexes_requires_update(true),
   _dis_dof_assembly_requires_update(false)
{
  operator= (uh_comp);
}
#endif // TO_CLEAN
// -------------------------------------------
// 4.1.2. field_component
// -------------------------------------------
// uh [i_comp] = expr;
template<class T, class M>
template<class Expr, class Sfinae>
inline
field_component<T,M>& 
field_component<T,M>::operator= (const Expr& expr) {
  check_macro (stamp() == expr.stamp(), "field [i_comp] = field_expression : incompatible spaces "
                << stamp() << " and " << expr.stamp());
  details::assign_with_operator (begin_dof(), end_dof(), expr.begin_dof(), details::assign_op());
  return *this;
}
#ifdef TO_CLEAN
template<class T, class M>
inline
field_component<T,M>&
field_component<T,M>::operator= (const field_component<T,M>& uh_comp)
{
  check_macro (stamp() == uh_comp.stamp(), "incompatible spaces "
        << stamp() << " and " << uh_comp.stamp()
        << " in field[i_comp] = field[j_comp]");
  std::copy (uh_comp.begin_dof(), uh_comp.end_dof(), begin_dof());
  return *this;
}
template<class T, class M>
inline
field_component<T,M>&
field_component<T,M>::operator= (const field_basic<T,M>& uh)
{
  check_macro (stamp() == uh.stamp(), "incompatible spaces "
        << stamp() << " and " << uh.stamp()
        << " in field[i_comp] = field");
  std::copy (uh.begin_dof(), uh.end_dof(), begin_dof());
  return *this;
}
#endif // TO_CLEAN
// -------------------------------------------
// 4.1.2. field_indirect
// -------------------------------------------
// uh [domain] = expr;
template<class T, class M>
template<class Expr, class Sfinae>
inline
field_indirect<T,M>& 
field_indirect<T,M>::operator= (const Expr& expr) {
#ifdef TO_CLEAN
  // is it still a restriction ? check it...
  check_macro (_V.valued_tag() == space_constant::scalar, 
               "field[domain] = expression: non-scalar valued field not yet supported -- sorry");
#endif // TO_CLEAN
  check_macro (stamp() == expr.stamp(), "field [domain] = field_expression : incompatible spaces "
                << stamp() << " and " << expr.stamp());
#ifdef TO_CLEAN
  algo::copy_n (expr.begin_dof(), expr.ndof(), begin_dof());
#endif // TO_CLEAN
  details::assign_with_operator (begin_dof(), end_dof(), expr.begin_dof(), details::assign_op());
  return *this;
}
#ifdef TO_CLEAN
template<class T, class M>
inline
field_indirect<T,M>&
field_indirect<T,M>::operator=  (const field_basic<T,M>& expr)
{
  check_macro (_V.valued_tag() == space_constant::scalar, "field[domain]: unsupported non-scalar field");
  check_macro (stamp() == expr.stamp(), "field[domain] = field : incompatible spaces "
		<< stamp() << " and " << expr.stamp());
  algo::copy_n (expr.begin_dof(), expr.ndof(), begin_dof());
  return *this;
}
template<class T, class M>
inline
field_indirect<T,M>&
field_indirect<T,M>::operator=  (const field_indirect<T,M>& expr)
{
  check_macro (_V.valued_tag() == space_constant::scalar, "field[domain]: unsupported non-scalar field");
  check_macro (stamp() == expr.stamp(), "field[domain] = field[domain]: incompatible spaces "
		<< stamp() << " and " << expr.stamp());
  algo::copy_n (expr.begin_dof(), expr.ndof(), begin_dof());
  return *this;
}
template<class T, class M>
inline
field_indirect<T,M>&
field_indirect<T,M>::operator=  (const field_indirect_const<T,M>& expr)
{
  check_macro (_V.valued_tag() == space_constant::scalar, "field[domain]: unsupported non-scalar field");
  check_macro (stamp() == expr.stamp(), "field[domain] = field[domain]: incompatible spaces "
		<< stamp() << " and " << expr.stamp());
  algo::copy_n (expr.begin_dof(), expr.ndof(), begin_dof());
  return *this;
}
#endif // TO_CLEAN
// ---------------------------------------------------------------------------
// 4.2. computed assignment
// ---------------------------------------------------------------------------

// uh -+= expr
// uh [i_comp] -+= expr; // note: requires a move &&
// uh [domain] -+= expr;
#define _RHEOLEF_field_expr_v2_op_assign_field(OP, FUNCTOR)			\
template<class T, class M, class Expr>						\
inline										\
typename std::enable_if<							\
  details::is_field_expr_v2_linear_arg<Expr>::value,				\
  field_basic<T,M>&								\
>::type										\
operator OP (field_basic<T,M>& uh, const Expr& expr)				\
{										\
  check_macro (uh.stamp() == expr.stamp(), "field " << #OP << " field_expression : incompatible spaces " \
            << uh.stamp() << " and " << expr.stamp());				\
  details::assign_with_operator (uh.begin_dof(), uh.end_dof(), expr.begin_dof(), FUNCTOR()); \
  return uh;									\
}

#define _RHEOLEF_field_expr_v2_op_assign_auxil(OP, FUNCTOR, NAME, IDX)		\
template<class T, class M, class Expr>						\
inline										\
typename std::enable_if<							\
  details::is_field_expr_v2_linear_arg<Expr>::value,				\
  NAME<T,M>&									\
>::type										\
operator OP (NAME<T,M>&& uh, const Expr& expr)					\
{										\
  check_macro (uh.stamp() == expr.stamp(), "field [" << #IDX << "] " << #OP << " field_expression : incompatible spaces " \
            << uh.stamp() << " and " << expr.stamp());				\
  details::assign_with_operator (uh.begin_dof(), uh.end_dof(), expr.begin_dof(), FUNCTOR()); \
  return uh;									\
}

#define _RHEOLEF_field_expr_v2_op_assign(OP, FUNCTOR) 				\
        _RHEOLEF_field_expr_v2_op_assign_field(OP, FUNCTOR) 			\
        _RHEOLEF_field_expr_v2_op_assign_auxil(OP, FUNCTOR, field_component, "i_comp")	\
        _RHEOLEF_field_expr_v2_op_assign_auxil(OP, FUNCTOR, field_indirect,  "domain")

_RHEOLEF_field_expr_v2_op_assign (+=, details::plus_assign)
_RHEOLEF_field_expr_v2_op_assign (-=, details::minus_assign)
#undef _RHEOLEF_field_expr_v2_op_assign_field
#undef _RHEOLEF_field_expr_v2_op_assign_auxil
#undef _RHEOLEF_field_expr_v2_op_assign

// uh -+*/= c
// uh [i_comp] -+*/= c; // requires a move &&
// uh [domain] -+*/= c; // TODO
#define _RHEOLEF_field_expr_v2_op_assign_constant_field(OP, FUNCTOR)		\
template<class T, class M, class Expr>						\
inline										\
typename std::enable_if<							\
  details::is_field_expr_v2_constant<Expr>::value				\
 ,field_basic<T,M>&								\
>::type										\
operator OP (field_basic<T,M>& uh, const Expr& expr)				\
{										\
  details::assign_with_operator (uh.begin_dof(), uh.end_dof(), details::iterator_on_constant<Expr>(expr), FUNCTOR()); \
  return uh;									\
}

#define _RHEOLEF_field_expr_v2_op_assign_constant_auxil(OP, FUNCTOR, NAME, IDX)	\
template<class T, class M, class Expr>						\
inline										\
typename std::enable_if<							\
  details::is_field_expr_v2_constant<Expr>::value				\
 ,NAME<T,M>&									\
>::type										\
operator OP (NAME<T,M>&& uh, const Expr& expr)					\
{										\
  details::assign_with_operator (uh.begin_dof(), uh.end_dof(), details::iterator_on_constant<Expr>(expr), FUNCTOR()); \
  return uh;									\
}

#define _RHEOLEF_field_expr_v2_op_assign_constant(OP, FUNCTOR) 				\
        _RHEOLEF_field_expr_v2_op_assign_constant_field(OP, FUNCTOR) 			\
        _RHEOLEF_field_expr_v2_op_assign_constant_auxil(OP, FUNCTOR, field_component, "i_comp")	\
        _RHEOLEF_field_expr_v2_op_assign_constant_auxil(OP, FUNCTOR, field_indirect,  "domain")

_RHEOLEF_field_expr_v2_op_assign_constant (+=, details::plus_assign)
_RHEOLEF_field_expr_v2_op_assign_constant (-=, details::minus_assign)
_RHEOLEF_field_expr_v2_op_assign_constant (*=, details::multiplies_assign)
_RHEOLEF_field_expr_v2_op_assign_constant (/=, details::divides_assign)
#undef _RHEOLEF_field_expr_v2_op_assign_constant_field
#undef _RHEOLEF_field_expr_v2_op_assign_constant_auxil
#undef _RHEOLEF_field_expr_v2_op_assign_constant

// ---------------------------------------------------------------------------
// 5. misc
// ---------------------------------------------------------------------------
// 5.1. duality product
// ---------------------------------------------------------------------------
// dual (uh,vh)
template <class Expr1, class Expr2>
inline
typename
std::enable_if<
  details::is_field_expr_v2_linear_arg<Expr1>::value &&
  details::is_field_expr_v2_linear_arg<Expr2>::value,
  typename promote<
    typename Expr1::float_type,
    typename Expr2::float_type>::type
>::type
dual (const Expr1& expr1, const Expr2& expr2)
{
  typedef typename Expr1::memory_type M;
  return dis_inner_product (expr1.begin_dof(), expr2.begin_dof(), expr1.get_space().ndof(), expr1.get_space().ownership().comm(), M());
}
// dual (c,uh)
template <class Expr1, class Expr2>
inline
typename
std::enable_if<
  details::is_field_expr_v2_constant  <Expr1>::value &&
  details::is_field_expr_v2_linear_arg<Expr2>::value
 ,typename Expr2::float_type
>::type
dual (const Expr1& expr1, const Expr2& expr2)
{
  typedef typename Expr2::memory_type M;
  return expr1*dis_accumulate (expr2.begin_dof(), expr2.get_space().ndof(), expr2.get_space().ownership().comm(), M());
}
// dual (uh,c) 
template <class Expr1, class Expr2>
inline
typename
std::enable_if<
  details::is_field_expr_v2_linear_arg<Expr1>::value &&
  details::is_field_expr_v2_constant  <Expr2>::value
 ,typename Expr1::float_type
>::type
dual (const Expr1& expr1, const Expr2& expr2)
{
  typedef typename Expr1::memory_type M;
  return dis_accumulate (expr1.begin_dof(), expr1.get_space().ndof(), expr1.get_space().ownership().comm(), M())*expr2;
}
// ---------------------------------------------------------------------------
// 5.2. form d = diag (expr)
// ---------------------------------------------------------------------------
template<class Expr>
inline
typename
std::enable_if<
       details::is_field_expr_v2_linear_arg<Expr>::value
  && ! details::is_field<Expr>::value
 ,form_basic <typename Expr::value_type, typename Expr::memory_type>
>::type
diag (const Expr& expr)
{
  typedef typename Expr::value_type  T;
  typedef typename Expr::memory_type M;
  return diag (field_basic<T,M>(expr));
}
// -------------------------------------------
// 5.3. output a linear expession
// -------------------------------------------
template <class Expr>
inline
typename
std::enable_if<
       details::is_field_expr_v2_linear_arg<Expr>::value
  && ! details::is_field<Expr>::value
 ,odiststream&
>::type
operator<< (odiststream& ops, const Expr& expr) {
  // distributed case: communications requires to store in memory
  // and create a temporary
  typedef typename Expr::value_type  T;
  typedef typename Expr::memory_type M;
  field_basic<T,M> tmp = expr;
  return tmp.put (ops);
}

} // namespace rheolef
#endif // _RHEOLEF_FIELD_EXPR_V2_LINEAR_H