3 #ifndef DUNE_PDELAB_LOCALOPERATOR_DEFAULTIMP_HH
4 #define DUNE_PDELAB_LOCALOPERATOR_DEFAULTIMP_HH
31 template<
typename Imp>
44 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
48 const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
51 typedef typename X::value_type D;
52 typedef typename Jacobian::value_type R;
54 typedef typename ResidualVector::WeightedAccumulationView ResidualView;
56 const int m=lfsv.
size();
57 const int n=lfsu.size();
62 ResidualVector down(mat.nrows(),0.),up(mat.nrows());
63 ResidualView downview = down.weightedAccumulationView(mat.weight());
64 ResidualView upview = up.weightedAccumulationView(mat.weight());
67 asImp().alpha_volume(eg,lfsu,u,lfsv,downview);
68 for (
int j=0; j<n; j++)
71 D delta = epsilon*(1.0+std::abs(u(lfsu,j)));
73 asImp().alpha_volume(eg,lfsu,u,lfsv,upview);
74 for (
int i=0; i<m; i++)
75 mat.rawAccumulate(lfsv,i,lfsu,j,(up(lfsv,i)-down(lfsv,i))/delta);
76 u(lfsu,j) = x(lfsu,j);
82 Imp& asImp () {
return static_cast<Imp &
> (*this); }
83 const Imp& asImp ()
const {
return static_cast<const Imp &
>(*this); }
95 template<
typename Imp>
108 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
112 const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
115 typedef typename X::value_type D;
116 typedef typename Jacobian::value_type R;
118 typedef typename ResidualVector::WeightedAccumulationView ResidualView;
120 const int m=lfsv.
size();
121 const int n=lfsu.size();
126 ResidualVector down(mat.nrows(),0.),up(mat.nrows());
127 ResidualView downview = down.weightedAccumulationView(mat.weight());
128 ResidualView upview = up.weightedAccumulationView(mat.weight());
130 asImp().alpha_volume_post_skeleton(eg,lfsu,u,lfsv,downview);
131 for (
int j=0; j<n; j++)
134 D delta = epsilon*(1.0+std::abs(u(lfsu,j)));
136 asImp().alpha_volume_post_skeleton(eg,lfsu,u,lfsv,upview);
137 for (
int i=0; i<m; i++)
138 mat.rawAccumulate(lfsv,i,lfsu,j,(up(lfsv,i)-down(lfsv,i))/delta);
139 u(lfsu,j) = x(lfsu,j);
144 const double epsilon;
145 Imp& asImp () {
return static_cast<Imp &
> (*this); }
146 const Imp& asImp ()
const {
return static_cast<const Imp &
>(*this); }
156 template<
typename Imp>
169 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
173 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
174 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
175 Jacobian& mat_ss, Jacobian& mat_sn,
176 Jacobian& mat_ns, Jacobian& mat_nn)
const
178 typedef typename X::value_type D;
179 typedef typename Jacobian::value_type R;
181 typedef typename ResidualVector::WeightedAccumulationView ResidualView;
183 const int m_s=lfsv_s.
size();
184 const int m_n=lfsv_n.size();
185 const int n_s=lfsu_s.size();
186 const int n_n=lfsu_n.size();
192 ResidualVector down_s(mat_ss.nrows()),up_s(mat_ss.nrows());
193 ResidualView downview_s = down_s.weightedAccumulationView(1.0);
194 ResidualView upview_s = up_s.weightedAccumulationView(1.0);
196 ResidualVector down_n(mat_nn.nrows()),up_n(mat_nn.nrows());
197 ResidualView downview_n = down_n.weightedAccumulationView(1.0);
198 ResidualView upview_n = up_n.weightedAccumulationView(1.0);
201 asImp().alpha_skeleton(ig,lfsu_s,u_s,lfsv_s,lfsu_n,u_n,lfsv_n,downview_s,
205 for (
int j=0; j<n_s; j++)
209 D delta = epsilon*(1.0+std::abs(u_s(lfsu_s,j)));
210 u_s(lfsu_s,j) += delta;
211 asImp().alpha_skeleton(ig,lfsu_s,u_s,lfsv_s,lfsu_n,u_n,lfsv_n,upview_s,
213 for (
int i=0; i<m_s; i++)
214 mat_ss.accumulate(lfsv_s,i,lfsu_s,j,(up_s(lfsv_s,i)-down_s(lfsv_s,i))/delta);
215 for (
int i=0; i<m_n; i++)
216 mat_ns.accumulate(lfsv_n,i,lfsu_s,j,(up_n(lfsv_n,i)-down_n(lfsv_n,i))/delta);
217 u_s(lfsu_s,j) = x_s(lfsu_s,j);
221 for (
int j=0; j<n_n; j++)
225 D delta = epsilon*(1.0+std::abs(u_n(lfsu_n,j)));
226 u_n(lfsu_n,j) += delta;
227 asImp().alpha_skeleton(ig,lfsu_s,u_s,lfsv_s,lfsu_n,u_n,lfsv_n,upview_s,
229 for (
int i=0; i<m_s; i++)
230 mat_sn.accumulate(lfsv_s,i,lfsu_n,j,(up_s(lfsv_s,i)-down_s(lfsv_s,i))/delta);
231 for (
int i=0; i<m_n; i++)
232 mat_nn.accumulate(lfsv_n,i,lfsu_n,j,(up_n(lfsv_n,i)-down_n(lfsv_n,i))/delta);
233 u_n(lfsu_n,j) = x_n(lfsu_n,j);
238 const double epsilon;
239 Imp& asImp () {
return static_cast<Imp &
> (*this); }
240 const Imp& asImp ()
const {
return static_cast<const Imp &
>(*this); }
250 template<
typename Imp>
263 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
267 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
268 Jacobian& mat_ss)
const
270 typedef typename X::value_type D;
271 typedef typename Jacobian::value_type R;
273 typedef typename ResidualVector::WeightedAccumulationView ResidualView;
275 const int m_s=lfsv_s.
size();
276 const int n_s=lfsu_s.size();
281 ResidualVector down_s(mat_ss.nrows()),up_s(mat_ss.nrows());
282 ResidualView downview_s = down_s.weightedAccumulationView(mat_ss.weight());
283 ResidualView upview_s = up_s.weightedAccumulationView(mat_ss.weight());;
287 asImp().alpha_boundary(ig,lfsu_s,u_s,lfsv_s,downview_s);
290 for (
int j=0; j<n_s; j++)
293 D delta = epsilon*(1.0+std::abs(u_s(lfsu_s,j)));
294 u_s(lfsu_s,j) += delta;
295 asImp().alpha_boundary(ig,lfsu_s,u_s,lfsv_s,upview_s);
296 for (
int i=0; i<m_s; i++)
297 mat_ss.rawAccumulate(lfsv_s,i,lfsu_s,j,(up_s(lfsv_s,i)-down_s(lfsv_s,i))/delta);
298 u_s(lfsu_s,j) = x_s(lfsu_s,j);
303 const double epsilon;
304 Imp& asImp () {
return static_cast<Imp &
> (*this); }
305 const Imp& asImp ()
const {
return static_cast<const Imp &
>(*this); }
320 template<
typename Imp>
333 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
337 const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
340 typedef typename X::value_type D;
341 typedef typename Y::value_type R;
343 typedef typename ResidualVector::WeightedAccumulationView ResidualView;
345 const int m=lfsv.
size();
346 const int n=lfsu.size();
351 ResidualVector down(y.size()),up(y.size());
352 ResidualView downview = down.weightedAccumulationView(y.weight());
353 ResidualView upview = up.weightedAccumulationView(y.weight());
355 asImp().alpha_volume(eg,lfsu,u,lfsv,downview);
356 for (
int j=0; j<n; j++)
359 D delta = epsilon*(1.0+std::abs(u(lfsu,j)));
361 asImp().alpha_volume(eg,lfsu,u,lfsv,upview);
362 for (
int i=0; i<m; i++)
363 y.rawAccumulate(lfsv,i,((up(lfsv,i)-down(lfsv,i))/delta)*x(lfsu,j));
364 u(lfsu,j) = x(lfsu,j);
369 const double epsilon;
370 Imp& asImp () {
return static_cast<Imp &
> (*this); }
371 const Imp& asImp ()
const {
return static_cast<const Imp &
>(*this); }
383 template<
typename Imp>
396 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
400 const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
403 typedef typename X::value_type D;
404 typedef typename Y::value_type R;
406 typedef typename ResidualVector::WeightedAccumulationView ResidualView;
408 const int m=lfsv.
size();
409 const int n=lfsu.size();
414 ResidualVector down(y.size()),up(y.size());
415 ResidualView downview = down.weightedAccumulationView(y.weight());
416 ResidualView upview = up.weightedAccumulationView(y.weight());
418 asImp().alpha_volume_post_skeleton(eg,lfsu,u,lfsv,downview);
419 for (
int j=0; j<n; j++)
422 D delta = epsilon*(1.0+std::abs(u(lfsu,j)));
424 asImp().alpha_volume_post_skeleton(eg,lfsu,u,lfsv,upview);
425 for (
int i=0; i<m; i++)
426 y.rawAccumulate(lfsv,i,((up(lfsv,i)-down(lfsv,i))/delta)*x(lfsu,j));
427 u(lfsu,j) = x(lfsu,j);
432 const double epsilon;
433 Imp& asImp () {
return static_cast<Imp &
> (*this);}
434 const Imp& asImp ()
const {
return static_cast<const Imp &
>(*this);}
444 template<
typename Imp>
457 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
461 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
462 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
463 Y& y_s, Y& y_n)
const
465 typedef typename X::value_type D;
466 typedef typename Y::value_type R;
468 typedef typename ResidualVector::WeightedAccumulationView ResidualView;
470 const int m_s=lfsv_s.
size();
471 const int m_n=lfsv_n.size();
472 const int n_s=lfsu_s.size();
473 const int n_n=lfsu_n.size();
479 ResidualVector down_s(y_s.size()),up_s(y_s.size());
480 ResidualView downview_s = down_s.weightedAccumulationView(1.0);
481 ResidualView upview_s = up_s.weightedAccumulationView(1.0);
483 ResidualVector down_n(y_n.size()),up_n(y_n.size());
484 ResidualView downview_n = down_n.weightedAccumulationView(1.0);
485 ResidualView upview_n = up_n.weightedAccumulationView(1.0);
488 asImp().alpha_skeleton(ig,lfsu_s,u_s,lfsv_s,lfsu_n,u_n,lfsv_n,downview_s,
492 for (
int j=0; j<n_s; j++)
496 D delta = epsilon*(1.0+std::abs(u_s(lfsu_s,j)));
497 u_s(lfsu_s,j) += delta;
498 asImp().alpha_skeleton(ig,lfsu_s,u_s,lfsv_s,lfsu_n,u_n,lfsv_n,upview_s,
500 for (
int i=0; i<m_s; i++)
501 y_s.accumulate(lfsv_s,i,((up_s(lfsv_s,i)-down_s(lfsv_s,i))/delta)*x_s(lfsu_s,j));
502 for (
int i=0; i<m_n; i++)
503 y_n.accumulate(lfsv_n,i,((up_n(lfsv_n,i)-down_n(lfsv_n,i))/delta)*x_s(lfsu_s,j));
504 u_s(lfsu_s,j) = x_s(lfsu_s,j);
508 for (
int j=0; j<n_n; j++)
512 D delta = epsilon*(1.0+std::abs(u_n(lfsu_n,j)));
513 u_n(lfsu_n,j) += delta;
514 asImp().alpha_skeleton(ig,lfsu_s,u_s,lfsv_s,lfsu_n,u_n,lfsv_n,upview_s,
516 for (
int i=0; i<m_s; i++)
517 y_s.accumulate(lfsv_s,i,((up_s(lfsv_s,i)-down_s(lfsv_s,i))/delta)*x_n(lfsu_n,j));
518 for (
int i=0; i<m_n; i++)
519 y_n.accumulate(lfsv_n,i,((up_n(lfsv_n,i)-down_n(lfsv_n,i))/delta)*x_n(lfsu_n,j));
520 u_n(lfsu_n,j) = x_n(lfsu_n,j);
525 const double epsilon;
526 Imp& asImp () {
return static_cast<Imp &
> (*this); }
527 const Imp& asImp ()
const {
return static_cast<const Imp &
>(*this); }
537 template<
typename Imp>
550 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
554 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
557 typedef typename X::value_type D;
558 typedef typename Y::value_type R;
560 typedef typename ResidualVector::WeightedAccumulationView ResidualView;
562 const int m_s=lfsv_s.
size();
563 const int n_s=lfsu_s.size();
568 ResidualVector down_s(y_s.size()),up_s(y_s.size());
569 ResidualView downview_s = down_s.weightedAccumulationView(1.0);
570 ResidualView upview_s = up_s.weightedAccumulationView(1.0);
573 asImp().alpha_boundary(ig,lfsu_s,u_s,lfsv_s,downview_s);
576 for (
int j=0; j<n_s; j++)
579 D delta = epsilon*(1.0+std::abs(u_s(lfsu_s,j)));
580 u_s(lfsu_s,j) += delta;
581 asImp().alpha_boundary(ig,lfsu_s,u_s,lfsv_s,upview_s);
582 for (
int i=0; i<m_s; i++)
583 y_s.rawAccumulate(lfsv_s,i,((up_s(lfsv_s,i)-down_s(lfsv_s,i))/delta)*x_s(lfsu_s,j));
584 u_s(lfsu_s,j) = x_s(lfsu_s,j);
589 const double epsilon;
590 Imp& asImp () {
return static_cast<Imp &
> (*this); }
591 const Imp& asImp ()
const {
return static_cast<const Imp &
>(*this); }
610 template<
typename Imp>
616 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
620 const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
624 typedef typename Jacobian::WeightedAccumulationView JacobianView;
626 Jacobian mat(r.size(),x.size(), 0);
628 asImp().jacobian_volume(eg, lfsu, x, lfsv, matview);
630 mat.usmv(r.weight(),x,r);
634 Imp& asImp () {
return static_cast<Imp &
> (*this); }
635 const Imp& asImp ()
const {
return static_cast<const Imp &
>(*this); }
649 template<
typename Imp>
655 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
659 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
660 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
661 R& r_s, R& r_n)
const
664 typedef typename Jacobian::WeightedAccumulationView JacobianView;
666 Jacobian mat_ss(r_s.size(),x_s.size(),0);
667 Jacobian mat_sn(r_s.size(),x_n.size(),0);
668 Jacobian mat_ns(r_n.size(),x_s.size(),0);
669 Jacobian mat_nn(r_n.size(),x_n.size(),0);
676 asImp().jacobian_skeleton(ig,
679 view_ss, view_sn, view_ns, view_nn);
681 mat_ss.usmv(r_s.weight(),x_s,r_s);
682 mat_ns.usmv(r_n.weight(),x_s,r_n);
683 mat_sn.usmv(r_s.weight(),x_n,r_s);
684 mat_nn.usmv(r_n.weight(),x_n,r_n);
688 Imp& asImp () {
return static_cast<Imp &
> (*this); }
689 const Imp& asImp ()
const {
return static_cast<const Imp &
>(*this); }
703 template<
typename Imp>
709 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
713 const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
717 typedef typename Jacobian::WeightedAccumulationView JacobianView;
719 Jacobian mat(x.size(),r.size(), 0);
721 asImp().jacobian_boundary(ig, lfsu, x, lfsv, view);
723 mat.usmv(r.weight(),x,r);
727 Imp& asImp () {
return static_cast<Imp &
> (*this); }
728 const Imp& asImp ()
const {
return static_cast<const Imp &
>(*this); }
735 #endif // DUNE_PDELAB_LOCALOPERATOR_DEFAULTIMP_HH
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: defaultimp.hh:172
NumericalJacobianVolumePostSkeleton()
Definition: defaultimp.hh:99
Implement jacobian_apply_skeleton() based on alpha_skeleton()
Definition: defaultimp.hh:445
NumericalJacobianSkeleton()
Definition: defaultimp.hh:160
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: defaultimp.hh:460
Implement jacobian_boundary() based on alpha_boundary()
Definition: defaultimp.hh:251
NumericalJacobianVolumePostSkeleton(double epsilon_)
Definition: defaultimp.hh:103
Definition: defaultimp.hh:96
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: defaultimp.hh:336
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: defaultimp.hh:266
Implement jacobian_apply_boundary() based on alpha_boundary()
Definition: defaultimp.hh:538
Implement alpha_volume() based on jacobian_volume()
Definition: defaultimp.hh:611
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: defaultimp.hh:553
NumericalJacobianApplyBoundary()
Definition: defaultimp.hh:541
NumericalJacobianApplyVolumePostSkeleton()
Definition: defaultimp.hh:387
NumericalJacobianApplyVolume()
Definition: defaultimp.hh:324
NumericalJacobianVolume()
Definition: defaultimp.hh:35
void alpha_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, R &r_s, R &r_n) const
compute
Definition: defaultimp.hh:658
NumericalJacobianBoundary()
Definition: defaultimp.hh:254
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
compute
Definition: defaultimp.hh:619
Implement jacobian_volume() based on alpha_volume()
Definition: defaultimp.hh:32
NumericalJacobianBoundary(double epsilon_)
Definition: defaultimp.hh:258
size_type size() const
The size of the container.
Definition: localvector.hh:245
A dense matrix for storing data associated with the degrees of freedom of a pair of LocalFunctionSpac...
Definition: common/localmatrix.hh:184
const IG & ig
Definition: common/constraints.hh:146
Definition: defaultimp.hh:384
NumericalJacobianApplySkeleton(double epsilon_)
Definition: defaultimp.hh:452
NumericalJacobianApplyVolumePostSkeleton(double epsilon_)
Definition: defaultimp.hh:391
WeightedAccumulationView weightedAccumulationView(weight_type weight)
Returns a weighted accumulate-only view of this matrix with the given weight.
Definition: common/localmatrix.hh:398
NumericalJacobianVolume(double epsilon_)
Definition: defaultimp.hh:39
WeightedAccumulationView weightedAccumulationView(weight_type weight)
Returns a WeighedAccumulationView with some weight in addition to this view's weight.
Definition: common/localmatrix.hh:47
NumericalJacobianSkeleton(double epsilon_)
Definition: defaultimp.hh:164
void alpha_boundary(const IG &ig, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
compute
Definition: defaultimp.hh:712
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: defaultimp.hh:111
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: defaultimp.hh:399
Implement alpha_skeleton() based on jacobian_skeleton()
Definition: defaultimp.hh:650
A container for storing data associated with the degrees of freedom of a LocalFunctionSpace.
Definition: localvector.hh:171
NumericalJacobianApplySkeleton()
Definition: defaultimp.hh:448
Implement jacobian_apply_volume() based on alpha_volume()
Definition: defaultimp.hh:321
const EG & eg
Definition: common/constraints.hh:277
Implement jacobian_skeleton() based on alpha_skeleton()
Definition: defaultimp.hh:157
NumericalJacobianApplyBoundary(double epsilon_)
Definition: defaultimp.hh:545
const E & e
Definition: interpolate.hh:172
Implement alpha_boundary() based on jacobian_boundary()
Definition: defaultimp.hh:704
NumericalJacobianApplyVolume(double epsilon_)
Definition: defaultimp.hh:328
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: defaultimp.hh:47