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

#ifndef _RHEO_FORM_VF_ASSEMBLY_H
#define _RHEO_FORM_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/form.h"
#include "rheolef/test.h"
#include "rheolef/quadrature.h"
#include "rheolef/field_vf_assembly.h" // for dg_dis_idof()
#include "rheolef/ublas-io.h" // for debug
namespace rheolef {
/*
   let:
     a(u,v) = int_domain expr(u,v) dx
  
   The integrals are evaluated over each element K of the domain
   by using a quadrature formulae given by fopt
  
   expr(u,v) is a bilinear expression with respect to the
   trial and test functions u and v

   The trial function u is replaced by each of the basis function of 
   the corresponding finite element space Xh: (phi_j), j=0..dim(Xh)-1
 
   The test function v is replaced by each of the basis function of 
   the corresponding finite element space Yh: (psi_i), i=0..dim(Yh)-1
 
   The integrals over the domain omega is the sum of integrals over K.

   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) = u*v
    int_K phi_j(x)*psi_i(x) dx 
     = int_{hat_K} hat_phi_j(hat_x)*hat_psi_i(hat_x) det(DF(hat_x)) d hat_x
     = sum_q hat_phi_j(hat_xq)*hat_psi_i(hat_xq) det(DF(hat_xq)) hat_wq

    The value(q,i,j) = (hat_phi_j(hat_xq)*hat_psi_i(hat_xq))
    refers to basis values on the reference element.
    There 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)

   3) expr(v) = dot(grad(u),grad(v)) dx
    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,j) = dot (value1(q,i), value2(q,j))

   This approch generalize for an expression tree.
*/

// external utilities:
template <class T> bool is_symmetric (ublas::matrix<T>& m, const T& tol);
template <class T> void local_lump   (ublas::matrix<T>& m);
template <class T> void local_invert (ublas::matrix<T>& m, bool is_diag);

// ====================================================================
// common implementation for integration on a band or an usual domain
// ====================================================================
template <class T, class M>
template <class Expr>
void
form_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 form_option_type&       fopt,
   bool                          is_on_band)
{
  // ----------------------------------------
  // 0) init assembly loop
  // ----------------------------------------
  _X = expr.get_trial_space();
  _Y = expr.get_test_space();
  check_macro (band.get_background_geo() == _X.get_geo().get_background_geo(),
	"assembly: incompatible integration domain "<<band.name() << " and trial function based domain "
	<< _X.get_geo().name());
  check_macro (band.get_background_geo() == _Y.get_geo().get_background_geo(),
	"assembly: incompatible integration domain "<<band.name() << " and test function based domain "
	<< _Y.get_geo().name());

  quadrature<T> quad;
  size_type n_derivative = expr.n_derivative();
  if (fopt.get_order() != std::numeric_limits<quadrature_option_type::size_type>::max()) {
    quad.set_order (fopt.get_order());
  } else {
    size_type k1 = _X.degree();
    size_type k2 = _Y.degree();
    size_type quad_order = k1 + k2 + 1;
    if (dom.get_background_geo().sizes().ownership_by_variant[reference_element::q].dis_size() != 0 ||
        dom.get_background_geo().sizes().ownership_by_variant[reference_element::P].dis_size() != 0 ||
        dom.get_background_geo().sizes().ownership_by_variant[reference_element::H].dis_size() != 0) {
          // integrate exactly ??
          quad_order += 2;
    }
    if (dom.coordinate_system() != space_constant::cartesian) quad_order++; // multiplies by a 'r' weight
    if (quad_order >= n_derivative) quad_order -= n_derivative;
    quad.set_order (quad_order);
  }
  quad.set_family (fopt.get_family());
  trace_macro ("quadrature : " << quad.get_family_name() << " order " << quad.get_order());

  if (fopt.invert) check_macro (
	_X.get_numbering().has_compact_support_inside_element() &&
	_Y.get_numbering().has_compact_support_inside_element(),
	"local inversion requires compact support inside elements (e.g. discontinuous or bubble)");
  if (fopt.lump) check_macro (n_derivative == 0,
	"local mass lumping requires no derivative operators");

  if (!is_on_band) {
    expr.initialize (dom, quad, fopt.ignore_sys_coord);
  } else {
    expr.initialize (gh,  quad, fopt.ignore_sys_coord);
  }
  expr.template valued_check<T>();
  basis_on_pointset<T> piola_on_quad;
  piola_on_quad.set (quad, _X.get_geo().get_piola_basis());
  asr<T,M> auu (_Y.iu_ownership(), _X.iu_ownership()),
           aub (_Y.iu_ownership(), _X.ib_ownership()),
           abu (_Y.ib_ownership(), _X.iu_ownership()),
           abb (_Y.ib_ownership(), _X.ib_ownership());
  std::vector<size_type> dis_idy, dis_jdx, dis_inod_K;
  ublas::matrix<T> ak, value;
  value.clear();
  size_type     d = dom.dimension();
  size_type map_d = dom.map_dimension();
  tensor_basic<T> DF;
  bool X_is_dg = (!_X.get_numbering().is_continuous() && map_d < _X.get_geo().map_dimension());
  bool Y_is_dg = (!_Y.get_numbering().is_continuous() && map_d < _Y.get_geo().map_dimension());
  if (X_is_dg) _X.get_geo().neighbour_guard();
  if (Y_is_dg) _Y.get_geo().neighbour_guard();
  bool is_sym = true;
  T sym_tol = 1e3*std::numeric_limits<T>::epsilon();
  for (size_type ie = 0, ne = dom.size(map_d); ie < ne; ie++) {
    // ----------------------------------------
    // 1) compute local form ak
    // ----------------------------------------
    const geo_element& K = dom.get_geo_element (map_d, ie);
    expr.element_initialize (K);
    if (! is_on_band) {
      assembly_dis_idof (_X, dom, K, dis_jdx);
      assembly_dis_idof (_Y, dom, K, dis_idy);
    } else {
      size_type L_ie = gh.sid_ie2bnd_ie (ie);
      const geo_element& L = band [L_ie];
      _X.dis_idof (L, dis_jdx);
      _Y.dis_idof (L, dis_idy);
    }
    ak.resize    (dis_idy.size(), dis_jdx.size());
    value.resize (dis_idy.size(), dis_jdx.size());
    ak.clear();
    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 (! fopt.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.size1(); loc_idof < loc_ndof; ++loc_idof) {
      for (size_type loc_jdof = 0, loc_mdof = value.size2(); loc_jdof < loc_mdof; ++loc_jdof) {
        ak (loc_idof,loc_jdof) += value (loc_idof,loc_jdof)*wq;
      }}
    }
#ifdef TO_CLEAN
    std::cerr << "ak=" << ak;
#endif // TO_CLEAN
    // ----------------------------------------
    // 2) optional local post-traitement
    // ----------------------------------------
    if (is_sym) is_sym = is_symmetric (ak, sym_tol);
    if (fopt.lump  ) local_lump   (ak); 
    if (fopt.invert) local_invert (ak, fopt.lump); 
    // ----------------------------------------
    // 3) assembly local ak in global form a
    // ----------------------------------------
    for (size_type loc_idof = 0, ny = ak.size1(); loc_idof < ny; loc_idof++) {
    for (size_type loc_jdof = 0, nx = ak.size2(); loc_jdof < nx; loc_jdof++) {

      const T& value = ak (loc_idof, loc_jdof);
      // reason to perform:
      // - efficient : lumped mass, structured meshes => sparsity increases
      // reason to avoid:
      // - conserve the sparsity pattern, even with some zero coefs
      //   usefull when dealing with solver::update_values()
      // - also solver_pastix: assume sparsity pattern symmetry
      //   and failed when a_ij=0 (skipped) and a_ji=1e-15 (conserved) i.e. non-sym pattern
      //   note: this actual pastix wrapper limitation could be suppressed
      if (1+value == 1) continue;

      size_type dis_idof = dis_idy [loc_idof];
      size_type dis_jdof = dis_jdx [loc_jdof];

      size_type dis_iub = _Y.dis_iub (dis_idof);
      size_type dis_jub = _X.dis_iub (dis_jdof);

      if   (_Y.dis_is_blocked(dis_idof))
        if (_X.dis_is_blocked(dis_jdof)) abb.dis_entry (dis_iub, dis_jub) += value;
        else                             abu.dis_entry (dis_iub, dis_jub) += value;
      else 
	if (_X.dis_is_blocked(dis_jdof)) aub.dis_entry (dis_iub, dis_jub) += value;
        else                             auu.dis_entry (dis_iub, dis_jub) += value;
    }}
  }
  // ----------------------------------------
  // 4) finalize the assembly process
  // ----------------------------------------
  //
  // since all is local, axx.dis_entry_assembly() compute only axx.dis_nnz
  //
  auu.dis_entry_assembly();
  aub.dis_entry_assembly();
  abu.dis_entry_assembly();
  abb.dis_entry_assembly();
  //
  // convert dynamic matrix asr to fixed-size one csr
  //
  _uu = csr<T,M>(auu);
  _ub = csr<T,M>(aub);
  _bu = csr<T,M>(abu);
  _bb = csr<T,M>(abb);
  //
  // set pattern dimension to uu:
  // => used by solvers, for efficiency: direct(d<3) or iterative(d=3)
  //
  _uu.set_pattern_dimension (map_d);
  _ub.set_pattern_dimension (map_d);
  _bu.set_pattern_dimension (map_d);
  _bb.set_pattern_dimension (map_d);
  //
  // symmetry is used by solvers, for efficiency: LDL^t or LU, CG or GMRES
  //
  // Implementation note: cannot be set at compile time
  // ex: expression=(eta*u)*v is structurally unsym, but numerical sym
  //     expression=(eta_h*grad(u))*(nu_h*grad(v)) is structurally sym,
  //    but numerical unsym when eta and nu are different tensors
  // So, test it numerically, at element level:
#ifdef _RHEOLEF_HAVE_MPI
  if (dom.comm().size() > 1 && is_distributed<M>::value) {
    is_sym = mpi::all_reduce (dom.comm(), size_type(is_sym), mpi::minimum<size_type>());
  }
#endif // _RHEOLEF_HAVE_MPI
  _uu.set_symmetry (is_sym);
  _bb.set_symmetry (is_sym);
  // when sym, the main matrix is set definite and positive by default
  _uu.set_definite_positive (is_sym);
  _bb.set_definite_positive (is_sym);
}

template <class T, class M>
template <class Expr>
inline
void
form_basic<T,M>::assembly (
    const geo_basic<T,M>&         dom,
    const Expr&                   expr,
    const form_option_type&       fopt)
{
  assembly_internal (dom, dom, band_basic<T,M>(), expr, fopt, false);
}
template <class T, class M>
template <class Expr>
inline
void
form_basic<T,M>::assembly (
    const band_basic<T,M>&        gh,
    const Expr&                   expr,
    const form_option_type&       fopt)
{
  assembly_internal (gh.level_set(), gh.band(), gh, expr, fopt, true);
}

}// namespace rheolef
#endif // _RHEO_FORM_VF_ASSEMBLY_H