///
/// 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
/// 
/// =========================================================================
#include "rheolef/geo_domain.h"
#include "rheolef/test.h"

namespace rheolef {

// --------------------------------------------------------------------------
// 1. cstors
// --------------------------------------------------------------------------
template <class T, class M>
test_rep<T,M>::test_rep (const space_type& V)
: _V(V),
  _fops(),
  _is_inside_on_local_sides(false),
  _is_on_band(false),
  _gh()
{
}
template <class T, class M>
test_rep<T,M>::test_rep (const test_rep<T,M>& x)
: _V                        (x._V),
  _fops                     (x._fops),
  _is_inside_on_local_sides (x._is_inside_on_local_sides),
  _is_on_band               (x._is_on_band),
  _gh                       (x._gh)
{
  trace_macro ("*** PHYSICAL COPY OF TEST_REP ***");
}
template <class T, class M>
test_rep<T,M>&
test_rep<T,M>::operator= (const test_rep<T,M>& x)
{
  trace_macro ("*** PHYSICAL ASSIGN OF TEST_REP ***");
  _V                        = x._V;
  _fops                     = x._fops;
  _is_inside_on_local_sides = x._is_inside_on_local_sides;
  _is_on_band               = x._is_on_band;
  _gh                       = x._gh;
  return *this;
}
// --------------------------------------------------------------------------
// 2. initializers
// --------------------------------------------------------------------------
template <class T, class M>
void
test_rep<T,M>::initialize (
  const piola_on_pointset<float_type>&             pops,
  const integrate_option&                          iopt)
{
  _is_on_band = false;
  _is_inside_on_local_sides = iopt._is_inside_on_local_sides;
  _fops.initialize (get_vf_space().get_basis(), pops);
}
template <class T, class M>
void
test_rep<T,M>::initialize (
  const space_basic<float_type,memory_type>&       Xh,
  const piola_on_pointset<float_type>&             pops,
  const integrate_option&                          iopt)
{
  _is_on_band = false;
  _fops.initialize (get_vf_space().get_basis(), pops);
}
template <class T, class M>
void
test_rep<T,M>::initialize (
  const band_basic<float_type,M>&                  gh,
  const piola_on_pointset<float_type>&             pops,
  const integrate_option&                          iopt)
{
  _is_on_band = true;
  _gh = gh;
  _fops.initialize (get_vf_space().get_basis(), pops);
}
// ----------------------------------------------------------------------------
// 3. band : specific cases
// ----------------------------------------------------------------------------
template <class T, class M, class Value, details::differentiate_option::type Diff>
struct evaluate_band_continued {
  void operator() (
    const test_rep<T,M>&                                obj,
    const geo_basic<T,M>&                               omega_K,
    const Eigen::Matrix<piola<T>,Eigen::Dynamic,1>&     piola,
    const reference_element&                            hat_K,
    const reference_element&                            tilde_L,
    const std::vector<size_t>&                          dis_inod_K,
    const std::vector<size_t>&                          dis_inod_L,
    const details::differentiate_option&                gopt,
    Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value) const
  {
    fatal_macro ("on band: unexpected call to differential operator (id="<<Diff<<")");
  }
};
template <class T, class M, class Value>
struct evaluate_band_continued<T,M,Value,details::differentiate_option::none> {
  void operator() (
    const test_rep<T,M>&                                obj,
    const geo_basic<T,M>&                               omega_K,
    const Eigen::Matrix<piola<T>,Eigen::Dynamic,1>&     piola,
    const reference_element&                            hat_K,
    const reference_element&                            tilde_L,
    const std::vector<size_t>&                          dis_inod_K,
    const std::vector<size_t>&                          dis_inod_L,
    const details::differentiate_option&                gopt,
    Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value) const
  {
    Eigen::Matrix<Value,Eigen::Dynamic,1> value_i;
    size_t loc_nnod = value.rows();
    for (size_t loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
      point_basic<T> tilde_xi = inverse_piola_transformation (obj._gh.band(), tilde_L, dis_inod_L, piola[loc_inod].F);
      obj._fops.get_basis_on_pointset().get_basis().evaluate (tilde_L, tilde_xi, value_i);
      value.row(loc_inod) = value_i.transpose();
    }
  }
};
// TODO: DVT_CLEAN_TEST_FEM_DG_BAND regrouper dans fem_on_pointset
template <class T, class GradValue>
static
void 
band_grad_post (
  const tensor_basic<T>&           invDF,
  const tensor_basic<T>&           P,
  const details::differentiate_option& gopt,
  const GradValue&                 hat_grad_u,
        GradValue&                     grad_u)
{
  fatal_macro("band_grad_post: unsupported gradient type: "<<typename_macro(GradValue));
}
template <class T>
static
void 
band_grad_post (
  const tensor_basic<T>&           invDF,
  const tensor_basic<T>&           P,
  const details::differentiate_option& gopt,
  const point_basic<T>&            tilde_grad_u,
        point_basic<T>&                  grad_u)
{
   grad_u = invDF.trans_mult (tilde_grad_u); // TODO: DVT_OPTIM_2D
   if (gopt.surfacic) {
     grad_u = P.trans_mult (grad_u); // TODO: DVT_OPTIM_2D
   }
}
template <class T>
static
void 
band_grad_post (
  const tensor_basic<T>&           invDF,
  const tensor_basic<T>&           P,
  const details::differentiate_option& gopt,
  const tensor_basic<T>&            tilde_grad_u,
        tensor_basic<T>&                  grad_u)
{
  grad_u = tilde_grad_u*invDF; // TODO: DVT_OPTIM_2D
  if (gopt.symmetrized) {
    grad_u = (grad_u + trans(grad_u))/2;
  }
  if (gopt.surfacic) {
    if (! gopt.symmetrized) {
      grad_u =   grad_u*P;
    } else {
      grad_u = P*grad_u*P;
    }
  }
}
// ----------------------------------------------    
// grad evaluate on band
// ----------------------------------------------    
/*
 Let F be the Piola transformation from the reference element:
   F : hat_K --> K
       hat_x --> x = F(hat_x)
 Then the gradient of a basis function u defined on K writes:
   grad(u)(xq) = DF^{-T}*hat_grad(hat_u)(hat_xq)

 When we are working with K as a side of a banded level set surface {phi(x)=0},
 then things are more complex.
 Let
   L in Lambda : a element of the bounding box Lambda
   K in Omega  : a side of the surface Omega, included in L, as K = L inter {phi(x)=0}.
 Let the two Piola transformations from the reference elements:
   F : tilde_L --> L
       tilde_x --> x = F(tilde_x)
   G : hat_K --> K
       hat_x --> x = G(hat_x)

 Let u is a basis function on L: it is defined over tilde_L as u_tilde:
   u(x) = tilde_u (F^{-1}(x))    for all x in K

 The quadrature formula is defined in hat_K with integration point hat_xq and weights hat_wq.
 Thus, integration point hat_xq may be transformed via F^{-1} o G.
 Then, for x=xq=G(hat_xq) :

   u(G(hat_xq)) = tilde_u(F^{-1}(G(hat_xq)))

 Its derivative expresses with a product of the jacobian DF^-T(tilde_x) :

   grad(u)(x) = DF^-T (F^{-1}(x)) grad(u)(F^{-1}(x))  for all x in K

 and then, for x=xq=G(hat_xq)

   grad(u)(G(hat_xq)) = DF^-T (F^{-1}(G(hat_xq))) grad(u)(F^{-1}(G(hat_xq)))

*/

template <class T, class M, class Value>
struct evaluate_band_continued<T,M,Value,details::differentiate_option::gradient> {
  void operator() (
    const test_rep<T,M>&                                obj,
    const geo_basic<T,M>&                               omega_K,
    const Eigen::Matrix<piola<T>,Eigen::Dynamic,1>&     piola,
    const reference_element&                            hat_K,
    const reference_element&                            tilde_L,
    const std::vector<size_t>&                          dis_inod_K,
    const std::vector<size_t>&                          dis_inod_L,
    const details::differentiate_option&                gopt,
    Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value) const
  {
    Eigen::Matrix<Value,Eigen::Dynamic,1> value_i;
    Eigen::Matrix<tensor_basic<T>,Eigen::Dynamic,1> DG;
    if (gopt.surfacic) {
      jacobian_piola_transformation (obj._gh.level_set(), obj._fops.get_piola_on_pointset().get_basis_on_pointset(), hat_K, dis_inod_K, DG);
    }
    size_t d = omega_K.dimension();
    Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> x;
    Eigen::Matrix<Value,Eigen::Dynamic,1> tilde_value_i;
    piola_transformation (obj._gh.level_set(), obj._fops.get_piola_on_pointset().get_basis_on_pointset(), hat_K, dis_inod_K, x);
    size_t loc_nnod = value.rows();
    size_t loc_ndof = value.cols();
    for (size_t loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
      // tilde_x = F_L^{-1}(G_K(hat_x))
      point_basic<T> tilde_xi = inverse_piola_transformation (obj._gh.band(), tilde_L, dis_inod_L, x[loc_inod]);
      // tilde_value_i: tilde_grad_u (tilde_xi)
      obj._fops.get_basis_on_pointset().get_basis().grad_evaluate (tilde_L, tilde_xi, tilde_value_i);
      // grad_u(xi) = DF^{-T}(tilde_xi)*tilde_grad_u(tilde_xi)
      tensor_basic<T> DF;
      tensor_basic<T> P;
      jacobian_piola_transformation (obj._gh.band(), obj._fops.get_piola_on_pointset().get_basis_on_pointset().get_basis(), tilde_L, dis_inod_L, tilde_xi, DF);
      tensor_basic<T> invDF = pseudo_inverse_jacobian_piola_transformation (DF, d, tilde_L.dimension());
      if (gopt.surfacic) {
        // apply also the tangential projection P=(I-n*n)
        // grad(u) = P*DF^{-T}*hat_grad(hat_u) = (DF^{-1}*P)^T*hat_grad(hat_u)
        map_projector (DG[loc_inod], d, hat_K.dimension(), P);
      }
      for (size_t loc_jdof = 0; loc_jdof < loc_ndof; ++loc_jdof) {
        const Value& hat_grad_u = tilde_value_i [loc_jdof];
              Value&     grad_u = value (loc_inod, loc_jdof);
        band_grad_post (invDF, P, gopt, hat_grad_u, grad_u);
      }
    }
  }
};
// --------------------------------------------------------------------------
// 4. evaluate
// --------------------------------------------------------------------------
// external:
// assembly loop on geo_elements K on a domain omega_K
// give its corresponding element K1 in space.geo
template <class T, class M>
const geo_element&
assembly2space_geo_element (
  const geo_basic<T,M>&                space_geo,
  const geo_basic<T,M>&                omega_K,
  const geo_element&                   K_in);

template <class T, class M>
template <class Value, details::differentiate_option::type Diff>
void
test_rep<T,M>::evaluate (
  const geo_basic<T,M>&                               omega_K,
  const geo_element&                                  K,
  const details::differentiate_option&                gopt,
  Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value) const
{
  // omega_K is the domaine associated to element K:
  // it gives the vertices table for K, required by piola
  // ----------------------------------------------------------------
  // identify the effective (omega_K,K) that belongs to V.get_geo
  // and then evaluate
  // ----------------------------------------------------------------
  // TODO DVT_CLEAN_TEST_FEM_DG_BAND 
  if (!_is_on_band) {
    bool use_dom2bgd = (omega_K.variant()      == geo_abstract_base_rep<T>::geo_domain &&
                        _V.get_geo().variant() != geo_abstract_base_rep<T>::geo_domain);
    bool use_bgd2dom = (omega_K.variant()      != geo_abstract_base_rep<T>::geo_domain &&
                        _V.get_geo().variant() == geo_abstract_base_rep<T>::geo_domain);
    if (!use_bgd2dom && !use_dom2bgd) {
      bool use_dg_on_sides 
        =  omega_K.variant() == geo_abstract_base_rep<T>::geo_domain_indirect &&
           1 + K.dimension() == omega_K.get_background_geo().map_dimension()  &&
  	  _V.get_basis().is_discontinuous();
      if (!use_dg_on_sides) {
        // usual case
        _fops.template evaluate<M,Value,Diff> (omega_K, K, gopt, value);
      } else {
        // DG on omega_K="sides" or "boundary" or other d-1 sides domain
        // but for a side K, will access to its neighbours L0 & L1 for e.g. dot(grad(u),normal())
        // omits "inner(u)" on a sides-domain => assume a boundary side
        geo_basic<T,M> omega_L_eff = omega_K.get_background_geo();
        size_type L_map_d = K.dimension() + 1;
        check_macro (L_map_d == omega_L_eff.map_dimension(),
           "unexpected dimension for side S="<<K.name()<<K.dis_ie()<<" in domain "<<omega_L_eff.name());
        check_macro (K.master(1) == std::numeric_limits<size_type>::max(),
          "unexpected non-boundary side S="<<K.name()<<K.dis_ie()<<" in domain "<<omega_L_eff.name());
        size_type L_dis_ie = K.master(0);
        check_macro (L_dis_ie != std::numeric_limits<size_type>::max(),
          "unexpected isolated side S="<<K.name()<<K.dis_ie());
        const geo_element& L_eff = omega_L_eff.dis_get_geo_element (L_map_d, L_dis_ie);
        side_information_type sid;
        L_eff.get_side_informations (K, sid);
        _fops.template evaluate_on_side<M,Value,Diff> (omega_L_eff, L_eff, sid, gopt, value);
      }
    } else if (use_dom2bgd) {
      bool use_dg_on_sides 
        =  omega_K.variant() != geo_abstract_base_rep<T>::geo &&
           1 + K.dimension() == _V.get_geo().get_background_geo().map_dimension()  &&
  	  _V.get_basis().is_discontinuous();
      if (!use_dg_on_sides) {
        // usual case
        _fops.template evaluate<M,Value,Diff> (omega_K, K, gopt, value);
      } else {
        // e.g. wh = interpolate(Wh,uh);
        // with uh in Vh=space("square",P1d) and Wh=space("square[left]",P1d);
        const geo_element& K_space = assembly2space_geo_element (_V.get_geo(), omega_K, K);
        size_type L_space_map_d = K_space.dimension() + 1;
        check_macro (L_space_map_d == _V.get_geo().map_dimension(),
           "unexpected dimension for side S="<<K_space.name()<<K_space.dis_ie()<<" in domain "<<_V.get_geo().name());
        check_macro (K_space.master(1) == std::numeric_limits<size_type>::max(),
          "unexpected non-boundary side S="<<K_space.name()<<K_space.dis_ie()<<" in domain "<<_V.get_geo().name());
        size_type L_space_dis_ie = K_space.master(0);
        check_macro (L_space_dis_ie != std::numeric_limits<size_type>::max(),
          "unexpected isolated side S="<<K_space.name()<<K_space.dis_ie());
        const geo_element& L_space = _V.get_geo().dis_get_geo_element (L_space_map_d, L_space_dis_ie);
        side_information_type sid;
        L_space.get_side_informations (K_space, sid);
        _fops.template evaluate_on_side<M,Value,Diff> (_V.get_geo(), L_space, sid, gopt, value);
      }
    } else { // use_bdg2dom
      // space X (omega);
      // space Y (omega["boundary"]);
      // trial u(X);
      // test  v(Y);
      // form m = integrate ("boundary", u*v);
      //  => v is on compacted geo_domain Y.geo but elements K 
      //    during integration process are on domain omega["boundary"]
      //  => element K need to be replaced by an element in Y.geo
      //     for computing correctly Y.dis_idof(K)
      check_macro (
         omega_K.get_background_geo().name() == _V.get_geo().name() ||
         omega_K.get_background_geo().name() == _V.get_geo().get_background_geo().name(),
  	"unexpected integration domain \"" << omega_K.name()
          << "\" based on geo \"" << omega_K.get_background_geo().name()
          << "\" together with a space on geo \"" << _V.get_geo().name()
          << "\" based on geo \"" << _V.get_geo().get_background_geo().name()<<"\"");
      geo_basic<T,M> omega_K_eff = _V.get_geo();
      const geo_element& K_eff = omega_K_eff.bgd2dom_geo_element (K);
      _fops.template evaluate<M,Value,Diff> (omega_K_eff, K_eff, gopt, value);
    }
  } else { // is_on_band
    // find L in the band such that K is the intersected side of L
    // as the zero level set isosurface
    std::vector<size_type> dis_inod_L;
    std::vector<size_type> dis_inod_K;
    size_type first_dis_ie = _gh.level_set().sizes().ownership_by_dimension [K.dimension()].first_index();
    check_macro (K.dis_ie() >= first_dis_ie, "unexpected intersected side element K.dis_ie="<<K.dis_ie());
    size_type K_ie = K.dis_ie() - first_dis_ie;
    size_type L_ie = _gh.sid_ie2bnd_ie (K_ie);
    const geo_element& L = _gh.band() [L_ie];
    reference_element hat_K   = K.variant();
    reference_element tilde_L = L.variant();
    // use the basis on L, the element from the band that contains the 
    // intersected side K on the zero level set surface
    const piola_fem<T>& pf = _fops.get_basis_on_pointset().get_basis().get_piola_fem();
    check_macro (! pf.transform_need_piola(), "band: unsupported piola-based basis");
    const Eigen::Matrix<piola<T>,Eigen::Dynamic,1>& piola = _fops.get_piola_on_pointset().get_piola (_gh.level_set(), K);
    size_type loc_nnod = piola.size();
    size_type loc_ndof = get_vf_space().get_basis().ndof (tilde_L);
    value.resize (loc_nnod, loc_ndof);
    _gh.level_set().dis_inod (K, dis_inod_K);
    _gh.band()     .dis_inod (L, dis_inod_L);
    evaluate_band_continued<T,M,Value,Diff> eval_band_cont;
    eval_band_cont (*this, omega_K, piola, hat_K, tilde_L, dis_inod_K, dis_inod_L, gopt, value);
  }
}
// --------------------------------------------------------------------------
// 5. evaluate on side
// --------------------------------------------------------------------------
template <class T, class M>
template <class Value, details::differentiate_option::type Diff>
void
test_rep<T,M>::evaluate_on_side (
  const geo_basic<T,M>&                               omega_K,
  const geo_element&                                  K,
  const side_information_type&                        sid,
  const details::differentiate_option&                gopt,
  Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value) const
{
  size_type space_map_d = get_vf_space().get_geo().map_dimension();
  // TODO DVT_CLEAN_TEST_FEM_DG_BAND 
  bool is_broken = get_vf_space().get_geo().is_broken() && space_map_d + 1 == K.dimension();
  if (!_is_on_band) {
    if (!is_broken) {
      _fops.template evaluate_on_side<M,Value,Diff> (omega_K, K, sid, gopt, value);
    } else { // is_broken: _is_inside_on_local_sides
      check_macro(_is_inside_on_local_sides, "unexpected broken case");
      size_type isid = sid.loc_isid;
      size_type dis_isid = (K.dimension() == 1) ? K[isid] : (K.dimension() == 2) ? K.edge(isid) : K.face(isid);
      const geo_element& S = omega_K.dis_get_geo_element (space_map_d, dis_isid);
      bool have_bgd_dom = get_vf_space().get_geo().variant() == geo_abstract_base_rep<T>::geo_domain && get_vf_space().get_geo().name() != omega_K.name();
      const geo_element* dom_S_ptr = (!have_bgd_dom) ? &S : &(get_vf_space().get_geo().bgd2dom_geo_element(S));
      const geo_element& dom_S = *dom_S_ptr;
      Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic> value_sid;
      _fops.template evaluate<M,Value,Diff> (omega_K, dom_S, gopt, value_sid);
      // next, locally assembly value_sid into value
      reference_element hat_K = K;
      size_type loc_nnod = value_sid.rows();
      size_type loc_ndof = 0;
      size_type loc_isid_ndof = 0;
      size_type start_loc_isid_idof = 0;
      for (size_type jsid = 0, nsid = hat_K.n_side(); jsid < nsid; ++jsid) {
        reference_element hat_Sj = hat_K.side(jsid);
        size_type loc_jsid_ndof = get_vf_space().get_basis().ndof (hat_Sj);
        if (jsid == isid) {
          start_loc_isid_idof = loc_ndof;
                loc_isid_ndof = loc_jsid_ndof;
        }
        loc_ndof += loc_jsid_ndof;
      }
      value.resize (loc_nnod, loc_ndof);
      value.fill (Value());
      for (size_type loc_isid_jdof = 0; loc_isid_jdof < loc_isid_ndof; ++loc_isid_jdof) {
        size_type loc_jdof = start_loc_isid_idof + loc_isid_jdof;
        for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
          value(loc_inod,loc_jdof) = value_sid (loc_inod, loc_isid_jdof); 
        }
      }
    }
  } else {
    check_macro (!_is_on_band,  "evaluate_on_side: not yet on band");
  }
}
// ----------------------------------------------    
// 6. local dg merge dofs on an internal side
// ----------------------------------------------    
template <class T, class M>
template <class Value>
void
test_rep<T,M>::local_dg_merge_on_side (
  const geo_basic<T,M>&                                     omega_K,
  const geo_element&                                        S,
  const geo_element&                                        K0,
  const geo_element&                                        K1,
  const Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value0,
  const Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value1,
  Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>&       value) const
{
  check_macro (value0.rows() == value1.rows(), "invalid node number");
  size_type loc_nnod  = value0.rows();
  size_type loc_ndof0 = value0.cols();
  size_type loc_ndof1 = value1.cols();
  value.resize (loc_nnod, loc_ndof0 + loc_ndof1);
  for (size_type loc_inod = 0; loc_inod < loc_nnod;  ++loc_inod) {
    for (size_type loc_idof = 0; loc_idof < loc_ndof0; ++loc_idof) {
      value (loc_inod,loc_idof) = value0 (loc_inod,loc_idof);
    }
    for (size_type loc_idof = 0; loc_idof < loc_ndof1; ++loc_idof) {
      value (loc_inod,loc_ndof0+loc_idof) = value1 (loc_inod,loc_idof); // TODO: DVT_EIGEN_BLAS2
    }
  }
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------

#define _RHEOLEF_instanciation_value_diff(T,M,Value,Diff)			\
template void test_rep<T,M>::evaluate<Value,Diff> (				\
  const geo_basic<T,M>&                               omega_K,			\
  const geo_element&                                  K,			\
  const details::differentiate_option&                gopt,			\
  Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value) const;		\
template void test_rep<T,M>::evaluate_on_side<Value,Diff> (			\
  const geo_basic<T,M>&                               omega_K,			\
  const geo_element&                                  K,			\
  const side_information_type&                        sid,			\
  const details::differentiate_option&                gopt,			\
  Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value) const;		\

#define _RHEOLEF_instanciation_value(T,M,Value)					\
template void test_rep<T,M>::local_dg_merge_on_side (				\
  const geo_basic<T,M>&                                     omega_K,		\
  const geo_element&                                        S,			\
  const geo_element&                                        K0,			\
  const geo_element&                                        K1,			\
  const Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value0,		\
  const Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value1,		\
  Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>&       value) const;	\
_RHEOLEF_instanciation_value_diff(T,M,Value,details::differentiate_option::none) \
_RHEOLEF_instanciation_value_diff(T,M,Value,details::differentiate_option::gradient) \
_RHEOLEF_instanciation_value_diff(T,M,Value,details::differentiate_option::divergence) \
_RHEOLEF_instanciation_value_diff(T,M,Value,details::differentiate_option::curl) \

#define _RHEOLEF_instanciation(T,M)	 					\
template class test_rep<T,M>;							\
_RHEOLEF_instanciation_value(T,M,T)			 			\
_RHEOLEF_instanciation_value(T,M,point_basic<T>)	 			\
_RHEOLEF_instanciation_value(T,M,tensor_basic<T>)	 			\
_RHEOLEF_instanciation_value(T,M,tensor3_basic<T>)	 			\
_RHEOLEF_instanciation_value(T,M,tensor4_basic<T>)	 			\

_RHEOLEF_instanciation(Float,sequential)
#ifdef _RHEOLEF_HAVE_MPI
_RHEOLEF_instanciation(Float,distributed)
#endif // _RHEOLEF_HAVE_MPI

} // namespace rheolef