1 #ifndef _RHEO_FORM_VF_ASSEMBLY_H
2 #define _RHEO_FORM_VF_ASSEMBLY_H
23 #include "rheolef/form.h"
24 #include "rheolef/test.h"
25 #include "rheolef/quadrature.h"
26 #include "rheolef/field_vf_assembly.h"
27 #include "rheolef/form_expr_quadrature.h"
77 template <
class T>
T norm_max (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
m);
78 template <
class T>
bool check_is_symmetric (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
m,
const T& tol_m_max);
79 template <
class T>
void local_lump (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
m);
80 template <
class T>
void local_invert (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
m,
bool is_diag);
86 template <
class T,
class M>
101 _X =
expr.get_trial_space();
102 _Y =
expr.get_test_space();
104 check_macro (
band.get_background_geo() == _X.get_geo().get_background_geo(),
105 "do_integratessembly: incompatible integration domain "<<
band.name() <<
" and trial function based domain "
106 << _X.get_geo().name());
107 check_macro (
band.get_background_geo() == _Y.get_geo().get_background_geo(),
108 "do_integratessembly: incompatible integration domain "<<
band.name() <<
" and test function based domain "
109 << _Y.get_geo().name());
114 _X.get_constitution().have_compact_support_inside_element() &&
115 _Y.get_constitution().have_compact_support_inside_element(),
116 "local inversion requires compact support inside elements (e.g. discontinuous or bubble)");
118 "local mass lumping requires no derivative operators");
123 expr.initialize (omega, iopt);
125 expr.initialize (
gh, iopt);
127 expr.template valued_check<T>();
128 asr<T,M> auu (_Y.iu_ownership(), _X.iu_ownership()),
129 aub (_Y.iu_ownership(), _X.ib_ownership()),
130 abu (_Y.ib_ownership(), _X.iu_ownership()),
131 abb (_Y.ib_ownership(), _X.ib_ownership());
132 std::vector<size_type> dis_idy, dis_jdx;
133 Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> ak;
135 if (_X.get_constitution().is_discontinuous()) _X.get_constitution().neighbour_guard();
136 if (_Y.get_constitution().is_discontinuous()) _Y.get_constitution().neighbour_guard();
139 for (
size_type ie = 0, ne = omega.size(map_d); ie < ne; ie++) {
143 const geo_element& K = omega.get_geo_element (map_d, ie);
145 _X.get_constitution().assembly_dis_idof (omega, K, dis_jdx);
146 _Y.get_constitution().assembly_dis_idof (omega, K, dis_idy);
150 _X.dis_idof (
L, dis_jdx);
151 _Y.dis_idof (
L, dis_idy);
153 expr.evaluate (omega, K, ak);
158 T eps_ak_max = eps*ak_max;
166 "invalid sizes ak("<<ak.rows()<<
","<<ak.cols()
167 <<
") with dis_idy("<<dis_idy.size()<<
") and dis_jdx("<<dis_jdx.size()<<
")");
168 for (
size_type loc_idof = 0,
ny = ak.rows(); loc_idof <
ny; loc_idof++) {
169 for (
size_type loc_jdof = 0,
nx = ak.cols(); loc_jdof <
nx; loc_jdof++) {
171 const T&
value = ak (loc_idof, loc_jdof);
180 if (fabs(
value) < eps_ak_max)
continue;
185 size_type dis_jub = _X.dis_idof2dis_iub (dis_jdof);
188 if (_X.dis_is_blocked(dis_jdof)) abb.
dis_entry (dis_iub, dis_jub) +=
value;
189 else abu.dis_entry (dis_iub, dis_jub) +=
value;
191 if (_X.dis_is_blocked(dis_jdof)) aub.dis_entry (dis_iub, dis_jub) +=
value;
192 else auu.dis_entry (dis_iub, dis_jub) +=
value;
201 auu.dis_entry_assembly();
202 aub.dis_entry_assembly();
203 abu.dis_entry_assembly();
216 _uu.set_pattern_dimension (map_d);
217 _ub.set_pattern_dimension (map_d);
218 _bu.set_pattern_dimension (map_d);
219 _bb.set_pattern_dimension (map_d);
228 #ifdef _RHEOLEF_HAVE_MPI
230 is_sym = mpi::all_reduce (omega.comm(),
size_type(is_sym), mpi::minimum<size_type>());
233 _uu.set_symmetry (is_sym);
234 _bb.set_symmetry (is_sym);
236 _uu.set_definite_positive (is_sym);
237 _bb.set_definite_positive (is_sym);
240 template <
class T,
class M>
241 template <
class Expr>
251 template <
class T,
class M>
252 template <
class Expr>
260 do_integrate_internal (
gh.level_set(),
gh.band(),
gh,
expr, iopt,
true);
field gh(Float epsilon, Float t, const field &uh, const test &v)
see the band page for the full documentation
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
see the integrate_option page for the full documentation
bool _is_inside_on_local_sides
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
void local_lump(Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m)
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 local_invert(Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m, bool is_diag)
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.