dune-pdelab  2.7-git
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*FESwitch::basis(lfsu.finiteElement()).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*FESwitch::basis(lfsu.finiteElement()).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 
171  } // namespace PDELab
172 } // namespace Dune
173 
174 #endif // DUNE_PDELAB_LOCALOPERATOR_L2_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::L2
Definition: l2.hh:149
Dune::PDELab::impl::ScalarL2::ScalarL2
ScalarL2(int intorderadd, double scaling)
Definition: l2.hh:38
defaultimp.hh
Dune::PDELab::LocalOperatorDefaultFlags
Default flags for all local operators.
Definition: flags.hh:18
Dune::PDELab::InstationaryLocalOperatorDefaultMethods
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
Dune
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
Dune::PDELab::impl::ScalarL2::doAlphaVolume
@ doAlphaVolume
Definition: l2.hh:36
idefault.hh
pattern.hh
Dune::PDELab::BlockDiagonalLocalOperatorFullCoupling
Block diagonal extension of scalar local operator.
Definition: blockdiagonal.hh:75
Dune::PDELab::FullVolumePattern
sparsity pattern generator
Definition: pattern.hh:13
blockdiagonal.hh
Dune::PDELab::impl::ScalarL2::doPatternVolume
@ doPatternVolume
Definition: l2.hh:33
Dune::PDELab::impl::ScalarL2
Definition: l2.hh:26
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::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
flags.hh
Dune::PDELab::L2::L2
L2(int intorderadd=0, double scaling=1.0)
Constructs a new L2 operator.
Definition: l2.hh:164
quadraturerules.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