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();
75 const T eps = 1e3*std::numeric_limits<T>::epsilon();
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;
120 size_type dis_iub = _Y.dis_idof2dis_iub (dis_idof);
121 size_type dis_jub = _X.dis_idof2dis_iub (dis_jdof);
123 if (_Y.dis_is_blocked(dis_idof))
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);