Loading...
Searching...
No Matches
form_lazy_convert.h
Go to the documentation of this file.
#ifndef _RHEO_FORM_LAZY_CONVERT_H
#define _RHEO_FORM_LAZY_CONVERT_H
// convert to form from an un-assembled lazy_form expression
// AUTHOR: Pierre.Saramito@imag.fr
// DATE: 30 march 2020
/*
a(u,v) = int_omega expr dx
= sum_K int_K expr dx
The integrals are evaluated over each element K in omega
and expr is a bilinear expression, as returned by form_lazy.
u & v are trial and test functions, involved in expr.
*/
#include "rheolef/form.h"
namespace rheolef {
namespace details {
// 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);
} // namespace details
// ====================================================================
// common implementation for integration on a band or an usual domain
// ====================================================================
template <class T, class M>
template <class Expr, class Sfinae>
void
form_basic<T,M>::convert_from_form_lazy (const Expr& expr0)
{
Expr expr = expr0;
// ----------------------------------------
// 0) init assembly loop
// ----------------------------------------
_X = expr.get_trial_space();
_Y = expr.get_test_space();
geo_basic<T,M> omega = expr.get_geo();
expr.initialize (omega);
bool is_on_band = expr.is_on_band();
band_basic<T,M> gh;
if (is_on_band) gh = expr.get_band();
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 = gh.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 = details::norm_max (ak);
T eps_ak_max = eps*ak_max;
if (is_sym) is_sym = details::check_is_symmetric (ak, eps_ak_max);
// ----------------------------------------
// 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);
// filter too small values
// * 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);
}
}// namespace rheolef
#endif // _RHEO_FORM_LAZY_CONVERT_H
size_type
field::size_type size_type
Definition branch.cc:430
gh
field gh(Float epsilon, Float t, const field &uh, const test &v)
Definition burgers_diffusion_operators.icc:37
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
see the csr page for the full documentation
Definition csr.h:317
rheolef::form_basic::convert_from_form_lazy
void convert_from_form_lazy(const Expr &expr)
Definition form_lazy_convert.h:52
rheolef::form_basic::size_type
csr< T, M >::size_type size_type
Definition form.h:150
rheolef::geo_basic
generic mesh with rerefence counting
Definition geo.h:1089
rheolef::geo_element
see the geo_element page for the full documentation
Definition geo_element.h:102
T
Expr1::float_type T
Definition field_expr.h:230
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
rheolef::details::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::details::norm_max
T norm_max(Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m)
Definition form_vf_assembly.cc:33
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
rheolef::is_distributed
Definition communicator.h:70