# ifndef _RHEOLEF_FIELD_LAZY_TERMINAL_H
# define _RHEOLEF_FIELD_LAZY_TERMINAL_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_lazy = un-assembled field
// as returned by lazy_integrate, lazy_interpolate or encapsulated field_basic
// AUTHOR: Pierre.Saramito@imag.fr
// DATE:   6 april 1920

// SUMMARY:
// 1. concept & base class : see field_lazy.h
// 2. terminal
//    2.1. field
//    2.2. integrate
//    2.3. integrate on a band
//    2.4. interpolate
// see also "field_lazy_node.h"
//
#include "rheolef/field_lazy.h"

namespace rheolef {

#ifdef TO_CLEAN
// -------------------------------------------------------------------
// 1. concept & base class
// -------------------------------------------------------------------
namespace details {
// Define a trait type for detecting field expression valid arguments
//     template<class T> struct is_field_lazy: std::false_type {};
// => should be defined in field.h for compilation reason,
//    otherwise Sfinae is always failed in field.h

// Define a base class
template<class Derived>
class field_lazy_base {
public:

// no common methods yet

protected:
  Derived&       derived()       { return *static_cast<      Derived*>(this); }
  const Derived& derived() const { return *static_cast<const Derived*>(this); }
};

} // namespace details
#endif // TO_CLEAN
// -------------------------------------------------------------------
// 2. terminal
// -------------------------------------------------------------------
// 2.1. field
// -------------------------------------------------------------------
// field_lazy_terminal_field encapsulates the field_basic class
// It then appears as an unassembled field, for combining with 
// others unassembled fields in expressions
//
namespace details {

template<class T, class M>
class field_lazy_terminal_field: public field_lazy_base<field_lazy_terminal_field<T,M>> {
public :
// definitions:

  using base        = field_lazy_base<field_lazy_terminal_field<T,M>>;
  using size_type   = geo_element::size_type;
  using memory_type = M;
  using scalar_type = T;
  using field_type  = field_basic<T,M>;
  using space_type  = typename field_type::space_type;
  using geo_type    = typename field_type::geo_type;
  using float_type  = typename float_traits<scalar_type>::type;
  using band_type   = band_basic<float_type,memory_type>;
  using vector_element_type = Eigen::Matrix<scalar_type,Eigen::Dynamic,1>;

// allocator:

  explicit field_lazy_terminal_field (const field_type& uh) : base(), _uh(uh) {}

// accessors:

  const geo_type&   get_geo()   const { return _uh.get_geo(); }
  const space_type& get_space() const { return _uh.get_space(); }
  bool  is_on_band()            const { return false; }
  band_type get_band()          const { return band_type(); }

  void initialize (const geo_type& omega_K);

  void evaluate (
        const geo_type&            omega_K,
        const geo_element&         K,
              vector_element_type& uk) const;
// data:
protected:
  field_type  _uh;
};
// concept;
template<class T, class M>
struct is_field_lazy <field_lazy_terminal_field<T,M> > : std::true_type {};

// inlined;
template<class T, class M>
void
field_lazy_terminal_field<T,M>::initialize (const geo_type& omega_K)
{
  _uh.dis_dof_update();
}
template<class T, class M>
void
field_lazy_terminal_field<T,M>::evaluate (
  const geo_type&            omega_K,
  const geo_element&         K,
        vector_element_type& uk) const
{
  std::vector<size_type> dis_idof;
  _uh.get_space().get_constitution().assembly_dis_idof (omega_K, K, dis_idof);
  size_type loc_ndof = dis_idof.size();
  uk.resize (loc_ndof);
  for (size_type loc_jdof = 0; loc_jdof < loc_ndof; ++loc_jdof) {
    size_type dis_jdof = dis_idof[loc_jdof];
    uk [loc_jdof] = _uh.dis_dof (dis_jdof);
  }
}

}// namespace details
// -------------------------------------------------------------------
// 2.2. integrate
// -------------------------------------------------------------------
// field_lazy_terminal_integrate is returned by the integrate(domain) function
// It contains the domain of integration and integration options (quadrature)
// The sumitted elements K during evaluation belongs to a superset omega_K of "domain"
// => the boolean is_on_domain[K_ie] behaves as an indicator function for "domain"
//
namespace details {

template<class Expr>
class field_lazy_terminal_integrate_rep {
public :
// definitions:

  using size_type   = geo_element::size_type;
  using memory_type = typename Expr::memory_type;
  using scalar_type = typename Expr::scalar_type;
  using space_type  = typename Expr::space_type;
  using geo_type    = typename Expr::geo_type;
  using float_type  = typename float_traits<scalar_type>::type;
  using band_type   = band_basic<float_type,memory_type>;
  using vector_element_type = Eigen::Matrix<scalar_type,Eigen::Dynamic,1>;

// allocators:

  field_lazy_terminal_integrate_rep (
    const Expr&             expr,
    const integrate_option& iopt);

  field_lazy_terminal_integrate_rep (
    const geo_type&         domain,
    const Expr&             expr,
    const integrate_option& iopt);

  field_lazy_terminal_integrate_rep (
    std::string             domname,
    const Expr&             expr,
    const integrate_option& iopt);

// accessors:

  const geo_type&   get_geo()   const { return _domain; }
  const space_type& get_space() const { return _expr.get_vf_space(); }
  bool  is_on_band()            const { return false; }
  band_type get_band()          const { return band_type(); }

  void initialize (const geo_type& omega_K) const;

  void evaluate (
      const geo_type&            omega_K,
      const geo_element&         K,
            vector_element_type& uk) const;
// data:
protected:
  geo_type                     _domain;
  mutable Expr                 _expr;
  mutable integrate_option     _iopt;
  mutable disarray<int,memory_type>  _is_on_domain;
  mutable geo_type             _prev_omega_K;
  mutable size_type            _prev_K_dis_ie;
  mutable vector_element_type  _prev_uk;
// internal:
  const geo_type& get_default_geo () const;
};
// inlined;
template<class Expr>
field_lazy_terminal_integrate_rep<Expr>::field_lazy_terminal_integrate_rep (
  const geo_type&         domain,
  const Expr&             expr,
  const integrate_option& iopt)
: _domain(domain),
  _expr(expr),
  _iopt(iopt),
  _is_on_domain(),
  _prev_omega_K(),
  _prev_K_dis_ie(std::numeric_limits<size_type>::max()), 
  _prev_uk()
{
}
// missing domain:
template<class Expr>
field_lazy_terminal_integrate_rep<Expr>::field_lazy_terminal_integrate_rep (
  const Expr&             expr,
  const integrate_option& iopt)
: _domain(),
  _expr(expr),
  _iopt(iopt),
  _is_on_domain(),
  _prev_omega_K(),
  _prev_K_dis_ie(std::numeric_limits<size_type>::max()), 
  _prev_uk()
{
  _domain = get_default_geo();
}
// missing domain:
template<class Expr>
field_lazy_terminal_integrate_rep<Expr>::field_lazy_terminal_integrate_rep (
  std::string             domname,
  const Expr&             expr,
  const integrate_option& iopt)
: _domain(),
  _expr(expr),
  _iopt(iopt),
  _is_on_domain(),
  _prev_omega_K(),
  _prev_K_dis_ie(std::numeric_limits<size_type>::max()), 
  _prev_uk()
{
  _domain = get_default_geo() [domname];
}
template<class Expr>
const typename field_lazy_terminal_integrate_rep<Expr>::geo_type&
field_lazy_terminal_integrate_rep<Expr>::get_default_geo () const
{
  // TODO: improve the automatic determination of the domain ?
  return get_space().get_constitution().get_geo();
}
template<class Expr>
void
field_lazy_terminal_integrate_rep<Expr>::initialize (const geo_type& omega_K) const
{
  using float_type = typename float_traits<scalar_type>::type;
  _iopt._is_on_interface          = false;
  _iopt._is_inside_on_local_sides = false;
  _expr.initialize (_domain, _iopt);
  _prev_omega_K  = omega_K;
  _prev_K_dis_ie = std::numeric_limits<size_type>::max();
  // --------------------------------
  // initialize the domain indicator:
  // --------------------------------
  // TODO: sub-domain
  // _domain could be a subdomain of omega_K:
  // in that case evaluate(omega_K,K) should return a zero uk vector
  // when K do not belongs to _domain
  size_type K_map_d = omega_K.dimension();
  _is_on_domain.resize (omega_K.geo_element_ownership(K_map_d));
  std::fill (_is_on_domain.begin(), _is_on_domain.end(), false);
  if (_domain.map_dimension() == omega_K.map_dimension()) {
    if (_domain == omega_K || _domain.variant() != geo_abstract_base_rep<float_type>::geo_domain) {
#ifdef TODO
      // => _domain has the same element numbering as omega_K
      //    _domain.variant() is a geo or a geo_domain_indirect
      trace_macro ("lazy_int(init): domain="<<_domain.name()<<": same numbering as "<<omega_K.name()<<"...");
      size_type first_dis_ie = omega_K.geo_element_ownership (omega_K.map_dimension()).first_index();
      trace_macro ("lazy_int(init): first_dis_ie="<<first_dis_ie);
      trace_macro ("lazy_int(init): _is_on_domain.size="<<_is_on_domain.size());
      for (size_type dom_ie = 0, dom_ne = _domain.size(); dom_ie < dom_ne; ++dom_ie) {
        trace_macro ("lazy_int(init): dom_ie="<<dom_ie);
        const geo_element& K = _domain [dom_ie];
        size_type dis_ie = K.dis_ie();
        trace_macro ("lazy_int(init): K.dis_ie="<<dis_ie);
        check_macro (first_dis_ie <= dis_ie, "invalid index"); // PARANO
        size_type ie = dis_ie - first_dis_ie;
        trace_macro ("lazy_int(init): ie="<<ie);
        check_macro (ie <= _is_on_domain.size(), "invalid index"); // PARANO
        _is_on_domain [ie] = true;
      }
#endif // TODO
      trace_macro ("lazy_int(init): domain="<<_domain.name()<<": same numbering as "<<omega_K.name()<<" done");
    } else {
      error_macro ("geo_domain="<<_domain.name()<<" subset of "<<omega_K.name()<<": not yet");
    }
  } else if (_domain.map_dimension()+1 == omega_K.map_dimension()) {
    check_macro (_domain.get_background_geo() == omega_K, "unexpected domain \""
        <<_domain.name()<<"\" as subset of \""<<omega_K.name()<<"\"");
    if (_domain.variant() != geo_abstract_base_rep<float_type>::geo_domain) {
      trace_macro ("lazy_int(init): domain(bdry or sides)="<<_domain.name()<<" subset of "<<omega_K.name()<<"...");
      for (size_type dom_ie = 0, dom_ne = _domain.size(); dom_ie < dom_ne; ++dom_ie) {
        const geo_element& S = _domain [dom_ie];
        size_type dis_is = S.dis_ie();
        // _domain="sides" or "boundary" or other d-1 sides domain
        // for a side S, will access to its neighbours K0 & K1
        size_type K_dis_ie = S.master(0); // TODO: domain = "interface" : orient could be -1 and then choose S.master(1) !!
        check_macro (K_dis_ie != std::numeric_limits<size_type>::max(),
          "unexpected isolated side S="<<S.name()<<S.dis_ie());
        const geo_element& K = omega_K.dis_get_geo_element (K_map_d, K_dis_ie);
        size_type dis_ie = K.dis_ie(); // not necessarily on the same proc as S
        _is_on_domain.dis_entry (dis_ie) = true;
      }
      _is_on_domain.dis_entry_assembly();
      trace_macro ("lazy_int(init): domain(bdry or sides)="<<_domain.name()<<" subset of "<<omega_K.name()<<" done");
    } else {
      error_macro ("geo_domain(bdry or sides)="<<_domain.name()<<" subset of "<<omega_K.name()<<": not yet");
    }
  } else {
    error_macro ("unsupported domain \"" << _domain.name()
       << "\" with map_dimension=" << _domain.map_dimension()
       << " when integrated in geometry \"" << omega_K.name()
       << "\" with map_dimension=" << omega_K.map_dimension());
  }
}
/*
 compute local idofs on S from the basis in K ?
 all dofs in K are numberd from 0 to ndof(K)-1
 for each S :
   is_in_S[0:ndof(K)-1] = false
   for each subgeo_dim=0 to S.dim :
     for each subgeo T of S with T.dim=subgeo_dim
       get T.first_idof and T.last_idof from basis infos
       for idof=T.first_idof to T.last_idof-1
         is_in_S[idof] = true
   ndof_S = 0
   for idof = 0 to ndof(K)-1
     if is_in_S[idof] then ndof_S++
   idof_S2idof.resize (ndof_S)
   idof_S = 0
   for idof = 0 to ndof(K)-1
     if is_in_S[idof] then
       idof_S2idof[idof_S++] = idof

   finally:
   for idof_S = 0 to ndof_S-1
     uk[idof_S2idof[idof_S] = us[idof_S] 
  The computation of idof_S2idof could be done on the basis(hat_K) 
  one time for all for all sides S at initialization or on the fly
  at the first request
  - first  step : we do it here, for all K
  - second step : we move it on the basis class in an internal mutable data
  For a space product eg Yh=Xh*Mh we have to concactenate the idof_S2idof arrays
  => will be a member function of space class ?
*/
template <class T>
void
compute_idof_S2idof_K (
  const basis_basic<T>&                                     b,
  const reference_element&                                  hat_K,
        std::array<std::vector<size_t>,
                   reference_element::max_side_by_variant>& idof_S2idof_K,
        std::array<size_t,
                   reference_element::max_side_by_variant>& ndof_S,
        size_t&                                             ndof_K)
{
  trace_macro ("compute_idof_S2idof_K: hat_K="<<hat_K.name()<<"...");
  using size_type = geo_element::size_type;
  check_macro (hat_K.dimension() > 0, "idof_S2idof_K: invalid K.dimension=0");
  size_type sid_dim = hat_K.dimension()-1;
  trace_macro ("compute_idof_S2idof_K: sid_dim="<<sid_dim);
  ndof_K = b.ndof (hat_K);
  // ------------------------------
  // 1) is_in_S [loc_isid] [idof_K]
  // ------------------------------
  std::array<std::vector<bool>, reference_element::max_side_by_variant> is_in_S;
  std::array<std::map<size_type,size_type>, reference_element::max_side_by_variant> idof_S2idof_K_map;
  for (size_type loc_isid = 0, loc_nsid = hat_K.n_subgeo(sid_dim); loc_isid < loc_nsid; ++loc_isid) {
    is_in_S[loc_isid].resize (ndof_K);
    std::fill (is_in_S[loc_isid].begin(), is_in_S[loc_isid].end(), false);
    reference_element hat_S = hat_K.subgeo (sid_dim, loc_isid);
    trace_macro ("compute_idof_S2idof_K: loc_isid="<<loc_isid<<", hat_S="<<hat_S.name());
    size_type idof_S = 0;
    for (size_type sub_sid_dim = 0; sub_sid_dim <= sid_dim; ++sub_sid_dim) {
      size_type first_idof_K = b.first_idof_by_dimension (hat_K, sub_sid_dim);
      trace_macro ("compute_idof_S2idof_K: first_idof_K(hat_K="<<hat_K.name()<<",sub_sid_dim="<<sub_sid_dim<<") = " <<first_idof_K);
      for (size_type loc_isub_sid = 0; loc_isub_sid < hat_S.n_subgeo(sub_sid_dim); ++loc_isub_sid) {
        reference_element hat_T = hat_S.subgeo (sub_sid_dim, loc_isub_sid);
        size_type ndof_K_on_T = b.ndof_on_subgeo (hat_K.dimension(), hat_T.variant());
        size_type loc_T_idx = 0;
        switch (sub_sid_dim) {
          case 0: {
                    loc_T_idx = hat_K.subgeo_local_vertex (sid_dim, loc_isid, loc_isub_sid);
                    break;
                  }
          case 1: {
                    if (sid_dim == 1) {
                      loc_T_idx = loc_isid; // the side=edge in 2D itself
                    } else {
                      error_macro("compute_idof_S2idof_K: dof on edge in 3D: not yet");
                    }
                    break;
                  }
          case 2:
          default:{
                    if (sid_dim == 2) {
                      loc_T_idx = loc_isid; // the side=face in 3D itself
                    } else {
                      error_macro("compute_idof_S2idof_K: dof on face in 3D: not yet");
                    }
                    break;
                  }
        }
        trace_macro ("compute_idof_S2idof_K: ndof_K_on_T(hat_K.dim="<<hat_K.dimension()<<",hat_T="<<hat_T.name()<<") = " <<ndof_K_on_T);
        trace_macro ("compute_idof_S2idof_K: T_idx="<<loc_T_idx);
        size_type start_idof_K = first_idof_K + ndof_K_on_T*loc_T_idx;
        size_type  last_idof_K = first_idof_K + ndof_K_on_T*(loc_T_idx+1);
        for (size_type idof_K = start_idof_K; idof_K < last_idof_K; ++idof_K) {
          if (is_in_S [loc_isid] [idof_K]) continue;
          is_in_S [loc_isid] [idof_K] = true;
          idof_S2idof_K_map [loc_isid] [idof_S] = idof_K;
          trace_macro ("idof_S2idof_K_map(isid="<<loc_isid<<",idof_S="<<idof_S<<") = "<<idof_K);
	  ++idof_S;
        }
      }
    }
  }
#ifdef TO_CLEAN
  // ------------------------------
  // 2) count ndof on each S
  // ------------------------------
  for (size_type loc_isid = 0, loc_nsid = hat_K.n_subgeo(sid_dim); loc_isid < loc_nsid; ++loc_isid) {
    ndof_S [loc_isid] = 0;
    for (size_type idof_K = 0; idof_K < ndof_K; ++idof_K) {
      if (is_in_S [loc_isid] [idof_K]) {
        ndof_S [loc_isid]++;
      }
    }
    trace_macro ("compute_idof_S2idof_K: ndof_S[isid="<<loc_isid<<"]="<<ndof_S[loc_isid]);
  }
#endif // TO_CLEAN
  // ------------------------------
  // 3) set idof_S2idof_K 
  // ------------------------------
  for (size_type loc_isid = 0, loc_nsid = hat_K.n_subgeo(sid_dim); loc_isid < loc_nsid; ++loc_isid) {
    ndof_S [loc_isid] = idof_S2idof_K_map [loc_isid].size();
    idof_S2idof_K[loc_isid].resize (ndof_S [loc_isid]);
    for (auto p : idof_S2idof_K_map [loc_isid]) {
      size_type idof_S = p.first;
      size_type idof_K = p.second;
      if (is_in_S [loc_isid] [idof_K]) {
        idof_S2idof_K [loc_isid] [idof_S] = idof_K;
        trace_macro ("idof_K(isid="<<loc_isid<<",idof_S="<<idof_S<<") = "<< idof_K);
        idof_S++;
      }
    }
  }
  trace_macro ("compute_idof_S2idof_K done");
}
template<class Expr>
void
field_lazy_terminal_integrate_rep<Expr>::evaluate (
  const geo_type&            omega_K,
  const geo_element&         K,
        vector_element_type& uk) const
{
  if (_prev_omega_K == omega_K && _prev_K_dis_ie == K.dis_ie()) {
    uk = _prev_uk;
    return;
  }
  if (_domain.map_dimension() == omega_K.map_dimension()) {
    trace_macro ("lazy_int(eval): domain="<<_domain.name()<<": same numbering as "<<omega_K.name()<<"...");
    _expr.evaluate (omega_K, K, uk);
    trace_macro ("lazy_int(eval): domain="<<_domain.name()<<": same numbering as "<<omega_K.name()<<" done");
  } else { // _domain.map_dimension()+1 == omega_K.map_dimension()
    trace_macro ("lazy_int(eval): domain(bdry or sides)="<<_domain.name()<<" subset of "<<omega_K.name()<<"...");
    std::array<std::vector<size_type>, reference_element::max_side_by_variant> idof_S2idof_K;
    std::array<size_type, reference_element::max_side_by_variant>              ndof_S;
    size_type                                                                  ndof_K = 0;
    check_macro (get_space().size() == 0, "lazy_int(eval): "<<get_space().size()<<"-component space: not yet supported");
    compute_idof_S2idof_K (get_space().get_basis(), K, idof_S2idof_K, ndof_S, ndof_K); // TODO: do it in space:: with concat
    uk = vector_element_type::Zero (ndof_K, 1);
    trace_macro ("lazy_int(eval): uk.size="<<uk.size());
    check_macro (K.dimension() > 0, "lazy_int(eval): invalid K.dimension=0");
    size_type sid_dim = K.dimension()-1;
    vector_element_type us;
    for (size_type loc_isid = 0, nsid = K.n_subgeo(sid_dim); loc_isid < nsid; ++loc_isid) {
      size_type dis_isid = K.subgeo_dis_index(sid_dim,loc_isid);
      const geo_element& S = omega_K.dis_get_geo_element (sid_dim, dis_isid);
      side_information_type sid;
      K.get_side_informations (S, sid);
      _expr.evaluate (_domain, S, us);
      trace_macro ("lazy_int(eval): us(loc_isid="<<loc_isid<<").size="<<us.size());
      check_macro (size_type(us.size()) == ndof_S [loc_isid], "invalid us size");
      for (size_type idof_S = 0; idof_S < ndof_S [loc_isid]; ++idof_S) {
        size_type idof_K = idof_S2idof_K[loc_isid][idof_S];
        trace_macro ("lazy_int(eval): uk[idof_K="<<idof_K<<"] += us[idof_S="<<idof_S<<"] = " << us[idof_S]);
        uk[idof_K] += us[idof_S];
      }
    }
    trace_macro ("lazy_int(eval): domain(bdry or sides)="<<_domain.name()<<" subset of "<<omega_K.name()<<" done");
  }
  _prev_uk       = uk; // expensive to compute, so memorize it for common subexpressions
  _prev_omega_K  = omega_K;
  _prev_K_dis_ie = K.dis_ie();
  trace_macro("lazy_int(K="<<K.name()<<K.dis_ie()<<",prev="<<_prev_K_dis_ie<<"): compute");
}

template<class Expr>
class field_lazy_terminal_integrate: public smart_pointer_nocopy<field_lazy_terminal_integrate_rep<Expr>>,
                                     public field_lazy_base     <field_lazy_terminal_integrate<Expr>> {
public :
// definitions:

  using rep    = field_lazy_terminal_integrate_rep<Expr>;
  using base1  = smart_pointer_nocopy<rep>;
  using base2  = field_lazy_base<field_lazy_terminal_integrate<Expr>>;
  using size_type   = typename rep::size_type;
  using memory_type = typename rep::memory_type;
  using scalar_type = typename rep::scalar_type;
  using space_type  = typename rep::space_type;
  using geo_type    = typename rep::geo_type;
  using band_type   = typename rep::band_type;
  using vector_element_type = typename rep::vector_element_type;

// allocator:

  field_lazy_terminal_integrate (
    const geo_type&         domain,
    const Expr&             expr,
    const integrate_option& iopt)
   : base1(new_macro(rep(domain,expr,iopt))),
     base2()
   {}

  field_lazy_terminal_integrate (
    const Expr&             expr,
    const integrate_option& iopt)
   : base1(new_macro(rep(expr,iopt))),
     base2()
   {}

  field_lazy_terminal_integrate (
    std::string             domname,
    const Expr&             expr,
    const integrate_option& iopt)
   : base1(new_macro(rep(domname,expr,iopt))),
     base2()
   {}

// accessors:

  const geo_type&   get_geo()   const { return base1::data().get_geo(); }
  const space_type& get_space() const { return base1::data().get_space(); }
  bool  is_on_band()            const { return base1::data().is_on_band(); }
  band_type get_band()          const { return base1::data().get_band(); }

  void initialize (const geo_type& omega_K) const { return base1::data().initialize (omega_K); }

  void evaluate (
      const geo_type&            omega_K,
      const geo_element&         K,
            vector_element_type& uk) const
      { base1::data().evaluate (omega_K, K, uk); }
};
// concept;
template<class Expr> struct is_field_lazy <field_lazy_terminal_integrate<Expr> > : std::true_type {};

}// namespace details

// 2.2.a) general call
template <class Expr>
inline
typename
std::enable_if<
  details::is_field_expr_quadrature_arg<Expr>::value
 ,details::field_lazy_terminal_integrate <Expr>
>::type
//! @brief see the @ref integrate_3 page for the full documentation
lazy_integrate (
  const typename Expr::geo_type& domain, 
  const Expr&                    expr,
  const integrate_option&        iopt = integrate_option())
{
  return details::field_lazy_terminal_integrate<Expr> (domain, expr, iopt);
}
template <class Expr>
inline
typename
std::enable_if<
  details::is_field_expr_v2_variational_arg<Expr>::value
 ,details::field_lazy_terminal_integrate <details::field_expr_quadrature_on_element<Expr>>
>::type
//! @brief see the @ref integrate_3 page for the full documentation
lazy_integrate (
  const geo_basic<typename Expr::scalar_type, typename Expr::memory_type>& domain, 
  const Expr& expr,
  const integrate_option& iopt = integrate_option())
{
  details::field_expr_quadrature_on_element<Expr> expr_quad(expr);
  return lazy_integrate (domain, expr_quad, iopt);
}
// 2.2.b) missing domain
template <class Expr>
inline
typename
std::enable_if<
  details::is_field_expr_quadrature_arg<Expr>::value
 ,details::field_lazy_terminal_integrate <Expr>
>::type
//! @brief see the @ref integrate_3 page for the full documentation
lazy_integrate (
  const Expr& expr,
  const integrate_option& iopt = integrate_option())
{
  return details::field_lazy_terminal_integrate<Expr> (expr, iopt);
}
template <class Expr>
inline
typename
std::enable_if<
  details::is_field_expr_v2_variational_arg<Expr>::value
 ,details::field_lazy_terminal_integrate <details::field_expr_quadrature_on_element<Expr>>
>::type
//! @brief see the @ref integrate_3 page for the full documentation
lazy_integrate (
  const Expr& expr,
  const integrate_option& iopt = integrate_option())
{
  details::field_expr_quadrature_on_element<Expr> expr_quad(expr);
  return lazy_integrate (expr_quad, iopt);
}
// 2.2.c) subdomain by its name
template <class Expr>
inline
typename
std::enable_if<
  details::is_field_expr_quadrature_arg<Expr>::value
 ,details::field_lazy_terminal_integrate <Expr>
>::type
//! @brief see the @ref integrate_3 page for the full documentation
lazy_integrate (
  const std::string& domname, 
  const Expr& expr,
  const integrate_option& iopt = integrate_option())
{
  return details::field_lazy_terminal_integrate<Expr> (domname, expr, iopt);
}
template <class Expr>
inline
typename
std::enable_if<
  details::is_field_expr_v2_variational_arg<Expr>::value
 ,details::field_lazy_terminal_integrate <details::field_expr_quadrature_on_element<Expr>>
>::type
//! @brief see the @ref integrate_3 page for the full documentation
lazy_integrate (
  const std::string& domname, 
  const Expr& expr,
  const integrate_option& iopt = integrate_option())
{
  details::field_expr_quadrature_on_element<Expr> expr_quad(expr);
  return lazy_integrate (domname, expr_quad, iopt);
}
// ----------------------------------------------
// 2.3. integrate on a band
// ----------------------------------------------
namespace details {

template<class Expr>
class field_lazy_terminal_integrate_band_rep {
public :
// definitions:

  using size_type   = geo_element::size_type;
  using memory_type = typename Expr::memory_type;
  using scalar_type = typename Expr::scalar_type;
  using space_type  = typename Expr::space_type;
  using geo_type    = typename Expr::geo_type;
  using float_type  = typename float_traits<scalar_type>::type;
  using band_type   = band_basic<float_type,memory_type>;
  using vector_element_type = Eigen::Matrix<scalar_type,Eigen::Dynamic,1>;

// allocators:

  field_lazy_terminal_integrate_band_rep (
    const band_type&        gh,
    const Expr&             expr,
    const integrate_option& iopt);

// accessors:

  const geo_type&   get_geo()   const { return _gh.level_set(); }
  const space_type& get_space() const { return _expr.get_vf_space(); }
  bool  is_on_band()            const { return true; }
  band_type get_band()          const { return _gh; }

  void initialize (const geo_type& omega_K) const;

  void evaluate (
      const geo_type&            omega_K,
      const geo_element&         K,
            vector_element_type& uk) const;
// data:
protected:
  band_type                    _gh;
  mutable Expr                 _expr;
  mutable integrate_option     _iopt;
  mutable geo_type             _prev_omega_K;
  mutable size_type            _prev_K_dis_ie;
  mutable vector_element_type  _prev_uk;
};
// inlined;
template<class Expr>
field_lazy_terminal_integrate_band_rep<Expr>::field_lazy_terminal_integrate_band_rep (
  const band_type&        gh,
  const Expr&             expr,
  const integrate_option& iopt)
: _gh(gh),
  _expr(expr),
  _iopt(iopt),
  _prev_omega_K(),
  _prev_K_dis_ie(std::numeric_limits<size_type>::max()), 
  _prev_uk()
{
}
template<class Expr>
void
field_lazy_terminal_integrate_band_rep<Expr>::initialize (const geo_type& omega_K) const
{
  const geo_type& band  = _gh.band(); // the 3D intersected tetras
  const geo_type& bgd_omega = get_space().get_constitution().get_background_geo();
  check_macro (omega_K == get_geo(),
        "integrate on band: unexpected integration domain");
  check_macro (band.get_background_geo() == bgd_omega,
        "do_integrate: incompatible integration domain "<<_gh.level_set().name() << " and test function based domain "
        << bgd_omega.name());
  _expr.initialize (_gh, _iopt);
  _prev_omega_K  = omega_K;
  _prev_K_dis_ie = std::numeric_limits<size_type>::max();
}
template<class Expr>
void
field_lazy_terminal_integrate_band_rep<Expr>::evaluate (
  const geo_type&            omega_K,
  const geo_element&         K,
        vector_element_type& uk) const
{
  if (_prev_omega_K == omega_K && _prev_K_dis_ie == K.dis_ie()) {
    uk = _prev_uk;
    return;
  }
  _expr.evaluate (omega_K, K, uk);
  _prev_uk       = uk; // expensive to compute, so memorize it for common subexpressions
  _prev_omega_K  = omega_K;
  _prev_K_dis_ie = K.dis_ie();
}
template<class Expr>
class field_lazy_terminal_integrate_band: public smart_pointer_nocopy<field_lazy_terminal_integrate_band_rep<Expr>>,
                                          public field_lazy_base     <field_lazy_terminal_integrate_band<Expr>> {
public :
// definitions:

  using rep    = field_lazy_terminal_integrate_band_rep<Expr>;
  using base1  = smart_pointer_nocopy<rep>;
  using base2  = field_lazy_base<field_lazy_terminal_integrate_band<Expr>>;
  using size_type   = typename rep::size_type;
  using memory_type = typename rep::memory_type;
  using scalar_type = typename rep::scalar_type;
  using space_type  = typename rep::space_type;
  using geo_type    = typename rep::geo_type;
  using band_type   = typename rep::band_type;
  using vector_element_type = typename rep::vector_element_type;

// allocator:

  field_lazy_terminal_integrate_band (
    const band_type&        gh, 
    const Expr&             expr,
    const integrate_option& iopt)
   : base1(new_macro(rep(gh,expr,iopt))),
     base2()
   {}

// accessors:

  const geo_type&   get_geo()   const { return base1::data().get_geo(); }
  const space_type& get_space() const { return base1::data().get_space(); }
  bool  is_on_band()            const { return base1::data().is_on_band(); }
  band_type get_band()          const { return base1::data().get_band(); }

  void initialize (const geo_type& omega_K) const { return base1::data().initialize (omega_K); }

  void evaluate (
      const geo_type&            omega_K,
      const geo_element&         K,
            vector_element_type& uk) const
      { base1::data().evaluate (omega_K, K, uk); }
};
// concept;
template<class Expr> struct is_field_lazy <field_lazy_terminal_integrate_band<Expr> > : std::true_type {};

}// namespace details

template <class Expr>
typename
std::enable_if<
  details::is_field_expr_quadrature_arg<Expr>::value
 ,details::field_lazy_terminal_integrate_band <Expr>
>::type
//! @brief see the @ref integrate_3 page for the full documentation
lazy_integrate (
  const band_basic<typename float_traits<typename Expr::scalar_type>::type,
                   typename Expr::memory_type>& gh, 
  const Expr& expr,
  const integrate_option& iopt = integrate_option())
{
  return details::field_lazy_terminal_integrate_band<Expr> (gh, expr, iopt);
}
template <class Expr>
typename
std::enable_if<
  details::is_field_expr_v2_variational_arg<Expr>::value
 ,details::field_lazy_terminal_integrate_band <details::field_expr_quadrature_on_element<Expr>>
>::type
//! @brief see the @ref integrate_3 page for the full documentation
lazy_integrate (
  const band_basic<typename float_traits<typename Expr::scalar_type>::type,
                   typename Expr::memory_type>& gh, 
  const Expr& expr,
  const integrate_option& iopt = integrate_option())
{

  details::field_expr_quadrature_on_element<Expr> expr_quad(expr);
  return lazy_integrate (gh, expr_quad, iopt);
}
// ----------------------------------------------
// 2.4. interpolate
// ----------------------------------------------
namespace details {

template<class Expr>
class field_lazy_terminal_interpolate_rep {
public :
// definitions:

  using size_type   = geo_element::size_type;
  using memory_type = typename Expr::memory_type;
  using value_type  = typename Expr::value_type; // TODO: undeterminated_type<T>
  using scalar_type = typename Expr::scalar_type;
  using float_type  = typename float_traits<scalar_type>::type;
  using geo_type    = geo_basic  <float_type,memory_type>;
  using space_type  = space_basic<float_type,memory_type>;
  using band_type   = band_basic <float_type,memory_type>;
  using vector_element_type = Eigen::Matrix<scalar_type,Eigen::Dynamic,1>;

// allocators:

  field_lazy_terminal_interpolate_rep (
    const space_type&       Xh,
    const Expr&             expr);

// accessors:

  const geo_type&   get_geo()   const { return _Xh.get_geo(); }
  const space_type& get_space() const { return _Xh; }
  bool  is_on_band()            const { return false; }
  band_type get_band()          const { return band_type(); }

  void initialize (const geo_type& omega_K) const;

  void evaluate (
      const geo_type&            omega_K,
      const geo_element&         K,
            vector_element_type& uk) const;
protected:
// internals:

  template <class Value>
  typename std::enable_if<
    is_undeterminated<Value>::value
   ,void
  >::type
  evaluate_internal (
    const geo_type&            omega_K,
    const geo_element&         K,
          vector_element_type& uk) const;

  template <class Value>
  typename std::enable_if<
    ! is_undeterminated<Value>::value
   ,void
  >::type
  evaluate_internal (
    const geo_type&            omega_K,
    const geo_element&         K,
          vector_element_type& uk) const;

// data:
  space_type                              _Xh;
  mutable Expr                            _expr;
  mutable piola_on_pointset<float_type>   _pops;
  mutable disarray<int,memory_type>       _is_marked;
  mutable geo_type                        _prev_omega_K;
  mutable size_type                       _prev_K_dis_ie;
  mutable vector_element_type             _prev_uk;
};
// inlined;
template<class Expr>
field_lazy_terminal_interpolate_rep<Expr>::field_lazy_terminal_interpolate_rep (
  const space_type&       Xh,
  const Expr&             expr)
: _Xh(Xh),
  _expr(expr),
  _pops(),
  _is_marked(),
  _prev_omega_K(),
  _prev_K_dis_ie(std::numeric_limits<size_type>::max()), 
  _prev_uk()
{
}
template<class Expr>
void
field_lazy_terminal_interpolate_rep<Expr>::initialize (const geo_type& omega_K) const
{
  check_macro (omega_K == _Xh.get_geo(),
    "interpolate: incompatible space \""<<_Xh.name()<<"\" with geometry \""
	<< omega_K.name());
  const basis_basic<float_type>& b  = _Xh.get_basis();
  const piola_fem<float_type>&   pf = b.get_piola_fem();
  integrate_option iopt;
  _pops.initialize (omega_K.get_piola_basis(), b, iopt);
  _expr.initialize (_Xh, _pops, iopt);
  _expr.template valued_check<value_type>();
  _prev_omega_K  = omega_K;
  _prev_K_dis_ie = std::numeric_limits<size_type>::max();
  if (_Xh.get_basis().is_continuous()) {
    std::set<size_type> ext_dis_idof;
    _Xh.get_parent_subgeo_owner_dis_indexes (ext_dis_idof);
    _is_marked.resize (_Xh.ownership());
    std::fill (_is_marked.begin(), _is_marked.end(), false);
    _is_marked.set_dis_indexes (ext_dis_idof);
  } // if continuous
}
template<class Expr>
template <class Value>
typename std::enable_if<
  ! is_undeterminated<Value>::value
 ,void
>::type
field_lazy_terminal_interpolate_rep<Expr>::evaluate_internal (
  const geo_type&            omega_K,
  const geo_element&         K,
        vector_element_type& uk) const
{
  // 1) evaluate as a Value
  Eigen::Matrix<Value,Eigen::Dynamic,1> value; // TODO: allocate it once on the mesh
  _expr.evaluate (omega_K, K, value);
  // 2) convert field "Values" at nodes of K to scalar "dofs"
  const piola_fem<float_type>& pf = _Xh.get_basis().get_piola_fem();
  if (pf.transform_need_piola()) {
    const Eigen::Matrix<piola<float_type>,Eigen::Dynamic,1>& piola = _pops.get_piola (omega_K, K);
    for (size_type loc_inod = 0, loc_nnod = value.size(); loc_inod < loc_nnod; ++loc_inod) {
      // be carefull: piola_fem::inv_transform should support inplace call in the "value" arg
      pf.inv_transform (piola[loc_inod], value[loc_inod], value[loc_inod]);
    }
  }
  _Xh.get_basis().compute_dofs (K, value, uk);
}
template<class Expr>
template <class Value>
typename std::enable_if<
  is_undeterminated<Value>::value
 ,void
>::type
field_lazy_terminal_interpolate_rep<Expr>::evaluate_internal (
  const geo_type&            omega_K,
  const geo_element&         K,
        vector_element_type& uk) const
{
  using T = scalar_type;
  // when valued_type is undeterminated at compile time e.g.
  // field uh = lazy_interpolate (Xh, tau_h*vh);
  // => undeterminated : could be scalar, vector or tensor valued
  // so, ask to the Xh space at run-time:
  switch (_Xh.valued_tag()) {
#define _RHEOLEF_switch(VALUED,VALUE)                                    \
    case space_constant::VALUED:                                         \
      (*this).template evaluate_internal<VALUE> (omega_K, K, uk); break;
_RHEOLEF_switch(scalar,T)
_RHEOLEF_switch(vector,point_basic<T>)
_RHEOLEF_switch(tensor,tensor_basic<T>)
_RHEOLEF_switch(unsymmetric_tensor,tensor_basic<T>)
#ifdef TODO
_RHEOLEF_switch(tensor3,tensor3_basic<T>)
_RHEOLEF_switch(tensor4,tensor4_basic<T>)
#endif // TODO
#undef _RHEOLEF_switch
    default: error_macro ("unexpected argument valued tag="<<_Xh.valued_tag());
  }
}
template<class Expr>
void
field_lazy_terminal_interpolate_rep<Expr>::evaluate (
  const geo_type&            omega_K,
  const geo_element&         K,
        vector_element_type& uk) const
{
  if (_prev_omega_K == omega_K && _prev_K_dis_ie == K.dis_ie()) {
    trace_macro("interpolate(K="<<K.name()<<K.dis_ie()<<",prev="<<_prev_K_dis_ie<<"): re-use");
    uk = _prev_uk;
    return;
  }
  trace_macro("interpolate(K="<<K.name()<<K.dis_ie()<<",prev="<<_prev_K_dis_ie<<"): compute");
  (*this).template evaluate_internal<value_type> (omega_K, K, uk);

  // filter for continuous element: each dof value is transmitted only once
  // so the result could be summed as any others field_lazy<Expr>
  if (_Xh.get_basis().is_continuous()) {
    std::vector<size_type> dis_idx;
    _Xh.dis_idof (K, dis_idx);
    check_macro (dis_idx.size() == size_type(uk.size()),
      "invalid sizes: dis_idx.size="<<dis_idx.size()<<", uk.size="<<uk.size());
    size_type my_proc = _Xh.comm().rank();
    for (size_type loc_idof = 0, loc_ndof = dis_idx.size(); loc_idof < loc_ndof; ++loc_idof) {
      size_type dis_idof = dis_idx [loc_idof];
      size_type iproc = _Xh.get_parent_subgeo_owner (dis_idof);
      if (iproc != my_proc || _is_marked.dis_at (dis_idof)) {
        uk [loc_idof] = 0;
        continue;
      }
      // conserve the uk[loc_idof] value and remember it:
      _is_marked.dis_entry (dis_idof) = true;
    }
  }
  _prev_uk       = uk; // expensive to compute, so memorize it for common subexpressions
  _prev_omega_K  = omega_K;
  _prev_K_dis_ie = K.dis_ie();
}
template<class Expr>
class field_lazy_terminal_interpolate: public smart_pointer_nocopy<field_lazy_terminal_interpolate_rep<Expr>>,
                                       public field_lazy_base     <field_lazy_terminal_interpolate<Expr>> {
public :
// definitions:

  using rep    = field_lazy_terminal_interpolate_rep<Expr>;
  using base1  = smart_pointer_nocopy<rep>;
  using base2  = field_lazy_base<field_lazy_terminal_interpolate<Expr>>;
  using size_type   = typename rep::size_type;
  using memory_type = typename rep::memory_type;
  using scalar_type = typename rep::scalar_type;
  using space_type  = typename rep::space_type;
  using geo_type    = typename rep::geo_type;
  using band_type   = typename rep::band_type;
  using vector_element_type = typename rep::vector_element_type;

// allocator:

  field_lazy_terminal_interpolate (
    const space_type&       Xh, 
    const Expr&             expr)
   : base1(new_macro(rep(Xh,expr))),
     base2()
   {}

// accessors:

  const geo_type&   get_geo()   const { return base1::data().get_geo(); }
  const space_type& get_space() const { return base1::data().get_space(); }
  bool  is_on_band()            const { return base1::data().is_on_band(); }
  band_type get_band()          const { return base1::data().get_band(); }

  void initialize (const geo_type& omega_K) const { return base1::data().initialize (omega_K); }

  void evaluate (
      const geo_type&            omega_K,
      const geo_element&         K,
            vector_element_type& uk) const
      { base1::data().evaluate (omega_K, K, uk); }
};
// concept;
template<class Expr> struct is_field_lazy <field_lazy_terminal_interpolate<Expr> > : std::true_type {};

}// namespace details

// 2.4.a. general nonlinear expression
//! @brief see the @ref interpolate_3 page for the full documentation
template<class Expr>
typename std::enable_if<
       details::is_field_expr_v2_nonlinear_arg<Expr>::value
  && ! details::has_field_rdof_interface<Expr>::value // TODO: interpolate also rdof with an extended "evaluate" interface
  && ! details::is_field_function<Expr>::value
 ,details::field_lazy_terminal_interpolate<
    typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr>::type
  >
>::type
lazy_interpolate (
  const space_basic<
    typename float_traits<typename Expr::scalar_type>::type
   ,typename Expr::memory_type>& Xh,
  const Expr& expr)
{
  using wrap_t = typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr>::type;
  return details::field_lazy_terminal_interpolate<wrap_t> (Xh, wrap_t(expr));
}
// 2.4.b. function & functor
//! @brief see the @ref interpolate_3 page for the full documentation
template <class Expr>
inline
typename std::enable_if<
  details::is_field_function<Expr>::value
 ,details::field_lazy_terminal_interpolate<
    details::field_expr_v2_nonlinear_terminal_function<Expr>
  >
>::type
lazy_interpolate (
  const space_basic<
    typename float_traits<typename details::field_expr_v2_nonlinear_terminal_function<Expr>::scalar_type>::type
   ,typename details::field_expr_v2_nonlinear_terminal_function<Expr>::memory_type>& Xh,
  const Expr& expr)
{
  using wrap_t = details::field_expr_v2_nonlinear_terminal_function<Expr>;
  return details::field_lazy_terminal_interpolate<wrap_t> (Xh, wrap_t(expr));
}
#ifdef TO_CLEAN
// 2.4.c. 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
//! @brief see the @ref interpolate_3 page for the full documentation
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);
}

//! @brief see the @ref interpolate_3 page for the full documentation
template <class T, class M, class Expr>
inline
typename std::enable_if<
     details::is_field_wdof<Expr>::value
  && ! details::is_field<Expr>::value
 ,field_basic<T,M>
>::type
lazy_interpolate (const space_basic<T,M>& Xh, const Expr& expr)
{
  return interpolate (Xh, field_basic<T,M>(expr));
}
#endif // TO_CLEAN

}// namespace rheolef
# endif /* _RHEOLEF_FIELD_LAZY_TERMINAL_H */