File: field_expr_recursive.h

package info (click to toggle)
rheolef 7.2-6
links: PTS, VCS
area: main
in suites: forky, sid
size: 88,076 kB
sloc: cpp: 110,259; sh: 16,733; makefile: 5,438; python: 1,391; yacc: 218; javascript: 203; xml: 191; awk: 61; sed: 5
file content (1507 lines) | stat: -rw-r--r-- 66,862 bytes
parent folder | download | duplicates (4)
#ifndef _RHEOLEF_FIELD_EXPR_RECURSIVE_H
#define _RHEOLEF_FIELD_EXPR_RECURSIVE_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 nonlinear expressions, as:
//
//	field wh = interpolate  (Xh, 1/uh - uh*(1-vh)):
//	Float Ih = integrate (omega, 1/uh - uh*(1-vh), qopt):
//
// author: Pierre.Saramito@imag.fr
//
// date: 15 september 2015
//
// Notes; use template expressions and SFINAE techiques
//   The interpolation operator is required, as in
//   1/uh and uh*vh that do not belong to Xh when uh, vh in Xh
//
// SUMMARY:
// 1. unary operations
//    1.1. unary node
//    1.2. unary calls
//    1.3. unary compose
// 2. binary operations
//    2.1. binary node
//    2.2. binary calls
//    2.3. binary compose
// 
#include "rheolef/field_expr_terminal.h"

namespace rheolef {

// -------------------------------------------
// 1. unary operations
// -------------------------------------------
// 1.1. unary node
// -------------------------------------------
namespace details {

template<class UnaryFunction, class Expr>
class field_expr_v2_nonlinear_node_unary {
public:
// typedefs:

  using size_type   = geo_element::size_type;
  using memory_type = typename Expr::memory_type;
  using result_type = typename details::generic_unary_traits<UnaryFunction>::template result_hint<typename Expr::result_type>::type;
  using value_type  = result_type;
  using scalar_type = typename scalar_traits<value_type>::type;
  using float_type  = typename  float_traits<scalar_type>::type;
  using geo_type    = geo_basic<float_type,memory_type>;
  using band_type   = band_basic<float_type,memory_type>;
  using space_type  = space_basic<float_type,memory_type>;
  using self_type   = field_expr_v2_nonlinear_node_unary<UnaryFunction,Expr>;

// alocators:

  field_expr_v2_nonlinear_node_unary (const UnaryFunction& f, const Expr& expr);

// --------------------------------------------
// accessors for the affine & homogeneous case:
// --------------------------------------------

  using is_affine_homogeneous 
                    = and_type<
                        or_type<
                          std::is_same<UnaryFunction,details::unary_plus>
                         ,std::is_same<UnaryFunction,details::negate>
                         ,std::is_same<UnaryFunction,details::binder_first <details::plus,      scalar_type>>
                         ,std::is_same<UnaryFunction,details::binder_second<details::plus,      scalar_type>>
                         ,std::is_same<UnaryFunction,details::binder_first <details::minus,     scalar_type>>
                         ,std::is_same<UnaryFunction,details::binder_second<details::minus,     scalar_type>>
                         ,std::is_same<UnaryFunction,details::binder_first <details::multiplies,scalar_type>>
                         ,std::is_same<UnaryFunction,details::binder_second<details::multiplies,scalar_type>>
                         ,std::is_same<UnaryFunction,details::binder_second<details::divides,   scalar_type>>
                        >
                       ,is_field_expr_affine_homogeneous<Expr>
                      >;

  bool have_homogeneous_space (space_basic<scalar_type,memory_type>& Vh) const {
    return is_affine_homogeneous::value
        && _expr.have_homogeneous_space (Vh);
  }
  // minimal forward iterator interface: 
  struct const_iterator {
    using iterator_category = std::forward_iterator_tag;
    using value_type        = typename Expr::scalar_type;
    using reference         = value_type&;
    using pointer           = value_type*;
    using difference_type   = std::ptrdiff_t;
    const_iterator (UnaryFunction f, typename Expr::const_iterator expr_iter)
     : _f(f), _expr_iter (expr_iter) {}
    const_iterator& operator++ () { ++_expr_iter; return *this; }
    value_type operator* () const {  return _f (*_expr_iter); }
  protected:
    const UnaryFunction            _f;
    typename Expr::const_iterator _expr_iter;
  };
  const_iterator begin_dof() const { return const_iterator (_f, _expr.begin_dof()); }

// --------------------------------------------
// interface for the general nonlinear case:
// --------------------------------------------
// accessors:

  static const space_constant::valued_type valued_hint = space_constant::valued_tag_traits<result_type>::value;
  const Expr&       expr()      const { return _expr; }

  space_constant::valued_type valued_tag() const {
    return details::generic_unary_traits<UnaryFunction>::valued_tag (_expr.valued_tag());
  }

#ifdef TO_CLEAN
// field_lazy interface:

  const geo_type&   get_geo()   const { return _expr.get_geo(); }
  const band_type&  get_band()  const { return _expr.get_band(); }
  const space_type& get_space() const { return _expr.get_space(); }
  bool is_on_band () const { return _expr.is_on_band(); }
  void initialize (const geo_type& omega_K) { _expr.initialize (omega_K); }
#endif // TO_CLEAN

// initializators:

  // field_lazy interface:
	
  void initialize (
    const piola_on_pointset<float_type>&             pops,
    const integrate_option&                          iopt)
  		{ _expr.initialize (pops, iopt); }

  void initialize (
    const space_basic<float_type,memory_type>&       Xh,
    const piola_on_pointset<float_type>&             pops,
    const integrate_option&                          iopt)
  		{ _expr.initialize (Xh, pops, iopt); }

// -------------------------
// run time check args types
// -------------------------
  template<class Result, class Arg, class Status>
  struct evaluate_call_check {
    template<class M>
    void operator() (
      const self_type&                         obj,
      const geo_basic<float_type,M>&           omega_K,
      const geo_element&                       K,
      Eigen::Matrix<Result,Eigen::Dynamic,1>&  value)
    {
      fatal_macro ("invalid type resolution: Result="<<typename_macro(Result)
          << ", Arg="<<typename_macro(Arg)
          << ", UnaryFunction="<<typename_macro(UnaryFunction));
    }
    template<class M>
    void operator() (
      const self_type&                         obj,
      const geo_basic<float_type,M>&           omega_K,
      const geo_element&                       K,
      const side_information_type&             sid,
      Eigen::Matrix<Result,Eigen::Dynamic,1>&  value) const
    {
      fatal_macro ("invalid type resolution: Result="<<typename_macro(Result)
          << ", Arg="<<typename_macro(Arg)
          << ", UnaryFunction="<<typename_macro(UnaryFunction));
    }
  };
  template<class Result, class Arg>
  struct evaluate_call_check<Result,Arg,std::true_type> {
    // in an element:
    template<class M>
    void operator() (
      const self_type&                         obj,
      const geo_basic<float_type,M>&           omega_K,
      const geo_element&                       K,
      Eigen::Matrix<Result,Eigen::Dynamic,1>&  value) const
    {
      Eigen::Matrix<Arg,Eigen::Dynamic,1> value1; 
      obj._expr.evaluate (omega_K, K, value1);
      value.resize(value1.size());
      for (size_type i = 0, ni = value.rows(); i < ni; ++i) {
        value[i] = obj._f (value1[i]);
      }
    }
    // on a side:
    template<class M>
    void operator() (
      const self_type&                         obj,
      const geo_basic<float_type,M>&           omega_K,
      const geo_element&                       K,
      const side_information_type&             sid,
      Eigen::Matrix<Result,Eigen::Dynamic,1>&  value) const
    { 
      Eigen::Matrix<Arg,Eigen::Dynamic,1> value1; 
      obj._expr.evaluate_on_side (omega_K, K, sid, value1);
      value.resize (value1.size());
      for (size_type i = 0, ni = value.rows(); i < ni; ++i) {
        value[i] = obj._f (value1[i]);
      }
    }
  };
// -------------------------------------------
// evaluate in an element, with known arg type
// -------------------------------------------
  template<class Result, class Arg, class M>
  void evaluate_call (
      const geo_basic<float_type,M>&           omega_K,
      const geo_element&                       K,
      Eigen::Matrix<Result,Eigen::Dynamic,1>&  value) const
  {
    typedef typename details::generic_unary_traits<UnaryFunction>::template hint<Arg,Result>::result_type result_type;
    typedef typename details::and_type<
		typename details::is_equal<Result,result_type>::type
	       ,typename details::not_type<typename details::is_error<Arg>::type>::type
	      >::type
      status_t;
    evaluate_call_check<Result,Arg,status_t> eval;
    eval (*this, omega_K, K, value);
  }
// -------------------------------------------
// evaluate on a side, with known arg type
// -------------------------------------------
  template<class Result, class Arg, class M>
  void evaluate_call (
      const geo_basic<float_type,M>&           omega_K,
      const geo_element&                       K,
      const side_information_type&             sid,
      Eigen::Matrix<Result,Eigen::Dynamic,1>&  value) const
  { 
    typedef typename details::generic_unary_traits<UnaryFunction>::template hint<Arg,Result>::result_type result_type;
    typedef typename details::and_type<
		typename details::is_equal<Result,result_type>::type
	       ,typename details::not_type<typename details::is_error<Arg>::type>::type
	      >::type
      status_t;
    evaluate_call_check<Result,Arg,status_t> eval;
    eval (*this, omega_K, K, sid, value);
  }
// -------------------------------------------
// determine args at compile-time or run-time:
// -------------------------------------------
  template<class This, class Result, class Arg, space_constant::valued_type ArgTag = space_constant::valued_tag_traits<Arg>::value>
  struct evaluate_switch {};

  // when arg is unknown at run-time:
  template<class This, class Result, class Arg>
  struct evaluate_switch<This, Result, Arg, space_constant::last_valued> {
    template<class M>
    void evaluate (
      const This&                              obj,
      const geo_basic<float_type,M>&           omega_K,
      const geo_element&                       K,
      Eigen::Matrix<Result,Eigen::Dynamic,1>&  value) const
    {
      typedef typename scalar_traits<Arg>::type T;
      space_constant::valued_type arg_valued_tag = obj._expr.valued_tag();
      switch (arg_valued_tag) {
        case space_constant::scalar:
	  obj.template evaluate_call<Result,T,M> (omega_K, K, value); break;
        case space_constant::vector:
	  obj.template evaluate_call<Result, point_basic<T> > (omega_K, K, value); break;
        case space_constant::tensor:
        case space_constant::unsymmetric_tensor:
	  obj.template evaluate_call<Result, tensor_basic<T> > (omega_K, K, value); break;
        default: { error_macro ("unexpected valued tag="<<arg_valued_tag); }
      }
    }
    template<class M>
    void evaluate_on_side (
      const This&                              obj,
      const geo_basic<float_type,M>&           omega_K,
      const geo_element&                       K,
      const side_information_type&             sid,
      Eigen::Matrix<Result,Eigen::Dynamic,1>&  value) const
    {
      typedef typename scalar_traits<Arg>::type T;
      space_constant::valued_type arg_valued_tag = obj._expr.valued_tag();
      switch (arg_valued_tag) {
        case space_constant::scalar:
	  obj.template evaluate_call<Result,T> (omega_K, K, sid, value); break;
        case space_constant::vector:
	  obj.template evaluate_call<Result, point_basic<T> > (omega_K, K, value); break;
        case space_constant::tensor:
        case space_constant::unsymmetric_tensor:
	  obj.template evaluate_call<Result, tensor_basic<T> > (omega_K, K, value); break;
        default: { error_macro ("unexpected valued tag="<<arg_valued_tag); }
      }
    }
  };
  // specializations when arg is known at compile-time:
#define _RHEOLEF_evaluate_switch_specialization(VALUED,VALUE)		\
  template<class This, class Result, class Arg>				\
  struct evaluate_switch <This, Result, Arg, VALUED> {			\
    typedef typename scalar_traits<Arg>::type T;			\
    typedef typename float_traits<Arg>::type  float_type; 		\
									\
    void evaluate (							\
      const This&                              obj,			\
      const geo_basic<float_type,memory_type>& omega_K,			\
      const geo_element&                       K,			\
      Eigen::Matrix<Result,Eigen::Dynamic,1>&  value) const		\
      	{ obj.template evaluate_call<Result, VALUE> (omega_K, K, value); } \
									\
    template<class M>							\
    void evaluate_on_side (						\
      const This&                              obj,			\
      const geo_basic<float_type,M>&           omega_K,			\
      const geo_element&                       K,			\
      const side_information_type&             sid,			\
      Eigen::Matrix<Result,Eigen::Dynamic,1>&  value) const		\
      	{ obj.template evaluate_call<Result, VALUE> (omega_K, K, sid, value); }	\
  };									\

_RHEOLEF_evaluate_switch_specialization(space_constant::scalar,T)
_RHEOLEF_evaluate_switch_specialization(space_constant::vector,point_basic<T>)
_RHEOLEF_evaluate_switch_specialization(space_constant::tensor,tensor_basic<T>)
_RHEOLEF_evaluate_switch_specialization(space_constant::tensor3,tensor3_basic<T>)
_RHEOLEF_evaluate_switch_specialization(space_constant::tensor4,tensor4_basic<T>)
#undef _RHEOLEF_evaluate_switch_specialization

// ----------------------
// evaluate in an element
// ----------------------
  template<class Result>
  void
  evaluate (
    const geo_basic<float_type,memory_type>& omega_K,
    const geo_element&                       K,
    Eigen::Matrix<Result,Eigen::Dynamic,1>&  value) const
  {
    typedef field_expr_v2_nonlinear_node_unary<UnaryFunction, Expr> This;
    typedef typename details::generic_unary_traits<UnaryFunction>::template hint<typename Expr::value_type,Result>::argument_type
            A1;
    evaluate_switch <This, Result, A1> helper;
    helper.evaluate (*this, omega_K, K, value);
  }
// -------------------
// evaluate on a side:
// -------------------
  template<class Result>
  void
  evaluate_on_side (
    const geo_basic<float_type,memory_type>& omega_K,
    const geo_element&                       K,
    const side_information_type&             sid,
    Eigen::Matrix<Result,Eigen::Dynamic,1>&  value) const
 { 
    typedef field_expr_v2_nonlinear_node_unary<UnaryFunction, Expr> This;
    typedef typename details::generic_unary_traits<UnaryFunction>::template hint<typename Expr::value_type,Result>::argument_type
            A1;
    evaluate_switch <This, Result, A1> helper;
    helper.evaluate_on_side (*this, omega_K, K, sid, value);
  }

  template<class Result>
  bool valued_check() const {
    typedef typename details::generic_unary_traits<UnaryFunction>::template hint<typename Expr::value_type,Result>::argument_type
                     A1;
    if (! is_undeterminated<A1>::value) return _expr.template valued_check<A1>();
    return true;
  }
protected:
// data:
  UnaryFunction   _f;
  Expr            _expr;
// working area:
public:
  mutable std::array<
    Eigen::Matrix<scalar_type,Eigen::Dynamic,1>
   ,reference_element::max_variant>                           	_scalar_val;
  mutable std::array<
    Eigen::Matrix<point_basic<scalar_type>,Eigen::Dynamic,1>
   ,reference_element::max_variant>          			_vector_val;
  mutable std::array<
    Eigen::Matrix<tensor_basic<scalar_type>,Eigen::Dynamic,1>
   ,reference_element::max_variant>          			_tensor_val;
  mutable std::array<
    Eigen::Matrix<tensor3_basic<scalar_type>,Eigen::Dynamic,1>
   ,reference_element::max_variant>          			_tensor3_val;
  mutable std::array<
    Eigen::Matrix<tensor4_basic<scalar_type>,Eigen::Dynamic,1>
   ,reference_element::max_variant>          			_tensor4_val;
};

template<class UnaryFunction, class Expr>
field_expr_v2_nonlinear_node_unary<UnaryFunction,Expr>::field_expr_v2_nonlinear_node_unary (
  const UnaryFunction& f,
  const Expr& expr) 
  : _f(f),
    _expr(expr),
    _scalar_val(),
    _vector_val(),
    _tensor_val(),
    _tensor3_val(),
    _tensor4_val()
{
}

template<class F, class Expr> struct is_field_expr_v2_nonlinear_arg  <field_expr_v2_nonlinear_node_unary<F,Expr>> : std::true_type {};
//template<class F, class Expr> struct has_field_lazy_interface        <field_expr_v2_nonlinear_node_unary<F,Expr>> : std::true_type {};
template<class F, class Expr> struct is_field_expr_affine_homogeneous<field_expr_v2_nonlinear_node_unary<F,Expr>, typename std::enable_if<
                                                                      field_expr_v2_nonlinear_node_unary<F,Expr>::is_affine_homogeneous::value>::type>: std::true_type {};


} // namespace details
// -------------------------------------------
// 1.2. unary calls
// -------------------------------------------
// unary operators +- and std::math

// ------------------------
// standard unary operators
// ------------------------
#define _RHEOLEF_make_field_expr_v2_nonlinear_unary_operator(FUNCTION,FUNCTOR)	\
template<class Expr>								\
inline										\
typename									\
std::enable_if<									\
       details::is_field_expr_v2_nonlinear_arg<Expr>::value			\
  && ! details::is_field_expr_v2_constant     <Expr>::value			\
  && ! details::has_field_rdof_interface      <Expr>::value			\
 ,details::field_expr_v2_nonlinear_node_unary<					\
    FUNCTOR									\
   ,typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr>::type \
  >										\
>::type										\
FUNCTION (const Expr& expr)							\
{										\
  typedef typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr>::type wrap_t;  \
  return details::field_expr_v2_nonlinear_node_unary <FUNCTOR,wrap_t> (FUNCTOR(), wrap_t(expr)); \
}

_RHEOLEF_make_field_expr_v2_nonlinear_unary_operator (operator+, details::unary_plus)
_RHEOLEF_make_field_expr_v2_nonlinear_unary_operator (operator-, details::negate)
#undef _RHEOLEF_make_field_expr_v2_nonlinear_unary_operator

// ------------------------
// std::cmath
// ------------------------
#define _RHEOLEF_make_field_expr_v2_nonlinear_unary_operator(FUNCTION,FUNCTOR)	\
template<class Expr>								\
inline										\
typename									\
std::enable_if<									\
  (    details::is_field_expr_v2_nonlinear_arg<Expr>::value			\
  ||   details::is_field<Expr>::value)						\
  && ! details::is_field_expr_v2_constant     <Expr>::value			\
 ,details::field_expr_v2_nonlinear_node_unary<					\
    FUNCTOR									\
   ,typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr>::type \
  >										\
>::type										\
FUNCTION (const Expr& expr)							\
{										\
  typedef typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr>::type wrap_t;  \
  return details::field_expr_v2_nonlinear_node_unary <FUNCTOR,wrap_t> (FUNCTOR(), wrap_t(expr)); \
}
#define _RHEOLEF_make_field_expr_v2_nonlinear_unary_function(FUNCTION)		\
        _RHEOLEF_make_field_expr_v2_nonlinear_unary_operator(FUNCTION, details::FUNCTION##_)

_RHEOLEF_make_field_expr_v2_nonlinear_unary_function (cos)
_RHEOLEF_make_field_expr_v2_nonlinear_unary_function (sin)
_RHEOLEF_make_field_expr_v2_nonlinear_unary_function (tan)
_RHEOLEF_make_field_expr_v2_nonlinear_unary_function (acos)
_RHEOLEF_make_field_expr_v2_nonlinear_unary_function (asin)
_RHEOLEF_make_field_expr_v2_nonlinear_unary_function (atan)
_RHEOLEF_make_field_expr_v2_nonlinear_unary_function (cosh)
_RHEOLEF_make_field_expr_v2_nonlinear_unary_function (sinh)
_RHEOLEF_make_field_expr_v2_nonlinear_unary_function (tanh)
_RHEOLEF_make_field_expr_v2_nonlinear_unary_function (exp)
_RHEOLEF_make_field_expr_v2_nonlinear_unary_function (log)
_RHEOLEF_make_field_expr_v2_nonlinear_unary_function (log10)
_RHEOLEF_make_field_expr_v2_nonlinear_unary_function (sqrt)
_RHEOLEF_make_field_expr_v2_nonlinear_unary_function (abs)
_RHEOLEF_make_field_expr_v2_nonlinear_unary_function (fabs)
_RHEOLEF_make_field_expr_v2_nonlinear_unary_function (floor)
_RHEOLEF_make_field_expr_v2_nonlinear_unary_function (ceil)
// rheolef extensions
_RHEOLEF_make_field_expr_v2_nonlinear_unary_function (sqr)
_RHEOLEF_make_field_expr_v2_nonlinear_unary_function (norm)
_RHEOLEF_make_field_expr_v2_nonlinear_unary_function (norm2)

// tr(sigma_h) & trans(sigma_h) : trace & transpose of a tensor-valued field
_RHEOLEF_make_field_expr_v2_nonlinear_unary_function (tr)
_RHEOLEF_make_field_expr_v2_nonlinear_unary_function (trans)

#undef _RHEOLEF_make_field_expr_v2_nonlinear_unary_function
#undef _RHEOLEF_make_field_expr_v2_nonlinear_unary_operator

// -------------------------------------------
// 1.3. unary compose
// -------------------------------------------

template<class Function, class Expr>
inline
typename
std::enable_if<
  details::is_field_expr_v2_nonlinear_arg<Expr>::value
 ,details::field_expr_v2_nonlinear_node_unary<
    typename details::function_traits<Function>::functor_type
   ,typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr>::type
  >
>::type
compose (const Function& f, const Expr& expr)
{
  typedef typename details::function_traits<Function>::functor_type                     fun_wrap_t;
  typedef typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr>::type expr_wrap_t;
  return details::field_expr_v2_nonlinear_node_unary <fun_wrap_t, expr_wrap_t> (fun_wrap_t(f), expr_wrap_t(expr));
}
// ---------------------------------------------------------------------------
// 2. binary operations
// ---------------------------------------------------------------------------
// 2.1. binary node
// -------------------------------------------
namespace details {

template<class BinaryFunction, class Expr1, class Expr2>
class field_expr_v2_nonlinear_node_binary {
public:
// typedefs:

  using size_type   = geo_element::size_type;
  using result_type = typename details::generic_binary_traits<BinaryFunction>::template result_hint<typename Expr1::result_type,typename Expr2::result_type>::type;
  using value_type  = result_type;
  using scalar_type = typename scalar_traits<value_type>::type;
  using float_type  = typename  float_traits<value_type>::type;
  using memory_type = typename Expr1::memory_type;

// alocators:

  field_expr_v2_nonlinear_node_binary (
	const BinaryFunction&  f, 
	const Expr1&           expr1,
        const Expr2&           expr2);

// --------------------------------------------
// accessors for the affine & homogeneous case:
// --------------------------------------------

  // the result expr is affine-homogeneous if and only if:
  //    binop expr1   expr2   
  //    +-    A|C     A|C
  //    */    A       C 
  //    *     C       A
  using is_affine_homogeneous 
                  = or_type<
                      and_type<
                        or_type<
                          std::is_same<BinaryFunction,details::plus>
                         ,std::is_same<BinaryFunction,details::minus>
                        >
                       ,is_field_expr_affine_homogeneous<Expr1>
                       ,is_field_expr_affine_homogeneous<Expr2>
                      >
                     ,and_type<
                        or_type<
                          std::is_same<BinaryFunction,details::multiplies>
                         ,std::is_same<BinaryFunction,details::divides>
                        >
                       ,is_field_expr_affine_homogeneous<Expr1>
                       ,is_field_expr_v2_constant       <Expr2>
                      >
                     ,and_type<
                        std::is_same<BinaryFunction,details::multiplies>
                       ,is_field_expr_v2_constant       <Expr1>
                       ,is_field_expr_affine_homogeneous<Expr2>
                      >
		    >;
  bool have_homogeneous_space (space_basic<scalar_type,memory_type>& Vh) const {
    space_basic<scalar_type,memory_type> Vh2;
    return is_affine_homogeneous::value
        && _expr1.have_homogeneous_space (Vh)
        && _expr2.have_homogeneous_space (Vh2)
        && Vh.name() == Vh2.name();
  }
  // minimal forward iterator interface:
  struct const_iterator {
    using iterator_category = std::forward_iterator_tag;
    using value_type        = typename promote<
                                typename Expr1::scalar_type,
                                typename Expr2::scalar_type>::type;
    using reference         = value_type&;
    using pointer           = value_type*;
    using difference_type   = std::ptrdiff_t;
    const_iterator (const BinaryFunction& f, typename Expr1::const_iterator iter1, typename Expr2::const_iterator iter2)
     : _f(f), _iter1 (iter1), _iter2 (iter2) {}
    const_iterator& operator++ () { ++_iter1; ++_iter2; return *this; }
    value_type operator* () const {  return _f (*_iter1, *_iter2); }
  protected:
    const BinaryFunction           _f;
    typename Expr1::const_iterator _iter1;
    typename Expr2::const_iterator _iter2;
  };
  const_iterator begin_dof() const { return const_iterator (_f, _expr1.begin_dof(), _expr2.begin_dof()); }

// --------------------------------------------
// interface for the general nonlinear case:
// --------------------------------------------

// accessors:
 
  static const space_constant::valued_type valued_hint = space_constant::valued_tag_traits<result_type>::value;

  space_constant::valued_type valued_tag() const {
    return details::generic_binary_traits<BinaryFunction>::valued_tag(_expr1.valued_tag(), _expr2.valued_tag());
  }

// initializers:

  void initialize (
    const piola_on_pointset<float_type>&             pops,
    const integrate_option&                          iopt)
  {
    _expr1.initialize (pops, iopt);
    _expr2.initialize (pops, iopt);
  }
  void initialize (
    const space_basic<float_type,memory_type>&       Xh,
    const piola_on_pointset<float_type>&             pops,
    const integrate_option&                          iopt)
  {
    _expr1.initialize (Xh, pops, iopt);
    _expr2.initialize (Xh, pops, iopt);
  }

// evaluators:

  template<class Result, class Arg1, class Arg2, class M>
  void evaluate_internal2 (
      const geo_basic<float_type,M>&          omega_K,
      const geo_element&                      K,
      Eigen::Matrix<Result,Eigen::Dynamic,1>& value) const
  {
    Eigen::Matrix<Arg1,Eigen::Dynamic,1> value1; _expr1.evaluate (omega_K, K, value1);
    Eigen::Matrix<Arg2,Eigen::Dynamic,1> value2; _expr2.evaluate (omega_K, K, value2);
    value.resize (value1.size());
    // TODO: DVT_EIGEN_BLAS1
    for (size_t i = 0, ni = value.rows(); i < ni; ++i) {
      value[i] = _f (value1[i], value2[i]);
    }
  }
  template<class Result, class Arg1, class Arg2, class M>
  void evaluate_internal2 (
      const geo_basic<float_type,M>&          omega_K,
      const geo_element&                      K,
      const side_information_type&            sid,
      Eigen::Matrix<Result,Eigen::Dynamic,1>& value) const
  {
    Eigen::Matrix<Arg1,Eigen::Dynamic,1> value1; _expr1.evaluate_on_side (omega_K, K, sid, value1);
    Eigen::Matrix<Arg2,Eigen::Dynamic,1> value2; _expr2.evaluate_on_side (omega_K, K, sid, value2);
    value.resize (value1.size());
    for (size_t i = 0, ni = value.rows(); i < ni; ++i) {
      value[i] = _f (value1[i], value2[i]);
    }
  }
  template<class This, class Result, class ReturnType, class Arg1, class Arg2>
  struct evaluate_internal {
    template<class M>
    void operator() (
      const This&                             obj,
      const geo_basic<float_type,M>&          omega_K,
      const geo_element&                      K,
      Eigen::Matrix<Result,Eigen::Dynamic,1>& value) const
    {
      fatal_macro ("unexpected return type "	
	<< pretty_typename_macro(ReturnType) << ": "
	<< pretty_typename_macro(Result) << " was expected for function "
	<< pretty_typename_macro(BinaryFunction) << "("
	<< pretty_typename_macro(Arg1) << ","
	<< pretty_typename_macro(Arg2) << ")");
    }
    template<class M>
    void operator() (
      const This&                             obj,
      const geo_basic<float_type,M>&          omega_K,
      const geo_element&                      K,
      const side_information_type&            sid,
      Eigen::Matrix<Result,Eigen::Dynamic,1>& value) const
    {
      fatal_macro ("unexpected return type "	
	<< pretty_typename_macro(ReturnType) << ": "
	<< pretty_typename_macro(Result) << " was expected for function "
	<< pretty_typename_macro(BinaryFunction) << "("
	<< pretty_typename_macro(Arg1) << ","
	<< pretty_typename_macro(Arg2) << ")");
    }
  };
  template<class This, class Result, class Arg1, class Arg2>
  struct evaluate_internal<This,Result,Result,Arg1,Arg2> {
    template<class M>
    void operator() (
      const This&                             obj,
      const geo_basic<float_type,M>&          omega_K,
      const geo_element&                      K,
      Eigen::Matrix<Result,Eigen::Dynamic,1>& value) const
    	{ obj.template evaluate_internal2<Result,Arg1,Arg2,M> (omega_K, K, value); }

    template<class M>
    void operator() (
      const This&                             obj,
      const geo_basic<float_type,M>&          omega_K,
      const geo_element&                      K,
      const side_information_type&            sid,
      Eigen::Matrix<Result,Eigen::Dynamic,1>& value) const
	{ obj.template evaluate_internal2 <Result,Arg1,Arg2,M> (omega_K, K, sid, value);
    }
  };
  template<class Result, class Arg1, class Arg2, class M>
  void evaluate_call (
      const geo_basic<float_type,M>&          omega_K,
      const geo_element&                      K,
      Eigen::Matrix<Result,Eigen::Dynamic,1>& value) const
  {
    typedef typename details::generic_binary_traits<BinaryFunction>::template result_hint<Arg1,Arg2>::type ReturnType;
    typedef field_expr_v2_nonlinear_node_binary<BinaryFunction, Expr1, Expr2> This;
    evaluate_internal<This,Result,ReturnType,Arg1,Arg2> eval_int;
    eval_int (*this, omega_K, K, value);
  }
  template<class Result, class Arg1, class Arg2, class M>
  void evaluate_call (
      const geo_basic<float_type,M>&          omega_K,
      const geo_element&                      K,
      const side_information_type&            sid,
      Eigen::Matrix<Result,Eigen::Dynamic,1>& value) const
  { 
    typedef typename details::generic_binary_traits<BinaryFunction>::template result_hint<Arg1,Arg2>::type ReturnType;
    typedef field_expr_v2_nonlinear_node_binary<BinaryFunction, Expr1, Expr2> This;
    evaluate_internal<This,Result,ReturnType,Arg1,Arg2> eval_int;
    eval_int (*this, omega_K, K, sid, value);
  }
  // when both args are defined at compile time:
  template<class This, class Result,
	class Arg1,        space_constant::valued_type Arg1Tag,
  	class Arg2,        space_constant::valued_type Arg2Tag>
  struct evaluate_switch {
    template<class M>
    void operator() (
      const This&                             obj,
      const geo_basic<float_type,M>&          omega_K,
      const geo_element&                      K,
      Eigen::Matrix<Result,Eigen::Dynamic,1>& value) const
    	{ obj.template evaluate_call<Result, Arg1, Arg2> (omega_K, K, value); }

    template<class M>
    void operator() (
      const This&                             obj,
      const geo_basic<float_type,M>&          omega_K,
      const geo_element&                      K,
      const side_information_type&            sid,
      Eigen::Matrix<Result,Eigen::Dynamic,1>& value) const
      	{ obj.template evaluate_call<Result, Arg1, Arg2> (omega_K, K, sid, value); }
  };
  // specialization when both args are undefined at compile time:
  template<class This, class Result,
	class Arg1,
  	class Arg2>
  struct evaluate_switch<This, Result,
	Arg1,   space_constant::last_valued,
        Arg2,   space_constant::last_valued> {
    template<class M>
    void operator() (
      const This&                             obj,
      const geo_basic<float_type,M>&          omega_K,
      const geo_element&                      K,
      Eigen::Matrix<Result,Eigen::Dynamic,1>& value) const
    {
      typedef typename scalar_traits<Arg1>::type T1;
      typedef typename scalar_traits<Arg2>::type T2;
      space_constant::valued_type arg1_valued_tag = obj._expr1.valued_tag();
      space_constant::valued_type arg2_valued_tag = obj._expr2.valued_tag();
      switch (arg1_valued_tag) {
        case space_constant::scalar: {
          switch (arg2_valued_tag) {
            case space_constant::scalar:
	      return obj.template evaluate_call<Result, T1, T2>                (omega_K, K, value); break;
            case space_constant::vector:
	      return obj.template evaluate_call<Result, T1, point_basic<T2> >  (omega_K, K, value); break;
            case space_constant::tensor:
            case space_constant::unsymmetric_tensor:
	      return obj.template evaluate_call<Result, T1, tensor_basic<T2> > (omega_K, K, value); break;
            case space_constant::tensor3:
	      return obj.template evaluate_call<Result, T1, tensor3_basic<T2> >(omega_K, K, value); break;
            default: error_macro ("unexpected second argument valued tag="<<arg2_valued_tag);
          }
          break;
        }
        case space_constant::vector: {
          switch (arg2_valued_tag) {
            case space_constant::scalar:
	      return obj.template evaluate_call<Result, point_basic<T1>, T2>                (omega_K, K, value); break;
            case space_constant::vector:
	      return obj.template evaluate_call<Result, point_basic<T1>, point_basic<T2> >  (omega_K, K, value); break;
            case space_constant::tensor:
            case space_constant::unsymmetric_tensor:
	      return obj.template evaluate_call<Result, point_basic<T1>, tensor_basic<T2> > (omega_K, K, value); break;
            case space_constant::tensor3:
	      return obj.template evaluate_call<Result, point_basic<T1>, tensor3_basic<T2> >(omega_K, K, value); break;
            default: error_macro ("unexpected second argument valued tag="<<arg2_valued_tag);
          }
          break;
        }
        case space_constant::tensor:
        case space_constant::unsymmetric_tensor: {
          switch (arg2_valued_tag) {
            case space_constant::scalar:
	      return obj.template evaluate_call<Result, tensor_basic<T1>, T2>                (omega_K, K, value); break;
            case space_constant::vector:
	      return obj.template evaluate_call<Result, tensor_basic<T1>, point_basic<T2> >  (omega_K, K, value); break;
            case space_constant::tensor:
            case space_constant::unsymmetric_tensor:
	      return obj.template evaluate_call<Result, tensor_basic<T1>, tensor_basic<T2> > (omega_K, K, value); break;
            case space_constant::tensor3:
	      return obj.template evaluate_call<Result, tensor_basic<T1>, tensor3_basic<T2> >(omega_K, K, value); break;
            default: error_macro ("unexpected second argument valued tag="<<arg2_valued_tag);
          }
          break;
        }
        case space_constant::tensor3: {
          switch (arg2_valued_tag) {
            case space_constant::scalar:
	      return obj.template evaluate_call<Result, tensor3_basic<T1>, T2>                (omega_K, K, value); break;
            case space_constant::vector:
	      return obj.template evaluate_call<Result, tensor3_basic<T1>, point_basic<T2> >  (omega_K, K, value); break;
            case space_constant::tensor:
            case space_constant::unsymmetric_tensor:
	      return obj.template evaluate_call<Result, tensor3_basic<T1>, tensor_basic<T2> > (omega_K, K, value); break;
            case space_constant::tensor3:
	      return obj.template evaluate_call<Result, tensor3_basic<T1>, tensor3_basic<T2> >(omega_K, K, value); break;
            default: error_macro ("unexpected second argument valued tag="<<arg2_valued_tag);
          }
          break;
        }
        default: error_macro ("unexpected first argument valued tag="<<arg1_valued_tag);
      }
    }
  };
  // specialization when only first arg is defined at compile time:
  template<class This, class Result,
	class Arg1,        space_constant::valued_type Arg1Tag,
  	class Arg2>
  struct evaluate_switch<This, Result,
	Arg1,   Arg1Tag,
        Arg2,   space_constant::last_valued> {
    template<class M>
    void operator() (
      const This&                             obj,
      const geo_basic<float_type,M>&          omega_K,
      const geo_element&                      K,
      Eigen::Matrix<Result,Eigen::Dynamic,1>& value) const
    {
      typedef typename scalar_traits<Arg2>::type T2;
      space_constant::valued_type arg2_valued_tag = obj._expr2.valued_tag();
      switch (arg2_valued_tag) {
        case space_constant::scalar:
	  return obj.template evaluate_call<Result, Arg1, T2>                (omega_K, K, value); break;
        case space_constant::vector:
	  return obj.template evaluate_call<Result, Arg1, point_basic<T2> >  (omega_K, K, value); break;
        case space_constant::tensor:
        case space_constant::unsymmetric_tensor:
	  return obj.template evaluate_call<Result, Arg1, tensor_basic<T2> > (omega_K, K, value); break;
        case space_constant::tensor3:
	  return obj.template evaluate_call<Result, Arg1, tensor3_basic<T2> > (omega_K, K, value); break;
        default: error_macro ("unexpected second argument valued tag="<<arg2_valued_tag);
      }
    }
  };
  // specialization when only second arg is defined at compile time:
  template<class This, class Result,
	class Arg1,     
  	class Arg2,        space_constant::valued_type Arg2Tag>
  struct evaluate_switch<This, Result,
	Arg1,              space_constant::last_valued,
        Arg2,              Arg2Tag> {
    template<class M>
    void operator() (
      const This&                             obj,
      const geo_basic<float_type,M>&          omega_K,
      const geo_element&                      K,
      Eigen::Matrix<Result,Eigen::Dynamic,1>& value) const
    {
      typedef typename scalar_traits<Arg1>::type T1;
      space_constant::valued_type arg1_valued_tag = obj._expr1.valued_tag();
      switch (arg1_valued_tag) {
        case space_constant::scalar:
	  return obj.template evaluate_call<Result, T1, Arg2>               (omega_K, K, value); break;
        case space_constant::vector:
	  return obj.template evaluate_call<Result, point_basic<T1>, Arg2>  (omega_K, K, value); break;
        case space_constant::tensor:
        case space_constant::unsymmetric_tensor:
	  return obj.template evaluate_call<Result, tensor_basic<T1>, Arg2> (omega_K, K, value); break;
        case space_constant::tensor3:
	  return obj.template evaluate_call<Result, tensor3_basic<T1>, Arg2>(omega_K, K, value); break;
        default: error_macro ("unexpected first argument valued tag="<<arg1_valued_tag);
      }
    }
  };
  template<class Result, class M>
  void evaluate (
      const geo_basic<float_type,M>&          omega_K,
      const geo_element&                      K,
      Eigen::Matrix<Result,Eigen::Dynamic,1>& value) const
  {
    typedef typename details::generic_binary_traits<BinaryFunction>::template hint<
          typename Expr1::value_type
         ,typename Expr2::value_type
         ,Result>::first_argument_type   A1;
    typedef typename details::generic_binary_traits<BinaryFunction>::template hint<
          typename Expr1::value_type
         ,typename Expr2::value_type
         ,Result>::second_argument_type A2;
    static const space_constant::valued_type  first_argument_tag = space_constant::valued_tag_traits<A1>::value;
    static const space_constant::valued_type second_argument_tag = space_constant::valued_tag_traits<A2>::value;
    typedef field_expr_v2_nonlinear_node_binary<BinaryFunction, Expr1, Expr2> This;
    evaluate_switch <This, Result, A1, first_argument_tag, A2, second_argument_tag>   eval;
    eval (*this, omega_K, K, value);
  }
  template<class Result, class M>
  void
  evaluate_on_side (
      const geo_basic<float_type,M>&          omega_K,
      const geo_element&                      K,
      const side_information_type&            sid,
      Eigen::Matrix<Result,Eigen::Dynamic,1>& value) const
  {
    typedef typename details::generic_binary_traits<BinaryFunction>::template hint<
          typename Expr1::value_type
         ,typename Expr2::value_type
         ,Result>::first_argument_type   A1;
    typedef typename details::generic_binary_traits<BinaryFunction>::template hint<
          typename Expr1::value_type
         ,typename Expr2::value_type
         ,Result>::second_argument_type A2;
    static const space_constant::valued_type  first_argument_tag = space_constant::valued_tag_traits<A1>::value;
    static const space_constant::valued_type second_argument_tag = space_constant::valued_tag_traits<A2>::value;
    typedef field_expr_v2_nonlinear_node_binary<BinaryFunction, Expr1, Expr2> This;
    evaluate_switch <This, Result, A1, first_argument_tag, A2, second_argument_tag>   eval;
    eval (*this, omega_K, K, sid, value);
  }
  template<class Result>
  bool valued_check() const {
    typedef typename details::generic_binary_traits<BinaryFunction>::template hint<
          typename Expr1::value_type
         ,typename Expr2::value_type
         ,Result>::first_argument_type   A1;
    typedef typename details::generic_binary_traits<BinaryFunction>::template hint<
          typename Expr1::value_type
         ,typename Expr2::value_type
         ,Result>::second_argument_type A2;
    bool status = true;
    if (! is_undeterminated<A1>::value)  status &= _expr1.template valued_check<A1>();
    if (! is_undeterminated<A2>::value)  status &= _expr2.template valued_check<A2>();
    return status;
  }

protected:
// data:
  BinaryFunction  _f;
  Expr1           _expr1;
  Expr2           _expr2;
};
template<class BinaryFunction, class Expr1, class Expr2>
field_expr_v2_nonlinear_node_binary<BinaryFunction,Expr1,Expr2>::field_expr_v2_nonlinear_node_binary (
  const BinaryFunction&  f, 
  const Expr1&           expr1,
  const Expr2&           expr2)
  : _f(f),
    _expr1(expr1),
    _expr2(expr2)
{
}

template<class F, class Expr1, class Expr2> struct is_field_expr_v2_nonlinear_arg  <field_expr_v2_nonlinear_node_binary<F,Expr1,Expr2>> : std::true_type {};
//template<class F, class Expr1, class Expr2> struct has_field_lazy_interface        <field_expr_v2_nonlinear_node_binary<F,Expr1,Expr2>> : std::true_type {};
template<class F, class Expr1, class Expr2> struct is_field_expr_affine_homogeneous<field_expr_v2_nonlinear_node_binary<F,Expr1,Expr2>, typename std::enable_if<
                                                                                    field_expr_v2_nonlinear_node_binary<F,Expr1,Expr2>::is_affine_homogeneous::value>::type>: std::true_type {};

} // namespace details
// -------------------------------------------
// 2.2. binary calls
// -------------------------------------------
/*
    combination table:

      +-   | c a n
        ---|-------
         c | C A N
         a | A A N
         n | N N N

      *    | c a n
        ---|-------
         c | C A N
         a | A N N
         n | N N N

      /    | c a n
        ---|-------
         c | C N N
         a | A N N
         n | N N N

    argument:
      c : constant, as scalar, point, tensor, ect
      l : affine homogeneous expr argument: as field, field_rdof_indirect or field_expr_node::is_affine_homogeneous
      n : function, functor or ! field_expr_node::is_affine_homogeneous
    result:
      C : constant : this combination is not implemented here
      A,N : are implemented here
    rules:
      at least one of the two args is not of type "c"
      when c: c value is embeded in bind_first or bind_second
              and the operation reduces to an unary one
      when a: if it is a field_convertible, it should be wrapped
              in field_expr_v2_nonlinear_terminal_field
      when c: no wrapper is need
    implementation:
      The a and n cases are grouped, thanks to the wrapper_traits
      and it remains to cases :
        1) both args are  field_expr_v2_nonlinear or a function 
        2) one  arg  is a field_expr_v2_nonlinear or a function and the second argument is a constant
*/

// -------------------------------------------
// binary+-
// -------------------------------------------
// TODO: remove rdof+-rdof when rdof_binary is available 
#define _RHEOLEF_make_field_expr_v2_nonlinear_binary(FUNCTION,FUNCTOR)	\
template<class Expr1, class Expr2>						\
inline										\
typename									\
std::enable_if<									\
     details::is_field_expr_v2_nonlinear_arg<Expr1>::value			\
  && details::is_field_expr_v2_nonlinear_arg<Expr2>::value			\
  && ! details::is_field_expr_v2_constant   <Expr1>::value			\
  && ! details::is_field_expr_v2_constant   <Expr2>::value			\
 ,details::field_expr_v2_nonlinear_node_binary<					\
    FUNCTOR									\
   ,typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr1>::type	\
   ,typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr2>::type 	\
  >										\
>::type										\
FUNCTION (const Expr1& expr1, const Expr2& expr2)				\
{										\
  typedef typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr1>::type wrap1_t; \
  typedef typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr2>::type wrap2_t; \
  return details::field_expr_v2_nonlinear_node_binary <FUNCTOR,wrap1_t,wrap2_t> \
	(FUNCTOR(), wrap1_t(expr1), wrap2_t(expr2)); 				\
}										\
template<class Expr1, class Expr2>						\
inline										\
typename 									\
std::enable_if<									\
       details::is_field_expr_v2_constant     <Expr1>::value			\
  &&   details::is_field_expr_v2_nonlinear_arg<Expr2>::value			\
  && ! details::is_field_expr_v2_constant     <Expr2>::value			\
  && !(details::has_field_rdof_interface      <Expr2>::value			\
       && details::is_rheolef_arithmetic      <Expr1>::value)			\
 ,details::field_expr_v2_nonlinear_node_unary<					\
    details::binder_first<							\
      FUNCTOR									\
     ,typename details::field_promote_first_argument<				\
        Expr1									\
       ,typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr2>::type::value_type \
      >::type									\
    >										\
   ,typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr2>::type \
  >										\
>::type										\
FUNCTION (const Expr1& expr1, const Expr2& expr2)				\
{										\
  typedef typename details::field_promote_first_argument<			\
        Expr1									\
       ,typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr2>::type::value_type \
      >::type									\
      value_type;								\
  typedef details::binder_first<FUNCTOR,value_type> fun_t;			\
  typedef typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr2>::type wrap2_t; \
  return details::field_expr_v2_nonlinear_node_unary<fun_t,wrap2_t>(fun_t(FUNCTOR(), expr1), wrap2_t(expr2)); \
}										\
template<class Expr1, class Expr2>						\
inline										\
typename 									\
std::enable_if<									\
     details::is_field_expr_v2_constant     <Expr2>::value			\
  && details::is_field_expr_v2_nonlinear_arg<Expr1>::value			\
  && ! details::is_field_expr_v2_constant   <Expr1>::value			\
  && !(details::has_field_rdof_interface    <Expr1>::value			\
       && details::is_rheolef_arithmetic    <Expr2>::value)			\
 ,details::field_expr_v2_nonlinear_node_unary<					\
    details::binder_second<							\
      FUNCTOR									\
     ,typename details::field_promote_second_argument<				\
        typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr1>::type::value_type \
       ,Expr2									\
      >::type									\
    >										\
   ,typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr1>::type \
  >										\
>::type										\
FUNCTION (const Expr1& expr1, const Expr2& expr2)				\
{										\
  typedef typename details::field_promote_second_argument<			\
        typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr1>::type::value_type \
       ,Expr2									\
      >::type									\
      value_type;								\
  typedef details::binder_second<FUNCTOR,value_type> fun_t;			\
  typedef typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr1>::type wrap1_t; \
  return details::field_expr_v2_nonlinear_node_unary<fun_t,wrap1_t>(fun_t(FUNCTOR(), expr2), wrap1_t(expr1)); \
}

_RHEOLEF_make_field_expr_v2_nonlinear_binary (operator+, details::plus)
_RHEOLEF_make_field_expr_v2_nonlinear_binary (operator-, details::minus)
#undef _RHEOLEF_make_field_expr_v2_nonlinear_binary

// -------------------------------------------
// binary*
// -------------------------------------------
#define _RHEOLEF_make_field_expr_v2_nonlinear_binary(FUNCTION,FUNCTOR)	\
template<class Expr1, class Expr2>						\
inline										\
typename									\
std::enable_if<									\
     details::is_field_expr_v2_nonlinear_arg<Expr1>::value			\
  && details::is_field_expr_v2_nonlinear_arg<Expr2>::value			\
  && ! details::is_field_expr_v2_constant   <Expr1>::value			\
  && ! details::is_field_expr_v2_constant   <Expr2>::value			\
 ,details::field_expr_v2_nonlinear_node_binary<					\
    FUNCTOR									\
   ,typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr1>::type	\
   ,typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr2>::type 	\
  >										\
>::type										\
FUNCTION (const Expr1& expr1, const Expr2& expr2)				\
{										\
  typedef typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr1>::type wrap1_t; \
  typedef typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr2>::type wrap2_t; \
  return details::field_expr_v2_nonlinear_node_binary <FUNCTOR,wrap1_t,wrap2_t> \
	(FUNCTOR(), wrap1_t(expr1), wrap2_t(expr2)); 				\
}										\
template<class Expr1, class Expr2>						\
inline										\
typename 									\
std::enable_if<									\
     details::is_field_expr_v2_constant     <Expr1>::value			\
  && details::is_field_expr_v2_nonlinear_arg<Expr2>::value			\
  && ! details::is_field_expr_v2_constant   <Expr2>::value			\
  && !(details::has_field_rdof_interface    <Expr2>::value			\
       && details::is_rheolef_arithmetic    <Expr1>::value)			\
 ,details::field_expr_v2_nonlinear_node_unary<					\
    details::binder_first<							\
      FUNCTOR									\
     ,typename details::field_promote_first_argument<				\
        Expr1									\
       ,typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr2>::type::value_type \
      >::type									\
    >										\
   ,typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr2>::type \
  >										\
>::type										\
FUNCTION (const Expr1& expr1, const Expr2& expr2)				\
{										\
  typedef typename details::field_promote_first_argument<			\
        Expr1									\
       ,typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr2>::type::value_type \
      >::type									\
      value_type;								\
  typedef details::binder_first<FUNCTOR,value_type> fun_t;			\
  typedef typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr2>::type wrap2_t; \
  return details::field_expr_v2_nonlinear_node_unary<fun_t,wrap2_t>(fun_t(FUNCTOR(), expr1), wrap2_t(expr2)); \
}										\
template<class Expr1, class Expr2>						\
inline										\
typename 									\
std::enable_if<									\
     details::is_field_expr_v2_constant     <Expr2>::value			\
  && details::is_field_expr_v2_nonlinear_arg<Expr1>::value			\
  && ! details::is_field_expr_v2_constant   <Expr1>::value			\
  && !(details::has_field_rdof_interface    <Expr1>::value			\
       && details::is_rheolef_arithmetic    <Expr2>::value)			\
 ,details::field_expr_v2_nonlinear_node_unary<					\
    details::binder_second<							\
      FUNCTOR									\
     ,typename details::field_promote_second_argument<				\
        typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr1>::type::value_type \
       ,Expr2									\
      >::type									\
    >										\
   ,typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr1>::type \
  >										\
>::type										\
FUNCTION (const Expr1& expr1, const Expr2& expr2)				\
{										\
  typedef typename details::field_promote_second_argument<			\
        typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr1>::type::value_type \
       ,Expr2									\
      >::type									\
      value_type;								\
  typedef details::binder_second<FUNCTOR,value_type> fun_t;			\
  typedef typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr1>::type wrap1_t; \
  return details::field_expr_v2_nonlinear_node_unary<fun_t,wrap1_t>(fun_t(FUNCTOR(), expr2), wrap1_t(expr1)); \
}
_RHEOLEF_make_field_expr_v2_nonlinear_binary (operator*, details::multiplies)
#undef _RHEOLEF_make_field_expr_v2_nonlinear_binary

// -------------------------------------------
// binary/
// -------------------------------------------
#define _RHEOLEF_make_field_expr_v2_nonlinear_binary(FUNCTION,FUNCTOR)	\
template<class Expr1, class Expr2>						\
inline										\
typename									\
std::enable_if<									\
     details::is_field_expr_v2_nonlinear_arg<Expr1>::value			\
  && details::is_field_expr_v2_nonlinear_arg<Expr2>::value			\
  && ! details::is_field_expr_v2_constant   <Expr1>::value			\
  && ! details::is_field_expr_v2_constant   <Expr2>::value			\
 ,details::field_expr_v2_nonlinear_node_binary<					\
    FUNCTOR									\
   ,typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr1>::type	\
   ,typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr2>::type 	\
  >										\
>::type										\
FUNCTION (const Expr1& expr1, const Expr2& expr2)				\
{										\
  typedef typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr1>::type wrap1_t; \
  typedef typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr2>::type wrap2_t; \
  return details::field_expr_v2_nonlinear_node_binary <FUNCTOR,wrap1_t,wrap2_t> \
	(FUNCTOR(), wrap1_t(expr1), wrap2_t(expr2)); 				\
}										\
template<class Expr1, class Expr2>						\
inline										\
typename 									\
std::enable_if<									\
     details::is_field_expr_v2_constant     <Expr1>::value			\
  && details::is_field_expr_v2_nonlinear_arg<Expr2>::value			\
  && ! details::is_field_expr_v2_constant   <Expr2>::value			\
 ,details::field_expr_v2_nonlinear_node_unary<					\
    details::binder_first<							\
      FUNCTOR									\
     ,typename details::field_promote_first_argument<				\
        Expr1									\
       ,typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr2>::type::value_type \
      >::type									\
    >										\
   ,typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr2>::type \
  >										\
>::type										\
FUNCTION (const Expr1& expr1, const Expr2& expr2)				\
{										\
  typedef typename details::field_promote_first_argument<			\
        Expr1									\
       ,typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr2>::type::value_type \
      >::type									\
      value_type;								\
  typedef details::binder_first<FUNCTOR,value_type> fun_t;			\
  typedef typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr2>::type wrap2_t; \
  return details::field_expr_v2_nonlinear_node_unary<fun_t,wrap2_t>(fun_t(FUNCTOR(), expr1), wrap2_t(expr2)); \
}										\
template<class Expr1, class Expr2>						\
inline										\
typename 									\
std::enable_if<									\
     details::is_field_expr_v2_constant     <Expr2>::value			\
  && details::is_field_expr_v2_nonlinear_arg<Expr1>::value			\
  && ! details::is_field_expr_v2_constant   <Expr1>::value			\
  && !(details::has_field_rdof_interface    <Expr1>::value			\
       && details::is_rheolef_arithmetic    <Expr2>::value)			\
 ,details::field_expr_v2_nonlinear_node_unary<					\
    details::binder_second<							\
      FUNCTOR									\
     ,typename details::field_promote_second_argument<				\
        typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr1>::type::value_type \
       ,Expr2									\
      >::type									\
    >										\
   ,typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr1>::type \
  >										\
>::type										\
FUNCTION (const Expr1& expr1, const Expr2& expr2)				\
{										\
  typedef typename details::field_promote_second_argument<			\
        typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr1>::type::value_type \
       ,Expr2									\
      >::type									\
      value_type;								\
  typedef details::binder_second<FUNCTOR,value_type> fun_t;			\
  typedef typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr1>::type wrap1_t; \
  return details::field_expr_v2_nonlinear_node_unary<fun_t,wrap1_t>(fun_t(FUNCTOR(), expr2), wrap1_t(expr1)); \
}
_RHEOLEF_make_field_expr_v2_nonlinear_binary (operator/, details::divides)
#undef _RHEOLEF_make_field_expr_v2_nonlinear_binary

// -------------------------------------------
// std::maths
// -------------------------------------------
#define _RHEOLEF_make_field_expr_v2_nonlinear_binary(FUNCTION,FUNCTOR)	\
template<class Expr1, class Expr2>						\
inline										\
typename									\
std::enable_if<									\
     details::is_field_expr_v2_nonlinear_arg<Expr1>::value			\
  && details::is_field_expr_v2_nonlinear_arg<Expr2>::value			\
  && ! details::is_field_expr_v2_constant   <Expr1>::value			\
  && ! details::is_field_expr_v2_constant   <Expr2>::value			\
 ,details::field_expr_v2_nonlinear_node_binary<					\
    FUNCTOR									\
   ,typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr1>::type	\
   ,typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr2>::type 	\
  >										\
>::type										\
FUNCTION (const Expr1& expr1, const Expr2& expr2)				\
{										\
  typedef typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr1>::type wrap1_t; \
  typedef typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr2>::type wrap2_t; \
  return details::field_expr_v2_nonlinear_node_binary <FUNCTOR,wrap1_t,wrap2_t> \
	(FUNCTOR(), wrap1_t(expr1), wrap2_t(expr2)); 				\
}										\
template<class Expr1, class Expr2>						\
inline										\
typename 									\
std::enable_if<									\
     details::is_field_expr_v2_constant     <Expr1>::value			\
  && details::is_field_expr_v2_nonlinear_arg<Expr2>::value			\
  && ! details::is_field_expr_v2_constant   <Expr2>::value			\
 ,details::field_expr_v2_nonlinear_node_unary<					\
    details::binder_first<							\
      FUNCTOR									\
     ,typename details::field_promote_first_argument<				\
        Expr1									\
       ,typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr2>::type::value_type \
      >::type									\
    >										\
   ,typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr2>::type \
  >										\
>::type										\
FUNCTION (const Expr1& expr1, const Expr2& expr2)				\
{										\
  typedef typename details::field_promote_first_argument<			\
        Expr1									\
       ,typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr2>::type::value_type \
      >::type									\
      value_type;								\
  typedef details::binder_first<FUNCTOR,value_type> fun_t;			\
  typedef typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr2>::type wrap2_t; \
  return details::field_expr_v2_nonlinear_node_unary<fun_t,wrap2_t>(fun_t(FUNCTOR(), expr1), wrap2_t(expr2)); \
}										\
template<class Expr1, class Expr2>						\
inline										\
typename 									\
std::enable_if<									\
     details::is_field_expr_v2_constant     <Expr2>::value			\
  && details::is_field_expr_v2_nonlinear_arg<Expr1>::value			\
  && ! details::is_field_expr_v2_constant   <Expr1>::value			\
 ,details::field_expr_v2_nonlinear_node_unary<					\
    details::binder_second<							\
      FUNCTOR									\
     ,typename details::field_promote_second_argument<				\
        typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr1>::type::value_type \
       ,Expr2									\
      >::type									\
    >										\
   ,typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr1>::type \
  >										\
>::type										\
FUNCTION (const Expr1& expr1, const Expr2& expr2)				\
{										\
  typedef typename details::field_promote_second_argument<			\
        typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr1>::type::value_type \
       ,Expr2									\
      >::type									\
      value_type;								\
  typedef details::binder_second<FUNCTOR,value_type> fun_t;			\
  typedef typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr1>::type wrap1_t; \
  return details::field_expr_v2_nonlinear_node_unary<fun_t,wrap1_t>(fun_t(FUNCTOR(), expr2), wrap1_t(expr1)); \
}

#define _RHEOLEF_make_field_expr_v2_nonlinear_binary_function(FUNCTION)				\
        _RHEOLEF_make_field_expr_v2_nonlinear_binary (FUNCTION, details::FUNCTION##_) 								\

_RHEOLEF_make_field_expr_v2_nonlinear_binary_function (atan2)
_RHEOLEF_make_field_expr_v2_nonlinear_binary_function (pow)
_RHEOLEF_make_field_expr_v2_nonlinear_binary_function (fmod)
_RHEOLEF_make_field_expr_v2_nonlinear_binary_function (min)
_RHEOLEF_make_field_expr_v2_nonlinear_binary_function (max)
_RHEOLEF_make_field_expr_v2_nonlinear_binary_function (dot)
_RHEOLEF_make_field_expr_v2_nonlinear_binary_function (ddot)
_RHEOLEF_make_field_expr_v2_nonlinear_binary_function (dddot)

#undef _RHEOLEF_make_field_expr_v2_nonlinear_binary_function
#undef _RHEOLEF_make_field_expr_v2_nonlinear_binary

// -------------------------------------------
// 2.3. binary compose
// -------------------------------------------
// note: compose 1 & 2 are not reductible to n-ary
//       as it uses deductible return types
// TODO: do not use deductible types => reduces to n-ary !!

// two args are field-expressions
template<class Function, class Expr1, class Expr2>
inline
typename
std::enable_if<
       details::is_field_expr_v2_nonlinear_arg<Expr1>::value
  && ! details::is_field_expr_v2_constant     <Expr1>::value
  &&   details::is_field_expr_v2_nonlinear_arg<Expr2>::value
  && ! details::is_field_expr_v2_constant     <Expr2>::value
 ,details::field_expr_v2_nonlinear_node_binary<
    typename details::function_traits<Function>::functor_type
   ,typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr1>::type
   ,typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr2>::type
  >
>::type
compose (const Function& f, const Expr1& expr1, const Expr2& expr2)
{
  typedef typename details::function_traits<Function>::functor_type                      fun_wrap_t;
  typedef typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr1>::type expr1_wrap_t;
  typedef typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr2>::type expr2_wrap_t;
  return details::field_expr_v2_nonlinear_node_binary 
	<fun_wrap_t,    expr1_wrap_t,        expr2_wrap_t> 
	(fun_wrap_t(f), expr1_wrap_t(expr1), expr2_wrap_t(expr2));
}
// left arg is a constant
template <class Function, class Expr1, class Expr2>
inline
typename
std::enable_if<
       details::is_field_expr_v2_constant     <Expr1>::value
  &&   details::is_field_expr_v2_nonlinear_arg<Expr2>::value
  && ! details::is_field_expr_v2_constant     <Expr2>::value
 ,details::field_expr_v2_nonlinear_node_unary<
    details::binder_first<
      typename details::function_traits<Function>::functor_type
     ,typename promote<
        Expr1
       ,typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr2>::type::value_type
      >::type
    >
   ,typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr2>::type
  >
>::type
compose (const Function& f, const Expr1& expr1, const Expr2& expr2)
{
  typedef typename promote<
    Expr1
   ,typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr2>::type::value_type
  >::type value_type;
  typedef typename details::function_traits<Function>::functor_type wrap_fun_t;
  typedef details::binder_first<wrap_fun_t,value_type> binded_fun_t;
  typedef typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr2>::type wrap2_t;
  return details::field_expr_v2_nonlinear_node_unary
	   <binded_fun_t,                       wrap2_t>
           (binded_fun_t(wrap_fun_t(f), expr1), wrap2_t(expr2));
}
// right arg is a constant
template <class Function, class Expr1, class Expr2>
inline
typename
std::enable_if<
       details::is_field_expr_v2_nonlinear_arg<Expr1>::value
  && ! details::is_field_expr_v2_constant     <Expr1>::value
  &&   details::is_field_expr_v2_constant     <Expr2>::value
 ,details::field_expr_v2_nonlinear_node_unary<
    details::binder_second<
      typename details::function_traits<Function>::functor_type
     ,typename promote<
        typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr1>::type::value_type
       ,Expr2
      >::type
    >
   ,typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr1>::type
  >
>::type
compose (const Function& f, const Expr1& expr1, const Expr2& expr2)
{
  typedef typename promote<
    typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr1>::type::value_type
   ,Expr2
  >::type value_type;
  typedef typename details::function_traits<Function>::functor_type wrap_fun_t;
  typedef details::binder_second<wrap_fun_t,value_type>             binded_fun_t;
  typedef typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr1>::type wrap1_t;
  return details::field_expr_v2_nonlinear_node_unary
	   <binded_fun_t,                       wrap1_t>
           (binded_fun_t(wrap_fun_t(f), expr2), wrap1_t(expr1));
}

} // namespace rheolef
#endif // _RHEOLEF_FIELD_EXPR_RECURSIVE_H