librheolef-dev 6.7-6 / usr / include / rheolef / field_expr_v2_variational_terminal.h

#ifndef _RHEOLEF_FIELD_EXPR_V2_VARIATIONAL_TERMINAL_H
#define _RHEOLEF_FIELD_EXPR_V2_VARIATIONAL_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
///
/// =========================================================================
//
// terminals variational expressions are used for field assembly
// e.g. for right-hand-side of linear systems
//
// author: Pierre.Saramito@imag.fr
//
// date: 21 september 2015 
//
// Notes: use template expressions and SFINAE techniques
//
// SUMMARY:
// 1. concept
// 2. grad, grad_s, D, etc
// 3. div, div_s
// 4. curl
// 5. jump, average, inner, outer
//
#include "rheolef/field_expr_v2_nonlinear.h"
#include "rheolef/test.h" // for grad_option_type

namespace rheolef {

// -------------------------------------------------------------------
// 1. concept
// -------------------------------------------------------------------
namespace details {

// Define a trait type for detecting field expression valid arguments
template<class T>                       struct is_field_expr_v2_variational_arg                         : std::false_type {};
template<class T, class M, class VfTag> struct is_field_expr_v2_variational_arg<test_basic<T,M,VfTag> > : std::true_type {};

} // namespace details
// ---------------------------------------------------------------------------
// 2. grad, grad_s, D, etc
// ---------------------------------------------------------------------------
namespace details {

template<class Expr>
class field_expr_v2_variational_grad {
public:
// typedefs:

  typedef geo_element::size_type                        size_type;
  typedef typename Expr::memory_type                    memory_type;
  typedef typename scalar_traits<typename Expr::value_type>::type 		
	  						scalar_type;
  typedef typename space_constant::rank_up<typename Expr::value_type>::type
                                                        value_type;
  typedef typename  float_traits<scalar_type>::type 	float_type;
  typedef space_basic<scalar_type,memory_type>		space_type;
  typedef typename Expr::vf_tag_type                    vf_tag_type;
  typedef typename details::dual_vf_tag<vf_tag_type>::type
                                                        vf_dual_tag_type;
  typedef field_expr_v2_variational_grad<Expr>          self_type;
  typedef field_expr_v2_variational_grad<typename Expr::dual_self_type>  
                                                        dual_self_type;

// alocators:

  field_expr_v2_variational_grad (const Expr& expr, const grad_option_type& opt = grad_option_type())
    : _expr(expr),
      _opt(opt)
  {
    check_macro (opt.broken
              || get_vf_space().get_numbering().is_continuous()
    	      || get_vf_space().get_numbering().name() == "bubble",
      "grad(.): unexpected " << get_vf_space().get_numbering().name() 
	<< " discontinuous approximation (HINT: consider grad_h(.))");
  }
  field_expr_v2_variational_grad (const field_expr_v2_variational_grad<Expr>& x)
    : _expr(x._expr),
      _opt(x._opt)
  {}

// accessors:

  const space_type&  get_vf_space()  const { return _expr.get_vf_space(); }
  static const space_constant::valued_type valued_hint = space_constant::valued_tag_traits<value_type>::value;
  space_constant::valued_type valued_tag() const {
    space_constant::valued_type v = _expr.valued_tag();
    switch (v) {
      case space_constant::scalar: return space_constant::vector;
      case space_constant::vector: return space_constant::unsymmetric_tensor;
      case space_constant::tensor: return space_constant::tensor3;
      default:
        fatal_macro ("unexpected " << space_constant::valued_name(v) << "-valued argument for grad() operator");
        return space_constant::last_valued;
    }
  }
  size_type n_derivative() const { return _expr.n_derivative() + 1; }

// mutable modifiers:

  void initialize (const geo_basic<float_type,memory_type>& dom, const quadrature<float_type>& quad, bool ignore_sys_coord) const { 
    _expr.initialize (dom, quad, ignore_sys_coord);
  }
  void initialize (const band_basic<float_type,memory_type>& gh, const quadrature<float_type>& quad, bool ignore_sys_coord) const {  
    _expr.initialize (gh, quad, ignore_sys_coord); 
  }
  void element_initialize (const geo_element& K) const {
    _expr.element_initialize (K);
  }
  void element_initialize_on_side (const geo_element& K, const side_information_type& sid) {
    _expr.element_initialize_on_side (K, sid);
  }
  template<class ValueType>
  void basis_evaluate (const reference_element& hat_K, size_type q, std::vector<ValueType>& value) const {
    _expr.grad_basis_evaluate (hat_K, q, _opt, value);
  }
  template<class ValueType>
  void valued_check() const {
    typedef typename space_constant::rank_down<ValueType>::type A1; // may be defined when ValueType is
    _expr.template valued_check<A1>();
  }
protected:
// data:
  Expr               _expr;
  grad_option_type   _opt;
};
template<class Expr> struct is_field_expr_v2_variational_arg    <field_expr_v2_variational_grad<Expr> > : std::true_type {};

} // namespace details

// grad(v)
template<class Expr>
inline
typename
std::enable_if<
  details::is_field_expr_v2_variational_arg<Expr>::value
 ,details::field_expr_v2_variational_grad<Expr>
>::type
grad (const Expr& expr)
{
  return details::field_expr_v2_variational_grad <Expr> (expr);
}
// grad_s(v) 
template<class Expr>
inline
typename
std::enable_if<
  details::is_field_expr_v2_variational_arg<Expr>::value
 ,details::field_expr_v2_variational_grad<Expr>
>::type
grad_s (const Expr& expr)
{
  details::grad_option_type   opt;
  opt.surfacic = true;
  return details::field_expr_v2_variational_grad <Expr> (expr, opt);
}
// grad_h(v) 
template<class Expr>
inline
typename
std::enable_if<
  details::is_field_expr_v2_variational_arg<Expr>::value
 ,details::field_expr_v2_variational_grad<Expr>
>::type
grad_h (const Expr& expr)
{
  details::grad_option_type   opt;
  opt.broken = true;
  return details::field_expr_v2_variational_grad <Expr> (expr, opt);
}
// D(v) 
template<class Expr>
inline
typename
std::enable_if<
  details::is_field_expr_v2_variational_arg<Expr>::value
 ,details::field_expr_v2_variational_grad<Expr>
>::type
D (const Expr& expr)
{
  details::grad_option_type   opt;
  opt.symmetrized = true;
  return details::field_expr_v2_variational_grad <Expr> (expr, opt);
}
// Ds(v) 
template<class Expr>
inline
typename
std::enable_if<
  details::is_field_expr_v2_variational_arg<Expr>::value
 ,details::field_expr_v2_variational_grad<Expr>
>::type
Ds (const Expr& expr)
{
  details::grad_option_type   opt;
  opt.symmetrized = true;
  opt.surfacic    = true;
  return details::field_expr_v2_variational_grad <Expr> (expr, opt);
}
// Dh(v) 
template<class Expr>
inline
typename
std::enable_if<
  details::is_field_expr_v2_variational_arg<Expr>::value
 ,details::field_expr_v2_variational_grad<Expr>
>::type
Dh (const Expr& expr)
{
  details::grad_option_type   opt;
  opt.symmetrized = true;
  opt.broken      = true;
  return details::field_expr_v2_variational_grad <Expr> (expr, opt);
}
// ---------------------------------------------------------------------------
// 3. div, div_s
// ---------------------------------------------------------------------------
namespace details {

template<class Expr>
class field_expr_v2_variational_div {
public:
// typedefs:

  typedef geo_element::size_type                        size_type;
  typedef typename Expr::memory_type                    memory_type;
  typedef typename space_constant::rank_down<typename Expr::value_type>::type
                                                        value_type;
  typedef typename scalar_traits<typename Expr::value_type>::type 		
	  						scalar_type;
  typedef typename  float_traits<scalar_type>::type 	float_type;
  typedef space_basic<scalar_type,memory_type>		space_type;
  typedef typename Expr::vf_tag_type                    vf_tag_type;
  typedef typename details::dual_vf_tag<vf_tag_type>::type
                                                        vf_dual_tag_type;
  typedef field_expr_v2_variational_div<Expr>                       self_type;
  typedef field_expr_v2_variational_div<typename Expr::dual_self_type>  
                                                        dual_self_type;

// alocators:

  field_expr_v2_variational_div (const Expr& expr, const grad_option_type& opt = grad_option_type())
    : _expr(expr),
      _opt(opt)
  {
    check_macro (opt.broken
              || get_vf_space().get_numbering().is_continuous()
    	      || get_vf_space().get_numbering().name() == "bubble",
      "div(.): unexpected " << get_vf_space().get_numbering().name() 
	<< " discontinuous approximation (HINT: consider div_h(.))");
  }

// accessors:

  const space_type&  get_vf_space()  const { return _expr.get_vf_space(); }
  static const space_constant::valued_type valued_hint = space_constant::valued_tag_traits<value_type>::value;
  space_constant::valued_type valued_tag() const {
    space_constant::valued_type v = _expr.valued_tag();
    switch (v) {
      case space_constant::vector:             return space_constant::scalar;
      case space_constant::tensor:
      case space_constant::unsymmetric_tensor: return space_constant::vector;
      default:
        fatal_macro ("unexpected " << space_constant::valued_name(v) << "-valued argument for div() operator");
        return space_constant::last_valued;
    }
  }
  size_type n_derivative() const { return _expr.n_derivative() + 1; }

// mutable modifiers:

  void initialize (const geo_basic<float_type,memory_type>& dom, const quadrature<float_type>& quad, bool ignore_sys_coord) const { 
    _expr.initialize (dom, quad, ignore_sys_coord);
  }
  void initialize (const band_basic<float_type,memory_type>& gh, const quadrature<float_type>& quad, bool ignore_sys_coord) const {  
    _expr.initialize (gh, quad, ignore_sys_coord); 
  }
  void element_initialize (const geo_element& K) const {
    _expr.element_initialize (K);
  }
  void element_initialize_on_side (const geo_element& K, const side_information_type& sid) {
    _expr.element_initialize_on_side (K, sid);
  }
  template<class ValueType>
  void basis_evaluate (const reference_element& hat_K, size_type q, std::vector<ValueType>& value) const {
    _expr.div_basis_evaluate (hat_K, q, _opt, value);
  }
  template<class ValueType>
  void valued_check() const {
    _expr.template div_valued_check<ValueType>();
  }
protected:
// data:
  Expr              _expr;
  grad_option_type  _opt;
};
template<class Expr> struct is_field_expr_v2_variational_arg    <field_expr_v2_variational_div<Expr> > : std::true_type {};

} // namespace details

// div(v)
template<class Expr>
inline
typename
std::enable_if<
  details::is_field_expr_v2_variational_arg<Expr>::value
 ,details::field_expr_v2_variational_div<Expr>
>::type
div (const Expr& expr)
{
  return details::field_expr_v2_variational_div <Expr> (expr);
}
// div_s(v)
template<class Expr>
inline
typename
std::enable_if<
  details::is_field_expr_v2_variational_arg<Expr>::value
 ,details::field_expr_v2_variational_div<Expr>
>::type
div_s (const Expr& expr)
{
  details::grad_option_type   opt;
  opt.surfacic = true;
  return details::field_expr_v2_variational_div <Expr> (expr, opt);
}
// div_h(v)
template<class Expr>
inline
typename
std::enable_if<
  details::is_field_expr_v2_variational_arg<Expr>::value
 ,details::field_expr_v2_variational_div<Expr>
>::type
div_h (const Expr& expr)
{
  details::grad_option_type   opt;
  opt.broken = true;
  return details::field_expr_v2_variational_div <Expr> (expr, opt);
}
// ---------------------------------------------------------------------------
// 4. curl
// ---------------------------------------------------------------------------
namespace details {

template<class Expr>
class field_expr_v2_variational_curl {
public:
// typedefs:

  typedef geo_element::size_type                        size_type;
  typedef typename Expr::memory_type                    memory_type;
  typedef typename scalar_traits<typename Expr::value_type>::type 		
	  						scalar_type;
  // value_type = vctor  when d=2 and Expr is scalar or when d=3
  //            = scalar when d=2 and Expr is vector 
  // thus is undeterminated at compile-time
  typedef undeterminated_basic<scalar_type>		value_type;
  typedef typename  float_traits<scalar_type>::type 	float_type;
  typedef space_basic<scalar_type,memory_type>		space_type;
  typedef typename Expr::vf_tag_type                    vf_tag_type;
  typedef typename details::dual_vf_tag<vf_tag_type>::type
                                                        vf_dual_tag_type;
  typedef field_expr_v2_variational_curl<Expr>          self_type;
  typedef field_expr_v2_variational_curl<typename Expr::dual_self_type>  
                                                        dual_self_type;

// alocators:

  field_expr_v2_variational_curl (const Expr& expr, const grad_option_type& opt = grad_option_type())
    : _expr(expr),
      _opt(opt)
  {
    check_macro (opt.broken
              || get_vf_space().get_numbering().is_continuous()
    	      || get_vf_space().get_numbering().name() == "bubble",
      "curl(.): unexpected " << get_vf_space().get_numbering().name() 
	<< " discontinuous approximation (HINT: consider curl_h(.))");
  }

// accessors:

  const space_type&  get_vf_space()  const { return _expr.get_vf_space(); }
  static const space_constant::valued_type valued_hint = space_constant::valued_tag_traits<value_type>::value;
  space_constant::valued_type valued_tag() const {
    space_constant::valued_type arg_v = _expr.valued_tag();
    switch (arg_v) {
      case space_constant::scalar:    return space_constant::vector;
      case space_constant::vector: {
        size_type d = get_vf_space().get_geo().dimension();
        return (d==2) ? space_constant::scalar : space_constant::vector;
      }
      default:
        fatal_macro ("unexpected " << space_constant::valued_name(arg_v) << "-valued argument for curl() operator");
        return space_constant::last_valued;
    }
  }
  size_type n_derivative() const { return _expr.n_derivative() + 1; }

// mutable modifiers:

  void initialize (const geo_basic<float_type,memory_type>& dom, const quadrature<float_type>& quad, bool ignore_sys_coord) const { 
    _expr.initialize (dom, quad, ignore_sys_coord);
  }
  void initialize (const band_basic<float_type,memory_type>& gh, const quadrature<float_type>& quad, bool ignore_sys_coord) const {  
    _expr.initialize (gh, quad, ignore_sys_coord); 
  }
  void element_initialize (const geo_element& K) const {
    _expr.element_initialize (K);
  }
  void element_initialize_on_side (const geo_element& K, const side_information_type& sid) {
    _expr.element_initialize_on_side (K, sid);
  }
  template<class ValueType>
  void basis_evaluate (const reference_element& hat_K, size_type q, std::vector<ValueType>& value) const {
    _expr.curl_basis_evaluate (hat_K, q, _opt, value);
  }
  void _valued_check_internal(const scalar_type&) const {
    _expr.template valued_check<point_basic<scalar_type> >();
    size_type d = get_vf_space().get_geo().dimension();
    check_macro (d==2, "unexpected "<<d<<"D physical dimension for the scalar-valued curl() operator");
  }
  void _valued_check_internal(const point_basic<scalar_type>&) const {
    size_type d = get_vf_space().get_geo().dimension();
    check_macro (d==2 || d==3, "unexpected "<<d<<"D physical dimension for the vector-valued curl() operator");
    if (d == 2) {
      _expr.template valued_check<scalar_type>();
    } else {
      _expr.template valued_check<point_basic<scalar_type> >();
    }
  }
  void _valued_check_internal(const tensor_basic<scalar_type>&) const {
    fatal_macro ("unexpected tensor-valued result for the curl() operator");
  }
  template<class ValueType>
  void valued_check() const {
    _valued_check_internal(ValueType());
  }
protected:
// data:
  Expr              _expr;
  grad_option_type  _opt;
};
template<class Expr> struct is_field_expr_v2_variational_arg    <field_expr_v2_variational_curl<Expr> > : std::true_type {};

} // namespace details

// curl(v)
template<class Expr>
inline
typename
std::enable_if<
  details::is_field_expr_v2_variational_arg<Expr>::value
 ,details::field_expr_v2_variational_curl<Expr>
>::type
curl (const Expr& expr)
{
  return details::field_expr_v2_variational_curl <Expr> (expr);
}
// bcurl(v) = Batchelor curl
template<class Expr>
inline
typename
std::enable_if<
  details::is_field_expr_v2_variational_arg<Expr>::value
 ,details::field_expr_v2_variational_curl<Expr>
>::type
bcurl (const Expr& expr)
{
  details::grad_option_type    opt;
  opt.batchelor_curl = true;
  return details::field_expr_v2_variational_curl <Expr> (expr, opt);
}
// ---------------------------------------------------------------------------
// 5. jump, average, inner, outer
// ---------------------------------------------------------------------------
// discontinuous Galerkin operators
// in expressions templates for variationnal formulations
//
/*

SPECIFICATION: a first example

  Let v be a function defined over Omega
        and discontinuous accross internal sides (e.g. Pkd).
  Let f be a function defined on Oemga.

  We want to assembly
    l(v) = int_{internal sides} f(x) [v](x) ds
  where [v] is the jump of v accross internal sides.

  => l(v) = sum_{K is internal sides} int_K f(x) [v](x) ds

  Let K be an internal side of the mesh of Omega.

    int_K f(x) [v](x) ds
    = int_{hat_K} f(F(hat_x))
                  [v](F(hat_x))
                  det(DF(hat_x)) d hat_s

  where F is the piola transformation from the reference element hat_K to K:
    F : hat_K ---> K
	hat_x |--> x = F(hat_x)

  The fonction v is not defined on a basis over internal sides K but over
  elements L of the mesh of Omega.
  Let L0 and L1 the two elements such that K is the common side of L0 and L1
  and K is oriented from L0 to L1:
 	[v] = v0 - v1 on K, where v0=v/L0 and v1=v/L1.
 
  Let G0 the piola transformation from the reference element tilde_L to L0:
    G0 : tilde_L ---> L0
	 tilde_x |--> x = G0(tilde_x)

  Conversely, let G1 the piola transformation from the reference element tilde_L to L1.

    int_K f(x) [v](x) ds
    = int_{hat_K} f(F(hat_x))
                  (v0-v1)(F(hat_x))
                  det(DF(hat_x)) d hat_s

  The the basis fonction v0 and v1 are defined by using tilde_v, on the reference element tilde_L:
                  v0(x) = tilde_v (G0^{-1}(x))
                  v1(x) = tilde_v (G1^{-1}(x))
  and with x=F(hat_x):
                  v0(F(hat_x)) = tilde_v (G0^{-1}(F(hat_x)))
                  v1(F(hat_x)) = tilde_v (G1^{-1}(F(hat_x)))

  Thus:
    int_K f(x) [v](x) ds
    = int_{hat_K} f(F(hat_x)) 
                  (   tilde_v (G0^{-1}(F(hat_x)))
                    - tilde_v (G1^{-1}(F(hat_x))) )
                  det(DF(hat_x)) ds
  
  Observe that H0=G0^{-1}oF is linear:
     H0 : hat_K ---> tilde0_K subset tilde_L
          hat_x ---> tilde0_x = H0(hat_x)
  Conversely:
     H1 : hat_K ---> tilde1_K subset tilde_L
          hat_x ---> tilde1_x = H1(hat_x)
  Thus, K linearly transforms by H0 into a side tilde0_K of the reference element tilde_L
  and, by H1, into another side tilde1_K of tilde_L.

    int_K f(x) [v](x) ds
    = int_{hat_K} f(F(hat_x)) 
                  (   tilde_v (H0(hat_x))
                    - tilde_v (H1(hat_x)) )
                  det(DF(hat_x)) ds

  Let (hat_xq, hat_wq)_{q=0...} a quadrature formulae over hat_K.  
  The integral becomes:

    int_K f(x) [v](x) ds
    = sum_q f(F(hat_xq)) 
            (   tilde_v (H0(hat_xq))
              - tilde_v (H1(hat_xq)) )
            det(DF(hat_xq)) hat_wq

  Then, the basis functions tilde_v can be computed one time for all
  over all the sides tilde(i)_K of the reference element tilde_L, i=0..nsides(tilde_L)
  at the quadratures nodes tilde(i)_xq = Hi(hat_xq):
            tilde_v (Hi(hat_xq)),  i=0..nsides(tilde_L),  q=0..nq(hat_K)

SPECIFICATION: a second example

  We want to assembly
    l(v) = int_{internal sides} f(x) [grad(v).n](x) ds
  where [grad(v).n] is the jump of the normal derivative of v accross internal sides.

  Let K be an internal side of the mesh of Omega.

    int_K f(x) [grad(v).n](x) ds
    = int_{hat_K} f(F(hat_x))
                 [grad(v).n](F(hat_x))
                 det(DF(hat_x)) d hat_s

    =   int_{hat_K} f(F(hat_x))
                 (grad(v0).n)(F(hat_x))
                  det(DF(hat_x)) d hat_s
      - int_{hat_K} f(F(hat_x))
                 (grad(v1).n)(F(hat_x))
                 det(DF(hat_x)) d hat_s

   where v0=v/L0 and v1=v/L1 and Li are the two elements containing the side K.
   Let us fix one of the Li and omits the i subscript.
   The computation reduces to evaluate:

    int_K f(x) grad(v).n(x) ds
    = sum_q f(F(hat_xq))
                 (grad(v).n)(F(hat_xq))
                 det(DF(hat_xq)) hat_wq
 
   From the gradient transformation:

     grad(v)(F(hat_xq)) = DG^{-T}(H(hat_xq)) * tilde_grad(tilde_v)(H(hat_xq))

   where H = G^{-1}oF is linear from hat_K to tilde_K subset tilde_L.

    int_K f(x) grad(v).n(x) ds
    = sum_q f(F(hat_xq))
                 DG^{-T}(H(hat_xq))*tilde_grad(tilde_v)(H(hat_xq))
                 .n(F(hat_xq))
                 det(DF(hat_xq)) hat_wq

   We can evaluate one time for all the gradients of basis functions tilde_v
   on the quadrature nodes of each sides tilde_K of tilde_L :
                 tilde_grad(tilde_v)(H(hat_xq))
   The piola basis functions and their derivatives are also evaluated one time for all on these nodes :
     		 DG^{-T}(H(hat_xq))

   The normal vector 
                 n(xq), xq=F(hat_xq), q=...
   should be evaluated on K, not on L that has no normal vector.

IMPLEMENTATION: bassis evaluation => test.cc

  The basis_on_pointset class extends to the case of an integration over a side of

     test_rep<T,M>::initialize (const geo_basic<float_type,M>& dom, const quadrature<T>& quad, bool ignore_sys_coord) const {
       _basis_on_quad.set (quad, get_vf_space().get_numbering().get_basis());
       _piola_on_quad.set (quad, get_vf_space().get_geo().get_piola_basis());
	 => inchange'
     }
     test_rep<T,M>::element_initialize (const geo_element& L, size_type loc_isid=-1) const {
        if (loc_isid != -1) {
	  basis_on_quad.restrict_on_side (tilde_L, loc_isid);
	  piola_on_quad.restrict_on_side (tilde_L, loc_isid);
        }
     }
     test_rep<T,M>::basis_evaluate (...) {
        // Then, a subsequent call to 
	basis_on_quad.evaluate (tilde_L, q);
        // will restrict to the loc_isid part.
     }

IMPLEMENTATION: normal vector => field_vf_expr.h & field_nl_expr_terminal.h
  on propage des vf_expr aux nl_expr le fait qu'on travaille sur une face :
    class nl_helper {
      void element_initialize (const This& obj, const geo_element& L, size_type loc_isid=-1) const {
        obj._nl_expr.evaluate (L, isid, obj._vector_nl_value_quad);
      }
    };
  pour la classe normal :
    field_expr_terminal_normal::evaluate (L, loc_isid, value) {
      if (loc_isid != -1) K=side(L,loc_isid); else K=L;
      puis inchange. 
    }
  pour la classe terminal_field: si on evalue un field uh qui est discontinu :
  on sait sur quelle face il se restreint :
    field_expr_terminal_field::evaluate (L, loc_isid, value) {
      if (loc_isid != -1) {
        _basis_on_quad.restrict_on_side (tilde_L, loc_isid);
      }
      for (q..) {
        general_field_evaluate (_uh, _basis_on_quad, tilde_L, _dis_idof, q, value[q]);
      }
    }
IMPLEMENTATION: bassis evaluation => basis_on_pointset.cc
  c'est la que se fait le coeur du travail :
    basis_on_pointset::restrict_on_side (tilde_L, loc_isid)
    => initialise

  a l'initialisation, on evalue une fois pour tte
  sur toutes les faces en transformant la quadrature via 
	tilde(i)_xq = Hi(hat_xq)
	tilde(i)_wq = ci*hat_wq
  avec
	ci = meas(tilde(i)_K)/meas(hat_K)

  puis :
    basis_on_pointset::evaluate (tilde_L, q)
  on se baladera dans la tranche [loc_isid*nq, (loc_isid+1)*nq[
  du coup, on positionne un pointeur de debut q_start  = loc_isid*nq
                               et une taille  q_size = nq
  si les faces sont differentes (tri,qua) dans un prisme, il faudra
  un tableau de pointeurs pour gerer cela :
		q_start [loc_nsid+1]
		q_size  [loc_isid] = q_start[loc_isid+1] - q_start[loc_isid] 

    basis_on_pointset::begin() { return _val[_curr_K_variant][_curr_q].begin() + q_start[_curr_K_variant][loc_isid]; }
    basis_on_pointset::begin() { return _val[_curr_K_variant][_curr_q].begin() + q_start[_curr_K_variant][loc_isid+1]; }

  et le tour est joue' !

PLAN DE DEVELOPPEMENT: 
 1) DG transport
    basis_on_pointset.cc 
    test.cc
    essais :
     lh = integrate(jump(v)*f);
     convect_dg2.cc
 2) DG diffusion : avec normale et gradient
    field_vf_expr.h 
    	class nl_helper 
    field_nl_expr_terminal.h
    	field_expr_terminal_normal::evaluate (L, loc_isid, value) 
    	field_expr_terminal_field ::evaluate (L, loc_isid, value)

*/ 

namespace details {
// ---------------------------------------------------------------------------
// 5.1 class dg
// ---------------------------------------------------------------------------
template<class Expr>
class field_expr_v2_variational_dg {
public:
// typedefs:

  typedef typename Expr::size_type            size_type;
  typedef typename Expr::memory_type          memory_type;
  typedef typename Expr::value_type           value_type;
  typedef typename Expr::scalar_type          scalar_type;
  typedef typename Expr::float_type           float_type;
  typedef typename Expr::space_type           space_type;
  typedef typename Expr::vf_tag_type          vf_tag_type;
  typedef typename details::dual_vf_tag<vf_tag_type>::type
                                              vf_dual_tag_type;
  typedef field_expr_v2_variational_dg<Expr>  self_type;
  typedef field_expr_v2_variational_dg<typename Expr::dual_self_type>  
                                              dual_self_type;

// alocators:

  field_expr_v2_variational_dg (const Expr& expr, const float_type& c0, const float_type& c1) 
    : _expr0(expr), 
      _expr1(expr),
      _c0(c0), 
      _c1(c1),
      _tilde0_L0(),
      _tilde1_L1(),
      _bgd_omega()
   {
   }

// accessors:

  const space_type& get_vf_space() const { return _expr0.get_vf_space(); }
  static const space_constant::valued_type valued_hint = Expr::valued_hint;
  space_constant::valued_type valued_tag() const { return _expr0.valued_tag(); } 
  size_type n_derivative() const { return _expr0.n_derivative(); }

// mutable modifiers:

  void initialize (const geo_basic<float_type,memory_type>& dom, const quadrature<float_type>& quad, bool ignore_sys_coord) const {
    _expr0.initialize (dom, quad, ignore_sys_coord);
    _expr1.initialize (dom, quad, ignore_sys_coord);
    _bgd_omega = get_vf_space().get_geo().get_background_geo();
    check_macro (_bgd_omega == dom.get_background_geo(),
        "discontinuous Galerkin: incompatible integration domain "<<dom.name() << " and test function based domain "
        << get_vf_space().get_geo().name());
  }
  void initialize (const band_basic<float_type,memory_type>& gh, const quadrature<float_type>& quad, bool ignore_sys_coord) const {  
    _expr0.initialize (gh, quad, ignore_sys_coord);
    _expr1.initialize (gh, quad, ignore_sys_coord);
    fatal_macro ("unsupported discontinuous Galerkin on a band"); // how to define background mesh _bgd_omega ?
  }
  void element_initialize (const geo_element& K) const;

  template<class ValueType>
  void basis_evaluate (const reference_element& hat_K, size_type q, std::vector<ValueType>& value) const;

  template<class ValueType>
  void valued_check() const {
    check_macro (get_vf_space().get_numbering().is_discontinuous(),
	"unexpected continuous test-function in space " << get_vf_space().stamp()
	<< " for jump or average operator (HINT: omit jump or average)");
    _expr0.template valued_check<ValueType>();
  }

protected:
// data:
  mutable Expr                              _expr0;
  mutable Expr                              _expr1;
  scalar_type                               _c0;
  scalar_type                               _c1;
  mutable reference_element                 _tilde0_L0;
  mutable reference_element                 _tilde1_L1;
  mutable geo_basic<float_type,memory_type> _bgd_omega;
};
template<class Expr> struct is_field_expr_v2_variational_arg    <field_expr_v2_variational_dg<Expr> > : std::true_type {};
// ---------------------------------------------------------------------------
// 5.2. basis_evaluate
// ---------------------------------------------------------------------------
template<class Expr>
template<class ValueType>
void
field_expr_v2_variational_dg<Expr>::basis_evaluate (
  const reference_element& hat_K, 
  size_type                q, 
  std::vector<ValueType>&  value) const
{
  size_type loc_ndof0 = _expr0.get_vf_space().get_constitution().loc_ndof (_tilde0_L0),
            loc_ndof1 = 0;
  if (_tilde1_L1.variant() != reference_element::max_variant) {
    loc_ndof1 = _expr1.get_vf_space().get_constitution().loc_ndof (_tilde1_L1);
  }
  check_macro (loc_ndof0+loc_ndof1 == value.size(),
	"unexpected value.size="<<value.size()<<": expect size="<<loc_ndof0+loc_ndof1
	<< " for ValueType="<<typename_macro(ValueType)<<" and valued=" << _expr1.get_vf_space().valued());
  std::vector<ValueType> value0 (loc_ndof0), 
                         value1 (loc_ndof1);
  _expr0.basis_evaluate (_tilde0_L0, q, value0);
  // average (i.e. _c0==0.5): fix it on the boundary where c0=1 : average(v)=v on the boundary
  Float c0 = (_tilde1_L1.variant() != reference_element::max_variant || _c0 != 0.5) ? _c0 : 1;
  for (size_type loc_idof = 0; loc_idof < loc_ndof0; ++loc_idof) {
    value[loc_idof] = c0*value0 [loc_idof];
  }
  if (_tilde1_L1.variant() != reference_element::max_variant) {
    _expr1.basis_evaluate (_tilde1_L1, q, value1);
    for (size_type loc_idof = 0; loc_idof < loc_ndof1; ++loc_idof) {
      value[loc_idof+loc_ndof0] = _c1*value1 [loc_idof];
    }
  }
} 
// ---------------------------------------------------------------------------
// 5.3. element_initialize
// ---------------------------------------------------------------------------
template<class Expr>
void
field_expr_v2_variational_dg<Expr>::element_initialize (const geo_element& K) const
{
  size_type L_map_d = K.dimension() + 1;
  size_type L_dis_ie0, L_dis_ie1;
  side_information_type sid0, sid1;

  L_dis_ie0 = K.master(0);
  L_dis_ie1 = K.master(1);
  check_macro (L_dis_ie0 != std::numeric_limits<size_type>::max(),
      "unexpected isolated mesh side K.dis_ie="<<K.dis_ie());
  if (L_dis_ie1 != std::numeric_limits<size_type>::max()) {
    // K is an internal side
    const geo_element& L0 = _bgd_omega.dis_get_geo_element (L_map_d, L_dis_ie0);
    const geo_element& L1 = _bgd_omega.dis_get_geo_element (L_map_d, L_dis_ie1);
    L0.get_side_informations (K, sid0);
    L1.get_side_informations (K, sid1);
    _tilde0_L0 = L0;
    _tilde1_L1 = L1;
    // methode "non-const": provoque une copie physique au 1er appel (c'est ce qu'il faut)
    _expr0.element_initialize_on_side (L0, sid0);
    _expr1.element_initialize_on_side (L1, sid1);
  } else {
    // K is a boundary side
    const geo_element& L0 = _bgd_omega.dis_get_geo_element (L_map_d, L_dis_ie0);
    L0.get_side_informations (K, sid0);
    _tilde0_L0 = L0;
    _tilde1_L1 = reference_element::max_variant;
    _expr0.element_initialize_on_side (L0, sid0);
  }
}

} // namespace details
// ---------------------------------------------------------------------------
// 5.4. user-level interface
// ---------------------------------------------------------------------------
#define _RHEOLEF_make_field_expr_v2_variational_dg(op,c0,c1) \
template<class Expr>					\
inline							\
typename						\
std::enable_if<						\
  details::is_field_expr_v2_variational_arg<Expr>::value \
 ,details::field_expr_v2_variational_dg<Expr>		\
>::type							\
op (const Expr& expr)					\
{							\
  return details::field_expr_v2_variational_dg <Expr> (expr, c0, c1); \
}

_RHEOLEF_make_field_expr_v2_variational_dg (jump,    1,    -1)
_RHEOLEF_make_field_expr_v2_variational_dg (average, 0.5,  0.5)
_RHEOLEF_make_field_expr_v2_variational_dg (inner,   1,    0)
_RHEOLEF_make_field_expr_v2_variational_dg (outer,   0,    1)
#undef _RHEOLEF_make_field_expr_v2_variational_dg

#ifdef TODO
#endif // TODO

} // namespace rheolef
#endif // _RHEOLEF_FIELD_EXPR_V2_VARIATIONAL_TERMINAL_H