24 #include "rheolef/limiter.h"
25 #include "rheolef/integrate.h"
40 if (fabs(
a) < tol)
return 0;
41 if (fabs(
b) < tol)
return 0;
42 T sgn_a = (
a >= 0) ? 1 : -1;
43 T sgn_b = (
b >= 0) ? 1 : -1;
44 T res = (sgn_a == sgn_b) ? sgn_a*min(fabs(
a),fabs(
b)) : 0;
53 trace_macro (
"minmod_tdb("<<yield_a<<
","<<
a<<
","<<
b<<
") = " << res);
63 template <
class T,
class M>
70 if (! opt.
active)
return uh;
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;
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) {
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));
252 #define _RHEOLEF_instanciation(T,M) \
253 template field_basic<T,M> limiter ( \
254 const field_basic<T,M>&, \
256 const limiter_option&);
259 #ifdef _RHEOLEF_HAVE_MPI
261 #endif // _RHEOLEF_HAVE_MPI