70 if (! opt.
active)
return uh;
72 T epsilon = std::numeric_limits<T>::epsilon();
75 omega.neighbour_guard();
77 "argument may be a discontinuous approximation: "
80 if (k == 0)
return uh;
81 check_macro (k == 1,
"order = " << k <<
" not yet supported");
83 check_macro (map_d == 1,
"dimension map_d = " << map_d <<
" not yet supported");
95 std::array<bool,ns_loc_max> is_on_boundary;
96 std::array<bool,ns_loc_max> is_upstream;
97 std::array<point_basic<T>,ns_loc_max> xK1;
98 size_type index [ns_loc_max][nss_loc_max];
99 T alpha [ns_loc_max][nss_loc_max];
100 std::array<T,ns_loc_max> tilde,
delta, Delta, hat_Delta;
101 std::array<geo_element::orientation_type,ns_loc_max> orient;
103 for (
size_type ie = 0, ne = omega.size(); ie < ne; ++ie) {
104 const geo_element& K = omega.get_geo_element (map_d, ie);
113 for (
size_type iv_loc = 0; iv_loc < nv_loc; ++iv_loc) {
115 xK += omega.dis_node (dis_iv);
117 xK /=
T(
int(nv_loc));
121 for (
size_type is_loc = 0, ns_loc = K.
n_subgeo(map_d-1); is_loc < ns_loc; ++is_loc) {
123 size_type dis_is = (map_d == 1) ? K[is_loc] : ((map_d == 2) ? K.
edge(is_loc) : K.
face(is_loc));
125 const geo_element& S = omega.dis_get_geo_element (map_d-1, dis_is);
132 is_on_boundary [is_loc] = (dis_ie1 == std::numeric_limits<size_type>::max());
133 if (is_on_boundary[is_loc]) {
134 index [is_loc][0] = is_loc;
135 alpha [is_loc][0] = 1;
138 const geo_element& K1 = omega.dis_get_geo_element (map_d, dis_ie1);
141 xK1 [is_loc] = {0,0,0};
143 for (
size_type iv1_loc = 0; iv1_loc < nv1_loc; ++iv1_loc) {
146 xK1 [is_loc] += omega.dis_node (dis_iv1);
148 xK1 [is_loc] /=
T(
int(nv1_loc));
151 index [is_loc][0] = is_loc;
152 alpha [is_loc][0] = meas_K.
dof(ie)/(meas_K.
dis_dof(dis_ie1) + meas_K.
dof(ie));
161 for (
size_type iv_loc = 0; iv_loc < nv_loc; ++iv_loc) {
163 bar_u_K += uh.
dof (idof);
167 T hK =
pow (meas_K.
dof(ie), 1./map_d);
169 T sum_Delta_plus = 0,
171 for (
size_type is_loc = 0, ns_loc = K.
n_subgeo(map_d-1); is_loc < ns_loc; ++is_loc) {
173 size_type dis_is = (map_d == 1) ? K[is_loc] : ((map_d == 2) ? K.
edge(is_loc) : K.
face(is_loc));
174 const geo_element& S = omega.dis_get_geo_element (map_d-1, dis_is);
179 is_on_boundary [is_loc] = (dis_ie1 == std::numeric_limits<size_type>::max());
180 is_upstream[is_loc] = is_on_boundary [is_loc] && (K.
ios_dis_ie() == 0);
183 if (
false && is_upstream [is_loc]) {
187 bar_u_S = uh.
dof (idof);
190 tilde [is_loc] = orient[is_loc]*(bar_u_S - bar_u_K);
194 if (is_on_boundary [is_loc]) {
195 if (is_upstream [is_loc]) {
202 const geo_element& K1 = omega.dis_get_geo_element (map_d, dis_ie1);
204 for (
size_type iv1_loc = 0; iv1_loc < nv1_loc; ++iv1_loc) {
208 bar_u_K1 += uh.
dis_dof (dis_idof);
214 delta [is_loc] = orient[is_loc]*alpha[is_loc][0]*(bar_u_K1 - bar_u_K);
218 sum_Delta_plus += max (
T(0.), Delta [is_loc]);
219 sum_Delta_minus += max (
T(0.), -Delta [is_loc]);
223 T r = (fabs(sum_Delta_plus) > tol) ? sum_Delta_minus/sum_Delta_plus : 0;
227 for (
size_type is_loc = 0, ns_loc = K.
n_subgeo(map_d-1); is_loc < ns_loc; ++is_loc) {
228 hat_Delta [is_loc] = (r == 0) ? Delta [is_loc]
229 : min(
T(1.), r)*max(
T(0.), Delta [is_loc])
230 - min(
T(1.),1./r)*max(
T(0.), -Delta [is_loc]);
234 for (
size_type iv_loc = 0, nv_loc = K.
size(); iv_loc < nv_loc; ++iv_loc) {
237 if (
false && is_on_boundary [is_loc] && ! is_upstream[is_loc]) {
238 vh.
dof (idof) = uh.
dof (idof);
240 vh.
dof (idof) = bar_u_K + orient[is_loc]*hat_Delta [is_loc];
243 <<
", u="<<uh.
dof (idof)<<
" -> " << vh.
dof (idof));
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
std::enable_if< details::is_field_expr_v2_nonlinear_arg< Expr >::value &&!is_undeterminated< Result >::value, Result >::type integrate(const geo_basic< T, M > &omega, const Expr &expr, const integrate_option &iopt, Result dummy=Result())
see the integrate page for the full documentation
field_basic< T, M > limiter(const field_basic< T, M > &uh, const T &bar_g_S, const limiter_option &opt)
see the limiter page for the full documentation