dune-pdelab  2.4-dev
l2.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil -*-
2 #ifndef DUNE_PDELAB_L2_HH
3 #define DUNE_PDELAB_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 
16 #include"defaultimp.hh"
17 #include"pattern.hh"
18 #include"flags.hh"
19 #include"idefault.hh"
20 
21 namespace Dune {
22  namespace PDELab {
26 
33  class L2 : public NumericalJacobianApplyVolume<L2>,
34  public FullVolumePattern,
37  {
38  public:
39  // pattern assembly flags
40  enum { doPatternVolume = true };
41 
42  // residual assembly flags
43  enum { doAlphaVolume = true };
44 
45  L2 (int intorder_=2,double scaling=1.0)
46  : intorder(intorder_)
47  , _scaling(scaling)
48  {}
49 
50  // volume integral depending on test and ansatz functions
51  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
52  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
53  {
54  // Switches between local and global interface
55  typedef FiniteElementInterfaceSwitch<
56  typename LFSU::Traits::FiniteElementType
57  > FESwitch;
58  typedef BasisInterfaceSwitch<
59  typename FESwitch::Basis
60  > BasisSwitch;
61 
62  // domain and range field type
63  typedef typename BasisSwitch::DomainField DF;
64  typedef typename BasisSwitch::RangeField RF;
65  typedef typename BasisSwitch::Range RangeType;
66 
67  typedef typename LFSU::Traits::SizeType size_type;
68 
69  // dimensions
70  const int dim = EG::Geometry::mydimension;
71 
72  // select quadrature rule
73  const auto& geometry = eg.geometry();
74  Dune::GeometryType gt = geometry.type();
75  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
76 
77  // loop over quadrature points
78  for (const auto& qp : rule)
79  {
80  // evaluate basis functions
81  std::vector<RangeType> phi(lfsu.size());
82  FESwitch::basis(lfsu.finiteElement()).evaluateFunction(qp.position(),phi);
83 
84  // evaluate u
85  RF u=0.0;
86  for (size_type i=0; i<lfsu.size(); i++)
87  u += RF(x(lfsu,i)*phi[i]);
88 
89  // u*phi_i
90  RF factor = _scaling * qp.weight() * geometry.integrationElement(qp.position());
91  for (size_type i=0; i<lfsu.size(); i++)
92  r.accumulate(lfsv,i, u*phi[i]*factor);
93  }
94  }
95 
96  // jacobian of volume term
97  template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
98  void jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv,
99  M & mat) const
100  {
101  // Switches between local and global interface
102  typedef FiniteElementInterfaceSwitch<
103  typename LFSU::Traits::FiniteElementType
104  > FESwitch;
105  typedef BasisInterfaceSwitch<
106  typename FESwitch::Basis
107  > BasisSwitch;
108 
109  // domain and range field type
110  typedef typename BasisSwitch::DomainField DF;
111  typedef typename BasisSwitch::RangeField RF;
112  typedef typename BasisSwitch::Range RangeType;
113  typedef typename LFSU::Traits::SizeType size_type;
114 
115  // dimensions
116  const int dim = EG::Geometry::mydimension;
117 
118  // select quadrature rule
119  const auto& geometry = eg.geometry();
120  Dune::GeometryType gt = geometry.type();
121  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
122 
123  // loop over quadrature points
124  for (const auto& qp : rule)
125  {
126  // evaluate basis functions
127  std::vector<RangeType> phi(lfsu.size());
128  FESwitch::basis(lfsu.finiteElement()).evaluateFunction(qp.position(),phi);
129 
130  // integrate phi_j*phi_i
131  RF factor = _scaling * qp.weight() * geometry.integrationElement(qp.position());
132  for (size_type j=0; j<lfsu.size(); j++)
133  for (size_type i=0; i<lfsu.size(); i++)
134  mat.accumulate(lfsv,i,lfsu,j, phi[j]*phi[i]*factor);
135  }
136  }
137 
138  private:
139  int intorder;
140  const double _scaling;
141  };
142 
149  class PowerL2 : public NumericalJacobianApplyVolume<PowerL2>,
150  public FullVolumePattern,
153  {
154  public:
155  // pattern assembly flags
156  enum { doPatternVolume = true };
157 
158  // residual assembly flags
159  enum { doAlphaVolume = true };
160 
161  PowerL2 (int intorder_=2)
162  : scalar_operator(intorder_)
163  {}
164 
165  // volume integral depending on test and ansatz functions
166  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
167  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
168  {
169  for(unsigned int i=0; i<LFSU::CHILDREN; ++i)
170  scalar_operator.alpha_volume(eg,lfsu.child(i),x,lfsv.child(i),r);
171  }
172 
173  // jacobian of volume term
174  template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
175  void jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv,
176  M& mat) const
177  {
178  for(unsigned int i=0; i<LFSU::CHILDREN; ++i)
179  scalar_operator.jacobian_volume(eg,lfsu.child(i),x,lfsv.child(i),mat);
180  }
181 
182  private:
183  L2 scalar_operator;
184  };
185 
187  } // namespace PDELab
188 } // namespace Dune
189 
190 #endif
Implement jacobian_apply_volume() based on alpha_volume()
Definition: defaultimp.hh:321
L2(int intorder_=2, double scaling=1.0)
Definition: l2.hh:45
PowerL2(int intorder_=2)
Definition: l2.hh:161
Definition: l2.hh:149
Default flags for all local operators.
Definition: flags.hh:18
sparsity pattern generator
Definition: pattern.hh:13
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: l2.hh:167
Definition: adaptivity.hh:27
static const int dim
Definition: adaptivity.hh:83
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: l2.hh:52
Definition: l2.hh:33
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: l2.hh:175
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: l2.hh:98
const EG & eg
Definition: constraints.hh:280
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89