1#ifndef _RHEO_FORM_LAZY_CONVERT_H
2#define _RHEO_FORM_LAZY_CONVERT_H
36#include "rheolef/form.h"
42template <
class T>
T norm_max (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
m);
43template <
class T>
bool check_is_symmetric (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
m,
const T& tol_m_max);
49template <
class T,
class M>
50template <
class Expr,
class Sfinae>
58 _X =
expr.get_trial_space();
59 _Y =
expr.get_test_space();
61 expr.initialize (omega);
62 bool is_on_band =
expr.is_on_band();
64 if (is_on_band)
gh =
expr.get_band();
65 asr<T,M> auu (_Y.iu_ownership(), _X.iu_ownership()),
66 aub (_Y.iu_ownership(), _X.ib_ownership()),
67 abu (_Y.ib_ownership(), _X.iu_ownership()),
68 abb (_Y.ib_ownership(), _X.ib_ownership());
69 std::vector<size_type> dis_idy, dis_jdx;
70 Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> ak;
72 if (_X.get_constitution().is_discontinuous()) _X.get_constitution().neighbour_guard();
73 if (_Y.get_constitution().is_discontinuous()) _Y.get_constitution().neighbour_guard();
76 for (
size_type ie = 0, ne = omega.size(map_d); ie < ne; ie++) {
80 const geo_element& K = omega.get_geo_element (map_d, ie);
82 _X.get_constitution().assembly_dis_idof (omega, K, dis_jdx);
83 _Y.get_constitution().assembly_dis_idof (omega, K, dis_idy);
87 _X.dis_idof (
L, dis_jdx);
88 _Y.dis_idof (
L, dis_idy);
90 expr.evaluate (omega, K, ak);
95 T eps_ak_max = eps*ak_max;
101 "invalid sizes ak("<<ak.rows()<<
","<<ak.cols()
102 <<
") with dis_idy("<<dis_idy.size()<<
") and dis_jdx("<<dis_jdx.size()<<
")");
103 for (
size_type loc_idof = 0,
ny = ak.rows(); loc_idof <
ny; loc_idof++) {
104 for (
size_type loc_jdof = 0,
nx = ak.cols(); loc_jdof <
nx; loc_jdof++) {
106 const T&
value = ak (loc_idof, loc_jdof);
116 if (fabs(
value) < eps_ak_max)
continue;
121 size_type dis_jub = _X.dis_idof2dis_iub (dis_jdof);
124 if (_X.dis_is_blocked(dis_jdof)) abb.
dis_entry (dis_iub, dis_jub) +=
value;
125 else abu.dis_entry (dis_iub, dis_jub) +=
value;
127 if (_X.dis_is_blocked(dis_jdof)) aub.dis_entry (dis_iub, dis_jub) +=
value;
128 else auu.dis_entry (dis_iub, dis_jub) +=
value;
137 auu.dis_entry_assembly();
138 aub.dis_entry_assembly();
139 abu.dis_entry_assembly();
152 _uu.set_pattern_dimension (map_d);
153 _ub.set_pattern_dimension (map_d);
154 _bu.set_pattern_dimension (map_d);
155 _bb.set_pattern_dimension (map_d);
164#ifdef _RHEOLEF_HAVE_MPI
166 is_sym = mpi::all_reduce (omega.comm(),
size_type(is_sym), mpi::minimum<size_type>());
169 _uu.set_symmetry (is_sym);
170 _bb.set_symmetry (is_sym);
172 _uu.set_definite_positive (is_sym);
173 _bb.set_definite_positive (is_sym);
field gh(Float epsilon, Float t, const field &uh, const test &v)
void dis_entry_assembly()
dis_reference dis_entry(size_type dis_i, size_type dis_j)
see the geo_element page for the full documentation
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
bool check_is_symmetric(Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m, const T &tol_m_max)
T norm_max(Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m)
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)
This file is part of Rheolef.