librheolef-dev 6.5-1build1 / usr / include / rheolef / field_vf_assembly.h

#ifndef _RHEO_FIELD_VF_ASSEMBLY_H
#define _RHEO_FIELD_VF_ASSEMBLY_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
/// 
/// =========================================================================
#include "rheolef/field.h"
#include "rheolef/test.h"
#include "rheolef/quadrature.h"
/*
   let:
     l(v) = int_domain expr(v) dx
  
   The integrals are evaluated over each element K of the domain
   by using a quadrature formulae given by qopt
  
   expr(v) is a linear expression with respect to the test-function v

   The test-function v is replaced by each of the basis function of 
   the corresponding finite element space Xh: (phi_i), i=0..dim(Xh)-1
 
   The integrals over K are transformed on the reference element with
   the piola transformation:
     F : hat_K ---> K
         hat_x |--> x = F(hat_x)
 
   exemples:
   1) expr(v) = v
    int_K phi_i(x) dx 
     = int_{hat_K} hat_phi_i(hat_x) det(DF(hat_x)) d hat_x
     = sum_q hat_phi_i(hat_xq) det(DF(hat_xq)) hat_wq

    The value(q,i) = (hat_phi_i(hat_xq)) are evaluated on time for all over the 
    reference element hat_K and for the given quadrature formulae by:
  	expr.initialize (dom, quad);
    This expression is represented by the 'test' class (see test.h)

   2) expr(v) = f*v
    int_K f(x)*phi_i(x) dx 
     = int_{hat_K} f(F(hat_x)*hat_phi_i(hat_x) det(DF(hat_x)) d hat_x
     = sum_q f(F(hat_xq))*hat_phi_i(hat_xq) det(DF(hat_xq)) hat_wq
 
    On each K, the computation needs two imbricated loops in (q,i).
    The expresion is represented by a tree at compile-time.
    The first subexpr 'f' is represented by field_vf_expr_terminal<f> : it evaluates : 
      value1(q,i) = f(F(hat_xq)) : the values depends only of q and are independent of i.
    They could be evaluated juste before the 'i' loop.
    As the field_vf_expr_terminal<> is general and can handle also more complex expressions
    involving fields with various approximations, all the values on K are evaluated 
    in a specific 'q' loop by
      subexpr1.element_initialize (K);
    The second subexpr 'phi_i' is represdented by the 'test' class (see before).
      value2(q,i) = (hat_phi_i(hat_xq))
    The '*' performs on the fly the product
      value(q,i) = value1(q,i)*value2(q,i)
  
   3) expr(v) = dot(f,grad(v)) dx
    The 'f' function is here vector-valued.
    The expresion is represented by a tree at compile-time :
      dot -> f 
          -> grad -> v
    The 'grad' node returns  
      value(q,i) = trans(inv(DF(hat_wq))*grad_phi_i(hat_xq) that is vector-valued
    The grad_phi values are obtained by a grad_value(q,i) method on the 'test' class.
    The 'dot' performs on the fly the product
      value(q,i) = ot (value1(q,i), value2(q,i))

   This approch generilize for an expression tree.
*/
namespace rheolef {

// ====================================================================
// compute dis_idof on K
// ====================================================================
template <class T, class M>
void
assembly_dis_idof (
  const space_basic<T,M>&              X,
  const geo_basic<T,M>&                dom,
  const geo_element&                   K,
  std::vector<geo_element::size_type>& dis_idx);

// ====================================================================
// common implementation for integration on a band or an usual domain
// ====================================================================
template <class T, class M>
template <class Expr>
void
field_basic<T,M>::assembly_internal (
    const geo_basic<T,M>&         dom,
    const geo_basic<T,M>&         band,
    const band_basic<T,M>&        gh,
    const Expr&                   expr,
    const quadrature_option_type& qopt,
    bool                          is_on_band)
{
  // ------------------------------------
  // 0) initialize the loop
  // ------------------------------------
  const bool ignore_sys_coord = false; // TODO: available as an option via fopt
  const space_basic<T,M>& Xh = expr.get_vf_space();
  const geo_basic<T,M>& omega     = Xh.get_geo();
  const geo_basic<T,M>& bgd_omega = Xh.get_geo().get_background_geo();
  check_macro (band.get_background_geo() == bgd_omega,
	"assembly: incompatible integration domain "<<dom.name() << " and test function based domain "
	<< Xh.get_geo().name());
  resize (Xh, 0);

  quadrature<T> quad;
  if (qopt.get_order() != 0 && qopt.get_order() != std::numeric_limits<quadrature_option_type::size_type>::max()) {
    quad.set_order (qopt.get_order());
  } else {
    size_type k = Xh.degree();
    if (dom.coordinate_system() != space_constant::cartesian) k++; // multiplies by a 'r' weight
    size_type quad_order = 2*k+1; // see Raviart & Thomas, Masson, 1988, page 130, theorem 5.3.1
    quad.set_order (quad_order);
  }
  quad.set_family (qopt.get_family());
  if (is_on_band) {
    expr.initialize (gh, quad, ignore_sys_coord);
  } else {
    expr.initialize (dom, quad, ignore_sys_coord);
  }
  expr.template valued_check<T>();
  basis_on_pointset<T> piola_on_quad;
  piola_on_quad.set (quad, Xh.get_geo().get_piola_basis());
  vec<T,M>& u = set_u();
  vec<T,M>& b = set_b();
  std::vector<size_type> dis_idx;
  std::vector<T> lk, value;
  tensor_basic<T> DF;
  std::vector<size_type> dis_inod_K;
  size_type d     = dom.dimension();
  size_type map_d = dom.map_dimension();
  bool is_dg = (!Xh.get_numbering().is_continuous() && map_d < omega.map_dimension());
  // ------------------------------------
  // 1) loop on elements
  // ------------------------------------
  for (size_type ie = 0, ne = dom.size(map_d); ie < ne; ie++) {
    const geo_element& K = dom.get_geo_element (map_d, ie);
    // ------------------------------------
    // 1.1) locally compute dofs
    // ------------------------------------
    if (!is_on_band) {
      assembly_dis_idof (Xh, dom, K, dis_idx);
    } else { // on a band
      size_type L_ie = gh.sid_ie2bnd_ie (ie);
      const geo_element& L = band [L_ie];
      Xh.dis_idof (L, dis_idx);
    }
    // ------------------------------------
    // 1.2) locally compute values
    // ------------------------------------
    // locally evaluate int_K expr dx
    // with loop on a quadrature formulae 
    // and accumulate in lk
    value.resize (dis_idx.size());
    lk.resize    (dis_idx.size());
    fill (lk.begin(), lk.end(), T(0));
    expr.element_initialize (K);
    reference_element hat_K = K.variant();
    dom.dis_inod (K, dis_inod_K);
    typename quadrature<T>::const_iterator
      first_quad = quad.begin(hat_K),
      last_quad  = quad.end  (hat_K);
    for (size_type q = 0; first_quad != last_quad; ++first_quad, ++q) {
      jacobian_piola_transformation (dom, piola_on_quad, K, dis_inod_K, q, DF);
      T det_DF = det_jacobian_piola_transformation (DF, d, map_d);
      T wq = 1;
      if (!ignore_sys_coord) {
        wq *= weight_coordinate_system (dom, piola_on_quad, K, dis_inod_K, q);
      }
      wq *= det_DF;
      wq *= (*first_quad).w;
      expr.basis_evaluate (hat_K, q, value);
      for (size_type loc_idof = 0, loc_ndof = value.size(); loc_idof < loc_ndof; ++loc_idof) {
        lk [loc_idof] += value[loc_idof]*wq;
      }
    }
    // -----------------------------------------
    // 1.3) assembly local lk in global field lh
    // -----------------------------------------
    check_macro (dis_idx.size() == lk.size(),
      "incompatible sizes dis_idx.size="<<dis_idx.size()<<", lk.size="<<lk.size());
    for (size_type loc_idof = 0, loc_ndof = lk.size(); loc_idof < loc_ndof; loc_idof++) {
      const T& value = lk [loc_idof];
      size_type dis_idof = dis_idx [loc_idof];
      size_type dis_iub  = Xh.dis_iub(dis_idof);
      if (Xh.dis_is_blocked(dis_idof)) b.dis_entry(dis_iub) += value;
      else                             u.dis_entry(dis_iub) += value;
    }
  }
  // -----------------------------------------
  // 2) finalize the assembly process
  // -----------------------------------------
  u.dis_entry_assembly(set_add_op<T,T>());
  b.dis_entry_assembly(set_add_op<T,T>());
}
template <class T, class M>
template <class Expr>
inline
void
field_basic<T,M>::assembly (
    const geo_basic<T,M>&         dom,
    const Expr&                   expr,
    const quadrature_option_type& qopt)
{
  assembly_internal (dom, dom, band_basic<T,M>(), expr, qopt, false);
}
template <class T, class M>
template <class Expr>
inline
void
field_basic<T,M>::assembly (
    const band_basic<T,M>&        gh,
    const Expr&                   expr,
    const quadrature_option_type& qopt)
{
  assembly_internal (gh.level_set(), gh.band(), gh, expr, qopt, true);
}

}// namespace rheolef
#endif // _RHEO_FIELD_VF_ASSEMBLY_H