dune-pdelab  2.5-dev
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 
16 #include <dune/geometry/quadraturerules.hh>
17 
21 
22 namespace Dune {
23  namespace PDELab {
27 
28  namespace ElectrodynamicImpl {
29 
30 
32  // Curl manipulation
33 
34  // dimension of the curl for a given space dimension
35  constexpr std::size_t dimOfCurl(std::size_t dimOfSpace)
36  {
37  return
38  dimOfSpace == 1 ? 2 :
39  dimOfSpace == 2 ? 1 :
40  dimOfSpace == 3 ? 3 :
41  // Dune exceptions are difficult to construct in constexpr functions
42  throw std::invalid_argument("Only applicable for dimensions 1-3");
43  }
44 
45  template<typename RF>
46  void jacobianToCurl(FieldVector<RF, 1> &curl,
47  const FieldMatrix<RF, 2, 2> &jacobian)
48  {
49  curl[0] = jacobian[1][0] - jacobian[0][1];
50  }
51  template<typename RF>
52  void jacobianToCurl(FieldVector<RF, 3> &curl,
53  const FieldMatrix<RF, 3, 3> &jacobian)
54  {
55  for(unsigned i = 0; i < 3; ++i)
56  curl[i] = jacobian[(i+2)%3][(i+1)%3] - jacobian[(i+1)%3][(i+2)%3];
57  }
58 
59  } // namespace ElectrodynamicImpl
60 
62 
71  template<typename Eps>
73  : public FullVolumePattern
75  , public JacobianBasedAlphaVolume<Electrodynamic_T<Eps> >
76  {
77 
78  public:
79 
80  // pattern assembly flags
81  static constexpr bool doPatternVolume = true;
82  static constexpr bool doAlphaVolume = true;
83 
85 
89  Electrodynamic_T(const Eps &eps, int qorder = 2) :
90  eps_(eps), qorder_(qorder)
91  {}
92 
93  Electrodynamic_T(Eps &&eps, int qorder = 2) :
94  eps_(std::move(eps)), qorder_(qorder)
95  {}
96 
100  template<typename EG, typename LFS, typename X, typename M>
101  void jacobian_volume (const EG& eg, const LFS& lfsu, const X& x,
102  const LFS& lfsv, M& mat) const
103  {
104  using BasisTraits =
105  typename LFS::Traits::FiniteElementType::Traits::Basis::Traits;
106 
107  // static checks
108  static constexpr unsigned dimR = BasisTraits::dimRange;
109  static_assert(dimR == 3 || dimR == 2, "Works only in 2D or 3D");
110 
111  using ctype = typename EG::Geometry::ctype;
112  using DF = typename BasisTraits::DomainField;
113  static_assert(std::is_same<ctype, DF>::value, "Grids ctype and "
114  "Finite Elements DomainFieldType must match");
115 
116  using Range = typename BasisTraits::Range;
117  std::vector<Range> phi(lfsu.size());
118 
119  // loop over quadrature points
120  for(const auto &qp : quadratureRule(eg.geometry(), qorder_))
121  {
122  // values of basefunctions
123  lfsu.finiteElement().basis().evaluateFunction(qp.position(),phi);
124 
125  // calculate T
126  auto factor = qp.weight()
127  * eg.geometry().integrationElement(qp.position())
128  * eps_(eg.entity(), qp.position());
129 
130  for(unsigned i = 0; i < lfsu.size(); ++i)
131  for(unsigned j = 0; j < lfsu.size(); ++j)
132  mat.accumulate(lfsv,i,lfsu,j,factor * (phi[i] * phi[j]));
133  }
134  }
135 
136  private:
137  Eps eps_;
138  int qorder_;
139  };
140 
142 
145  template<class Eps>
146  Electrodynamic_T<std::decay_t<Eps> >
147  makeLocalOperatorEdynT(Eps &&eps, int qorder = 2)
148  {
149  return { std::forward<Eps>(eps), qorder };
150  }
151 
153 
163  template<typename Mu>
165  : public FullVolumePattern
167  , public JacobianBasedAlphaVolume<Electrodynamic_S<Mu> >
168  {
169 
170  public:
171 
172  // pattern assembly flags
173  static constexpr bool doPatternVolume = true;
174  static constexpr bool doAlphaVolume = true;
175 
177 
181  Electrodynamic_S(const Mu &mu, int qorder = 2) :
182  mu_(mu), qorder_(qorder)
183  {}
184 
185  Electrodynamic_S(Mu &&mu, int qorder = 2) :
186  mu_(std::move(mu)), qorder_(qorder)
187  {}
188 
192  template<typename EG, typename LFS, typename X, typename M>
193  void jacobian_volume (const EG& eg, const LFS& lfsu, const X& x,
194  const LFS& lfsv, M& mat) const
195  {
198 
199  using BasisTraits =
200  typename LFS::Traits::FiniteElementType::Traits::Basis::Traits;
201 
202  // static checks
203  static constexpr unsigned dimR = BasisTraits::dimRange;
204  static_assert(dimR == 3 || dimR == 2, "Works only in 2D or 3D");
205 
206  using ctype = typename EG::Geometry::ctype;
207  using DF = typename BasisTraits::DomainField;
208  static_assert(std::is_same<ctype, DF>::value, "Grids ctype and "
209  "Finite Elements DomainFieldType must match");
210 
211  using Jacobian = typename BasisTraits::Jacobian;
212  std::vector<Jacobian> J(lfsu.size());
213 
214  using RF = typename BasisTraits::RangeField;
215  using Curl = FieldVector<RF, dimOfCurl(dimR)>;
216  std::vector<Curl> rotphi(lfsu.size());
217 
218  // loop over quadrature points
219  for(const auto &qp : quadratureRule(eg.geometry(), qorder_))
220  {
221  // curl of the basefunctions
222  lfsu.finiteElement().basis().evaluateJacobian(qp.position(),J);
223  for(unsigned i = 0; i < lfsu.size(); ++i)
224  jacobianToCurl(rotphi[i], J[i]);
225 
226  // calculate S
227  auto factor = qp.weight()
228  * eg.geometry().integrationElement(qp.position())
229  / mu_(eg.entity(), qp.position());
230 
231  for(unsigned i = 0; i < lfsu.size(); ++i)
232  for(unsigned j = 0; j < lfsu.size(); ++j)
233  mat.accumulate(lfsv,i,lfsu,j,factor * (rotphi[i] * rotphi[j]));
234 
235  }
236  }
237 
238  private:
239  Mu mu_;
240  int qorder_;
241  };
242 
244 
247  template<class Mu>
248  Electrodynamic_S<std::decay_t<Mu> >
249  makeLocalOperatorEdynS(Mu &&mu, int qorder = 2)
250  {
251  return { std::forward<Mu>(mu), qorder };
252  }
253 
255  } // namespace PDELab
256 } // namespace Dune
257 
258 #endif // DUNE_PDELAB_LOCALOPERATOR_ELECTRODYNAMIC_HH
Dune::PDELab::Electrodynamic_S::Electrodynamic_S
Electrodynamic_S(const Mu &mu, int qorder=2)
Construct an Electrodynamic_S localoperator.
Definition: electrodynamic.hh:181
flags.hh
Dune::PDELab::JacobianBasedAlphaVolume
Implement alpha_volume() based on jacobian_volume()
Definition: numericalresidual.hh:35
Dune::PDELab::Electrodynamic_S::Electrodynamic_S
Electrodynamic_S(Mu &&mu, int qorder=2)
Definition: electrodynamic.hh:185
Dune::PDELab::LocalOperatorDefaultFlags
Default flags for all local operators.
Definition: flags.hh:18
Dune
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
Dune::PDELab::ElectrodynamicImpl::jacobianToCurl
void jacobianToCurl(FieldVector< RF, 3 > &curl, const FieldMatrix< RF, 3, 3 > &jacobian)
Definition: electrodynamic.hh:52
numericalresidual.hh
Dune::PDELab::FullVolumePattern
sparsity pattern generator
Definition: pattern.hh:13
Dune::PDELab::makeLocalOperatorEdynS
Electrodynamic_S< std::decay_t< Mu > > makeLocalOperatorEdynS(Mu &&mu, int qorder=2)
construct an Electrodynamic_S operator
Definition: electrodynamic.hh:249
Dune::PDELab::Electrodynamic_T
Construct matrix T for the Electrodynamic operator.
Definition: electrodynamic.hh:72
Dune::PDELab::Electrodynamic_T::Electrodynamic_T
Electrodynamic_T(Eps &&eps, int qorder=2)
Definition: electrodynamic.hh:93
Dune::PDELab::ElectrodynamicImpl::jacobianToCurl
void jacobianToCurl(FieldVector< RF, 1 > &curl, const FieldMatrix< RF, 2, 2 > &jacobian)
Definition: electrodynamic.hh:46
Dune::PDELab::makeLocalOperatorEdynT
Electrodynamic_T< std::decay_t< Eps > > makeLocalOperatorEdynT(Eps &&eps, int qorder=2)
construct an Electrodynamic_T operator
Definition: electrodynamic.hh:147
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
Dune::PDELab::Electrodynamic_S::doPatternVolume
static constexpr bool doPatternVolume
Definition: electrodynamic.hh:173
pattern.hh
Dune::PDELab::ElectrodynamicImpl::dimOfCurl
constexpr std::size_t dimOfCurl(std::size_t dimOfSpace)
Definition: electrodynamic.hh:35
value
static const unsigned int value
Definition: gridfunctionspace/tags.hh:139
Dune::PDELab::Electrodynamic_T::jacobian_volume
void jacobian_volume(const EG &eg, const LFS &lfsu, const X &x, const LFS &lfsv, M &mat) const
Definition: electrodynamic.hh:101
Dune::PDELab::Electrodynamic_S
Contruct matrix S for the Electrodynamic operator.
Definition: electrodynamic.hh:164
Dune::PDELab::Electrodynamic_T::doPatternVolume
static constexpr bool doPatternVolume
Definition: electrodynamic.hh:81
Dune::PDELab::Electrodynamic_T::Electrodynamic_T
Electrodynamic_T(const Eps &eps, int qorder=2)
Construct an Electrodynamic_T localoperator.
Definition: electrodynamic.hh:89
Dune::PDELab::Electrodynamic_S::jacobian_volume
void jacobian_volume(const EG &eg, const LFS &lfsu, const X &x, const LFS &lfsv, M &mat) const
Definition: electrodynamic.hh:193
Dune::PDELab::Electrodynamic_S::doAlphaVolume
static constexpr bool doAlphaVolume
Definition: electrodynamic.hh:174
Dune::PDELab::Electrodynamic_T::doAlphaVolume
static constexpr bool doAlphaVolume
Definition: electrodynamic.hh:82