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

#ifndef _RHEOLEF_COMPOSE_H
#define _RHEOLEF_COMPOSE_H
//
// This file is part of Rheolef.
//
// Copyright (C) 2000-2009 Pierre Saramito 
//
// 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
// 
// ==========================================================================
// 
// compose a n-ary function with n fields : compose(f,uh1...uhN)
//
// author: Pierre.Saramito@imag.fr
//
// date: 4 september 2015
//
//<compose:  
/*Class:compose
NAME:  @code{compose} - a n-ary function with n fields
DESCRIPTION:       
  @noindent
  Compose a n-ary function f with n fields.
  @example
        geo omega ("circle");
        space Xh (omega, "P1");
        field uh (Xh, 1.0);
        field vh = interpolate (compose(f,uh));
  @end example
  The @code{compose} operator could be used in all non-linear expressions
  involved in either the @code{interpolate} or the @code{integrate} functions
  (see @ref{interpolate algorithm} and @ref{integrate algorithm}).
  The @code{f} function could be either a usual function or a functor.
CHARACTERISTIC:
  The @code{compose} function supports also the characteristic algorithm
  (see @ref{characteristic class})
  used for convection.
  @example
        characteristic X (-delta_t*uh);
        test v (Xh);
        field lh = integrate (compose(uh,X)*v);
  @end example
IMPLEMENTATION:
  The n-arity bases on the variadic template feature of the 2011 c++ normalisation.
  When this feature is not available, only unary and binary functions are supported.
AUTHOR: Pierre.Saramito@imag.fr
DATE:   3 september 2015
End:
*/
#include "rheolef/field_expr_v2_nonlinear.h"
namespace rheolef {

// ---------------------------------------------------------------------------
// N-ary function call: (f expr1...exprN) , N >= 3 only
// ---------------------------------------------------------------------------
namespace details {

template<class NaryFunctor, class... Exprs>
class field_expr_v2_nonlinear_node_nary {
public:
// constants:

  static const size_t N = sizeof...(Exprs);
  typedef typename range_builder<0,N>::type IndexRange;

// typedefs:

  typedef geo_element::size_type                   size_type;
  using nary_functor_traits = functor_traits<typename std::decay<NaryFunctor>::type>;
  using result_type = typename nary_functor_traits::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 rheo_default_memory_model                memory_type;
#ifdef TODO
  // TODO: extract first type Expr1 from Exprs (HOWTO extract ?) ; 
  //       also, check that all args have the same memory model
  typedef typename Expr1::memory_type              memory_type;
#endif // TODO

// alocators:

  field_expr_v2_nonlinear_node_nary (const NaryFunctor& f, const Exprs&... exprs)
    : _f(f), _exprs(exprs...) {}

// 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 valued_hint; // when N >= 3 : return type should be solved at compile time
#ifdef TODO
    // TODO: when N=1,2 : possible unsolved return type until run-time:
    return details::generic_binary_traits<NaryFunctor>::valued_tag(_expr1.valued_tag(), _expr2.valued_tag());
#endif // TODO
  }

// initializers:

  template<size_t ...Is>
  bool _initialize_internal (const geo_basic<float_type,memory_type>& omega, const quadrature<float_type>& hat_x, index_list<Is...>) const {
    bool status_list[] = {std::get<Is>(_exprs).initialize (omega, hat_x)...};
    bool status = true;
    for (bool status_i : status_list) {
      status &= status_i;
    }
    return status;
  }
  bool initialize (const geo_basic<float_type,memory_type>& omega, const quadrature<float_type>& hat_x) const {
    return _initialize_internal (omega, hat_x, IndexRange());
  }
  template<size_t ...Is>
  bool _initialize_internal (const space_basic<float_type,memory_type>& Xh, index_list<Is...>) const {
    bool status_list[] = {std::get<Is>(_exprs).initialize (Xh)...};
    bool status = true;
    for (bool status_i : status_list) { status &= status_i; }
    return status;
  }
  bool initialize (const space_basic<float_type,memory_type>& Xh) const {
    return _initialize_internal (Xh, IndexRange());
  }

// evaluators:

  template<class Result, size_t ...Is>
  bool _evaluate_internal (const geo_element& K, std::vector<Result>& value, index_list<Is...>) const {
    using traits = functor_traits<typename std::decay<NaryFunctor>::type>;
    typedef std::tuple<std::vector<typename std::decay<typename traits::template arg<Is>::type>::type>...> vec_args_type;
    vec_args_type tmps;
    bool status_list[] = {std::get<Is>(_exprs).evaluate (K, std::get<Is>(tmps))...};
    value.resize (std::get<0>(tmps).size());
    for (size_type i = 0, n_value = value.size(); i < n_value; ++i) {
      value[i] = _f (std::get<Is>(tmps)[i]...);
    }
    bool status = true;
    for (bool status_i : status_list) { status &= status_i; }
    return status;
  }
  template<class Result>
  bool evaluate (const geo_element& K, std::vector<Result>& value) const {
    return _evaluate_internal (K, value, IndexRange());
  }
  template<class Result, size_t ...Is>
  bool _valued_check_internal (Result, index_list<Is...>) const {
    bool are_equivalent = (decay_is_same<Result,result_type>::value);
    check_macro (are_equivalent,
	"compose; incompatible function " << typename_macro(NaryFunctor) 
        << " return value " << typename_macro(result_type)
	<< " and expected value " << typename_macro(Result));
    // check function argument type vs Exprs return types via recursive calls
    bool status_list[] = { std::get<Is>(_exprs).template valued_check<
    	typename nary_functor_traits::template arg<Is>::decay_type>()... };
    bool status = true;
    for (bool status_i : status_list) { status &= status_i; }
    return status;
  }
  template<class Result>
  bool valued_check() const {
    return _valued_check_internal (Result(), IndexRange());
  }
protected:
// data:
  NaryFunctor           _f;
  std::tuple<Exprs...>  _exprs;
};
template<class F, class... Exprs> struct is_field_expr_v2_nonlinear_arg     <field_expr_v2_nonlinear_node_nary<F,Exprs...> > : std::true_type {};

} // namespace details
// ------------------------------------------------------
// compose(f,u1...uN)
// ------------------------------------------------------
// TODO: check that Function is a valid n-ary function or functor
// TODO: check that args are valid field_expr or field_constant
//	details::and_type<details::is_field_expr_v2_nonlinear_arg<Exprs>...>::value
// TODO: when i-th arg is field_constant, use bind to support it
// TODO: possibly undetermined return type until run-time, when N=1,2
// TODO: then, unary-compose and binary-compose can be removed
template<class Function, class... Exprs>
inline
typename std::enable_if <
  sizeof...(Exprs) >= 3,
  details::field_expr_v2_nonlinear_node_nary<
    typename details::function_traits<Function>::functor_type
   ,typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Exprs>::type...
  >
>::type
compose (const Function& f, const Exprs&... exprs) {
  typedef typename details::function_traits<Function>::functor_type   fun_t;
  return details::field_expr_v2_nonlinear_node_nary
	<fun_t,    typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Exprs>::type...>
        (fun_t(f), typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Exprs>::type(exprs)...);
}

} // namespace rheolef
#endif // _RHEOLEF_COMPOSE_H