field.h

Go to the documentation of this file.

# ifndef _RHEOLEF_FIELD_H
# define _RHEOLEF_FIELD_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
// 
// =========================================================================
// AUTHOR: Pierre.Saramito@imag.fr
// DATE:   23 february 2011
 
namespace rheolef {
} // namespace rheolef
 
#include "rheolef/field_rdof.icc"
#include "rheolef/field_wdof.icc"
#include "rheolef/integrate_option.h" // TO_CLEAN with old interface
 
namespace rheolef {
 
namespace details {
 
// forward declarations:
template <class T, class M> class field_concat_value;
template <class FieldWdof> class field_wdof_sliced;
template <class FieldRdof> class field_rdof_sliced_const;
 
} // namespace details
 
// ----------------------------------------------------------------------------
// main class
// ----------------------------------------------------------------------------
// [verbatim_field_basic]
template <class T, class M = rheo_default_memory_model>
class field_basic: public details::field_wdof_base<field_basic<T,M>> {
public :
// typedefs:
 
  using wdof_base     = details::field_wdof_base<field_basic<T,M>>;
  using rdof_base     = details::field_rdof_base<field_basic<T,M>>;
  using size_type     = std::size_t;
  using scalar_type   = T;
  using memory_type   = M;
  using float_type    = typename float_traits<T>::type;
  using geo_type      = geo_basic  <float_type,memory_type>;
  using space_type    = space_basic<float_type,memory_type>;
  using dis_reference = typename vec<scalar_type,memory_type>::dis_reference;
  using valued_type   = space_constant::valued_type;
  using value_type    = T; // TODO: run-time dependent, set it to undeterminated<T>
  using result_type   = T; // TODO: run-time dependent, set it to undeterminated<T>
  class iterator;
  class const_iterator;
 
// allocator/deallocator:
  
  field_basic();
 
  explicit field_basic (
    const space_type& V, 
    const T& init_value = std::numeric_limits<T>::max());
   
  void resize (
    const space_type& V, 
    const T& init_value = std::numeric_limits<T>::max());
 
  // expressions: field_expr or field_lazy are accepted here
  // TODO: merge with FieldLazy piecewise_poly
  template <class Expr, class Sfinae =
    typename std::enable_if<
      (      details::is_field_expr_affine_homogeneous<Expr>::value
        && ! details::is_field<Expr>::value
      )
      ||   details::is_field_lazy<Expr>::value
    >::type>
  field_basic (const Expr& expr);
 
  field_basic (const std::initializer_list<details::field_concat_value<T,M> >& init_list);
 
// assignments:
 
  field_basic<T,M>& operator= (const field_basic<T,M>&);
  field_basic<T,M>& operator= (const std::initializer_list<details::field_concat_value<T,M> >& init_list);
 
// accessors:
 
  const space_type&  get_space()  const { return _V; }
  const geo_type&    get_geo()    const { return _V.get_geo(); }
  std::string        get_approx() const { return _V.get_approx(); }
  valued_type        valued_tag() const { return _V.valued_tag(); }
  const std::string& valued()     const { return _V.valued(); }
#ifdef TO_CLEAN
  std::string        name()      const { return _V.name(); }
  bool have_homogeneous_space (space_basic<T,M>& Xh) const { Xh = get_space(); return true; }
#endif // TO_CLEAN
 
// accessors & modifiers to unknown & blocked parts:
 
  const vec<T,M>&     u() const { return _u; }
  const vec<T,M>&     b() const { return _b; }
        vec<T,M>& set_u()       { dis_dof_indexes_requires_update(); return _u; }
        vec<T,M>& set_b()       { dis_dof_indexes_requires_update(); return _b; }
 
// accessors to extremas:
 
  T min() const;
  T max() const;
  T max_abs() const;
  T min_abs() const;
 
// accessors by degrees-of-freedom (dof):
 
  const distributor& ownership() const { return get_space().ownership(); }
  const communicator& comm() const { return ownership().comm(); }
  size_type     ndof() const { return ownership().size(); }
  size_type dis_ndof() const { return ownership().dis_size(); }
        T& dof (size_type idof);
  const T& dof (size_type idof) const;
  const T& dis_dof (size_type dis_idof) const;
  // write access to non-local dis_idof changes to others procs
  dis_reference dis_dof_entry (size_type dis_idof);
 
  iterator begin_dof();
  iterator end_dof();
  const_iterator begin_dof() const;
  const_iterator end_dof() const;
 
// input/output:
 
  idiststream& get (idiststream& ips);
  odiststream& put (odiststream& ops) const;
  odiststream& put_field (odiststream& ops) const;
 
// evaluate uh(x) where x is given locally as hat_x in K:
 
  T dis_evaluate (const point_basic<T>& x, size_type i_comp = 0) const;
  T operator()   (const point_basic<T>& x) const { return dis_evaluate (x,0); }
  point_basic<T> dis_vector_evaluate (const point_basic<T>& x) const;
 
// [verbatim_field_basic]
// internals:
  int constraint_process_rank() const { return u().constraint_process_rank(); }
public:
 
  // evaluate uh(x) where x is given locally as hat_x in K:
  // requiers to call field::dis_dof_upgrade() before.
  T evaluate (const geo_element& K, const point_basic<T>& hat_xq, size_type i_comp = 0) const;
 
  // propagate changed values shared at partition boundaries to others procs
  template <class SetOp = details::generic_set_op>
  void dis_dof_update (const SetOp& = SetOp()) const;
 
  // assignments from expressions
  template <class Value>
  typename std::enable_if<
    details::is_rheolef_arithmetic<Value>::value
   ,field_basic<T,M>&
  >::type
  operator=  (const Value& value)
  {
trace_macro("op=(value)...");
    check_macro (get_space().name() != "", "field=constant : uninitialized field in affectation");
    wdof_base::operator= (value);
trace_macro("op=(value) done");
    return *this;
  }
  // TODO: merge with FieldRdof/Lazy && piecewise_polynomial
  template <class Expr>
  typename std::enable_if<
         details::is_field_expr_affine_homogeneous<Expr>::value
    && ! details::has_field_rdof_interface        <Expr>::value
    && ! details::is_field_expr_v2_constant       <Expr>::value
    && ! details::is_field                        <Expr>::value
   ,field_basic<T, M>&
  >::type
  operator=  (const Expr& expr)
  {
trace_macro("op=(expr)...");
    space_type Xh;
    check_macro (expr.have_homogeneous_space (Xh),
      "field = expr; expr should have homogeneous space. HINT: use field = interpolate(Xh, expr)");
    if (get_space().name() != Xh.name()) resize (Xh, 0);
    details::assign_with_operator (begin_dof(), end_dof(), expr.begin_dof(), details::assign_op());
trace_macro("op=(expr) done");
    return *this;
  }
  // TODO: FieldRdof && piecewise_polynomial only: filter
  template<class FieldRdof>
  typename std::enable_if<
         details::has_field_rdof_interface<FieldRdof>::value
    && ! details::is_field                <FieldRdof>::value
   ,field_basic<T, M>&
  >::type
  operator= (const FieldRdof& rdof)
  {
trace_macro("op=(rdof)...");
    const space_type& Xh = rdof.get_space();
    if (get_space().name() != Xh.name()) resize (Xh, 0);
    wdof_base::operator= (rdof);
trace_macro("op=(rdof) done");
    return *this;
  }
  // TODO: FieldLazy && piecewise_polynomial only: filter
  template<class FieldLazy>
  typename std::enable_if<
         details::has_field_lazy_interface<FieldLazy>::value
    && ! details::has_field_rdof_interface<FieldLazy>::value
   ,field_basic<T, M>&
  >::type
  operator= (const FieldLazy& lazy)
  {
trace_macro("op=(lazy)...");
    const space_type& Xh = lazy.get_space();
    if (get_space().name() != Xh.name()) resize (Xh, 0);
    wdof_base::operator= (lazy);
trace_macro("op=(lazy) done");
    return *this;
  }
 
#ifdef TO_CLEAN
  // TODO: move in field_wdof
  template<class Expr>
  typename std::enable_if<
    details::is_field_lazy<Expr>::value
   ,field_basic<T, M>&
  >::type
  operator+= (const Expr&);
 
  // TODO: move in field_wdof
  template<class Expr>
  typename std::enable_if<
    details::is_field_lazy<Expr>::value
   ,field_basic<T, M>&
  >::type
  operator-= (const Expr&);
 
  template<class FieldLazy, class SetPlusOp>
  typename std::enable_if<
    details::is_field_lazy<FieldLazy>::value
   ,void
  >::type
  convert_from_field_lazy (const FieldLazy& expr, const SetPlusOp& set_plus_op);
#endif // TO_CLEAN
 
// accessors by domains and indexes:
 
  // TODO: why is it necessary to explicit them while its done implicitely in field_wdof_indirect ?
  details::field_wdof_indirect<field_basic<T,M>>       operator[] (const geo_type& dom)
                                                       { return wdof_base::operator[] (dom); }
  details::field_wdof_indirect<field_basic<T,M>>       operator[] (std::string dom_name)
                                                       { return wdof_base::operator[] (dom_name); }
  details::field_rdof_indirect_const<field_basic<T,M>> operator[] (const geo_type& dom) const
                                                       { return rdof_base::operator[] (dom); }
  details::field_rdof_indirect_const<field_basic<T,M>> operator[] (std::string dom_name) const
                                                       { return rdof_base::operator[] (dom_name); }
 
 
  size_type size() const { return _V.size(); }
  details::field_wdof_sliced<field_basic<T,M>>       operator[] (size_type i_comp)
                                                     { return wdof_base::operator[] (i_comp); }
  details::field_wdof_sliced<field_basic<T,M>>       operator() (size_type i_comp, size_type j_comp)
                                                     { return wdof_base::operator() (i_comp, j_comp); }
  details::field_rdof_sliced_const<field_basic<T,M>> operator[] (size_type i_comp) const
                                                     { return rdof_base::operator[] (i_comp); }
  details::field_rdof_sliced_const<field_basic<T,M>> operator() (size_type i_comp, size_type j_comp) const
                                                     { return rdof_base::operator() (i_comp, j_comp); }
 
    // old interface
    template <class Expr>
    void do_integrate_internal (
        const geo_basic<T,M>&         dom,
        const geo_basic<T,M>&         band,
        const band_basic<T,M>&        gh,
        const Expr&                   expr,
        const integrate_option&       qopt,
        bool                          is_on_band);
    template <class Expr>
    void do_integrate (
        const geo_basic<T,M>&         dom,
        const Expr&                   expr,
        const integrate_option&       iopt);
    template <class Expr>
    void do_integrate (
        const band_basic<T,M>&        gh,
        const Expr&                   expr,
        const integrate_option&       iopt);
 
protected:
    void dis_dof_indexes_requires_update() const;
    void dis_dof_assembly_requires_update() const;
 
// data:
            space_type   _V;
    mutable vec<T,M>     _u;
    mutable vec<T,M>     _b;
    mutable bool         _dis_dof_indexes_requires_update;
    mutable bool         _dis_dof_assembly_requires_update;
// [verbatim_field_basic_cont]
};
 
namespace details {
 
// concepts:
template<class T, class M>
struct is_field<field_basic<T,M>> : std::true_type {};
 
// field class is not an ordinary function/functor : for compose(field,args...) filtering
template <typename T, typename M>  struct function_traits <field_basic<T,M> > {};
 
template<class T, class M> 
struct field_traits<field_basic<T,M>> {
  using size_type     = std::size_t;
  using scalar_type   = T;
  using memory_type   = M;
};
 
} // namespace details
 
// [verbatim_field_basic_cont]
template <class T, class M>
idiststream& operator >> (odiststream& ips, field_basic<T,M>& u);
 
template <class T, class M>
odiststream& operator << (odiststream& ops, const field_basic<T,M>& uh);
 
// [verbatim_field]
typedef field_basic<Float> field;
// [verbatim_field]
 
typedef field_basic<Float,sequential> field_sequential;
 
// =================================================================================
// field::iterator
// =================================================================================
// TODO: use field_wdof_sliced iterators
template <class T, class M>
class field_basic<T,M>::iterator {
protected:
  typedef typename vec<T,M>::iterator                          data_t; // random acess
  typedef typename disarray<space_pair_type,M>::const_iterator iter_t; // forward access
public:
 
// typedefs:
 
  typedef std::forward_iterator_tag         iterator_category; // TODO: not fully random yet
  typedef typename vec<T,M>::size_type      size_type;
  typedef T                                 value_type;
  typedef T&                                reference;
  typedef T*                                pointer;
  typedef std::ptrdiff_t                    difference_type;
 
// allocators:
 
  iterator ()
    : _blk_dis_iub_iter(),
      _blk_dis_iub_incr(1),
      _u(),
      _b(),
      _first_iu(0),
      _first_ib(0)
    {}
  iterator (iter_t blk_dis_iub_iter, data_t u, data_t b, size_type first_iu, size_type first_ib)
    : _blk_dis_iub_iter(blk_dis_iub_iter),
      _blk_dis_iub_incr(1),
      _u(u),
      _b(b),
      _first_iu(first_iu),
      _first_ib(first_ib)
    {}
  iterator (iter_t blk_dis_iub_iter, size_type blk_dis_iub_incr, data_t u, data_t b, size_type first_iu, size_type first_ib)
    : _blk_dis_iub_iter(blk_dis_iub_iter),
      _blk_dis_iub_incr(blk_dis_iub_incr),
      _u(u),
      _b(b),
      _first_iu(first_iu),
      _first_ib(first_ib)
    {}
 
  void set_increment (size_type incr) { _blk_dis_iub_incr = incr; }
 
// accessors & modifiers:
 
  reference operator* () const { 
      bool      blk     = (*_blk_dis_iub_iter).is_blocked();
      size_type dis_iub = (*_blk_dis_iub_iter).dis_iub();
      size_type     iub = (!blk) ? dis_iub - _first_iu : dis_iub - _first_ib;
      return (!blk) ? _u[iub] : _b[iub];
  }
  iterator& operator++ () { _blk_dis_iub_iter += _blk_dis_iub_incr; return *this; }
  iterator  operator++ (int) { iterator tmp = *this; operator++(); return tmp; }
  iterator& operator+= (difference_type n) { _blk_dis_iub_iter += n*_blk_dis_iub_incr; return *this; }
  iterator  operator+  (difference_type n) const { iterator tmp = *this; return tmp += n; }
  reference operator[] (size_type n) const { return *(*this + n); }
 
// comparators:
 
  bool operator== (const iterator& j) const { return _blk_dis_iub_iter == j._blk_dis_iub_iter; }
  bool operator!= (const iterator& j) const { return ! operator== (j); }
//protected:
  friend class const_iterator;
// data:
  iter_t    _blk_dis_iub_iter;
  size_type _blk_dis_iub_incr;
  data_t _u;
  data_t _b;
  size_type _first_iu, _first_ib;
};
template <class T, class M>
inline
typename field_basic<T,M>::iterator
field_basic<T,M>::begin_dof ()
{
  dis_dof_indexes_requires_update();
  return iterator (_V.data()._idof2blk_dis_iub.begin(), _u.begin(), _b.begin(), _u.ownership().first_index(), _b.ownership().first_index());
}
template <class T, class M>
inline
typename field_basic<T,M>::iterator
field_basic<T,M>::end_dof ()
{
  dis_dof_indexes_requires_update();
  return iterator (_V.data()._idof2blk_dis_iub.end(),   _u.begin(), _b.begin(), _u.ownership().first_index(), _b.ownership().first_index());
}
// =================================================================================
// field::const_iterator
// =================================================================================
template <class T, class M>
class field_basic<T,M>::const_iterator {
protected:
  typedef typename vec<T,M>::const_iterator                    data_t;
  typedef typename disarray<space_pair_type,M>::const_iterator iter_t;
public:
 
// typedefs:
 
  typedef std::forward_iterator_tag         iterator_category;
  typedef typename vec<T,M>::size_type      size_type;
  typedef T                                 value_type;
  typedef const T&                          reference;
  typedef const T*                          pointer;
  typedef std::ptrdiff_t                    difference_type;
 
// allocators:
 
  const_iterator ()
    : _blk_dis_iub_iter(),
      _blk_dis_iub_incr(1),
      _u(),
      _b(),
      _first_iu(0),
      _first_ib(0)
    {}
  const_iterator (iter_t blk_dis_iub_iter, data_t u, data_t b, size_type first_iu, size_type first_ib)
    : _blk_dis_iub_iter(blk_dis_iub_iter),
      _blk_dis_iub_incr(1),
      _u(u), _b(b),
      _first_iu(first_iu),
      _first_ib(first_ib)
    {}
  const_iterator (iterator i)
    : _blk_dis_iub_iter(i._blk_dis_iub_iter),
      _blk_dis_iub_incr(i._blk_dis_iub_incr),
      _u(i._u),
      _b(i._b),
      _first_iu(i._first_iu),
      _first_ib(i._first_ib)
    {}
 
  void set_increment (size_type incr) { _blk_dis_iub_incr = incr; }
 
// accessors & modifiers:
 
  reference operator* () const { 
      bool      blk     = (*_blk_dis_iub_iter).is_blocked();
      size_type dis_iub = (*_blk_dis_iub_iter).dis_iub();
      size_type     iub = (!blk) ? dis_iub - _first_iu : dis_iub - _first_ib;
      return (!blk) ? _u[iub] : _b[iub];
  }
  const_iterator& operator++ () { _blk_dis_iub_iter += _blk_dis_iub_incr; return *this; }
  const_iterator  operator++ (int) { const_iterator tmp = *this; operator++(); return tmp; }
  const_iterator& operator+= (difference_type n) { _blk_dis_iub_iter += n*_blk_dis_iub_incr; return *this; }
  const_iterator  operator+  (difference_type n) const { const_iterator tmp = *this; return tmp += n; }
  reference operator[] (size_type n) const { return *(*this + n); }
 
// comparators:
 
  bool operator== (const const_iterator& j) const { return _blk_dis_iub_iter == j._blk_dis_iub_iter; }
  bool operator!= (const const_iterator& j) const { return ! operator== (j); }
//protected:
// data:
  iter_t    _blk_dis_iub_iter;
  size_type _blk_dis_iub_incr;
  data_t _u;
  data_t _b;
  size_type _first_iu, _first_ib;
};
template <class T, class M>
inline
typename field_basic<T,M>::const_iterator
field_basic<T,M>::begin_dof () const
{
  dis_dof_indexes_requires_update();
  return const_iterator (_V.data()._idof2blk_dis_iub.begin(), _u.begin(), _b.begin(), _u.ownership().first_index(), _b.ownership().first_index());
}
template <class T, class M>
inline
typename field_basic<T,M>::const_iterator
field_basic<T,M>::end_dof () const
{
  dis_dof_indexes_requires_update();
  return const_iterator (_V.data()._idof2blk_dis_iub.end(),   _u.begin(), _b.begin(), _u.ownership().first_index(), _b.ownership().first_index());
}
// =================================================================================
// field: inlined
// =================================================================================
template <class T, class M>
inline
field_basic<T,M>::field_basic ()
 : _V(),
   _u(),
   _b(),
   _dis_dof_indexes_requires_update(true),
   _dis_dof_assembly_requires_update(false)
{
}
// TODO: merge Expr with Fieldlazy
template<class T, class M>
template<class Expr, class Sfinae>
inline
field_basic<T,M>::field_basic (const Expr& expr)
 : _V (),
   _u (),
   _b (),
   _dis_dof_indexes_requires_update(true),
   _dis_dof_assembly_requires_update(false)
{
    operator= (expr);
}
template <class T, class M>
void
field_basic<T,M>::dis_dof_indexes_requires_update() const
{
  _dis_dof_indexes_requires_update = true;
}
template <class T, class M>
void
field_basic<T,M>::dis_dof_assembly_requires_update() const
{
  _dis_dof_assembly_requires_update = true;
}
template <class T, class M>
inline
T& 
field_basic<T,M>::dof (size_type idof)
{
  dis_dof_indexes_requires_update();
  bool      blk     = _V.is_blocked (idof);
  size_type dis_iub = _V.dis_iub    (idof);
  if (!blk) {
    size_type iub = dis_iub - _u.ownership().first_index();
    return _u [iub]; 
  } else {
    size_type iub = dis_iub - _b.ownership().first_index();
    return _b [iub];
  }
}
template <class T, class M>
inline
const T&
field_basic<T,M>::dof (size_type idof) const
{
  bool      blk     = _V.is_blocked (idof);
  size_type dis_iub = _V.dis_iub    (idof);
  return (!blk) ? _u.dis_at(dis_iub) : _b.dis_at(dis_iub);
}
template <class T, class M>
template <class SetOp>
void
field_basic<T,M>::dis_dof_update (const SetOp& set_op) const
{
#ifdef _RHEOLEF_HAVE_MPI
  std::size_t nproc = ownership().comm().size();
  if (is_distributed<M>::value && nproc > 1) {
    std::size_t do_assembly = mpi::all_reduce (ownership().comm(), size_t(_dis_dof_assembly_requires_update), std::plus<std::size_t>());
    std::size_t do_indexes  = mpi::all_reduce (ownership().comm(), size_t(_dis_dof_indexes_requires_update),  std::plus<std::size_t>());
    if (do_assembly) {
      _u.dis_entry_assembly (set_op);
      _b.dis_entry_assembly (set_op);
    }
    if (do_indexes) {
      _u.set_dis_indexes (_V.ext_iu_set());
      _b.set_dis_indexes (_V.ext_ib_set());
    }
  }
#endif // _RHEOLEF_HAVE_MPI
  _dis_dof_indexes_requires_update = false;
  _dis_dof_assembly_requires_update = false;
}
template <class T, class M>
inline
T
field_basic<T,M>::min () const
{
  T val = std::numeric_limits<T>::max();
  for (const_iterator iter = begin_dof(), last = end_dof(); iter != last; iter++) {
    val = std::min(val, *iter);
  }
#ifdef _RHEOLEF_HAVE_MPI
  if (is_distributed<M>::value) {
    val = mpi::all_reduce (comm(), val, mpi::minimum<T>());
  }
#endif // _RHEOLEF_HAVE_MPI
  return val;
}
template <class T, class M>
inline
T
field_basic<T,M>::max () const
{
  T val = std::numeric_limits<T>::min();
  for (const_iterator iter = begin_dof(), last = end_dof(); iter != last; iter++) {
    val = std::max(val, *iter);
  }
#ifdef _RHEOLEF_HAVE_MPI
  if (is_distributed<M>::value) {
    val = mpi::all_reduce (comm(), val, mpi::maximum<T>());
  }
#endif // _RHEOLEF_HAVE_MPI
  return val;
}
template <class T, class M>
inline
T
field_basic<T,M>::min_abs () const
{
  T val = std::numeric_limits<T>::max();
  for (const_iterator iter = begin_dof(), last = end_dof(); iter != last; iter++) {
    val = std::min(val, abs(*iter));
  }
#ifdef _RHEOLEF_HAVE_MPI
  if (is_distributed<M>::value) {
    val = mpi::all_reduce (comm(), val, mpi::minimum<T>());
  }
#endif // _RHEOLEF_HAVE_MPI
  return val;
}
template <class T, class M>
inline
T
field_basic<T,M>::max_abs () const
{
  T val = 0;
  for (const_iterator iter = begin_dof(), last = end_dof(); iter != last; iter++) {
    val = std::max(val, abs(*iter));
  }
#ifdef _RHEOLEF_HAVE_MPI
  if (is_distributed<M>::value) {
    val = mpi::all_reduce (comm(), val, mpi::maximum<T>());
  }
#endif // _RHEOLEF_HAVE_MPI
  return val;
}
template <class T, class M>
inline
idiststream&
operator >> (idiststream& ips, field_basic<T,M>& uh)
{
  return uh.get (ips);
}
template <class T, class M>
inline
odiststream&
operator << (odiststream& ops, const field_basic<T,M>& uh)
{
  return uh.put (ops);
}
// ----------------------------------
// special interpolate with field arg
// ----------------------------------
// re-interpolation of fields and linear field expressions
//      for change of mesh, of approx, ect
// not truly lazy, but here for the completeness of the interpolate interface
template<class T, class M>
inline
field_basic<T,M>
lazy_interpolate (const space_basic<T,M>& X2h, const field_basic<T,M>& u1h)
{
  return interpolate (X2h, u1h);
}
// interpolate e.g. uh[0] or uh["boundary"], i.e. a FieldRdof
template <class T, class M, class FieldRdof>
inline
typename std::enable_if<
       details::has_field_rdof_interface<FieldRdof>::value
  && ! details::is_field<FieldRdof>::value
 ,field_basic<T,M>
>::type
lazy_interpolate (const space_basic<T,M>& Xh, const FieldRdof& uh)
{
  return interpolate (Xh, field_basic<T,M>(uh));
}
// --------------------------
// operator=
// --------------------------
template <class T, class M>
inline
field_basic<T,M>&
field_basic<T,M>::operator=  (const field_basic<T,M>& x)
{
  _V.operator= (x._V);
  _u.operator= (x._u);
  _b.operator= (x._b);
  _dis_dof_indexes_requires_update =  x._dis_dof_indexes_requires_update;
  _dis_dof_assembly_requires_update = x._dis_dof_assembly_requires_update;
  return *this;
}
#ifdef TO_CLEAN
template<class T, class M>
template <class Value>
inline
typename std::enable_if<
  details::is_rheolef_arithmetic<Value>::value
 ,field_basic<T,M>&
>::type
field_basic<T,M>::operator=  (const Value& value)
{
  check_macro (name() != "", "field=constant : uninitialized field in affectation");
  std::fill (begin_dof(), end_dof(), value);
  return *this;
}
template<class T, class M>
template<class FieldLazy>
typename std::enable_if<
  details::is_field_lazy<FieldLazy>::value
 ,field_basic<T, M>&
>::type
field_basic<T,M>::operator+= (const FieldLazy& expr)
{
  const space_basic<T,M>& Xh = expr.get_space();
  convert_from_field_lazy (expr, details::generic_set_plus_op());
  return *this;
}
template<class T, class M>
template<class FieldLazy>
typename std::enable_if<
  details::is_field_lazy<FieldLazy>::value
 ,field_basic<T, M>&
>::type
field_basic<T,M>::operator-= (const FieldLazy& expr)
{
  const space_basic<T,M>& Xh = expr.get_space();
  convert_from_field_lazy (expr, details::generic_set_minus_op());
  return *this;
}
#endif // TO_CLEAN
// --------------------------------------
// io
// --------------------------------------
// TODO: one output for all rdof or lazy & one input for all wdof
template<class FieldLazy>
typename std::enable_if<
  details::has_field_lazy_interface<FieldLazy>::value
  && ! details::is_field<FieldLazy>::value
  // details::is_field_lazy<FieldLazy>::value
 ,odiststream&
>::type
operator<< (odiststream& out, const FieldLazy& expr)
{
  using T = typename FieldLazy::scalar_type;
  using M = typename FieldLazy::memory_type;
  field_basic<T,M> uh = expr;
  return out << uh;
}
#ifdef TODO
// TODO: as: field vh; din >> vh ; uh[i_comp] = vh;
template<class FieldWdof>
inline
idiststream&
operator>> (odiststream& ids, details::field_wdof_sliced<FieldWdof>& u);
#endif // TODO
// --------------------------------------
// interaction with field_indirect
// --------------------------------------
// TODO: use a template FieldRdof and a dof-based convert function
// => interaction=convert from field_rdof
#ifdef TO_CLEAN
template<class T, class M>
inline
details::field_wdof_indirect<field_basic<T,M>>
field_basic<T,M>::operator[] (const geo_type& dom)
{
  dis_dof_indexes_requires_update();
  return details::field_wdof_indirect<field_basic<T,M>> (*this, dom);
}
template<class T, class M>
inline
details::field_wdof_indirect<field_basic<T,M>>
field_basic<T,M>::operator[] (std::string dom_name)
{
  dis_dof_indexes_requires_update();
  return operator[] (get_space().get_geo().operator[] (dom_name));
}
template<class T, class M>
inline
details::field_rdof_indirect_const<field_basic<T,M>>
field_basic<T,M>::operator[] (const geo_type& dom) const
{
  return details::field_rdof_indirect_const<field_basic<T,M>> (*this, dom);
}
template<class T, class M>
inline
details::field_rdof_indirect_const<field_basic<T,M>>
field_basic<T,M>::operator[] (std::string dom_name) const
{
  return operator[] (get_space().get_geo().operator[] (dom_name));
}
#endif // TO_CLEAN
// --------------------------------------
// interaction with field_wdof_sliced
// --------------------------------------
#ifdef TO_CLEAN
template<class T, class M>
inline
details::field_rdof_sliced_const<field_basic<T,M>>
field_basic<T,M>::operator[] (size_type i_comp) const
{
  return details::field_rdof_sliced_const<field_basic<T,M>> (*this, i_comp);
}
template<class T, class M>
inline
details::field_wdof_sliced<field_basic<T,M>>
field_basic<T,M>::operator[] (size_type i_comp)
{
  dis_dof_indexes_requires_update();
  return details::field_wdof_sliced<field_basic<T,M>> (*this, i_comp);
}
// --------------------------------------
// convert from field_lazy
// --------------------------------------
template<class T, class M>
template<class FieldLazy, class SetPlusOp>
inline
typename std::enable_if<
  details::is_field_lazy<FieldLazy>::value
 ,void
>::type
field_basic<T,M>::convert_from_field_lazy (const FieldLazy& expr, const SetPlusOp& set_plus_op)
{
  details::convert_lazy2wdof (expr, set_plus_op, *this);
}
#endif // TO_CLEAN
 
}// namespace rheolef
# endif /* _RHEOLEF_FIELD_H */