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());
111 size_type n_derivative = expr.n_derivative();
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();
138 const T eps = 1e3*std::numeric_limits<T>::epsilon();
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;
184 size_type dis_iub = _Y.dis_idof2dis_iub (dis_idof);
185 size_type dis_jub = _X.dis_idof2dis_iub (dis_jdof);
187 if (_Y.dis_is_blocked(dis_idof))
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);