DUNE PDELab (git)

numericalnonlinearjacobianapply.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_NUMERICALNONLINEARJACOBIANAPPLY_HH
4#define DUNE_PDELAB_LOCALOPERATOR_NUMERICALNONLINEARJACOBIANAPPLY_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 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
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 nonlinear version of jacobian_apply_boundary() based on alpha_boundary()
Definition: numericalnonlinearjacobianapply.hh:249
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
Implements nonlinear version of jacobian_apply_skeleton() based on alpha_skeleton()
Definition: numericalnonlinearjacobianapply.hh:158
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
Definition: numericalnonlinearjacobianapply.hh:97
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
Implements nonlinear version of jacobian_apply_volume() based on alpha_volume()
Definition: numericalnonlinearjacobianapply.hh:34
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
Dune namespace.
Definition: alignedallocator.hh:13
STL namespace.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 23, 23:29, 2024)