dune-pdelab  2.5-dev
l2.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; 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 : public NumericalJacobianApplyVolume<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  // 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,
95  M & mat) const
96  {
97  // Switches between local and global interface
98  using FESwitch = FiniteElementInterfaceSwitch<
99  typename LFSU::Traits::FiniteElementType>;
100  using BasisSwitch = BasisInterfaceSwitch<
101  typename FESwitch::Basis>;
102 
103  // Define types
104  using RangeType = typename BasisSwitch::Range;
105  using size_type = typename LFSU::Traits::SizeType;
106 
107  // Get geometry
108  auto geo = eg.geometry();
109 
110  // Inititialize vectors outside for loop
111  std::vector<RangeType> phi(lfsu.size());
112 
113  // Loop over quadrature points
114  for (const auto& qp : quadratureRule(geo,intorder))
115  {
116  // Evaluate basis functions
117  FESwitch::basis(lfsu.finiteElement()).evaluateFunction(qp.position(),phi);
118 
119  // Integrate phi_j*phi_i
120  auto factor = _scaling * qp.weight() * geo.integrationElement(qp.position());
121  for (size_type j=0; j<lfsu.size(); j++)
122  for (size_type i=0; i<lfsu.size(); i++)
123  mat.accumulate(lfsv,i,lfsu,j, phi[j]*phi[i]*factor);
124  }
125  }
126 
127  private:
128  int intorder;
129  const double _scaling;
130  };
131 
138  class PowerL2 : public NumericalJacobianApplyVolume<PowerL2>,
139  public FullVolumePattern,
142  {
143  public:
144  // Pattern assembly flags
145  enum { doPatternVolume = true };
146 
147  // Residual assembly flags
148  enum { doAlphaVolume = true };
149 
150  PowerL2 (int intorder_=2)
151  : scalar_operator(intorder_)
152  {}
153 
154  // Volume integral depending on test and ansatz functions
155  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
156  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
157  {
158  for(unsigned int i=0; i<TypeTree::degree(lfsu); ++i)
159  scalar_operator.alpha_volume(eg,lfsu.child(i),x,lfsv.child(i),r);
160  }
161 
162  // Jacobian of volume term
163  template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
164  void jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv,
165  M& mat) const
166  {
167  for(unsigned int i=0; i<TypeTree::degree(lfsu); ++i)
168  scalar_operator.jacobian_volume(eg,lfsu.child(i),x,lfsv.child(i),mat);
169  }
170 
171  private:
172  L2 scalar_operator;
173  };
174 
176  } // namespace PDELab
177 } // namespace Dune
178 
179 #endif // DUNE_PDELAB_LOCALOPERATOR_L2_HH
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
Definition: l2.hh:138
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
Definition: l2.hh:36
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: l2.hh:164
Implements linear and nonlinear versions of jacobian_apply_volume() based on alpha_volume() ...
Definition: numericaljacobianapply.hh:32
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: l2.hh:94
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: l2.hh:55
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
sparsity pattern generator
Definition: pattern.hh:13
L2(int intorder_=2, double scaling=1.0)
Definition: l2.hh:48
PowerL2(int intorder_=2)
Definition: l2.hh:150
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:156