dune-pdelab  2.5-dev
l2.hh
Go to the documentation of this file.
1 // -*- tab-width: 2; indent-tabs-mode: nil -*-
2 #ifndef DUNE_PDELAB_LOCALOPERATOR_L2_HH
3 #define DUNE_PDELAB_LOCALOPERATOR_L2_HH
4 
5 #include <vector>
6 
7 #include <dune/common/fvector.hh>
8 
9 #include <dune/localfunctions/common/interfaceswitch.hh>
10 
12 
18 
19 namespace Dune {
20  namespace PDELab {
21 
22  namespace impl {
23 
24  // Scalar L2 operator. Only for internal use! Use the L2 class instead,
25  // as that will also work for non-scalar spaces.
26  class ScalarL2 :
27  public FullVolumePattern,
30  {
31  public:
32  // Pattern assembly flags
33  enum { doPatternVolume = true };
34 
35  // Residual assembly flags
36  enum { doAlphaVolume = true };
37 
38  ScalarL2 (int intorderadd, double scaling)
39  : _intorderadd(intorderadd)
40  , _scaling(scaling)
41  {}
42 
43  // Volume integral depending on test and ansatz functions
44  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
45  void alpha_volume(const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
46  {
47  // Switches between local and global interface
48  using FESwitch = FiniteElementInterfaceSwitch<
49  typename LFSU::Traits::FiniteElementType>;
50  using BasisSwitch = BasisInterfaceSwitch<
51  typename FESwitch::Basis>;
52 
53  // Define types
54  using RF = typename BasisSwitch::RangeField;
55  using RangeType = typename BasisSwitch::Range;
56  using size_type = typename LFSU::Traits::SizeType;
57 
58  // Get geometry
59  auto geo = eg.geometry();
60 
61  // Initialize vectors outside for loop
62  std::vector<RangeType> phi(lfsu.size());
63 
64  // determine integration order
65  auto intorder = 2*lfsu.finiteElement().localBasis().order() + _intorderadd;
66 
67  // Loop over quadrature points
68  for (const auto& qp : quadratureRule(geo,intorder))
69  {
70  // Evaluate basis functions
71  FESwitch::basis(lfsu.finiteElement()).evaluateFunction(qp.position(),phi);
72 
73  // Evaluate u
74  RF u=0.0;
75  for (size_type i=0; i<lfsu.size(); i++)
76  u += RF(x(lfsu,i)*phi[i]);
77 
78  // u*phi_i
79  auto factor = _scaling * qp.weight() * geo.integrationElement(qp.position());
80  for (size_type i=0; i<lfsu.size(); i++)
81  r.accumulate(lfsv,i, u*phi[i]*factor);
82  }
83  }
84 
85  // apply jacobian of volume term
86  template<typename EG, typename LFSU, typename X, typename LFSV, typename Y>
87  void jacobian_apply_volume(const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, Y& y) const
88  {
89  alpha_volume(eg,lfsu,x,lfsv,y);
90  }
91 
92  // Jacobian of volume term
93  template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
94  void jacobian_volume(const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, M & mat) const
95  {
96  // Switches between local and global interface
97  using FESwitch = FiniteElementInterfaceSwitch<
98  typename LFSU::Traits::FiniteElementType>;
99  using BasisSwitch = BasisInterfaceSwitch<
100  typename FESwitch::Basis>;
101 
102  // Define types
103  using RangeType = typename BasisSwitch::Range;
104  using size_type = typename LFSU::Traits::SizeType;
105 
106  // Get geometry
107  auto geo = eg.geometry();
108 
109  // Inititialize vectors outside for loop
110  std::vector<RangeType> phi(lfsu.size());
111 
112  // determine integration order
113  auto intorder = 2*lfsu.finiteElement().localBasis().order() + _intorderadd;
114 
115  // Loop over quadrature points
116  for (const auto& qp : quadratureRule(geo,intorder))
117  {
118  // Evaluate basis functions
119  FESwitch::basis(lfsu.finiteElement()).evaluateFunction(qp.position(),phi);
120 
121  // Integrate phi_j*phi_i
122  auto factor = _scaling * qp.weight() * geo.integrationElement(qp.position());
123  for (size_type j=0; j<lfsu.size(); j++)
124  for (size_type i=0; i<lfsu.size(); i++)
125  mat.accumulate(lfsv,i,lfsu,j, phi[j]*phi[i]*factor);
126  }
127  }
128 
129  private:
130  int _intorderadd;
131  double _scaling;
132  };
133 
134  } // namespace impl
135 
139 
149  class L2 :
150  public BlockDiagonalLocalOperatorFullCoupling<impl::ScalarL2>
151  {
152 
153  public:
154 
156 
164  L2 (int intorderadd = 0, double scaling = 1.0)
165  : BlockDiagonalLocalOperatorFullCoupling<impl::ScalarL2>(intorderadd,scaling)
166  {}
167 
168  };
169 
176  class PowerL2 :
177  public FullVolumePattern,
180  {
181  public:
182  // Pattern assembly flags
183  enum { doPatternVolume = true };
184 
185  // Residual assembly flags
186  enum { doAlphaVolume = true };
187 
188  DUNE_DEPRECATED_MSG("PowerL2 is deprecated in PDELab 2.6 and will be removed after this release. The standard L2 operator now works on nested spaces, please use that instead.")
189  PowerL2 (int intorderadd = 2)
190  : scalar_operator(intorderadd)
191  {}
192 
193  // Volume integral depending on test and ansatz functions
194  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
195  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
196  {
197  for(unsigned int i=0; i<TypeTree::degree(lfsu); ++i)
198  scalar_operator.alpha_volume(eg,lfsu.child(i),x,lfsv.child(i),r);
199  }
200 
201  // apply jacobian of volume term
202  template<typename EG, typename LFSU, typename X, typename LFSV, typename Y>
203  void jacobian_apply_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, Y& y) const
204  {
205  alpha_volume(eg,lfsu,x,lfsv,y);
206  }
207 
208  // Jacobian of volume term
209  template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
210  void jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv,
211  M& mat) const
212  {
213  for(unsigned int i=0; i<TypeTree::degree(lfsu); ++i)
214  scalar_operator.jacobian_volume(eg,lfsu.child(i),x,lfsv.child(i),mat);
215  }
216 
217  private:
218  L2 scalar_operator;
219  };
220 
222  } // namespace PDELab
223 } // namespace Dune
224 
225 #endif // DUNE_PDELAB_LOCALOPERATOR_L2_HH
Dune::PDELab::impl::ScalarL2::jacobian_apply_volume
void jacobian_apply_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, Y &y) const
Definition: l2.hh:87
flags.hh
Dune::PDELab::L2
Definition: l2.hh:149
Dune::PDELab::PowerL2::DUNE_DEPRECATED_MSG
DUNE_DEPRECATED_MSG("PowerL2 is deprecated in PDELab 2.6 and will be removed after this release. The standard L2 operator now works on nested spaces, please use that instead.") PowerL2(int intorderadd
Dune::PDELab::InstationaryLocalOperatorDefaultMethods
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
Dune::PDELab::BlockDiagonalLocalOperatorFullCoupling::jacobian_volume
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: blockdiagonal.hh:128
Dune::PDELab::LocalOperatorDefaultFlags
Default flags for all local operators.
Definition: flags.hh:18
Dune
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
Dune::PDELab::BlockDiagonalLocalOperatorFullCoupling::alpha_volume
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: blockdiagonal.hh:104
Dune::PDELab::impl::ScalarL2
Definition: l2.hh:26
Dune::PDELab::impl::ScalarL2::doAlphaVolume
@ doAlphaVolume
Definition: l2.hh:36
defaultimp.hh
idefault.hh
Dune::PDELab::PowerL2::doPatternVolume
@ doPatternVolume
Definition: l2.hh:183
Dune::PDELab::FullVolumePattern
sparsity pattern generator
Definition: pattern.hh:13
Dune::PDELab::PowerL2::X
X
Definition: l2.hh:194
Dune::PDELab::L2::L2
L2(int intorderadd=0, double scaling=1.0)
Constructs a new L2 operator.
Definition: l2.hh:164
Dune::PDELab::PowerL2::doAlphaVolume
@ doAlphaVolume
Definition: l2.hh:186
Dune::PDELab::quadratureRule
QuadratureRuleWrapper< QuadratureRule< typename Geometry::ctype, Geometry::mydimension > > quadratureRule(const Geometry &geo, std::size_t order, QuadratureType::Enum quadrature_type=QuadratureType::GaussLegendre)
Returns a quadrature rule for the given geometry.
Definition: quadraturerules.hh:117
Dune::PDELab::PowerL2
Definition: l2.hh:176
Dune::PDELab::impl::ScalarL2::jacobian_volume
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: l2.hh:94
Dune::PDELab::PowerL2::LFSV
LFSV
Definition: l2.hh:194
Dune::PDELab::PowerL2::jacobian_volume
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: l2.hh:210
pattern.hh
Dune::PDELab::impl::ScalarL2::alpha_volume
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: l2.hh:45
Dune::PDELab::impl::ScalarL2::doPatternVolume
@ doPatternVolume
Definition: l2.hh:33
Dune::PDELab::PowerL2::alpha_volume
R void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: l2.hh:195
Dune::PDELab::PowerL2::jacobian_apply_volume
void jacobian_apply_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, Y &y) const
Definition: l2.hh:203
Dune::PDELab::BlockDiagonalLocalOperatorFullCoupling
Block diagonal extension of scalar local operator.
Definition: blockdiagonal.hh:75
quadraturerules.hh
Dune::PDELab::PowerL2::LFSU
LFSU
Definition: l2.hh:194
Dune::PDELab::impl::ScalarL2::ScalarL2
ScalarL2(int intorderadd, double scaling)
Definition: l2.hh:38
blockdiagonal.hh