dune-pdelab  2.5-dev
numericalnonlinearjacobianapply.hh
Go to the documentation of this file.
1 // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=8 sw=2 sts=2:
3 #ifndef DUNE_PDELAB_LOCALOPERATOR_NUMERICALNONLINEARJACOBIANAPPLY_HH
4 #define DUNE_PDELAB_LOCALOPERATOR_NUMERICALNONLINEARJACOBIANAPPLY_HH
5 
6 #include <cmath>
7 
10 
11 namespace Dune {
12  namespace PDELab {
13 
17 
19  //
20  // Numerical implementation of nonlinear jacobian_apply_*() in terms of alpha_*()
21  //
22 
24 
32  template<typename Imp>
34  {
35  public:
37  : epsilon(1e-7)
38  {}
39 
41  : epsilon(epsilon_)
42  {}
43 
45  template<typename EG, typename LFSU, typename X, typename LFSV, typename Y>
47  const EG& eg,
48  const LFSU& lfsu, const X& x, const X& z,
49  const LFSV& lfsv, Y& y) const
50  {
51  using R = typename Y::value_type;
53 
54  auto m = lfsv.size();
55  auto n = lfsu.size();
56 
57  X u(x);
58 
59  // Notice that in general lfsv.size() != y.size()
60  ResidualVector down(y.size()),up(y.size());
61  auto downview = down.weightedAccumulationView(y.weight());
62  auto upview = up.weightedAccumulationView(y.weight());
63 
64  asImp().alpha_volume(eg,lfsu,u,lfsv,downview);
65  for (decltype(n) j = 0; j < n; ++j) // loop over columns
66  {
67  using namespace std;
68  up = 0.0;
69  R delta = epsilon*(1.0 + abs(u(lfsu,j)));
70  u(lfsu,j) += delta;
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);
75  }
76  }
77 
78  private:
79  const double epsilon; // problem: this depends on data type R!
80  Imp& asImp () { return static_cast<Imp &> (*this); }
81  const Imp& asImp () const { return static_cast<const Imp &>(*this); }
82  };
83 
86 
95  template<typename Imp>
97  {
98  public:
100  : epsilon(1e-7)
101  {}
102 
104  : epsilon(epsilon_)
105  {}
106 
108  template<typename EG, typename LFSU, typename X, typename LFSV, typename Y>
110  const EG& eg,
111  const LFSU& lfsu, const X& x, const X& z,
112  const LFSV& lfsv, Y& y) const
113  {
114  using R = typename Y::value_type;
116 
117  auto m = lfsv.size();
118  auto n = lfsu.size();
119 
120  X u(x);
121 
122  // Notice that in general lfsv.size() != y.size()
123  ResidualVector down(y.size()),up(y.size());
124  auto downview = down.weightedAccumulationView(y.weight());
125  auto upview = up.weightedAccumulationView(y.weight());
126 
127  asImp().alpha_volume_post_skeleton(eg,lfsu,u,lfsv,downview);
128  for (decltype(n) j = 0; j < n; ++j) // loop over columns
129  {
130  using namespace std;
131  up = 0.0;
132  R delta = epsilon*(1.0 + abs(u(lfsu,j)));
133  u(lfsu,j) += delta;
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);
138  }
139  }
140 
141  private:
142  const double epsilon; // problem: this depends on data type R!
143  Imp& asImp () {return static_cast<Imp &> (*this);}
144  const Imp& asImp () const {return static_cast<const Imp &>(*this);}
145  };
146 
148 
156  template<typename Imp>
158  {
159  public:
161  : epsilon(1e-7)
162  {}
163 
165  : epsilon(epsilon_)
166  {}
167 
169  template<typename IG, typename LFSU, typename X, typename LFSV, typename Y>
171  const IG& ig,
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
175  {
176  using R = typename X::value_type;
178 
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();
183 
184  X u_s(x_s);
185  X u_n(x_n);
186 
187  // Notice that in general lfsv_s.size() != y_s.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);
191 
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);
195 
196  // base line
197  asImp().alpha_skeleton(ig,lfsu_s,u_s,lfsv_s,lfsu_n,u_n,lfsv_n,downview_s,downview_n);
198 
199  // jiggle in self
200  for (decltype(n_s) j = 0; j < n_s; ++j)
201  {
202  using namespace std;
203  up_s = 0.0;
204  up_n = 0.0;
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);
213  }
214 
215  // jiggle in neighbor
216  for (decltype(n_n) j = 0; j < n_n; ++j)
217  {
218  using namespace std;
219  up_s = 0.0;
220  up_n = 0.0;
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);
229  }
230  }
231 
232  private:
233  const double epsilon; // problem: this depends on data type R!
234  Imp& asImp () { return static_cast<Imp &> (*this); }
235  const Imp& asImp () const { return static_cast<const Imp &>(*this); }
236  };
237 
239 
247  template<typename Imp>
249  {
250  public:
252  : epsilon(1e-7)
253  {}
254 
256  : epsilon(epsilon_)
257  {}
258 
260  template<typename IG, typename LFSU, typename X, typename LFSV, typename Y>
262  const IG& ig,
263  const LFSU& lfsu_s, const X& x_s, const X& z_s,
264  const LFSV& lfsv_s, Y& y_s) const
265  {
266  using R = typename Y::value_type;
268 
269  auto m_s=lfsv_s.size();
270  auto n_s=lfsu_s.size();
271 
272  X u_s(x_s);
273 
274  // Notice that in general lfsv_s.size() != y_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);
278 
279  // base line
280  asImp().alpha_boundary(ig,lfsu_s,u_s,lfsv_s,downview_s);
281 
282  // jiggle in self
283  for (decltype(n_s) j = 0; j < n_s; ++j)
284  {
285  up_s = 0.0;
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);
292  }
293  }
294 
295  private:
296  const double epsilon; // problem: this depends on data type R!
297  Imp& asImp () { return static_cast<Imp &> (*this); }
298  const Imp& asImp () const { return static_cast<const Imp &>(*this); }
299  };
300 
302  }
303 }
304 
305 #endif // DUNE_PDELAB_LOCALOPERATOR_NUMERICALNONLINEARJACOBIANAPPLY_HH
const IG & ig
Definition: constraints.hh:148
Implements nonlinear version of jacobian_apply_volume() based on alpha_volume()
Definition: numericalnonlinearjacobianapply.hh:33
const Entity & e
Definition: localfunctionspace.hh:111
void jacobian_apply_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const X &z_s, const LFSV &lfsv_s, Y &y_s) const
apply local jacobian of the boundaryterm
Definition: numericalnonlinearjacobianapply.hh:261
void jacobian_apply_volume(const EG &eg, const LFSU &lfsu, const X &x, const X &z, const LFSV &lfsv, Y &y) const
apply local jacobian of the volume term
Definition: numericalnonlinearjacobianapply.hh:46
void jacobian_apply_volume_post_skeleton(const EG &eg, const LFSU &lfsu, const X &x, const X &z, const LFSV &lfsv, Y &y) const
apply local jacobian of the volume term (post skeleton part)
Definition: numericalnonlinearjacobianapply.hh:109
A container for storing data associated with the degrees of freedom of a LocalFunctionSpace.
Definition: localvector.hh:171
Definition: numericalnonlinearjacobianapply.hh:96
size_type size() const
The size of the container.
Definition: localvector.hh:245
NumericalNonlinearJacobianApplyVolumePostSkeleton(double epsilon_)
Definition: numericalnonlinearjacobianapply.hh:103
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
NumericalNonlinearJacobianApplyBoundary(double epsilon_)
Definition: numericalnonlinearjacobianapply.hh:255
Implements nonlinear version of jacobian_apply_boundary() based on alpha_boundary() ...
Definition: numericalnonlinearjacobianapply.hh:248
NumericalNonlinearJacobianApplySkeleton(double epsilon_)
Definition: numericalnonlinearjacobianapply.hh:164
STL namespace.
Implements nonlinear version of jacobian_apply_skeleton() based on alpha_skeleton() ...
Definition: numericalnonlinearjacobianapply.hh:157
NumericalNonlinearJacobianApplyVolume(double epsilon_)
Definition: numericalnonlinearjacobianapply.hh:40
void jacobian_apply_skeleton(const IG &ig, const LFSU &lfsu_s, const X &x_s, const X &z_s, const LFSV &lfsv_s, const LFSU &lfsu_n, const X &x_n, const X &z_n, const LFSV &lfsv_n, Y &y_s, Y &y_n) const
apply local jacobian of the skeleton term
Definition: numericalnonlinearjacobianapply.hh:170
NumericalNonlinearJacobianApplyBoundary()
Definition: numericalnonlinearjacobianapply.hh:251
NumericalNonlinearJacobianApplyVolumePostSkeleton()
Definition: numericalnonlinearjacobianapply.hh:99
NumericalNonlinearJacobianApplyVolume()
Definition: numericalnonlinearjacobianapply.hh:36
NumericalNonlinearJacobianApplySkeleton()
Definition: numericalnonlinearjacobianapply.hh:160