dune-pdelab 2.7-git
Loading...
Searching...
No Matches
electrodynamic.hh
Go to the documentation of this file.
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3#ifndef DUNE_PDELAB_LOCALOPERATOR_ELECTRODYNAMIC_HH
4#define DUNE_PDELAB_LOCALOPERATOR_ELECTRODYNAMIC_HH
5
6#include <cstddef>
7#include <stdexcept>
8#include <type_traits>
9#include <utility>
10#include <vector>
11
12#include <dune/common/deprecated.hh>
13#include <dune/common/fmatrix.hh>
14#include <dune/common/fvector.hh>
15
20
21namespace Dune {
22 namespace PDELab {
26
27 namespace ElectrodynamicImpl {
28
29
31 // Curl manipulation
32
33 // dimension of the curl for a given space dimension
34 constexpr std::size_t dimOfCurl(std::size_t dimOfSpace)
35 {
36 return
37 dimOfSpace == 1 ? 2 :
38 dimOfSpace == 2 ? 1 :
39 dimOfSpace == 3 ? 3 :
40 // Dune exceptions are difficult to construct in constexpr functions
41 throw std::invalid_argument("Only applicable for dimensions 1-3");
42 }
43
44 template<typename RF>
45 void jacobianToCurl(FieldVector<RF, 1> &curl,
46 const FieldMatrix<RF, 2, 2> &jacobian)
47 {
48 curl[0] = jacobian[1][0] - jacobian[0][1];
49 }
50 template<typename RF>
51 void jacobianToCurl(FieldVector<RF, 3> &curl,
52 const FieldMatrix<RF, 3, 3> &jacobian)
53 {
54 for(unsigned i = 0; i < 3; ++i)
55 curl[i] = jacobian[(i+2)%3][(i+1)%3] - jacobian[(i+1)%3][(i+2)%3];
56 }
57
58 } // namespace ElectrodynamicImpl
59
61
70 template<typename Eps>
72 : public FullVolumePattern
74 , public JacobianBasedAlphaVolume<Electrodynamic_T<Eps> >
75 {
76
77 public:
78
79 // pattern assembly flags
80 static constexpr bool doPatternVolume = true;
81 static constexpr bool doAlphaVolume = true;
82
84
88 Electrodynamic_T(const Eps &eps, int qorder = 2) :
89 eps_(eps), qorder_(qorder)
90 {}
91
92 Electrodynamic_T(Eps &&eps, int qorder = 2) :
93 eps_(std::move(eps)), qorder_(qorder)
94 {}
95
99 template<typename EG, typename LFS, typename X, typename M>
100 void jacobian_volume (const EG& eg, const LFS& lfsu, const X& x,
101 const LFS& lfsv, M& mat) const
102 {
103 using BasisTraits =
104 typename LFS::Traits::FiniteElementType::Traits::Basis::Traits;
105
106 // static checks
107 static constexpr unsigned dimR = BasisTraits::dimRange;
108 static_assert(dimR == 3 || dimR == 2, "Works only in 2D or 3D");
109
110 using ctype = typename EG::Geometry::ctype;
111 using DF = typename BasisTraits::DomainField;
112 static_assert(std::is_same<ctype, DF>::value, "Grids ctype and "
113 "Finite Elements DomainFieldType must match");
114
115 using Range = typename BasisTraits::Range;
116 std::vector<Range> phi(lfsu.size());
117
118 // loop over quadrature points
119 for(const auto &qp : quadratureRule(eg.geometry(), qorder_))
120 {
121 // values of basefunctions
122 lfsu.finiteElement().basis().evaluateFunction(qp.position(),phi);
123
124 // calculate T
125 auto factor = qp.weight()
126 * eg.geometry().integrationElement(qp.position())
127 * eps_(eg.entity(), qp.position());
128
129 for(unsigned i = 0; i < lfsu.size(); ++i)
130 for(unsigned j = 0; j < lfsu.size(); ++j)
131 mat.accumulate(lfsv,i,lfsu,j,factor * (phi[i] * phi[j]));
132 }
133 }
134
135 private:
136 Eps eps_;
137 int qorder_;
138 };
139
141
144 template<class Eps>
145 Electrodynamic_T<std::decay_t<Eps> >
146 makeLocalOperatorEdynT(Eps &&eps, int qorder = 2)
147 {
148 return { std::forward<Eps>(eps), qorder };
149 }
150
152
162 template<typename Mu>
164 : public FullVolumePattern
166 , public JacobianBasedAlphaVolume<Electrodynamic_S<Mu> >
167 {
168
169 public:
170
171 // pattern assembly flags
172 static constexpr bool doPatternVolume = true;
173 static constexpr bool doAlphaVolume = true;
174
176
180 Electrodynamic_S(const Mu &mu, int qorder = 2) :
181 mu_(mu), qorder_(qorder)
182 {}
183
184 Electrodynamic_S(Mu &&mu, int qorder = 2) :
185 mu_(std::move(mu)), qorder_(qorder)
186 {}
187
191 template<typename EG, typename LFS, typename X, typename M>
192 void jacobian_volume (const EG& eg, const LFS& lfsu, const X& x,
193 const LFS& lfsv, M& mat) const
194 {
197
198 using BasisTraits =
199 typename LFS::Traits::FiniteElementType::Traits::Basis::Traits;
200
201 // static checks
202 static constexpr unsigned dimR = BasisTraits::dimRange;
203 static_assert(dimR == 3 || dimR == 2, "Works only in 2D or 3D");
204
205 using ctype = typename EG::Geometry::ctype;
206 using DF = typename BasisTraits::DomainField;
207 static_assert(std::is_same<ctype, DF>::value, "Grids ctype and "
208 "Finite Elements DomainFieldType must match");
209
210 using Jacobian = typename BasisTraits::Jacobian;
211 std::vector<Jacobian> J(lfsu.size());
212
213 using RF = typename BasisTraits::RangeField;
214 using Curl = FieldVector<RF, dimOfCurl(dimR)>;
215 std::vector<Curl> rotphi(lfsu.size());
216
217 // loop over quadrature points
218 for(const auto &qp : quadratureRule(eg.geometry(), qorder_))
219 {
220 // curl of the basefunctions
221 lfsu.finiteElement().basis().evaluateJacobian(qp.position(),J);
222 for(unsigned i = 0; i < lfsu.size(); ++i)
223 jacobianToCurl(rotphi[i], J[i]);
224
225 // calculate S
226 auto factor = qp.weight()
227 * eg.geometry().integrationElement(qp.position())
228 / mu_(eg.entity(), qp.position());
229
230 for(unsigned i = 0; i < lfsu.size(); ++i)
231 for(unsigned j = 0; j < lfsu.size(); ++j)
232 mat.accumulate(lfsv,i,lfsu,j,factor * (rotphi[i] * rotphi[j]));
233
234 }
235 }
236
237 private:
238 Mu mu_;
239 int qorder_;
240 };
241
243
246 template<class Mu>
247 Electrodynamic_S<std::decay_t<Mu> >
248 makeLocalOperatorEdynS(Mu &&mu, int qorder = 2)
249 {
250 return { std::forward<Mu>(mu), qorder };
251 }
252
254 } // namespace PDELab
255} // namespace Dune
256
257#endif // DUNE_PDELAB_LOCALOPERATOR_ELECTRODYNAMIC_HH
Electrodynamic_S< std::decay_t< Mu > > makeLocalOperatorEdynS(Mu &&mu, int qorder=2)
construct an Electrodynamic_S operator
Definition: electrodynamic.hh:248
Electrodynamic_T< std::decay_t< Eps > > makeLocalOperatorEdynT(Eps &&eps, int qorder=2)
construct an Electrodynamic_T operator
Definition: electrodynamic.hh:146
STL namespace.
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
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
constexpr std::size_t dimOfCurl(std::size_t dimOfSpace)
Definition: electrodynamic.hh:34
void jacobianToCurl(FieldVector< RF, 1 > &curl, const FieldMatrix< RF, 2, 2 > &jacobian)
Definition: electrodynamic.hh:45
Construct matrix T for the Electrodynamic operator.
Definition: electrodynamic.hh:75
static constexpr bool doPatternVolume
Definition: electrodynamic.hh:80
void jacobian_volume(const EG &eg, const LFS &lfsu, const X &x, const LFS &lfsv, M &mat) const
Definition: electrodynamic.hh:100
static constexpr bool doAlphaVolume
Definition: electrodynamic.hh:81
Electrodynamic_T(Eps &&eps, int qorder=2)
Definition: electrodynamic.hh:92
Electrodynamic_T(const Eps &eps, int qorder=2)
Construct an Electrodynamic_T localoperator.
Definition: electrodynamic.hh:88
Contruct matrix S for the Electrodynamic operator.
Definition: electrodynamic.hh:167
static constexpr bool doAlphaVolume
Definition: electrodynamic.hh:173
static constexpr bool doPatternVolume
Definition: electrodynamic.hh:172
void jacobian_volume(const EG &eg, const LFS &lfsu, const X &x, const LFS &lfsv, M &mat) const
Definition: electrodynamic.hh:192
Electrodynamic_S(Mu &&mu, int qorder=2)
Definition: electrodynamic.hh:184
Electrodynamic_S(const Mu &mu, int qorder=2)
Construct an Electrodynamic_S localoperator.
Definition: electrodynamic.hh:180
Default flags for all local operators.
Definition: flags.hh:19
Implement alpha_volume() based on jacobian_volume()
Definition: numericalresidual.hh:36
sparsity pattern generator
Definition: pattern.hh:14