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 
31  // Deprecated interface handling
32 
33  // This is an adaptor that turns the old-style policy-object based
34  // parameters into a function-object based one.
35  //
36  // This is here for backward-compatibility only.
37  template<class Params>
38  class Functor
39  {
40  const Params *params_;
41 
42  public:
43  DUNE_DEPRECATED_MSG("You are using Dune::PDELab::Electrodynamic_T or "
44  "Dune::PDELab::Electrodynamic_S with an old-style "
45  "parameter class for eps/mu. Please use a "
46  "callable (function object, (generic) lambda, or "
47  "a function pointer) instead.")
48  Functor(const Params &params) : params_(&params) {}
49 
50  template<class Element, class Pos>
51  typename Params::Traits::RangeType
52  operator()(const Element &e, const Pos &xl) const
53  {
54  typename Params::Traits::RangeType result;
55  params_->evaluate(e, xl, result);
56  return result;
57  };
58  };
59 
60  template<class Params, class = void>
61  struct IsOldstyleParams : std::false_type {};
62 
63  template<class Params>
65  <Params,
66  std::enable_if_t<sizeof(typename Params::Traits::RangeType)> > :
67  std::true_type
68  {};
69 
70  template<class Params>
71  using ConstRefOrFunctor =
73  Functor<Params>, Params>;
74 
76  // Curl manipulation
77 
78  // dimension of the curl for a given space dimension
79  constexpr std::size_t dimOfCurl(std::size_t dimOfSpace)
80  {
81  return
82  dimOfSpace == 1 ? 2 :
83  dimOfSpace == 2 ? 1 :
84  dimOfSpace == 3 ? 3 :
85  // Dune exceptions are difficult to construct in constexpr functions
86  throw std::invalid_argument("Only applicable for dimensions 1-3");
87  }
88 
89  template<typename RF>
90  void jacobianToCurl(FieldVector<RF, 1> &curl,
91  const FieldMatrix<RF, 2, 2> &jacobian)
92  {
93  curl[0] = jacobian[1][0] - jacobian[0][1];
94  }
95  template<typename RF>
96  void jacobianToCurl(FieldVector<RF, 3> &curl,
97  const FieldMatrix<RF, 3, 3> &jacobian)
98  {
99  for(unsigned i = 0; i < 3; ++i)
100  curl[i] = jacobian[(i+2)%3][(i+1)%3] - jacobian[(i+1)%3][(i+2)%3];
101  }
102 
103  } // namespace ElectrodynamicImpl
104 
106 
115  template<typename Eps>
117  : public FullVolumePattern
119  , public JacobianBasedAlphaVolume<Electrodynamic_T<Eps> >
120  {
121  static constexpr bool oldstyle =
123 
124  public:
125 
126  // pattern assembly flags
127  static constexpr bool doPatternVolume = true;
128  static constexpr bool doAlphaVolume = true;
129 
131 
143  Electrodynamic_T(const Eps &eps, int qorder = 2) :
144  eps_(eps), qorder_(qorder)
145  {}
146 
147  Electrodynamic_T(Eps &&eps, int qorder = 2) :
148  eps_(std::move(eps)), qorder_(qorder)
149  {}
150 
154  template<typename EG, typename LFS, typename X, typename M>
155  void jacobian_volume (const EG& eg, const LFS& lfsu, const X& x,
156  const LFS& lfsv, M& mat) const
157  {
158  using BasisTraits =
159  typename LFS::Traits::FiniteElementType::Traits::Basis::Traits;
160 
161  // static checks
162  static constexpr unsigned dimR = BasisTraits::dimRange;
163  static_assert(dimR == 3 || dimR == 2, "Works only in 2D or 3D");
164 
165  using ctype = typename EG::Geometry::ctype;
166  using DF = typename BasisTraits::DomainField;
167  static_assert(std::is_same<ctype, DF>::value, "Grids ctype and "
168  "Finite Elements DomainFieldType must match");
169 
170  using Range = typename BasisTraits::Range;
171  std::vector<Range> phi(lfsu.size());
172 
173  // loop over quadrature points
174  for(const auto &qp : quadratureRule(eg.geometry(), qorder_))
175  {
176  // values of basefunctions
177  lfsu.finiteElement().basis().evaluateFunction(qp.position(),phi);
178 
179  // calculate T
180  auto factor = qp.weight()
181  * eg.geometry().integrationElement(qp.position())
182  * eps_(eg.entity(), qp.position());
183 
184  for(unsigned i = 0; i < lfsu.size(); ++i)
185  for(unsigned j = 0; j < lfsu.size(); ++j)
186  mat.accumulate(lfsv,i,lfsu,j,factor * (phi[i] * phi[j]));
187  }
188  }
189 
190  private:
192  int qorder_;
193  };
194 
196 
199  template<class Eps>
201  makeLocalOperatorEdynT(Eps &&eps, int qorder = 2)
202  {
203  return { std::forward<Eps>(eps), qorder };
204  }
205 
207 
217  template<typename Mu>
219  : public FullVolumePattern
221  , public JacobianBasedAlphaVolume<Electrodynamic_S<Mu> >
222  {
223  static constexpr bool oldstyle =
225 
226  public:
227 
228  // pattern assembly flags
229  static constexpr bool doPatternVolume = true;
230  static constexpr bool doAlphaVolume = true;
231 
233 
245  Electrodynamic_S(const Mu &mu, int qorder = 2) :
246  mu_(mu), qorder_(qorder)
247  {}
248 
249  Electrodynamic_S(Mu &&mu, int qorder = 2) :
250  mu_(std::move(mu)), qorder_(qorder)
251  {}
252 
256  template<typename EG, typename LFS, typename X, typename M>
257  void jacobian_volume (const EG& eg, const LFS& lfsu, const X& x,
258  const LFS& lfsv, M& mat) const
259  {
262 
263  using BasisTraits =
264  typename LFS::Traits::FiniteElementType::Traits::Basis::Traits;
265 
266  // static checks
267  static constexpr unsigned dimR = BasisTraits::dimRange;
268  static_assert(dimR == 3 || dimR == 2, "Works only in 2D or 3D");
269 
270  using ctype = typename EG::Geometry::ctype;
271  using DF = typename BasisTraits::DomainField;
272  static_assert(std::is_same<ctype, DF>::value, "Grids ctype and "
273  "Finite Elements DomainFieldType must match");
274 
275  using Jacobian = typename BasisTraits::Jacobian;
276  std::vector<Jacobian> J(lfsu.size());
277 
278  using RF = typename BasisTraits::RangeField;
279  using Curl = FieldVector<RF, dimOfCurl(dimR)>;
280  std::vector<Curl> rotphi(lfsu.size());
281 
282  // loop over quadrature points
283  for(const auto &qp : quadratureRule(eg.geometry(), qorder_))
284  {
285  // curl of the basefunctions
286  lfsu.finiteElement().basis().evaluateJacobian(qp.position(),J);
287  for(unsigned i = 0; i < lfsu.size(); ++i)
288  jacobianToCurl(rotphi[i], J[i]);
289 
290  // calculate S
291  auto factor = qp.weight()
292  * eg.geometry().integrationElement(qp.position())
293  / mu_(eg.entity(), qp.position());
294 
295  for(unsigned i = 0; i < lfsu.size(); ++i)
296  for(unsigned j = 0; j < lfsu.size(); ++j)
297  mat.accumulate(lfsv,i,lfsu,j,factor * (rotphi[i] * rotphi[j]));
298 
299  }
300  }
301 
302  private:
304  int qorder_;
305  };
306 
308 
311  template<class Mu>
313  makeLocalOperatorEdynS(Mu &&mu, int qorder = 2)
314  {
315  return { std::forward<Mu>(mu), qorder };
316  }
317 
319  } // namespace PDELab
320 } // namespace Dune
321 
322 #endif // DUNE_PDELAB_LOCALOPERATOR_ELECTRODYNAMIC_HH
const Entity & e
Definition: localfunctionspace.hh:111
Electrodynamic_S< std::decay_t< Mu > > makeLocalOperatorEdynS(Mu &&mu, int qorder=2)
construct an Electrodynamic_S operator
Definition: electrodynamic.hh:313
void jacobian_volume(const EG &eg, const LFS &lfsu, const X &x, const LFS &lfsv, M &mat) const
Definition: electrodynamic.hh:257
void jacobian_volume(const EG &eg, const LFS &lfsu, const X &x, const LFS &lfsv, M &mat) const
Definition: electrodynamic.hh:155
Electrodynamic_S(Mu &&mu, int qorder=2)
Definition: electrodynamic.hh:249
Construct matrix T for the Electrodynamic operator.
Definition: electrodynamic.hh:116
Params::Traits::RangeType operator()(const Element &e, const Pos &xl) const
Definition: electrodynamic.hh:52
Electrodynamic_S(const Mu &mu, int qorder=2)
Construct an Electrodynamic_S localoperator.
Definition: electrodynamic.hh:245
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
void jacobianToCurl(FieldVector< RF, 1 > &curl, const FieldMatrix< RF, 2, 2 > &jacobian)
Definition: electrodynamic.hh:90
Definition: electrodynamic.hh:38
STL namespace.
Electrodynamic_T< std::decay_t< Eps > > makeLocalOperatorEdynT(Eps &&eps, int qorder=2)
construct an Electrodynamic_T operator
Definition: electrodynamic.hh:201
Electrodynamic_T(Eps &&eps, int qorder=2)
Definition: electrodynamic.hh:147
Contruct matrix S for the Electrodynamic operator.
Definition: electrodynamic.hh:218
sparsity pattern generator
Definition: pattern.hh:13
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
Implement alpha_volume() based on jacobian_volume()
Definition: numericalresidual.hh:35
Electrodynamic_T(const Eps &eps, int qorder=2)
Construct an Electrodynamic_T localoperator.
Definition: electrodynamic.hh:143
constexpr std::size_t dimOfCurl(std::size_t dimOfSpace)
Definition: electrodynamic.hh:79
Default flags for all local operators.
Definition: flags.hh:18
static const unsigned int value
Definition: gridfunctionspace/tags.hh:139
std::conditional_t< IsOldstyleParams< Params >::value, Functor< Params >, Params > ConstRefOrFunctor
Definition: electrodynamic.hh:73