librheolef-dev 6.5-1build1 / usr / include / rheolef / field_nonlinear_expr.h

#ifndef _RHEOLEF_FIELD_NONLINEAR_EXPR_H
#define _RHEOLEF_FIELD_NONLINEAR_EXPR_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_nonlinear_expr: separate from field_expr
//
// it is a prototype: only compose(f,uh) is still handled
//
// OUTLINE:
// 1) main wrapper class
// 2) unary function call: (f expr)
// 3) binary function call: (f expr1 expr2)
// 4) jump of a field
#include "rheolef/field_nonlinear_expr_terminal.h"
#include "rheolef/field_evaluate.h"
#include "rheolef/tensor4.h"
#include "rheolef/field_functor.h"
#include "rheolef/operators.h"

#include <boost/functional.hpp>

namespace rheolef {

// ---------------------------------------------------------------------------
// 1) main wrapper class
// ---------------------------------------------------------------------------
template<class RawExpr>
class field_nonlinear_expr {
public:
// typedefs:

  typedef geo_element::size_type         size_type;
  typedef typename RawExpr::memory_type  memory_type;
  typedef typename RawExpr::result_type  result_type;
  typedef typename RawExpr::value_type   value_type;
  typedef typename RawExpr::float_type   float_type;
  typedef typename RawExpr::scalar_type  scalar_type;

// alocators:

  field_nonlinear_expr (const RawExpr& raw_expr) 
    : _raw_expr(raw_expr) {}

// accessors:

  static const space_constant::valued_type valued_hint = RawExpr::valued_hint;

  space_constant::valued_type valued_tag() const { return _raw_expr.valued_tag(); }

  void initialize (const geo_basic<float_type,memory_type>& omega, const quadrature<float_type>& hat_x) const {
    _raw_expr.initialize (omega, hat_x);
  }
  bool initialize (const space_basic<float_type,memory_type>& Xh) const {
    return _raw_expr.initialize (Xh);
  }
  template<class Result>
  void evaluate (const geo_element& K, std::vector<Result>& value) const {
    _raw_expr.evaluate (K, value);
  }
  void evaluate_on_side (const geo_element& K, const side_information_type& sid, std::vector<result_type>& value) const {
    _raw_expr.evaluate_on_side (K, sid, value);
  }
  template<class Result>
  void valued_check() const {
    _raw_expr.valued_check<Result>();
  }

protected:
// data:
  RawExpr   _raw_expr;
};
// ---------------------------------------------------------------------------
// 2) unary function call: (f expr)
// ---------------------------------------------------------------------------
template<class UnaryFunction, class Expr>
class field_nonlinear_expr_uf {
public:
// typedefs:

  typedef geo_element::size_type                      size_type;
  typedef typename Expr::memory_type                  memory_type;
  typedef typename details::generic_unary_traits<UnaryFunction>::template result_hint<typename Expr::result_type>::type
						      result_type;
  typedef result_type                                 value_type;
  typedef typename scalar_traits<value_type>::type    scalar_type;
  typedef typename  float_traits<value_type>::type    float_type;
  typedef field_nonlinear_expr_uf<UnaryFunction,Expr> self_type;

// alocators:

  field_nonlinear_expr_uf (const UnaryFunction& f, const field_nonlinear_expr<Expr>& expr) 
    : _f(f), _expr(expr) {}

// 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_unary_traits<UnaryFunction>::valued_tag (_expr.valued_tag());
  }

// initializators:

  void initialize (const geo_basic<float_type,memory_type>& omega, const quadrature<float_type>& hat_x) const {
    _expr.initialize (omega, hat_x);
  }
  bool initialize (const space_basic<float_type,memory_type>& Xh) const {
    return _expr.initialize (Xh);
  }
// evaluator:

  template<class Result, class Arg, class Status>
  struct evaluate_call_check {
    void operator() (const self_type& obj, const geo_element& K, std::vector<Result>& 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,mpl::true_> {
    void operator() (const self_type& obj, const geo_element& K, std::vector<Result>& value) const {
      std::vector<Arg> tmp_value;
      obj._expr.evaluate (K, tmp_value);
      value.resize(tmp_value.size());
      typename std::vector<Arg>::const_iterator tmp = tmp_value.begin();
      for (typename std::vector<Result>::iterator
	iter = value.begin(),
	last = value.end(); iter != last; ++iter, ++tmp) {
          *iter = obj._f(*tmp);
      }
    }
  };
  template<class Result, class Arg>
  void evaluate_call (const geo_element& K, std::vector<Result>& value) const {
    typedef typename details::generic_unary_traits<UnaryFunction>::template hint<Arg,Result>::result_type result_type;
    typedef typename mpl::and_<
		typename details::is_equal<Result,result_type>::type
	       ,typename mpl::not_<typename details::is_error<Arg>::type>::type
	      >::type
      status_t;
    evaluate_call_check<Result,Arg,status_t> eval;
    eval (*this, K, value);
  }
  // when arg is known at run-time:
  template<class This, class Result, class Arg, space_constant::valued_type ArgTag>
  struct evaluate_switch {
    void evaluate (const This& obj, const geo_element& K, std::vector<Result>& 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> (K, value); break;
        case space_constant::vector:
	  obj.template evaluate_call<Result, point_basic<T> > (K, value); break;
        case space_constant::tensor:
        case space_constant::unsymmetric_tensor:
	  obj.template evaluate_call<Result, tensor_basic<T> > (K, value); break;
        default: { error_macro ("unexpected valued tag="<<arg_valued_tag); }
      }
    }
  };
  // when arg is known at compile-time:
  template<class This, class Result, class Arg>
  struct evaluate_switch <This, Result, Arg, space_constant::scalar> {
    typedef typename scalar_traits<Arg>::type T;
    void evaluate (const This& obj, const geo_element& K, std::vector<Result>& value) const {
      obj.template evaluate_call<Result,T> (K, value);
    }
  };
  template<class This, class Result, class Arg>
  struct evaluate_switch <This, Result, Arg, space_constant::vector> {
    typedef typename scalar_traits<Arg>::type T;
    void evaluate (const This& obj, const geo_element& K, std::vector<Result>& value) const {
      obj.template evaluate_call<Result, point_basic<T> > (K, value);
    }
  };
  template<class This, class Result, class Arg>
  struct evaluate_switch <This, Result, Arg, space_constant::tensor> {
    typedef typename scalar_traits<Arg>::type T;
    void evaluate (const This& obj, const geo_element& K, std::vector<Result>& value) const {
      obj.template evaluate_call<Result, tensor_basic<T> > (K, value);
    }
  };
  template<class Result>
  void evaluate (const geo_element& K, std::vector<Result>& value) const
  {
   typedef typename details::generic_unary_traits<UnaryFunction>::template hint<typename Expr::value_type,Result>::argument_type
                     A1;
    typedef field_nonlinear_expr_uf<UnaryFunction, Expr> This;
    static const space_constant::valued_type argument_tag = space_constant::valued_tag_traits<A1>::value;
    evaluate_switch <This, Result, A1, argument_tag> helper;
    helper.evaluate (*this, K, value);
  }
  template<class Result>
  void evaluate_on_side (const geo_element& K, const side_information_type& sid, std::vector<Result>& value) const { 
    // TODO: group with evaluate
    fatal_macro ("uf::evaluate_on_side: not yet");
  }

  template<class Result>
  void 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) _expr.valued_check<A1>();
  }

protected:
// data:
  UnaryFunction                 _f;
  field_nonlinear_expr<Expr>    _expr;
};
// ---------------------------------------------------------------------------
// 3) binary function call: (f expr1 expr2)
// ---------------------------------------------------------------------------
template<class BinaryFunction, class Expr1, class Expr2>
class field_nonlinear_expr_bf {
public:
// typedefs:

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

// alocators:

  field_nonlinear_expr_bf (const BinaryFunction&             f, 
		    const Expr1& expr1,
                    const Expr2& expr2)
    : _f(f), _expr1(expr1), _expr2(expr2)
  {
  }

// 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 geo_basic<float_type,memory_type>& omega, const quadrature<float_type>& hat_x) const
  {
    _expr1.initialize (omega, hat_x); 
    _expr2.initialize (omega, hat_x);
  }
  bool initialize (const space_basic<float_type,memory_type>& Xh) const {
    bool is_homogeneous1 = _expr1.initialize (Xh);
    bool is_homogeneous2 = _expr2.initialize (Xh);
    return is_homogeneous1 && is_homogeneous2;
  }

// evaluators:

  template<class Result, class Arg1, class Arg2>
  void evaluate_internal2 (const geo_element& K, std::vector<Result>& value) const {
#ifdef TO_CLEAN
    _check<Result,Arg1,Arg2> ();
#endif // TO_CLEAN
    std::vector<Arg1> tmp1_value;
    std::vector<Arg2> tmp2_value;
    _expr1.evaluate (K, tmp1_value);
    _expr2.evaluate (K, tmp2_value);
    value.resize(tmp1_value.size());
    typename std::vector<Arg1>::const_iterator tmp1 = tmp1_value.begin();
    typename std::vector<Arg2>::const_iterator tmp2 = tmp2_value.begin();
    for (typename std::vector<Result>::iterator
        iter = value.begin(),
	last = value.end(); iter != last; ++iter, ++tmp1, ++tmp2) {
      *iter = _f (*tmp1, *tmp2);
    }
  }
  template<class This, class Result, class ReturnType, class Arg1, class Arg2>
  struct evaluate_internal {
    void operator() (const This& obj, const geo_element& K, std::vector<Result>& 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> {
    void operator() (const This& obj, const geo_element& K, std::vector<Result>& value) const {
      obj.template evaluate_internal2 <Result,Arg1,Arg2> (K, value);
    }
  };
  template<class Result, class Arg1, class Arg2>
  void evaluate_call (const geo_element& K, std::vector<Result>& value) const {
    typedef typename details::generic_binary_traits<BinaryFunction>::template result_hint<Arg1,Arg2>::type ReturnType;
    typedef field_nonlinear_expr_bf<BinaryFunction, Expr1, Expr2> This;
    evaluate_internal<This,Result,ReturnType,Arg1,Arg2> eval_int;
    eval_int (*this, K, 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 {
    void operator() (const This& obj, const geo_element& K, std::vector<Result>& value) const {
      obj.template evaluate_call<Result, Arg1, Arg2> (K, value);
    }
  };
  // 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> {
    void operator() (const This& obj, const geo_element& K, std::vector<Result>& 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:
	      obj.template evaluate_call<Result, T1, T2>                (K, value); break;
            case space_constant::vector:
	      obj.template evaluate_call<Result, T1, point_basic<T2> >  (K, value); break;
            case space_constant::tensor:
            case space_constant::unsymmetric_tensor:
	      obj.template evaluate_call<Result, T1, tensor_basic<T2> > (K, value); break;
            case space_constant::tensor3:
	      obj.template evaluate_call<Result, T1, tensor3_basic<T2> > (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:
	      obj.template evaluate_call<Result, point_basic<T1>, T2>                (K, value); break;
            case space_constant::vector:
	      obj.template evaluate_call<Result, point_basic<T1>, point_basic<T2> >  (K, value); break;
            case space_constant::tensor:
            case space_constant::unsymmetric_tensor:
	      obj.template evaluate_call<Result, point_basic<T1>, tensor_basic<T2> > (K, value); break;
            case space_constant::tensor3:
	      obj.template evaluate_call<Result, point_basic<T1>, tensor3_basic<T2> > (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:
	      obj.template evaluate_call<Result, tensor_basic<T1>, T2>                (K, value); break;
            case space_constant::vector:
	      obj.template evaluate_call<Result, tensor_basic<T1>, point_basic<T2> >  (K, value); break;
            case space_constant::tensor:
            case space_constant::unsymmetric_tensor:
	      obj.template evaluate_call<Result, tensor_basic<T1>, tensor_basic<T2> > (K, value); break;
            case space_constant::tensor3:
	      obj.template evaluate_call<Result, tensor_basic<T1>, tensor3_basic<T2> > (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:
	      obj.template evaluate_call<Result, tensor3_basic<T1>, T2>                (K, value); break;
            case space_constant::vector:
	      obj.template evaluate_call<Result, tensor3_basic<T1>, point_basic<T2> >  (K, value); break;
            case space_constant::tensor:
            case space_constant::unsymmetric_tensor:
	      obj.template evaluate_call<Result, tensor3_basic<T1>, tensor_basic<T2> > (K, value); break;
            case space_constant::tensor3:
	      obj.template evaluate_call<Result, tensor3_basic<T1>, tensor3_basic<T2> > (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);
      }
    }
  };
  // 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> {
    void operator() (const This& obj, const geo_element& K, std::vector<Result>& 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:
	  obj.template evaluate_call<Result, Arg1, T2>                (K, value); break;
        case space_constant::vector:
	  obj.template evaluate_call<Result, Arg1, point_basic<T2> >  (K, value); break;
        case space_constant::tensor:
        case space_constant::unsymmetric_tensor:
	  obj.template evaluate_call<Result, Arg1, tensor_basic<T2> > (K, value); break;
        case space_constant::tensor3:
	  obj.template evaluate_call<Result, Arg1, tensor3_basic<T2> > (K, value); break;
        default: error_macro ("unexpected second argument valued tag="<<arg2_valued_tag);
      }
    }
  };
  // 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> {
    void operator() (const This& obj, const geo_element& K, std::vector<Result>& 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:
	  obj.template evaluate_call<Result, T1, Arg2>               (K, value); break;
        case space_constant::vector:
	  obj.template evaluate_call<Result, point_basic<T1>, Arg2>  (K, value); break;
        case space_constant::tensor:
        case space_constant::unsymmetric_tensor:
	  obj.template evaluate_call<Result, tensor_basic<T1>, Arg2> (K, value); break;
        case space_constant::tensor3:
	  obj.template evaluate_call<Result, tensor3_basic<T1>, Arg2> (K, value); break;
        default: error_macro ("unexpected first argument valued tag="<<arg1_valued_tag);
      }
    }
  };
  template<class Result>
  void evaluate (const geo_element& K, std::vector<Result>& 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_nonlinear_expr_bf<BinaryFunction, Expr1, Expr2> This;
    evaluate_switch <This, Result, A1, first_argument_tag, A2, second_argument_tag>   eval;
    eval (*this, K, value);
  }
  template<class Result>
  void 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;
    if (! is_undeterminated<A1>::value) _expr1.valued_check<A1>();
    if (! is_undeterminated<A2>::value) _expr2.valued_check<A2>();
  }

protected:
// data:
  BinaryFunction  _f;
  Expr1           _expr1;
  Expr2           _expr2;
};
// ----------------------------------------------------------------------------
// 4) jump of a field
// ----------------------------------------------------------------------------
template<class Expr>
class field_nonlinear_expr_dg {
public:
// typedefs:

  typedef geo_element::size_type          size_type;
  typedef typename Expr::memory_type      memory_type;
  typedef typename Expr::result_type      result_type;
  typedef typename Expr::value_type       value_type;
  typedef typename Expr::scalar_type      scalar_type;
  typedef typename Expr::float_type       float_type;

// alocators:

  field_nonlinear_expr_dg (const Expr& expr, const float_type& c0, const float_type& c1);

// accessors:

  static const space_constant::valued_type valued_hint = Expr::valued_hint;

  space_constant::valued_type valued_tag() const { return _expr0.valued_tag(); }

  void initialize (const geo_basic<float_type,memory_type>& omega, const quadrature<float_type>& quad) const;
  bool initialize (const space_basic<float_type,memory_type>& Xh) const;

  template<class ValueType>
  void evaluate (const geo_element& K, std::vector<ValueType>& value) const;

  template<class ValueType>
  void valued_check() const;

// data:
protected:
  Expr                                      _expr0, _expr1;
  float_type                                _c0,    _c1;
  mutable geo_basic<float_type,memory_type> _omega;
};

template<class Expr>
inline
field_nonlinear_expr_dg<Expr>::field_nonlinear_expr_dg (
  const Expr&       expr, 
  const float_type& c0,
  const float_type& c1)
 : _expr0(expr), _expr1(expr),
   _c0(c0),      _c1(c1),
   _omega()
{
}
template<class Expr>
inline
void
field_nonlinear_expr_dg<Expr>::initialize (
    const geo_basic<float_type,memory_type>& side_dom,
    const quadrature<float_type>&            quad) const
{
  _expr0.initialize (side_dom, quad);
  _expr1.initialize (side_dom, quad);
  _omega = side_dom.get_background_geo();
}
template<class Expr>
inline
bool
field_nonlinear_expr_dg<Expr>::initialize (
    const space_basic<float_type,memory_type>& Xh) const
{
  _expr0.initialize (Xh);
  _expr1.initialize (Xh);
  _omega = Xh.get_geo().get_background_geo();
  return true;
}
template<class Expr>
template<class ValueType>
inline
void 
field_nonlinear_expr_dg<Expr>::valued_check() const 
{
  space_constant::valued_type valued_tag = space_constant::valued_tag_traits<ValueType>::value;
  check_macro (_expr0.valued_tag() == valued_tag,
 	              "unexpected "<<  space_constant::valued_name(_expr0.valued_tag())
        << "-valued field while a " << space_constant::valued_name(valued_tag)
        << "-valued one is expected in expression");
}
template<class Expr>
template<class Result>
void
field_nonlinear_expr_dg<Expr>::evaluate (
    const geo_element&   K,
    std::vector<Result>& value) const
{
  size_type L_map_d = K.dimension() + 1;
  size_type L_dis_ie0, L_dis_ie1;
  side_information_type sid0, sid1;

  L_dis_ie0 = K.master(0);
  L_dis_ie1 = K.master(1);
  check_macro (L_dis_ie0 != std::numeric_limits<size_type>::max(),
      "unexpected isolated mesh side K="<<K);
  if (L_dis_ie1 == std::numeric_limits<size_type>::max()) {
    // K is a boundary side
    const geo_element& L0 = _omega.dis_get_geo_element (L_map_d, L_dis_ie0);
    L0.get_side_informations (K, sid0);
    _expr0.evaluate_on_side (L0, sid0, value);
    // average (i.e. _c0==0.5): fix it on the boundary where c0=1 : average(v)=v on the boundary
    Float c0 = (_c0 != 0.5) ? _c0 : 1;
    for (size_type loc_idof = 0, loc_ndof = value.size(); loc_idof < loc_ndof; ++loc_idof) {
      value[loc_idof] = c0*value[loc_idof];
    }
    return;
  }
  // K is an internal side
  std::vector<Result> value1;
  const geo_element& L0 = _omega.dis_get_geo_element (L_map_d, L_dis_ie0);
  const geo_element& L1 = _omega.dis_get_geo_element (L_map_d, L_dis_ie1);
  L0.get_side_informations (K, sid0);
  L1.get_side_informations (K, sid1);
  _expr0.evaluate_on_side (L0, sid0, value);
  _expr1.evaluate_on_side (L1, sid1, value1);
  for (size_type loc_idof = 0, loc_ndof = value.size(); loc_idof < loc_ndof; ++loc_idof) {
    value[loc_idof] = _c0*value[loc_idof] + _c1*value1[loc_idof];
  }
}

} // namespace rheolef
#endif // _RHEOLEF_FIELD_nONLINEAR_EXPR_H