librheolef-dev 6.6-1build2 / 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/promote.h"
#include "rheolef/field_expr.h"

#include <boost/function_types/function_type.hpp>
#include <boost/function_types/parameter_types.hpp>
#include <boost/function_types/result_type.hpp>
#include <boost/function_types/function_arity.hpp>

namespace rheolef {

// ------------------------------------------
// tools for creating index lists (TODO: move in utilities)
// ------------------------------------------
// TODO: C++2014 introduced index_sequence : test configure, etc
namespace details {

// the structure that encapsulates index lists
template <size_t... Is>
struct index_list {};

// Collects internal details for generating index ranges [MIN, MAX)
    // declare primary template for index range builder
    template <size_t MIN, size_t N, size_t... Is>
    struct range_builder;

    // base step
    template <size_t MIN, size_t... Is>
    struct range_builder<MIN, MIN, Is...>
    {
        typedef index_list<Is...> type;
    };

    // induction step
    template <size_t MIN, size_t N, size_t... Is>
    struct range_builder: public range_builder<MIN, N - 1, N - 1, Is...>
    {
    };

// Meta-function that returns a [MIN, MAX[ index range
template<size_t MIN, size_t MAX>
using index_range = typename range_builder<MIN, MAX>::type;

} // namespace details
// ---------------------------------------------------------------------------
// functor traits (TODO: move in utilities)
// ---------------------------------------------------------------------------
namespace details {

template <typename T>
struct functor_traits : public functor_traits<decltype(&T::operator())> {};

template <typename C, typename Ret, typename... Args>
struct functor_traits<Ret(C::*)(Args...) const> {
    using result_type =  Ret;
    template <std::size_t i>
    struct arg {
        using type = typename std::tuple_element<i, std::tuple<Args...> >::type;
        using decay_type = typename std::decay<type>::type;
    };
};

template <typename T1, typename T2>
struct decay_is_equivalent
  : std::is_same<
	typename std::decay<T1>::type, 
	typename std::decay<T2>::type
    >::type 
{};


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

template<class NaryFunctor, class... Exprs>
class field_nonlinear_expr_nf {
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;
#ifdef TO_CLEAN
  typedef typename NaryFunctor::result_type       result_type;
#endif // TO_CLEAN

  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_nonlinear_expr_nf (const NaryFunctor& f, const Exprs&... exprs)
    : _f(f), _exprs(exprs...) {}

#ifdef TO_CLEAN
  template<class TrueFunction, class Check = typename std::enable_if<std::is_function<TrueFunction>::value, TrueFunction>::type>
  explicit field_nonlinear_expr_nf (TrueFunction f, const Exprs&... exprs)
    : _f(std::function<TrueFunction>(f)), _exprs(exprs...) {}
#endif // TO_CLEAN

// 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 {
#ifdef TO_CLEAN
    using traits = functor_traits<typename std::decay<NaryFunctor>::type>;
    // check function return type vs Result
    using return_type = typename nary_functor_traits::result_type;
#endif // TO_CLEAN
    bool are_equivalent = (decay_is_equivalent<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).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;
};
// ------------------------------------------------------
// field_expr helpers, for filtering arguments
// TODO: move for all operators in exprs
// ------------------------------------------------------

template <class T>
struct is_field {
  static const bool value = false;
};
template <class T, class M>
struct is_field <field_basic<T,M> > {
  static const bool value = true;
};


template <class T>
struct is_field_constant {
  static const bool value = false;
};
template <>
struct is_field_constant<int> {
  static const bool value = true;
  typedef int type;
};

template <class F, class Check = F>
struct field_expr_traits {
  static const bool is_valid = false;
};
template <class RawExpr>
struct field_expr_traits <rheolef::field_nonlinear_expr<RawExpr> > {
   static const bool is_valid = true;
   typedef rheolef::field_nonlinear_expr<RawExpr> checked_type;
   typedef rheolef::field_nonlinear_expr<RawExpr> wrapped_type;
};
template <class T, class M>
struct field_expr_traits <field_basic<T,M> > {
   static const bool is_valid = true;
   typedef field_basic<T,M>               checked_type;
   typedef field_expr_terminal_field<T,M> wrapped_type;
};
template <class T, class M>
struct field_expr_traits <field_indirect<T,M> > {
   static const bool is_valid = true;
   typedef field_indirect<T,M>            checked_type;
   typedef field_expr_terminal_field<T,M> wrapped_type;
};
template <class T, class M>
struct field_expr_traits <field_indirect_const<T,M> > {
   static const bool is_valid = true;
   typedef field_indirect_const<T,M>      checked_type;
   typedef field_expr_terminal_field<T,M> wrapped_type;
};
template <class T, class M>
struct field_expr_traits <field_component<T,M> > {
   static const bool is_valid = true;
   typedef field_component<T,M>           checked_type;
   typedef field_expr_terminal_field<T,M> wrapped_type;
};
template <class T, class M>
struct field_expr_traits <field_component_const<T,M> > {
   static const bool is_valid = true;
   typedef field_component_const<T,M>     checked_type;
   typedef field_expr_terminal_field<T,M> wrapped_type;
};
template <class Expr>
struct field_expr_traits <field_expr<Expr> > {
   static const bool is_valid = true;
   typedef field_expr<Expr>           checked_type;
   typedef field_expr_terminal_field<typename field_expr<Expr>::scalar_type, typename field_expr<Expr>::memory_type> wrapped_type;
};
template <class Result>
struct field_expr_traits <Result(const point_basic<typename float_traits<Result>::type>&)> {
   static const bool is_valid = true;
   typedef Result(&checked_type)(const point_basic<typename float_traits<Result>::type>&);
   typedef std::pointer_to_unary_function<const point_basic<typename float_traits<Result>::type>&,Result> function_type;
   typedef field_expr_terminal_function<function_type>  wrapped_type;
};
template <class UnaryFunctor>
struct field_expr_traits
   <
       UnaryFunctor,
       typename std::enable_if
       <
	   // SFINAE filter: arity==1 and first argument_type==point
           is_point<typename UnaryFunctor::argument_type>::value
           && ! is_field<UnaryFunctor>::value // field class is also an unary_function but with different treatment in expr
          ,UnaryFunctor
       >::type
   > {
   static const bool is_valid = true;
   typedef UnaryFunctor                               checked_type;
   typedef field_expr_terminal_function<UnaryFunctor> wrapped_type;
};
template <typename Constant>
struct field_expr_traits
   <
       Constant,
       typename std::enable_if <std::is_arithmetic<Constant>::value, Constant>::type
   > {
   static const bool is_valid = true;
   typedef typename upgrade_integral_to_float<Constant>::type constant_type;
   typedef constant_type                                      checked_type;
   typedef typename float_traits<constant_type>::type         float_type;
   typedef f_constant<point_basic<float_type>,constant_type>  function_type;
   typedef field_expr_terminal_function<function_type>        wrapped_type;
};

template <class Functor, class Check = Functor>
struct wrap_function_traits {
  typedef Functor type;
};
template <class TrueFunction>
struct wrap_function_traits<TrueFunction, typename std::enable_if<std::is_function<TrueFunction>::value, TrueFunction>::type> {
  typedef std::function<TrueFunction> type;
};

} // namespace details
// ------------------------------------------------------
// compose(f,u1...uN)
// ------------------------------------------------------
template<class Function, class... Exprs>
inline
typename std::enable_if <
  sizeof...(Exprs) >= 3,
  field_nonlinear_expr<
    details::field_nonlinear_expr_nf<
      typename details::wrap_function_traits<Function>::type,
      typename details::field_expr_traits<Exprs>::wrapped_type...
    >
  >
>::type
compose (const Function& f, const Exprs&... exprs) {
  typedef typename details::wrap_function_traits<Function>::type   fun_t;
  typedef details::field_nonlinear_expr_nf<fun_t, typename details::field_expr_traits<Exprs>::wrapped_type...>   exprs_t;
  return  field_nonlinear_expr<exprs_t> (exprs_t(fun_t(f), typename details::field_expr_traits<Exprs>::wrapped_type(exprs)...));
}

} // namespace rheolef
#endif // _RHEOLEF_COMPOSE_H