librheolef-dev 6.5-1+b1 / usr / include / rheolef / field_expr.h

#ifndef _RHEOLEF_FIELD_EXPR_H
#define _RHEOLEF_FIELD_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
/// 
/// =========================================================================
//
// expressions of field without temporaries, based on boost::proto
// and that uses an iterator for the assignment
//
// author: Pierre.Saramito@imag.fr
//
// date: 23 february 2011
//
#include "rheolef/field.h"
#include "rheolef/field_component.h"
#include "rheolef/promote.h"
#include "rheolef/dis_accumulate.h"
#include "rheolef/dis_inner_product.h"
#include "rheolef/pretty_name.h"
#include "rheolef/form.h"

#include <boost/mpl/bool.hpp>
#include <boost/proto/core.hpp>
#include <boost/proto/debug.hpp>
#include <boost/proto/context.hpp>
#include <boost/proto/transform.hpp>
#include <boost/utility/enable_if.hpp>
#include <boost/typeof/std/vector.hpp>
#include <boost/typeof/std/complex.hpp>
#include <boost/type_traits/remove_reference.hpp>

namespace rheolef {
namespace mpl   = boost::mpl;
namespace proto = boost::proto;
using proto::_;

// forward declaration:
template <class Expr> struct field_expr;
template <class Expr> class field_nonlinear_expr;
template <class T, class M, class VfTag> class test_basic;
template <class RawExpr, class VfTag> class field_vf_expr;
template <class RawExpr> class form_vf_expr;

// -----------------------------------------------------------
// iterator:
// 	expr.begin_dof()
// -----------------------------------------------------------
template<class Iterator>
struct field_iterator_wrapper {
    typedef Iterator iterator;
    explicit field_iterator_wrapper (Iterator iter1) : iter(iter1) {} 
    mutable Iterator iter;
};
struct field_begin : proto::callable {
    template<class Sig>
    struct result;

    template<class This, class Container>
    struct result<This(Container)> : proto::result_of::as_expr
       <
        field_iterator_wrapper<
	  typename boost::remove_reference<Container>::type::const_iterator>
       > {};

    template<class Container>
    typename result<field_begin(const Container&)>::type
    operator ()(const Container& cont) const {
        field_iterator_wrapper<typename Container::const_iterator> iter (cont.begin_dof());
        return proto::as_expr (iter);
    }
};
// Here is the grammar that replaces field terminals with their begin iterators
struct field_grammar_begin
  : proto::or_<
        proto::when <proto::terminal<field_basic<_, _> >,                field_begin(proto::_value)>
      , proto::when <proto::terminal<field_indirect<_, _> >,             field_begin(proto::_value)>
      , proto::when <proto::terminal<field_indirect_const<_, _> >,       field_begin(proto::_value)>
      , proto::when <proto::terminal<field_component<_, _> >,            field_begin(proto::_value)>
      , proto::when <proto::terminal<field_component_const<_, _> >,      field_begin(proto::_value)>
      , proto::when <proto::terminal<_> >
      , proto::when <proto::nary_expr<_, proto::vararg<field_grammar_begin> > >
    >
{};

// Here is an evaluation context that dereferences iterator terminals.
struct field_dereference_context {
    // Unless this is an iterator terminal, use the default evaluation context
    template<class Expr, class EnableIf = void>
    struct eval : proto::default_eval<Expr, const field_dereference_context> {};

    // Dereference iterator terminals.
    template<class Expr>
    struct eval<Expr,
       typename boost::enable_if<
         proto::matches<Expr, proto::terminal<field_iterator_wrapper<_> > >
        >::type
       > {
        typedef typename proto::result_of::value<Expr>::type IteratorWrapper;
        typedef typename IteratorWrapper::iterator iterator;
        typedef typename std::iterator_traits<iterator>::reference result_type;

        result_type operator() (Expr &expr, const field_dereference_context&) const {
            return *proto::value(expr).iter;
        }
    };
};
// Here is an evaluation context that increments iterator
// terminals.
struct field_increment_context {
    // Unless this is an iterator terminal, use the default evaluation context
    template<class Expr, class EnableIf = void>
    struct eval : proto::null_eval<Expr, const field_increment_context> {};

    // advance iterator terminals.
    template<class Expr> struct eval<Expr,
	typename boost::enable_if<
            proto::matches<Expr, proto::terminal<field_iterator_wrapper<_> > >
        >::type
      > {
        typedef void result_type;

        result_type operator() (Expr &expr, const field_increment_context&) const {
            ++proto::value(expr).iter;
        }
    };
};
// -----------------------------------------------------------
// Here is an evaluation context that indexes into a field
// expression and combines the result:
// 	expr.dof(i)
// -----------------------------------------------------------
struct field_subscript_context {
    typedef std::size_t size_type;
    field_subscript_context (size_type i) : _i(i) {}
    // Unless this is a field terminal, use the default evaluation context
    template<class Expr, class EnableIf = void>
    struct eval : proto::default_eval<Expr, const field_subscript_context> {};
    // Index field terminals with our subscript.
    template<class Expr>
    struct eval<Expr, typename boost::enable_if<
        mpl::or_<
            proto::matches<Expr, proto::terminal<field_basic<_, _> > >
           ,proto::matches<Expr, proto::terminal<field_indirect<_, _> > >
           ,proto::matches<Expr, proto::terminal<field_indirect_const<_, _> > >
           ,proto::matches<Expr, proto::terminal<field_component<_, _> > >
	   ,proto::matches<Expr, proto::terminal<field_component_const<_, _> > >
        >
      >::type>
    {
             typedef typename proto::result_of::value<Expr>::type::value_type result_type;
             result_type operator() (Expr& expr, const field_subscript_context& ctx) const
             {
                 return proto::value(expr).dof(ctx._i);
             }
    };
    // Index iterator on field terminals
    template<class Expr>
    struct eval<Expr, typename boost::enable_if<
         proto::matches<Expr, proto::terminal<field_iterator_wrapper<_> > > >::type > {
        typedef typename proto::result_of::value<Expr>::type IteratorWrapper;
        typedef typename IteratorWrapper::iterator iterator;
        typedef typename std::iterator_traits<iterator>::reference result_type;
        result_type operator() (Expr &expr, const field_subscript_context& ctx) const {
            return proto::value(expr).iter [ctx._i];
        }
    };
    mutable size_type _i;
};
// -----------------------------------------------------------
// get size
// 	expr.ndof()
// -----------------------------------------------------------
struct field_get_size_context {
    typedef std::size_t size_type;
// allocator:
    field_get_size_context () {}

// eval:
    template<class Expr, class EnableIf = void>
    struct eval : proto::null_eval<Expr, const field_get_size_context> {};

    template<class Expr>
    struct eval<Expr, typename boost::enable_if<
        mpl::or_<
            proto::matches<Expr, proto::terminal<field_basic<_, _> > >
           ,proto::matches<Expr, proto::terminal<field_indirect<_, _> > >
           ,proto::matches<Expr, proto::terminal<field_indirect_const<_, _> > >
           ,proto::matches<Expr, proto::terminal<field_component<_, _> > >
           ,proto::matches<Expr, proto::terminal<field_component_const<_, _> > >
        >
      >::type>
    {
        typedef void result_type;
        result_type operator() (Expr& expr, const field_get_size_context& ctx) const
        {
            ctx._ownership = proto::value(expr).ownership();
        }
    };
// accessors:
    size_type        ndof() const { return _ownership.size(); };
    distributor ownership() const { return _ownership; };
// data:
    mutable distributor _ownership;
};
// -----------------------------------------------------------
// get space
// 	expr.get_space()
// -----------------------------------------------------------
template<class T, class M>
struct field_get_space_context {
    typedef typename field_basic<T,M>::space_type space_type;
// allocator:
    field_get_space_context () {}

// eval:
    template<class Expr, class EnableIf = void>
    struct eval : proto::null_eval<Expr, const field_get_space_context<T,M> > {};

    template<class Expr>
    struct eval<Expr, typename boost::enable_if<
        mpl::or_<
            proto::matches<Expr, proto::terminal<field_basic<_, _> > >
           ,proto::matches<Expr, proto::terminal<field_indirect<_, _> > >
           ,proto::matches<Expr, proto::terminal<field_indirect_const<_, _> > >
           ,proto::matches<Expr, proto::terminal<field_component<_, _> > >
	   ,proto::matches<Expr, proto::terminal<field_component_const<_, _> > >
        >
      >::type>
    {
        typedef void result_type;
        result_type operator() (Expr &expr, const field_get_space_context<T,M>& ctx) const
        {
            ctx._V = proto::value(expr).get_space();
        }
    };
// accessors:
    space_type get_space() const { return _V; };
// data:
    mutable space_type _V;
};
// -----------------------------------------------------------
// check that stamp match
// 	expr.stamp()
// -----------------------------------------------------------
struct field_check_stamp_context {
    field_check_stamp_context (std::string stamp) : _stamp(stamp) {}

    template<class Expr, class EnableIf = void>
    struct eval : proto::null_eval<Expr, const field_check_stamp_context> {};

    template<class Expr>
    struct eval<Expr, typename boost::enable_if<
        mpl::or_<
            proto::matches<Expr, proto::terminal<field_basic<_, _> > >
           ,proto::matches<Expr, proto::terminal<field_indirect<_, _> > >
           ,proto::matches<Expr, proto::terminal<field_indirect_const<_, _> > >
           ,proto::matches<Expr, proto::terminal<field_component<_, _> > >
           ,proto::matches<Expr, proto::terminal<field_component_const<_, _> > >
        >
      >::type>
    {
        typedef void result_type;
        result_type operator() (Expr &expr, const field_check_stamp_context& ctx) const
        {
            check_macro (ctx._stamp == proto::value(expr).stamp(),
                "incompatible spaces " << ctx._stamp << " and " 
		  << proto::value(expr).stamp() << " in field expression ");
        }
    };
    mutable std::string _stamp;
};
// -----------------------------------------------------------
// A grammar which matches all the assignment operators,
// so we can easily disable them.
// -----------------------------------------------------------
struct field_assign_operators : proto::switch_<struct field_assign_operators_cases> {};

// Here are the cases used by the switch_ above.
struct field_assign_operators_cases {
    template<class Tag, int D = 0> struct case_  : proto::not_<_> {};

    template<int D> struct case_< proto::tag::plus_assign, D >         : _ {};
    template<int D> struct case_< proto::tag::minus_assign, D >        : _ {};
    template<int D> struct case_< proto::tag::multiplies_assign, D >   : _ {};
    template<int D> struct case_< proto::tag::divides_assign, D >      : _ {};
};
// An expression conforms to the field_grammar if it is a terminal or some
// op that is not an assignment op. (Assignment will be handled specially.)
// expr : +expr | -expr | expr + expr | expr - expr
//      | c*expr | expr*c | expr/c
//      | f(expr) | f(expr,expr)
//      | += expr | -= expr
struct field_terminal :
          proto::or_<
            proto::terminal<field_basic<_, _> >
           ,proto::terminal<field_indirect<_, _> >
           ,proto::terminal<field_indirect_const<_, _> >
           ,proto::terminal<field_component<_, _> >
           ,proto::terminal<field_component_const<_, _> >
          >
{};

struct field_constant_expr :
        proto::and_<
        proto::terminal<_>
       ,proto::not_<
          proto::or_< // "proto::or" limited to 6 or 8 args
            proto::or_<
              proto::terminal<field_basic<_, _> >
             ,proto::terminal<field_indirect<_, _> >
             ,proto::terminal<field_indirect_const<_, _> >
             ,proto::terminal<field_component<_, _> >
             ,proto::terminal<field_component_const<_, _> >
            >
           ,proto::or_<
              proto::terminal<form_basic<_, _> >
             ,proto::terminal<field_nonlinear_expr<_> >
             ,proto::terminal<test_basic<_,_,_> >
             ,proto::terminal<field_vf_expr<_,_> >
             ,proto::terminal<form_vf_expr<_> >
            >
          >
        >
      >
{};
struct field_grammar
  : proto::or_< // "proto::or" limited to 6 or 8 args
      proto::or_<
        field_terminal
       ,field_constant_expr
       ,proto::unary_plus <field_grammar>	// -expr
       ,proto::negate     <field_grammar>	// +expr
       ,proto::plus       <field_grammar, field_grammar>
       ,proto::minus      <field_grammar, field_grammar>
      >
     ,proto::or_<
        proto::multiplies <field_constant_expr, field_grammar> // as 2*uh
       ,proto::multiplies <field_grammar, field_constant_expr> // as uh*2
       ,proto::divides    <field_grammar, field_constant_expr> // as uh/2
      >
    >
{};

// Expressions in the field_domain will be wrapped in field_expr<>
// and must conform to the field_grammar
struct field_domain : proto::domain<proto::generator<field_expr>, field_grammar> {};

// -------------------------------------------------------------------
// Here is field_expr, a wrapper for expression types in the field_domain.
// -------------------------------------------------------------------
template<class Expr>
struct field_expr : proto::extends<Expr, field_expr<Expr>, field_domain> {

// typedefs :

    typedef field_get_size_context::size_type size_type;
    typedef typename boost::result_of<field_grammar_begin(const Expr&)>::type raw_iterator;
    typedef rheo_default_memory_model memory_type; // TODO: compute it at compile time

// allocators:

    explicit field_expr (const Expr& expr)
      : proto::extends<Expr, field_expr<Expr>, field_domain>(expr)
    {}

// accessors:

    // TODO: space_type get_space() const; 
    // pb: space_type depends upon Expr
    distributor ownership() const {
        const field_get_size_context get_size;
        proto::eval (*this, get_size);
        return get_size.ownership();
    }
    size_type ndof() const {
        return ownership().size();
    }
    size_type dis_ndof() const {
        return ownership().dis_size();
    }
    struct const_iterator_begin : raw_iterator {

       typedef std::input_iterator_tag        iterator_category;
       typedef typename proto::result_of::eval<raw_iterator,field_dereference_context>::type value_type;
       typedef value_type reference;
       typedef value_type* pointer;
       typedef std::ptrdiff_t difference_type;

#if BOOST_VERSION < 104601
       const_iterator_begin (const field_expr<Expr>& expr)
         : raw_iterator(field_grammar_begin() (expr))
       {}
#else // BOOST_VERSION < 104601
       const_iterator_begin (const field_expr<Expr>& expr)
         : raw_iterator(field_grammar_begin() (expr.proto_expr_))
           // ICI: il faut recuperer base_type(expr) de type Expr au lieu field_expr<Expr>
       {}
#endif // BOOST_VERSION < 104601
       const_iterator_begin& operator++ () { 	
         const field_increment_context inc = {};
         proto::eval (*this, inc);
	 return *this;
       }
       const_iterator_begin operator++ (int) { 
         const_iterator_begin tmp = *this;
	 operator++();
	 return tmp;
       }
       reference operator* () {
         const field_dereference_context deref = {};
         typedef typename proto::result_of::eval<raw_iterator,field_dereference_context>::type my_ref;
         reference value = proto::eval(*this, deref);
         return value;
       }
    };
    typedef const_iterator_begin const_iterator;
    const_iterator begin_dof() const { return const_iterator(*this); }

// typedefs:

    typedef typename proto::result_of::eval<const Expr, const field_subscript_context>::type value_type;
    typedef value_type result_type;
    typedef typename scalar_traits<value_type>::type scalar_type;
    typedef typename float_traits<value_type>::type  float_type;

// random access:

    scalar_type dof (size_type i) const {
        const field_subscript_context get_subscript(i);
        return proto::eval(*this, get_subscript);
    }

#ifdef TODO
    scalar_type max () const {
      scalar_type val = std::numeric_limits<scalar_type>::min();
      for (const_iterator iter = begin_dof(), last = end_dof(); iter != last; iter++) {
        val = std::max(val, *iter);
      }
      return val;
    }
    scalar_type max_abs () const { return abs(*this).max(); }
    scalar_type min_abs () const { return abs(*this).min(); }
#endif // TODO
private:
    // hide this:
    using proto::extends<Expr, field_expr<Expr>, field_domain>::operator[];
};
template <class Expr>
inline
odiststream& operator<< (odiststream& ops, const field_expr<Expr>& uh)
{
    typedef typename field_expr<Expr>::scalar_type T;
    field_basic<T> tmp = uh; // TODO: add M as <T,M>
    return tmp.put (ops);
}
// Define a trait type for detecting field terminals, to
// be used by the BOOST_PROTO_DEFINE_OPERATORS macro below.
template<class T>
struct is_field : mpl::false_ {};

template<class T, class M>
struct is_field<field_basic<T, M> > : mpl::true_ {};

template<class T, class M>
struct is_field<field_indirect<T, M> > : mpl::true_ {};

template<class T, class M>
struct is_field<field_indirect_const<T, M> > : mpl::true_ {};

template<class T, class M>
struct is_field<field_component<T, M> > : mpl::true_ {};

template<class T, class M>
struct is_field<field_component_const<T, M> > : mpl::true_ {};

// -------------------------------------------
// x = expr
// -------------------------------------------
namespace field_detail {
    template<class FwdIter, class Expr, class Op>
    void evaluate (FwdIter begin, FwdIter end, const Expr& expr, Op op) {
        const field_increment_context   inc   = {};
        const field_dereference_context deref = {};
        typename boost::result_of<field_grammar_begin(const Expr&)>::type expr_begin
	  = field_grammar_begin() (expr);
        for (; begin != end; ++begin) {
            op (*begin, proto::eval(expr_begin, deref));
            proto::eval (expr_begin, inc);
        }
    }
    struct assign_op {
        template<class T, class U>
        void operator() (T &t, const U &u) const { t = u; }
    };
} // namespace field_detail
template<class T, class M>
template<class Expr>
inline
field_basic<T,M>& 
field_basic<T,M>::operator= (const field_expr<Expr>& expr) {
  dis_dof_update_needed();
  // get and check sizes:
  const field_get_size_context get_size;
  proto::eval (proto::as_expr<field_domain>(expr), get_size);
  size_type   expr_size      = get_size.ndof();
  if (stamp() == "") {
    const field_get_space_context<T,M> get_space;
    proto::eval (proto::as_expr<field_domain>(expr), get_space);
    space_type V = get_space.get_space();
    resize (V);
  } else {
    const field_check_stamp_context check_stamp (stamp());
    proto::eval (proto::as_expr<field_domain>(expr), check_stamp); // error if the spaces arent compatibles
  }
  // perform all computations here:
  field_detail::evaluate (begin_dof(), end_dof(), proto::as_expr<field_domain>(expr), field_detail::assign_op());
  return *this;
}
template<class T, class M>
template<class Expr>
inline
field_component<T,M>& 
field_component<T,M>::operator= (const field_expr<Expr>& expr) {
    // get and check sizes:
    const field_get_size_context get_size;
    proto::eval (proto::as_expr<field_domain>(expr), get_size);
    size_type   expr_size      = get_size.ndof();
    const field_check_stamp_context check_stamp (stamp());
    proto::eval (proto::as_expr<field_domain>(expr), check_stamp); // error if the spaces arent compatibles
    field_detail::evaluate (begin_dof(), end_dof(), proto::as_expr<field_domain>(expr), field_detail::assign_op());
    return *this;
}
template<class T, class M>
template<class Expr>
inline
field_basic<T,M>::field_basic (const field_expr<Expr>& expr)
 : _V (),
   _u (),
   _b (),
   _dis_dof_update_needed(true)
{
    operator= (expr);
}
// ---------------------------------------------------------------------------
// TODO: form := diag(field_expr) => how to deduce M from field_expr<Expr> ?
// ---------------------------------------------------------------------------
template<class Expr>
inline
form_basic<typename field_expr<Expr>::scalar_type, rheo_default_memory_model>
diag (const field_expr<Expr>& expr)
{
  typedef typename field_expr<Expr>::scalar_type T;
  return diag(field_basic<T,rheo_default_memory_model>(expr));
}

} // namespace rheolef
#endif // _RHEOLEF_FIELD_EXPR_H