3#ifndef DUNE_PDELAB_LOCALOPERATOR_NUMERICALJACOBIAN_HH
4#define DUNE_PDELAB_LOCALOPERATOR_NUMERICALJACOBIAN_HH
9#include <dune/pdelab/gridoperator/common/localmatrix.hh>
30 template<
typename Imp>
43 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
47 const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
50 typedef typename X::value_type D;
51 typedef typename Jacobian::value_type R;
53 typedef typename ResidualVector::WeightedAccumulationView ResidualView;
55 const int m=lfsv.
size();
56 const int n=lfsu.size();
61 ResidualVector down(mat.nrows(),0.),up(mat.nrows());
62 ResidualView downview = down.weightedAccumulationView(mat.weight());
63 ResidualView upview = up.weightedAccumulationView(mat.weight());
66 asImp().alpha_volume(eg,lfsu,u,lfsv,downview);
67 for (
int j=0; j<n; j++)
70 D delta = epsilon*(1.0+std::abs(u(lfsu,j)));
72 asImp().alpha_volume(eg,lfsu,u,lfsv,upview);
73 for (
int i=0; i<m; i++)
74 mat.rawAccumulate(lfsv,i,lfsu,j,(up(lfsv,i)-down(lfsv,i))/delta);
75 u(lfsu,j) = x(lfsu,j);
81 Imp& asImp () {
return static_cast<Imp &
> (*this); }
82 const Imp& asImp ()
const {
return static_cast<const Imp &
>(*this); }
94 template<
typename Imp>
107 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
111 const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
114 typedef typename X::value_type D;
115 typedef typename Jacobian::value_type R;
117 typedef typename ResidualVector::WeightedAccumulationView ResidualView;
119 const int m=lfsv.
size();
120 const int n=lfsu.size();
125 ResidualVector down(mat.nrows(),0.),up(mat.nrows());
126 ResidualView downview = down.weightedAccumulationView(mat.weight());
127 ResidualView upview = up.weightedAccumulationView(mat.weight());
129 asImp().alpha_volume_post_skeleton(eg,lfsu,u,lfsv,downview);
130 for (
int j=0; j<n; j++)
133 D delta = epsilon*(1.0+std::abs(u(lfsu,j)));
135 asImp().alpha_volume_post_skeleton(eg,lfsu,u,lfsv,upview);
136 for (
int i=0; i<m; i++)
137 mat.rawAccumulate(lfsv,i,lfsu,j,(up(lfsv,i)-down(lfsv,i))/delta);
138 u(lfsu,j) = x(lfsu,j);
143 const double epsilon;
144 Imp& asImp () {
return static_cast<Imp &
> (*this); }
145 const Imp& asImp ()
const {
return static_cast<const Imp &
>(*this); }
155 template<
typename Imp>
168 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
172 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
173 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
174 Jacobian& mat_ss, Jacobian& mat_sn,
175 Jacobian& mat_ns, Jacobian& mat_nn)
const
177 typedef typename X::value_type D;
178 typedef typename Jacobian::value_type R;
180 typedef typename ResidualVector::WeightedAccumulationView ResidualView;
182 const int m_s=lfsv_s.
size();
183 const int m_n=lfsv_n.size();
184 const int n_s=lfsu_s.size();
185 const int n_n=lfsu_n.size();
191 ResidualVector down_s(mat_ss.nrows()),up_s(mat_ss.nrows());
192 ResidualView downview_s = down_s.weightedAccumulationView(1.0);
193 ResidualView upview_s = up_s.weightedAccumulationView(1.0);
195 ResidualVector down_n(mat_nn.nrows()),up_n(mat_nn.nrows());
196 ResidualView downview_n = down_n.weightedAccumulationView(1.0);
197 ResidualView upview_n = up_n.weightedAccumulationView(1.0);
200 asImp().alpha_skeleton(ig,lfsu_s,u_s,lfsv_s,lfsu_n,u_n,lfsv_n,downview_s,
204 for (
int j=0; j<n_s; j++)
208 D delta = epsilon*(1.0+std::abs(u_s(lfsu_s,j)));
209 u_s(lfsu_s,j) += delta;
210 asImp().alpha_skeleton(ig,lfsu_s,u_s,lfsv_s,lfsu_n,u_n,lfsv_n,upview_s,
212 for (
int i=0; i<m_s; i++)
213 mat_ss.accumulate(lfsv_s,i,lfsu_s,j,(up_s(lfsv_s,i)-down_s(lfsv_s,i))/delta);
214 for (
int i=0; i<m_n; i++)
215 mat_ns.accumulate(lfsv_n,i,lfsu_s,j,(up_n(lfsv_n,i)-down_n(lfsv_n,i))/delta);
216 u_s(lfsu_s,j) = x_s(lfsu_s,j);
220 for (
int j=0; j<n_n; j++)
224 D delta = epsilon*(1.0+std::abs(u_n(lfsu_n,j)));
225 u_n(lfsu_n,j) += delta;
226 asImp().alpha_skeleton(ig,lfsu_s,u_s,lfsv_s,lfsu_n,u_n,lfsv_n,upview_s,
228 for (
int i=0; i<m_s; i++)
229 mat_sn.accumulate(lfsv_s,i,lfsu_n,j,(up_s(lfsv_s,i)-down_s(lfsv_s,i))/delta);
230 for (
int i=0; i<m_n; i++)
231 mat_nn.accumulate(lfsv_n,i,lfsu_n,j,(up_n(lfsv_n,i)-down_n(lfsv_n,i))/delta);
232 u_n(lfsu_n,j) = x_n(lfsu_n,j);
237 const double epsilon;
238 Imp& asImp () {
return static_cast<Imp &
> (*this); }
239 const Imp& asImp ()
const {
return static_cast<const Imp &
>(*this); }
249 template<
typename Imp>
262 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
266 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
267 Jacobian& mat_ss)
const
269 typedef typename X::value_type D;
270 typedef typename Jacobian::value_type R;
272 typedef typename ResidualVector::WeightedAccumulationView ResidualView;
274 const int m_s=lfsv_s.
size();
275 const int n_s=lfsu_s.size();
280 ResidualVector down_s(mat_ss.nrows()),up_s(mat_ss.nrows());
281 ResidualView downview_s = down_s.weightedAccumulationView(mat_ss.weight());
282 ResidualView upview_s = up_s.weightedAccumulationView(mat_ss.weight());;
286 asImp().alpha_boundary(ig,lfsu_s,u_s,lfsv_s,downview_s);
289 for (
int j=0; j<n_s; j++)
292 D delta = epsilon*(1.0+std::abs(u_s(lfsu_s,j)));
293 u_s(lfsu_s,j) += delta;
294 asImp().alpha_boundary(ig,lfsu_s,u_s,lfsv_s,upview_s);
295 for (
int i=0; i<m_s; i++)
296 mat_ss.rawAccumulate(lfsv_s,i,lfsu_s,j,(up_s(lfsv_s,i)-down_s(lfsv_s,i))/delta);
297 u_s(lfsu_s,j) = x_s(lfsu_s,j);
302 const double epsilon;
303 Imp& asImp () {
return static_cast<Imp &
> (*this); }
304 const Imp& asImp ()
const {
return static_cast<const Imp &
>(*this); }
A container for storing data associated with the degrees of freedom of a LocalFunctionSpace.
Definition: localvector.hh:184
size_type size() const
The size of the container.
Definition: localvector.hh:318
Implement jacobian_boundary() based on alpha_boundary()
Definition: numericaljacobian.hh:251
void jacobian_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, Jacobian &mat_ss) const
compute local jacobian of the boundary term
Definition: numericaljacobian.hh:265
Implement jacobian_skeleton() based on alpha_skeleton()
Definition: numericaljacobian.hh:157
void jacobian_skeleton(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, const LFSU &lfsu_n, const X &x_n, const LFSV &lfsv_n, Jacobian &mat_ss, Jacobian &mat_sn, Jacobian &mat_ns, Jacobian &mat_nn) const
compute local jacobian of the skeleton term
Definition: numericaljacobian.hh:171
Definition: numericaljacobian.hh:96
void jacobian_volume_post_skeleton(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, Jacobian &mat) const
compute local post-skeleton jacobian of the volume term
Definition: numericaljacobian.hh:110
Implement jacobian_volume() based on alpha_volume()
Definition: numericaljacobian.hh:32
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, Jacobian &mat) const
compute local jacobian of the volume term
Definition: numericaljacobian.hh:46
Dune namespace.
Definition: alignedallocator.hh:11