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

#ifndef _RHEOLEF_FIELD_VF_EXPR_DG_H
#define _RHEOLEF_FIELD_VF_EXPR_DG_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_vf_expr_dg: 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)

*/ 
#include "rheolef/field_vf_expr.h"
namespace rheolef {

// ---------------------------------------------------------------------------
// class dg
// ---------------------------------------------------------------------------
template<class Expr, class VfTag = typename Expr::vf_tag_type>
class field_vf_expr_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          VfTag                      vf_tag_type;
  typedef typename details::dual_vf_tag<vf_tag_type>::type
                                              vf_dual_tag_type;
  typedef field_vf_expr<Expr,VfTag>           self_type;
  typedef field_vf_expr<typename Expr::dual_self_type,vf_dual_tag_type>
                                              dual_self_type;

// alocators:

  field_vf_expr_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.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;
};
// ---------------------------------------------------------------------------
// basis_evaluate
// ---------------------------------------------------------------------------
template<class Expr, class VfTag>
template<class ValueType>
void
field_vf_expr_dg<Expr,VfTag>::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];
    }
  }
} 
// ---------------------------------------------------------------------------
// element_initialize
// ---------------------------------------------------------------------------
template<class Expr, class VfTag>
void
field_vf_expr_dg<Expr,VfTag>::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="<<K);
  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 rheolef
#endif // _RHEOLEF_FIELD_VF_EXPR_DG_H