DUNE PDELab (git)

numericaljacobianapply.hh
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_NUMERICALJACOBIANAPPLY_HH
4#define DUNE_PDELAB_LOCALOPERATOR_NUMERICALJACOBIANAPPLY_HH
5
6#include <dune/pdelab/localoperator/numericalnonlinearjacobianapply.hh>
7
8namespace Dune {
9 namespace PDELab {
10
14
16 //
17 // Numerical implementation of jacobian_apply_*() in terms of alpha_*()
18 //
19
21
31 template<typename Imp>
34 {
35 public:
36
37 // We need to reimport the name from the base class; otherwise clang plays stupid and refuses to
38 // find the overload in the base class
40
43 , epsilon(1e-7)
44 {}
45
46 NumericalJacobianApplyVolume (double epsilon_)
48 , epsilon(epsilon_)
49 {}
50
52 template<typename EG, typename LFSU, typename X, typename LFSV,
53 typename Y>
55 ( const EG& eg,
56 const LFSU& lfsu, const X& x, const LFSV& lfsv,
57 Y& y) const
58 {
59 typedef typename X::value_type D;
60 typedef typename Y::value_type R;
62 typedef typename ResidualVector::WeightedAccumulationView ResidualView;
63
64 const int m=lfsv.size();
65 const int n=lfsu.size();
66
67 X u(x);
68
69 // Notice that in general lfsv.size() != y.size()
70 ResidualVector down(y.size()),up(y.size());
71 ResidualView downview = down.weightedAccumulationView(y.weight());
72 ResidualView upview = up.weightedAccumulationView(y.weight());
73
74 asImp().alpha_volume(eg,lfsu,u,lfsv,downview);
75 for (int j=0; j<n; j++) // loop over columns
76 {
77 up = 0.0;
78 D delta = epsilon*(1.0+std::abs(u(lfsu,j)));
79 u(lfsu,j) += delta;
80 asImp().alpha_volume(eg,lfsu,u,lfsv,upview);
81 for (int i=0; i<m; i++)
82 y.rawAccumulate(lfsv,i,((up(lfsv,i)-down(lfsv,i))/delta)*x(lfsu,j));
83 u(lfsu,j) = x(lfsu,j);
84 }
85 }
86
87 private:
88 const double epsilon; // problem: this depends on data type R!
89 Imp& asImp () { return static_cast<Imp &> (*this); }
90 const Imp& asImp () const { return static_cast<const Imp &>(*this); }
91 };
92
95
106 template<typename Imp>
109 {
110 public:
111
112 // We need to reimport the name from the base class; otherwise clang plays stupid and refuses to
113 // find the overload in the base class
115
118 , epsilon(1e-7)
119 {}
120
123 , epsilon(epsilon_)
124 {}
125
127 template<typename EG, typename LFSU, typename X, typename LFSV,
128 typename Y>
130 ( const EG& eg,
131 const LFSU& lfsu, const X& x, const LFSV& lfsv,
132 Y& y) const
133 {
134 typedef typename X::value_type D;
135 typedef typename Y::value_type R;
137 typedef typename ResidualVector::WeightedAccumulationView ResidualView;
138
139 const int m=lfsv.size();
140 const int n=lfsu.size();
141
142 X u(x);
143
144 // Notice that in general lfsv.size() != y.size()
145 ResidualVector down(y.size()),up(y.size());
146 ResidualView downview = down.weightedAccumulationView(y.weight());
147 ResidualView upview = up.weightedAccumulationView(y.weight());
148
149 asImp().alpha_volume_post_skeleton(eg,lfsu,u,lfsv,downview);
150 for (int j=0; j<n; j++) // loop over columns
151 {
152 up = 0.0;
153 D delta = epsilon*(1.0+std::abs(u(lfsu,j)));
154 u(lfsu,j) += delta;
155 asImp().alpha_volume_post_skeleton(eg,lfsu,u,lfsv,upview);
156 for (int i=0; i<m; i++)
157 y.rawAccumulate(lfsv,i,((up(lfsv,i)-down(lfsv,i))/delta)*x(lfsu,j));
158 u(lfsu,j) = x(lfsu,j);
159 }
160 }
161
162 private:
163 const double epsilon; // problem: this depends on data type R!
164 Imp& asImp () {return static_cast<Imp &> (*this);}
165 const Imp& asImp () const {return static_cast<const Imp &>(*this);}
166 };
167
169
179 template<typename Imp>
182 {
183 public:
184
185 // We need to reimport the name from the base class; otherwise clang plays stupid and refuses to
186 // find the overload in the base class
188
191 , epsilon(1e-7)
192 {}
193
194 NumericalJacobianApplySkeleton (double epsilon_)
196 , epsilon(epsilon_)
197 {}
198
200 template<typename IG, typename LFSU, typename X, typename LFSV,
201 typename Y>
203 ( const IG& ig,
204 const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
205 const LFSU& lfsu_n, const X& x_n, const LFSV& lfsv_n,
206 Y& y_s, Y& y_n) const
207 {
208 typedef typename X::value_type D;
209 typedef typename Y::value_type R;
211 typedef typename ResidualVector::WeightedAccumulationView ResidualView;
212
213 const int m_s=lfsv_s.size();
214 const int m_n=lfsv_n.size();
215 const int n_s=lfsu_s.size();
216 const int n_n=lfsu_n.size();
217
218 X u_s(x_s);
219 X u_n(x_n);
220
221 // Notice that in general lfsv_s.size() != y_s.size()
222 ResidualVector down_s(y_s.size()),up_s(y_s.size());
223 ResidualView downview_s = down_s.weightedAccumulationView(1.0);
224 ResidualView upview_s = up_s.weightedAccumulationView(1.0);
225
226 ResidualVector down_n(y_n.size()),up_n(y_n.size());
227 ResidualView downview_n = down_n.weightedAccumulationView(1.0);
228 ResidualView upview_n = up_n.weightedAccumulationView(1.0);
229
230 // base line
231 asImp().alpha_skeleton(ig,lfsu_s,u_s,lfsv_s,lfsu_n,u_n,lfsv_n,downview_s,
232 downview_n);
233
234 // jiggle in self
235 for (int j=0; j<n_s; j++)
236 {
237 up_s = 0.0;
238 up_n = 0.0;
239 D delta = epsilon*(1.0+std::abs(u_s(lfsu_s,j)));
240 u_s(lfsu_s,j) += delta;
241 asImp().alpha_skeleton(ig,lfsu_s,u_s,lfsv_s,lfsu_n,u_n,lfsv_n,upview_s,
242 upview_n);
243 for (int i=0; i<m_s; i++)
244 y_s.accumulate(lfsv_s,i,((up_s(lfsv_s,i)-down_s(lfsv_s,i))/delta)*x_s(lfsu_s,j));
245 for (int i=0; i<m_n; i++)
246 y_n.accumulate(lfsv_n,i,((up_n(lfsv_n,i)-down_n(lfsv_n,i))/delta)*x_s(lfsu_s,j));
247 u_s(lfsu_s,j) = x_s(lfsu_s,j);
248 }
249
250 // jiggle in neighbor
251 for (int j=0; j<n_n; j++)
252 {
253 up_s = 0.0;
254 up_n = 0.0;
255 D delta = epsilon*(1.0+std::abs(u_n(lfsu_n,j)));
256 u_n(lfsu_n,j) += delta;
257 asImp().alpha_skeleton(ig,lfsu_s,u_s,lfsv_s,lfsu_n,u_n,lfsv_n,upview_s,
258 upview_n);
259 for (int i=0; i<m_s; i++)
260 y_s.accumulate(lfsv_s,i,((up_s(lfsv_s,i)-down_s(lfsv_s,i))/delta)*x_n(lfsu_n,j));
261 for (int i=0; i<m_n; i++)
262 y_n.accumulate(lfsv_n,i,((up_n(lfsv_n,i)-down_n(lfsv_n,i))/delta)*x_n(lfsu_n,j));
263 u_n(lfsu_n,j) = x_n(lfsu_n,j);
264 }
265 }
266
267 private:
268 const double epsilon; // problem: this depends on data type R!
269 Imp& asImp () { return static_cast<Imp &> (*this); }
270 const Imp& asImp () const { return static_cast<const Imp &>(*this); }
271 };
272
274
284 template<typename Imp>
287 {
288 public:
289
290 // We need to reimport the name from the base class; otherwise clang plays stupid and refuses to
291 // find the overload in the base class
293
296 , epsilon(1e-7)
297 {}
298
299 NumericalJacobianApplyBoundary (double epsilon_)
301 , epsilon(epsilon_)
302 {}
303
305 template<typename IG, typename LFSU, typename X, typename LFSV,
306 typename Y>
308 ( const IG& ig,
309 const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
310 Y& y_s) const
311 {
312 typedef typename X::value_type D;
313 typedef typename Y::value_type R;
315 typedef typename ResidualVector::WeightedAccumulationView ResidualView;
316
317 const int m_s=lfsv_s.size();
318 const int n_s=lfsu_s.size();
319
320 X u_s(x_s);
321
322 // Notice that in general lfsv_s.size() != y_s.size()
323 ResidualVector down_s(y_s.size()),up_s(y_s.size());
324 ResidualView downview_s = down_s.weightedAccumulationView(1.0);
325 ResidualView upview_s = up_s.weightedAccumulationView(1.0);
326
327 // base line
328 asImp().alpha_boundary(ig,lfsu_s,u_s,lfsv_s,downview_s);
329
330 // jiggle in self
331 for (int j=0; j<n_s; j++)
332 {
333 up_s = 0.0;
334 D delta = epsilon*(1.0+std::abs(u_s(lfsu_s,j)));
335 u_s(lfsu_s,j) += delta;
336 asImp().alpha_boundary(ig,lfsu_s,u_s,lfsv_s,upview_s);
337 for (int i=0; i<m_s; i++)
338 y_s.rawAccumulate(lfsv_s,i,((up_s(lfsv_s,i)-down_s(lfsv_s,i))/delta)*x_s(lfsu_s,j));
339 u_s(lfsu_s,j) = x_s(lfsu_s,j);
340 }
341 }
342
343 private:
344 const double epsilon; // problem: this depends on data type R!
345 Imp& asImp () { return static_cast<Imp &> (*this); }
346 const Imp& asImp () const { return static_cast<const Imp &>(*this); }
347 };
348
350 }
351}
352
353#endif // DUNE_PDELAB_LOCALOPERATOR_NUMERICALJACOBIANAPPLY_HH
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
Implements linear and nonlinear versions of jacobian_apply_boundary() based on alpha_boundary()
Definition: numericaljacobianapply.hh:287
void jacobian_apply_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, Y &y_s) const
apply local jacobian of the boundaryterm
Definition: numericaljacobianapply.hh:308
Implements linear and nonlinear versions of jacobian_apply_skeleton() based on alpha_skeleton()
Definition: numericaljacobianapply.hh:182
void jacobian_apply_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, Y &y_s, Y &y_n) const
apply local jacobian of the skeleton term
Definition: numericaljacobianapply.hh:203
Definition: numericaljacobianapply.hh:109
void jacobian_apply_volume_post_skeleton(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, Y &y) const
apply local jacobian of the volume term (post skeleton part)
Definition: numericaljacobianapply.hh:130
Implements linear and nonlinear versions of jacobian_apply_volume() based on alpha_volume()
Definition: numericaljacobianapply.hh:34
void jacobian_apply_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, Y &y) const
apply local jacobian of the volume term
Definition: numericaljacobianapply.hh:55
Implements nonlinear version of jacobian_apply_boundary() based on alpha_boundary()
Definition: numericalnonlinearjacobianapply.hh:249
Implements nonlinear version of jacobian_apply_skeleton() based on alpha_skeleton()
Definition: numericalnonlinearjacobianapply.hh:158
Definition: numericalnonlinearjacobianapply.hh:97
Implements nonlinear version of jacobian_apply_volume() based on alpha_volume()
Definition: numericalnonlinearjacobianapply.hh:34
Dune namespace.
Definition: alignedallocator.hh:13
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)