dune-pdelab  2.7-git
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 
21 namespace 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_T< std::decay_t< Eps > > makeLocalOperatorEdynT(Eps &&eps, int qorder=2)
construct an Electrodynamic_T operator
Definition: electrodynamic.hh:146
Electrodynamic_S< std::decay_t< Mu > > makeLocalOperatorEdynS(Mu &&mu, int qorder=2)
construct an Electrodynamic_S operator
Definition: electrodynamic.hh:248
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
void jacobianToCurl(FieldVector< RF, 3 > &curl, const FieldMatrix< RF, 3, 3 > &jacobian)
Definition: electrodynamic.hh:51
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
static const unsigned int value
Definition: gridfunctionspace/tags.hh:139