DUNE PDELab (git)

l2.hh
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
8
9#include <dune/localfunctions/common/interfaceswitch.hh>
10
11#include <dune/pdelab/common/quadraturerules.hh>
12
13#include <dune/pdelab/localoperator/defaultimp.hh>
14#include <dune/pdelab/localoperator/pattern.hh>
15#include <dune/pdelab/localoperator/flags.hh>
16#include <dune/pdelab/localoperator/idefault.hh>
17#include <dune/pdelab/localoperator/blockdiagonal.hh>
18
19namespace 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,
28 public LocalOperatorDefaultFlags,
29 public InstationaryLocalOperatorDefaultMethods<double>
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
Block diagonal extension of scalar local operator.
Definition: blockdiagonal.hh:78
Definition: l2.hh:151
L2(int intorderadd=0, double scaling=1.0)
Constructs a new L2 operator.
Definition: l2.hh:164
Implements a vector constructed from a given type representing a field and a compile-time given size.
Dune namespace.
Definition: alignedallocator.hh:13
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 23, 23:29, 2024)