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/exceptions.hh>
8 #include<dune/common/fvector.hh>
9 
10 #include<dune/geometry/type.hh>
11 #include<dune/geometry/referenceelements.hh>
12 #include<dune/geometry/quadraturerules.hh>
13 
14 #include<dune/localfunctions/common/interfaceswitch.hh>
15 
18 
19 #include"defaultimp.hh"
20 #include"pattern.hh"
21 #include"flags.hh"
22 #include"idefault.hh"
23 
24 namespace Dune {
25  namespace PDELab {
29 
36  class L2 :
37  public FullVolumePattern,
40  {
41  public:
42  // Pattern assembly flags
43  enum { doPatternVolume = true };
44 
45  // Residual assembly flags
46  enum { doAlphaVolume = true };
47 
48  L2 (int intorder_=2,double scaling=1.0)
49  : intorder(intorder_)
50  , _scaling(scaling)
51  {}
52 
53  // Volume integral depending on test and ansatz functions
54  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
55  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
56  {
57  // Switches between local and global interface
58  using FESwitch = FiniteElementInterfaceSwitch<
59  typename LFSU::Traits::FiniteElementType>;
60  using BasisSwitch = BasisInterfaceSwitch<
61  typename FESwitch::Basis>;
62 
63  // Define types
64  using RF = typename BasisSwitch::RangeField;
65  using RangeType = typename BasisSwitch::Range;
66  using size_type = typename LFSU::Traits::SizeType;
67 
68  // Get geometry
69  auto geo = eg.geometry();
70 
71  // Initialize vectors outside for loop
72  std::vector<RangeType> phi(lfsu.size());
73 
74  // Loop over quadrature points
75  for (const auto& qp : quadratureRule(geo,intorder))
76  {
77  // Evaluate basis functions
78  FESwitch::basis(lfsu.finiteElement()).evaluateFunction(qp.position(),phi);
79 
80  // Evaluate u
81  RF u=0.0;
82  for (size_type i=0; i<lfsu.size(); i++)
83  u += RF(x(lfsu,i)*phi[i]);
84 
85  // u*phi_i
86  auto factor = _scaling * qp.weight() * geo.integrationElement(qp.position());
87  for (size_type i=0; i<lfsu.size(); i++)
88  r.accumulate(lfsv,i, u*phi[i]*factor);
89  }
90  }
91 
92  // apply jacobian of volume term
93  template<typename EG, typename LFSU, typename X, typename LFSV, typename Y>
94  void jacobian_apply_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, Y& y) const
95  {
96  alpha_volume(eg,lfsu,x,lfsv,y);
97  }
98 
99  // Jacobian of volume term
100  template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
101  void jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv,
102  M & mat) const
103  {
104  // Switches between local and global interface
105  using FESwitch = FiniteElementInterfaceSwitch<
106  typename LFSU::Traits::FiniteElementType>;
107  using BasisSwitch = BasisInterfaceSwitch<
108  typename FESwitch::Basis>;
109 
110  // Define types
111  using RangeType = typename BasisSwitch::Range;
112  using size_type = typename LFSU::Traits::SizeType;
113 
114  // Get geometry
115  auto geo = eg.geometry();
116 
117  // Inititialize vectors outside for loop
118  std::vector<RangeType> phi(lfsu.size());
119 
120  // Loop over quadrature points
121  for (const auto& qp : quadratureRule(geo,intorder))
122  {
123  // Evaluate basis functions
124  FESwitch::basis(lfsu.finiteElement()).evaluateFunction(qp.position(),phi);
125 
126  // Integrate phi_j*phi_i
127  auto factor = _scaling * qp.weight() * geo.integrationElement(qp.position());
128  for (size_type j=0; j<lfsu.size(); j++)
129  for (size_type i=0; i<lfsu.size(); i++)
130  mat.accumulate(lfsv,i,lfsu,j, phi[j]*phi[i]*factor);
131  }
132  }
133 
134  private:
135  int intorder;
136  const double _scaling;
137  };
138 
145  class PowerL2 :
146  public FullVolumePattern,
149  {
150  public:
151  // Pattern assembly flags
152  enum { doPatternVolume = true };
153 
154  // Residual assembly flags
155  enum { doAlphaVolume = true };
156 
157  PowerL2 (int intorder_=2)
158  : scalar_operator(intorder_)
159  {}
160 
161  // Volume integral depending on test and ansatz functions
162  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
163  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
164  {
165  for(unsigned int i=0; i<TypeTree::degree(lfsu); ++i)
166  scalar_operator.alpha_volume(eg,lfsu.child(i),x,lfsv.child(i),r);
167  }
168 
169  // apply jacobian of volume term
170  template<typename EG, typename LFSU, typename X, typename LFSV, typename Y>
171  void jacobian_apply_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, Y& y) const
172  {
173  alpha_volume(eg,lfsu,x,lfsv,y);
174  }
175 
176  // Jacobian of volume term
177  template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
178  void jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv,
179  M& mat) const
180  {
181  for(unsigned int i=0; i<TypeTree::degree(lfsu); ++i)
182  scalar_operator.jacobian_volume(eg,lfsu.child(i),x,lfsv.child(i),mat);
183  }
184 
185  private:
186  L2 scalar_operator;
187  };
188 
190  } // namespace PDELab
191 } // namespace Dune
192 
193 #endif // DUNE_PDELAB_LOCALOPERATOR_L2_HH
void jacobian_apply_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, Y &y) const
Definition: l2.hh:94
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: l2.hh:101
void jacobian_apply_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, Y &y) const
Definition: l2.hh:171
L2(int intorder_=2, double scaling=1.0)
Definition: l2.hh:48
Definition: l2.hh:36
sparsity pattern generator
Definition: pattern.hh:13
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
PowerL2(int intorder_=2)
Definition: l2.hh:157
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: l2.hh:178
Default flags for all local operators.
Definition: flags.hh:18
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: l2.hh:163
Definition: l2.hh:145
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: l2.hh:55