1 #ifndef _RHEO_FORM_VF_ASSEMBLY_H
2 #define _RHEO_FORM_VF_ASSEMBLY_H
3 #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"
76 template <
class T>
T norm_max (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
m);
77 template <
class T>
bool check_is_symmetric (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
m,
const T& tol_m_max);
78 template <
class T>
void local_lump (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
m);
79 template <
class T>
void local_invert (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
m,
bool is_diag);
84 template <
class T,
class M>
98 _X =
expr.get_trial_space();
99 _Y =
expr.get_test_space();
101 check_macro (
band.get_background_geo() == _X.get_geo().get_background_geo(),
102 "assembly: incompatible integration domain "<<
band.name() <<
" and trial function based domain "
103 << _X.get_geo().name());
104 check_macro (
band.get_background_geo() == _Y.get_geo().get_background_geo(),
105 "assembly: incompatible integration domain "<<
band.name() <<
" and test function based domain "
106 << _Y.get_geo().name());
111 _X.get_constitution().have_compact_support_inside_element() &&
112 _Y.get_constitution().have_compact_support_inside_element(),
113 "local inversion requires compact support inside elements (e.g. discontinuous or bubble)");
115 "local mass lumping requires no derivative operators");
120 expr.initialize (omega, iopt);
122 expr.initialize (
gh, iopt);
124 expr.template valued_check<T>();
125 asr<T,M> auu (_Y.iu_ownership(), _X.iu_ownership()),
126 aub (_Y.iu_ownership(), _X.ib_ownership()),
127 abu (_Y.ib_ownership(), _X.iu_ownership()),
128 abb (_Y.ib_ownership(), _X.ib_ownership());
129 std::vector<size_type> dis_idy, dis_jdx;
130 Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> ak;
132 if (_X.get_constitution().is_discontinuous()) _X.get_constitution().neighbour_guard();
133 if (_Y.get_constitution().is_discontinuous()) _Y.get_constitution().neighbour_guard();
136 for (
size_type ie = 0, ne = omega.size(map_d); ie < ne; ie++) {
140 const geo_element& K = omega.get_geo_element (map_d, ie);
142 _X.get_constitution().assembly_dis_idof (omega, K, dis_jdx);
143 _Y.get_constitution().assembly_dis_idof (omega, K, dis_idy);
147 _X.dis_idof (
L, dis_jdx);
148 _Y.dis_idof (
L, dis_idy);
150 expr.evaluate (omega, K, ak);
155 T eps_ak_max = eps*ak_max;
163 "invalid sizes ak("<<ak.rows()<<
","<<ak.cols()
164 <<
") with dis_idy("<<dis_idy.size()<<
") and dis_jdx("<<dis_jdx.size()<<
")");
165 for (
size_type loc_idof = 0,
ny = ak.rows(); loc_idof <
ny; loc_idof++) {
166 for (
size_type loc_jdof = 0,
nx = ak.cols(); loc_jdof <
nx; loc_jdof++) {
168 const T&
value = ak (loc_idof, loc_jdof);
177 if (fabs(
value) < eps_ak_max)
continue;
182 size_type dis_jub = _X.dis_idof2dis_iub (dis_jdof);
185 if (_X.dis_is_blocked(dis_jdof)) abb.
dis_entry (dis_iub, dis_jub) +=
value;
186 else abu.dis_entry (dis_iub, dis_jub) +=
value;
188 if (_X.dis_is_blocked(dis_jdof)) aub.dis_entry (dis_iub, dis_jub) +=
value;
189 else auu.dis_entry (dis_iub, dis_jub) +=
value;
198 auu.dis_entry_assembly();
199 aub.dis_entry_assembly();
200 abu.dis_entry_assembly();
213 _uu.set_pattern_dimension (map_d);
214 _ub.set_pattern_dimension (map_d);
215 _bu.set_pattern_dimension (map_d);
216 _bb.set_pattern_dimension (map_d);
225 #ifdef _RHEOLEF_HAVE_MPI
227 is_sym = mpi::all_reduce (omega.comm(),
size_type(is_sym), mpi::minimum<size_type>());
229 #endif // _RHEOLEF_HAVE_MPI
230 _uu.set_symmetry (is_sym);
231 _bb.set_symmetry (is_sym);
233 _uu.set_definite_positive (is_sym);
234 _bb.set_definite_positive (is_sym);
237 template <
class T,
class M>
238 template <
class Expr>
248 template <
class T,
class M>
249 template <
class Expr>
257 assembly_internal (
gh.level_set(),
gh.band(),
gh,
expr, iopt,
true);
261 #endif // _RHEO_FORM_VF_ASSEMBLY_H