3 #ifndef DUNE_PDELAB_LOCALOPERATOR_NUMERICALNONLINEARJACOBIANAPPLY_HH
4 #define DUNE_PDELAB_LOCALOPERATOR_NUMERICALNONLINEARJACOBIANAPPLY_HH
32 template<
typename Imp>
45 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename Y>
48 const LFSU& lfsu,
const X& x,
const X& z,
49 const LFSV& lfsv, Y& y)
const
51 using R =
typename Y::value_type;
60 ResidualVector down(y.size()),up(y.size());
61 auto downview = down.weightedAccumulationView(y.weight());
62 auto upview = up.weightedAccumulationView(y.weight());
64 asImp().alpha_volume(eg,lfsu,u,lfsv,downview);
65 for (decltype(n) j = 0; j < n; ++j)
69 R delta = epsilon*(1.0 + abs(u(lfsu,j)));
71 asImp().alpha_volume(eg,lfsu,u,lfsv,upview);
72 for (decltype(m) i = 0; i < m; ++i)
73 y.rawAccumulate(lfsv,i,((up(lfsv,i)-down(lfsv,i))/delta)*z(lfsu,j));
74 u(lfsu,j) = x(lfsu,j);
80 Imp& asImp () {
return static_cast<Imp &
> (*this); }
81 const Imp& asImp ()
const {
return static_cast<const Imp &
>(*this); }
95 template<
typename Imp>
108 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename Y>
111 const LFSU& lfsu,
const X& x,
const X& z,
112 const LFSV& lfsv, Y& y)
const
114 using R =
typename Y::value_type;
117 auto m = lfsv.
size();
118 auto n = lfsu.size();
123 ResidualVector down(y.size()),up(y.size());
124 auto downview = down.weightedAccumulationView(y.weight());
125 auto upview = up.weightedAccumulationView(y.weight());
127 asImp().alpha_volume_post_skeleton(eg,lfsu,u,lfsv,downview);
128 for (decltype(n) j = 0; j < n; ++j)
132 R delta = epsilon*(1.0 + abs(u(lfsu,j)));
134 asImp().alpha_volume_post_skeleton(eg,lfsu,u,lfsv,upview);
135 for (decltype(m) i = 0; i < m; ++i)
136 y.rawAccumulate(lfsv,i,((up(lfsv,i)-down(lfsv,i))/delta)*z(lfsu,j));
137 u(lfsu,j) = x(lfsu,j);
142 const double epsilon;
143 Imp& asImp () {
return static_cast<Imp &
> (*this);}
144 const Imp& asImp ()
const {
return static_cast<const Imp &
>(*this);}
156 template<
typename Imp>
169 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename Y>
172 const LFSU& lfsu_s,
const X& x_s,
const X& z_s,
const LFSV& lfsv_s,
173 const LFSU& lfsu_n,
const X& x_n,
const X& z_n,
const LFSV& lfsv_n,
174 Y& y_s, Y& y_n)
const
176 using R =
typename X::value_type;
179 auto m_s=lfsv_s.
size();
180 auto m_n=lfsv_n.size();
181 auto n_s=lfsu_s.size();
182 auto n_n=lfsu_n.size();
188 ResidualVector down_s(y_s.size()),up_s(y_s.size());
189 auto downview_s = down_s.weightedAccumulationView(1.0);
190 auto upview_s = up_s.weightedAccumulationView(1.0);
192 ResidualVector down_n(y_n.size()),up_n(y_n.size());
193 auto downview_n = down_n.weightedAccumulationView(1.0);
194 auto upview_n = up_n.weightedAccumulationView(1.0);
197 asImp().alpha_skeleton(
ig,lfsu_s,u_s,lfsv_s,lfsu_n,u_n,lfsv_n,downview_s,downview_n);
200 for (decltype(n_s) j = 0; j < n_s; ++j)
205 R delta = epsilon*(1.0 + abs(u_s(lfsu_s,j)));
206 u_s(lfsu_s,j) += delta;
207 asImp().alpha_skeleton(
ig,lfsu_s,u_s,lfsv_s,lfsu_n,u_n,lfsv_n,upview_s,upview_n);
208 for (decltype(m_s) i = 0; i < m_s; ++i)
209 y_s.accumulate(lfsv_s,i,((up_s(lfsv_s,i)-down_s(lfsv_s,i))/delta)*z_s(lfsu_s,j));
210 for (decltype(m_n) i = 0; i < m_n; i++)
211 y_n.accumulate(lfsv_n,i,((up_n(lfsv_n,i)-down_n(lfsv_n,i))/delta)*x_s(lfsu_s,j));
212 u_s(lfsu_s,j) = x_s(lfsu_s,j);
216 for (decltype(n_n) j = 0; j < n_n; ++j)
221 R delta = epsilon*(1.0+abs(u_n(lfsu_n,j)));
222 u_n(lfsu_n,j) += delta;
223 asImp().alpha_skeleton(
ig,lfsu_s,u_s,lfsv_s,lfsu_n,u_n,lfsv_n,upview_s,upview_n);
224 for (decltype(m_s) i = 0; i < m_s; ++i)
225 y_s.accumulate(lfsv_s,i,((up_s(lfsv_s,i)-down_s(lfsv_s,i))/delta)*z_n(lfsu_n,j));
226 for (decltype(m_n) i = 0; i < m_n; ++i)
227 y_n.accumulate(lfsv_n,i,((up_n(lfsv_n,i)-down_n(lfsv_n,i))/delta)*x_n(lfsu_n,j));
228 u_n(lfsu_n,j) = x_n(lfsu_n,j);
233 const double epsilon;
234 Imp& asImp () {
return static_cast<Imp &
> (*this); }
235 const Imp& asImp ()
const {
return static_cast<const Imp &
>(*this); }
247 template<
typename Imp>
260 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename Y>
263 const LFSU& lfsu_s,
const X& x_s,
const X& z_s,
264 const LFSV& lfsv_s, Y& y_s)
const
266 using R =
typename Y::value_type;
269 auto m_s=lfsv_s.
size();
270 auto n_s=lfsu_s.size();
275 ResidualVector down_s(y_s.size()),up_s(y_s.size());
276 auto downview_s = down_s.weightedAccumulationView(1.0);
277 auto upview_s = up_s.weightedAccumulationView(1.0);
280 asImp().alpha_boundary(
ig,lfsu_s,u_s,lfsv_s,downview_s);
283 for (decltype(n_s) j = 0; j < n_s; ++j)
286 R delta = epsilon*(1.0+std::abs(u_s(lfsu_s,j)));
287 u_s(lfsu_s,j) += delta;
288 asImp().alpha_boundary(
ig,lfsu_s,u_s,lfsv_s,upview_s);
289 for (decltype(m_s) i = 0; i < m_s; ++i)
290 y_s.rawAccumulate(lfsv_s,i,((up_s(lfsv_s,i)-down_s(lfsv_s,i))/delta)*x_s(lfsu_s,j));
291 u_s(lfsu_s,j) = x_s(lfsu_s,j);
296 const double epsilon;
297 Imp& asImp () {
return static_cast<Imp &
> (*this); }
298 const Imp& asImp ()
const {
return static_cast<const Imp &
>(*this); }
305 #endif // DUNE_PDELAB_LOCALOPERATOR_NUMERICALNONLINEARJACOBIANAPPLY_HH