dune-pdelab  2.0.0
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 #include<dune/common/static_assert.hh>
10 
11 #include<dune/geometry/type.hh>
12 #include<dune/geometry/referenceelements.hh>
13 #include<dune/geometry/quadraturerules.hh>
14 
15 #include <dune/localfunctions/common/interfaceswitch.hh>
16 
17 #include"defaultimp.hh"
18 #include"pattern.hh"
19 #include"flags.hh"
20 #include"idefault.hh"
21 
22 namespace Dune {
23  namespace PDELab {
27 
34  class L2 : public NumericalJacobianApplyVolume<L2>,
35  public FullVolumePattern,
38  {
39  public:
40  // pattern assembly flags
41  enum { doPatternVolume = true };
42 
43  // residual assembly flags
44  enum { doAlphaVolume = true };
45 
46  L2 (int intorder_=2,double scaling=1.0)
47  : intorder(intorder_)
48  , _scaling(scaling)
49  {}
50 
51  // volume integral depending on test and ansatz functions
52  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
53  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
54  {
55  // Switches between local and global interface
56  typedef FiniteElementInterfaceSwitch<
57  typename LFSU::Traits::FiniteElementType
58  > FESwitch;
59  typedef BasisInterfaceSwitch<
60  typename FESwitch::Basis
61  > BasisSwitch;
62 
63  // domain and range field type
64  typedef typename BasisSwitch::DomainField DF;
65  typedef typename BasisSwitch::RangeField RF;
66  typedef typename BasisSwitch::Range RangeType;
67 
68  typedef typename LFSU::Traits::SizeType size_type;
69 
70  // dimensions
71  const int dim = EG::Geometry::dimension;
72 
73  // select quadrature rule
74  Dune::GeometryType gt = eg.geometry().type();
75  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
76 
77  // loop over quadrature points
78  for (typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
79  {
80  // evaluate basis functions
81  std::vector<RangeType> phi(lfsu.size());
82  FESwitch::basis(lfsu.finiteElement()).evaluateFunction(it->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 * it->weight() * eg.geometry().integrationElement(it->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::dimension;
117 
118  // select quadrature rule
119  Dune::GeometryType gt = eg.geometry().type();
120  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
121 
122  // loop over quadrature points
123  for (typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
124  {
125  // evaluate basis functions
126  std::vector<RangeType> phi(lfsu.size());
127  FESwitch::basis(lfsu.finiteElement()).evaluateFunction(it->position(),phi);
128 
129  // integrate phi_j*phi_i
130  RF factor = _scaling * it->weight() * eg.geometry().integrationElement(it->position());
131  for (size_type j=0; j<lfsu.size(); j++)
132  for (size_type i=0; i<lfsu.size(); i++)
133  mat.accumulate(lfsv,i,lfsu,j, phi[j]*phi[i]*factor);
134  }
135  }
136 
137  private:
138  int intorder;
139  const double _scaling;
140  };
141 
148  class PowerL2 : public NumericalJacobianApplyVolume<PowerL2>,
149  public FullVolumePattern,
152  {
153  public:
154  // pattern assembly flags
155  enum { doPatternVolume = true };
156 
157  // residual assembly flags
158  enum { doAlphaVolume = true };
159 
160  PowerL2 (int intorder_=2)
161  : scalar_operator(intorder_)
162  {}
163 
164  // volume integral depending on test and ansatz functions
165  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
166  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
167  {
168  for(unsigned int i=0; i<LFSU::CHILDREN; ++i)
169  scalar_operator.alpha_volume(eg,lfsu.child(i),x,lfsv.child(i),r);
170  }
171 
172  // jacobian of volume term
173  template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
174  void jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv,
175  M& mat) const
176  {
177  for(unsigned int i=0; i<LFSU::CHILDREN; ++i)
178  scalar_operator.jacobian_volume(eg,lfsu.child(i),x,lfsv.child(i),mat);
179  }
180 
181  private:
182  L2 scalar_operator;
183  };
184 
186  } // namespace PDELab
187 } // namespace Dune
188 
189 #endif
L2(int intorder_=2, double scaling=1.0)
Definition: l2.hh:46
Default flags for all local operators.
Definition: flags.hh:18
static const int dim
Definition: adaptivity.hh:82
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: l2.hh:174
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: l2.hh:53
PowerL2(int intorder_=2)
Definition: l2.hh:160
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
sparsity pattern generator
Definition: pattern.hh:14
Implement jacobian_apply_volume() based on alpha_volume()
Definition: defaultimp.hh:321
const EG & eg
Definition: common/constraints.hh:277
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: l2.hh:166
Definition: l2.hh:148
Definition: l2.hh:34
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: l2.hh:98