form_vf_assembly.h
Go to the documentation of this file.
#ifndef _RHEO_FORM_VF_ASSEMBLY_H
#define _RHEO_FORM_VF_ASSEMBLY_H
#include "rheolef/form.h"
#include "rheolef/test.h"
#include "rheolef/quadrature.h"
#include "rheolef/field_vf_assembly.h" // for dg_dis_idof()
#include "rheolef/form_expr_quadrature.h"
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 iopt
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 (omega, 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 generalizes for an expression tree.
*/
// external utilities:
template <class T> T norm_max (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m);
template <class T> bool check_is_symmetric (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m, const T& tol_m_max);
template <class T> void local_lump (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m);
template <class T> void local_invert (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& 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>& omega,
const geo_basic<T,M>& band,
const band_basic<T,M>& gh,
const Expr& expr,
const integrate_option& iopt,
bool is_on_band)
{
// ----------------------------------------
// 0) init assembly loop
// ----------------------------------------
_X = expr.get_trial_space();
_Y = expr.get_test_space();
if (is_on_band) {
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());
}
size_type n_derivative = expr.n_derivative();
if (iopt.invert) check_macro (
_X.get_constitution().have_compact_support_inside_element() &&
_Y.get_constitution().have_compact_support_inside_element(),
"local inversion requires compact support inside elements (e.g. discontinuous or bubble)");
if (iopt.lump) check_macro (n_derivative == 0,
"local mass lumping requires no derivative operators");
iopt._is_on_interface = false;
iopt._is_inside_on_local_sides = false;
if (!is_on_band) {
expr.initialize (omega, iopt);
} else {
expr.initialize (gh, iopt);
}
expr.template valued_check<T>();
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;
Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> ak;
size_type map_d = omega.map_dimension();
if (_X.get_constitution().is_discontinuous()) _X.get_constitution().neighbour_guard();
if (_Y.get_constitution().is_discontinuous()) _Y.get_constitution().neighbour_guard();
bool is_sym = true;
const T eps = 1e3*std::numeric_limits<T>::epsilon();
for (size_type ie = 0, ne = omega.size(map_d); ie < ne; ie++) {
// ----------------------------------------
// 1) compute local form ak
// ----------------------------------------
const geo_element& K = omega.get_geo_element (map_d, ie);
if (! is_on_band) {
_X.get_constitution().assembly_dis_idof (omega, K, dis_jdx);
_Y.get_constitution().assembly_dis_idof (omega, 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);
}
expr.evaluate (omega, K, ak);
// ----------------------------------------
// 2) optional local post-traitement
// ----------------------------------------
T ak_max = norm_max (ak);
T eps_ak_max = eps*ak_max;
if (is_sym) is_sym = check_is_symmetric (ak, eps_ak_max);
if (iopt.lump ) local_lump (ak);
if (iopt.invert) local_invert (ak, iopt.lump);
// ----------------------------------------
// 3) assembly local ak in global form a
// ----------------------------------------
check_macro (size_type(ak.rows()) == dis_idy.size() && size_type(ak.cols()) == dis_jdx.size(),
"invalid sizes ak("<<ak.rows()<<","<<ak.cols()
<<") with dis_idy("<<dis_idy.size()<<") and dis_jdx("<<dis_jdx.size()<<")");
for (size_type loc_idof = 0, ny = ak.rows(); loc_idof < ny; loc_idof++) {
for (size_type loc_jdof = 0, nx = ak.cols(); 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
// useful 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 (fabs(value) < eps_ak_max) continue;
size_type dis_idof = dis_idy [loc_idof];
size_type dis_jdof = dis_jdx [loc_jdof];
size_type dis_iub = _Y.dis_idof2dis_iub (dis_idof);
size_type dis_jub = _X.dis_idof2dis_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 (omega.comm().size() > 1 && is_distributed<M>::value) {
is_sym = mpi::all_reduce (omega.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>& omega,
const Expr& expr,
const integrate_option& iopt)
{
assembly_internal (omega, omega, band_basic<T,M>(), expr, iopt, 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 integrate_option& iopt)
{
assembly_internal (gh.level_set(), gh.band(), gh, expr, iopt, true);
}
}// namespace rheolef
#endif // _RHEO_FORM_VF_ASSEMBLY_H
gh
field gh(Float epsilon, Float t, const field &uh, const test &v)
Definition: burgers_diffusion_operators.icc:37
band
see the band page for the full documentation
rheolef::asr
Definition: asr.h:49
rheolef::asr::dis_entry_assembly
void dis_entry_assembly()
Definition: asr.h:112
rheolef::asr::dis_entry
dis_reference dis_entry(size_type dis_i, size_type dis_j)
Definition: asr.h:193
rheolef::band_basic
Definition: band.h:95
rheolef::csr< T, M >
rheolef::form_basic::size_type
csr< T, M >::size_type size_type
Definition: form.h:148
rheolef::form_basic::assembly_internal
void assembly_internal(const geo_basic< T, M > &dom, const geo_basic< T, M > &band, const band_basic< T, M > &gh, const Expr &expr, const integrate_option &fopt, bool is_on_band)
Definition: form_vf_assembly.h:87
rheolef::form_basic::assembly
void assembly(const geo_basic< T, M > &domain, const Expr &expr, const integrate_option &fopt)
Definition: form_vf_assembly.h:241
rheolef::geo_basic< T, M >
rheolef::geo_element
see the geo_element page for the full documentation
Definition: geo_element.h:102
rheolef::integrate_option
see the integrate_option page for the full documentation
Definition: integrate_option.h:125
rheolef::integrate_option::lump
bool lump
Definition: integrate_option.h:168
rheolef::integrate_option::invert
bool invert
Definition: integrate_option.h:168
rheolef::integrate_option::_is_on_interface
bool _is_on_interface
Definition: integrate_option.h:180
rheolef::integrate_option::_is_inside_on_local_sides
bool _is_inside_on_local_sides
Definition: integrate_option.h:180
rheolef::size_type
size_t size_type
Definition: basis_get.cc:76
rheolef::value
rheolef::std value
T
Expr1::float_type T
Definition: field_expr.h:261
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
mkgeo_ball.expr
expr
Definition: mkgeo_ball.sh:361
mkgeo_contraction.ny
ny
Definition: mkgeo_contraction.sh:150
mkgeo_contraction.nx
string nx
Definition: mkgeo_contraction.sh:124
mkgeo_obstacle.L
string L
Definition: mkgeo_obstacle.sh:133
mkgeo_sector.m
m
Definition: mkgeo_sector.sh:118
rheolef::space_numbering::dis_idof
void dis_idof(const basis_basic< T > &b, const geo_size &gs, const geo_element &K, typename std::vector< size_type >::iterator dis_idof_tab)
Definition: space_numbering.cc:187
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::local_lump
void local_lump(Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m)
Definition: form_vf_assembly.cc:61
rheolef::local_invert
void local_invert(Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m, bool is_diag)
Definition: form_vf_assembly.cc:77
rheolef::check_is_symmetric
bool check_is_symmetric(Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m, const T &tol_m_max)
Definition: form_vf_assembly.cc:46
rheolef::norm_max
T norm_max(Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m)
Definition: form_vf_assembly.cc:33
rheolef::is_distributed
Definition: distributed.h:36
epsilon
Float epsilon
Definition: transmission_error.cc:25