DUNE PDELab (git)

numericaljacobian.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_NUMERICALJACOBIAN_HH
4#define DUNE_PDELAB_LOCALOPERATOR_NUMERICALJACOBIAN_HH
5
6#include <cmath>
7
9#include <dune/pdelab/gridoperator/common/localmatrix.hh>
10
11namespace Dune {
12 namespace PDELab {
13
17
19 //
20 // Numerical implementation of jacobian_*() in terms of alpha_*()
21 //
22
24
30 template<typename Imp>
32 {
33 public:
35 : epsilon(1e-7)
36 {}
37
38 NumericalJacobianVolume (double epsilon_)
39 : epsilon(epsilon_)
40 {}
41
43 template<typename EG, typename LFSU, typename X, typename LFSV,
44 typename Jacobian>
46 ( const EG& eg,
47 const LFSU& lfsu, const X& x, const LFSV& lfsv,
48 Jacobian& mat) const
49 {
50 typedef typename X::value_type D;
51 typedef typename Jacobian::value_type R;
53 typedef typename ResidualVector::WeightedAccumulationView ResidualView;
54
55 const int m=lfsv.size();
56 const int n=lfsu.size();
57
58 X u(x);
59
60 // Notice that in general lfsv.size() != mat.nrows()
61 ResidualVector down(mat.nrows(),0.),up(mat.nrows());
62 ResidualView downview = down.weightedAccumulationView(mat.weight());
63 ResidualView upview = up.weightedAccumulationView(mat.weight());
64
65
66 asImp().alpha_volume(eg,lfsu,u,lfsv,downview);
67 for (int j=0; j<n; j++) // loop over columns
68 {
69 up = 0.0;
70 D delta = epsilon*(1.0+std::abs(u(lfsu,j)));
71 u(lfsu,j) += delta;
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);
76 }
77 }
78
79 private:
80 const double epsilon; // problem: this depends on data type R!
81 Imp& asImp () { return static_cast<Imp &> (*this); }
82 const Imp& asImp () const { return static_cast<const Imp &>(*this); }
83 };
84
87
94 template<typename Imp>
96 {
97 public:
99 : epsilon(1e-7)
100 {}
101
103 : epsilon(epsilon_)
104 {}
105
107 template<typename EG, typename LFSU, typename X, typename LFSV,
108 typename Jacobian>
110 ( const EG& eg,
111 const LFSU& lfsu, const X& x, const LFSV& lfsv,
112 Jacobian& mat) const
113 {
114 typedef typename X::value_type D;
115 typedef typename Jacobian::value_type R;
117 typedef typename ResidualVector::WeightedAccumulationView ResidualView;
118
119 const int m=lfsv.size();
120 const int n=lfsu.size();
121
122 X u(x);
123
124 // Notice that in general lfsv.size() != mat.nrows()
125 ResidualVector down(mat.nrows(),0.),up(mat.nrows());
126 ResidualView downview = down.weightedAccumulationView(mat.weight());
127 ResidualView upview = up.weightedAccumulationView(mat.weight());
128
129 asImp().alpha_volume_post_skeleton(eg,lfsu,u,lfsv,downview);
130 for (int j=0; j<n; j++) // loop over columns
131 {
132 up = 0.0;
133 D delta = epsilon*(1.0+std::abs(u(lfsu,j)));
134 u(lfsu,j) += delta;
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);
139 }
140 }
141
142 private:
143 const double epsilon; // problem: this depends on data type R!
144 Imp& asImp () { return static_cast<Imp &> (*this); }
145 const Imp& asImp () const { return static_cast<const Imp &>(*this); }
146 };
147
149
155 template<typename Imp>
157 {
158 public:
160 : epsilon(1e-7)
161 {}
162
163 NumericalJacobianSkeleton (double epsilon_)
164 : epsilon(epsilon_)
165 {}
166
168 template<typename IG, typename LFSU, typename X, typename LFSV,
169 typename Jacobian>
171 ( const IG& ig,
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
176 {
177 typedef typename X::value_type D;
178 typedef typename Jacobian::value_type R;
180 typedef typename ResidualVector::WeightedAccumulationView ResidualView;
181
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();
186
187 X u_s(x_s);
188 X u_n(x_n);
189
190 // Notice that in general lfsv.size() != mat.nrows()
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);
194
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);
198
199 // base line
200 asImp().alpha_skeleton(ig,lfsu_s,u_s,lfsv_s,lfsu_n,u_n,lfsv_n,downview_s,
201 downview_n);
202
203 // jiggle in self
204 for (int j=0; j<n_s; j++)
205 {
206 up_s = 0.0;
207 up_n = 0.0;
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,
211 upview_n);
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);
217 }
218
219 // jiggle in neighbor
220 for (int j=0; j<n_n; j++)
221 {
222 up_s = 0.0;
223 up_n = 0.0;
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,
227 upview_n);
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);
233 }
234 }
235
236 private:
237 const double epsilon; // problem: this depends on data type R!
238 Imp& asImp () { return static_cast<Imp &> (*this); }
239 const Imp& asImp () const { return static_cast<const Imp &>(*this); }
240 };
241
243
249 template<typename Imp>
251 {
252 public:
254 : epsilon(1e-7)
255 {}
256
257 NumericalJacobianBoundary (double epsilon_)
258 : epsilon(epsilon_)
259 {}
260
262 template<typename IG, typename LFSU, typename X, typename LFSV,
263 typename Jacobian>
265 ( const IG& ig,
266 const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
267 Jacobian& mat_ss) const
268 {
269 typedef typename X::value_type D;
270 typedef typename Jacobian::value_type R;
272 typedef typename ResidualVector::WeightedAccumulationView ResidualView;
273
274 const int m_s=lfsv_s.size();
275 const int n_s=lfsu_s.size();
276
277 X u_s(x_s);
278
279 // Notice that in general lfsv.size() != mat.nrows()
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());;
283
284
285 // base line
286 asImp().alpha_boundary(ig,lfsu_s,u_s,lfsv_s,downview_s);
287
288 // jiggle in self
289 for (int j=0; j<n_s; j++)
290 {
291 up_s = 0.0;
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);
298 }
299 }
300
301 private:
302 const double epsilon; // problem: this depends on data type R!
303 Imp& asImp () { return static_cast<Imp &> (*this); }
304 const Imp& asImp () const { return static_cast<const Imp &>(*this); }
305 };
306
308 }
309}
310
311#endif // DUNE_PDELAB_LOCALOPERATOR_NUMERICALJACOBIAN_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
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:13
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 23, 23:29, 2024)