dune-pdelab  2.4-dev
vectorwave.hh
Go to the documentation of this file.
1 // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=8 sw=2 sts=2:
3 
4 #ifndef DUNE_PDELAB_LOCALOPERATOR_VECTORWAVE_HH
5 #define DUNE_PDELAB_LOCALOPERATOR_VECTORWAVE_HH
6 
7 #include <cstddef>
8 #include <vector>
9 
10 #include <dune/common/fvector.hh>
11 
12 #include <dune/geometry/type.hh>
13 #include <dune/geometry/quadraturerules.hh>
14 
20 
21 namespace Dune {
22  namespace PDELab {
26 
27  namespace VectorWave {
28 
30 
43  template<class Imp, class GV, class RF = double, class Time_ = double>
44  struct Parameters {
46  typedef Time_ Time;
47 
49  static const std::size_t dimension = GV::dimension;
50 
52  typedef typename GV::ctype DomainField;
54  typedef FieldVector<DomainField, dimension> Domain;
55 
57  typedef RF RangeField;
59  typedef FieldVector<RangeField, dimension> Range;
60 
62  typedef typename GV::template Codim<0>::Entity Element;
63 
66 
73  RangeField epsilon(const Element& e, const Domain& xl) const
74  { return asImp().epsilonGlobal(e.geometry().global(xl)); }
76 
81  bool epsilonChanged(Time t1, Time t2) const { return true; }
83 
90  RangeField mu(const Element& e, const Domain& xl) const
91  { return asImp().muGlobal(e.geometry().global(xl)); }
93 
98  bool muChanged(Time t1, Time t2) const { return true; }
99 
101 
104  void setTime(const Time &time) { }
105 
106  private:
107  const Imp &asImp() const { return *static_cast<const Imp*>(this); }
108  };
109 
111 
120  template<class GV, class RF = double, class Time = double>
122  public Parameters<ConstantParameters<GV, RF, Time>, GV, RF, Time>
123  {
124  RF epsilon_;
125  RF mu_;
126 
127  public:
128  ConstantParameters(RF epsilon, RF mu) : epsilon_(epsilon), mu_(mu) { }
129 
130  template<typename Domain>
131  RF epsilonGlobal(const Domain &) const { return epsilon_; }
132  bool epsilonChanged(Time t1, Time t2) const { return false; }
133  template<typename Domain>
134  RF muGlobal(const Domain &) const { return mu_; }
135  bool muChanged(Time t1, Time t2) const { return false; }
136  };
137 
140 
163  template<class Params>
164  class R0 :
165  public FullVolumePattern,
167  public InstationaryLocalOperatorDefaultMethods<typename Params::Time>,
168  public JacobianBasedAlphaVolume<R0<Params> >
169  {
171  IBase;
172 
173  Params &params;
174  std::size_t qorder;
175 
176  public:
177  // pattern assembly flags
178  enum { doPatternVolume = true };
179  enum { doAlphaVolume = true };
180 
182 
189  R0(Params &params_, std::size_t qorder_ = 0) :
190  params(params_), qorder(qorder_)
191  {}
192 
193  template<typename EG, typename LFSU, typename X, typename LFSV,
194  typename M>
195  void jacobian_volume(const EG& eg, const LFSU& lfsu, const X& x,
196  const LFSV& lfsv, M& mat) const
197  {
198  // domain and range field type
199  typedef typename LFSU::Traits::FiniteElementType FEU;
200  typedef typename LFSV::Traits::FiniteElementType FEV;
201 
202  const FEU &feU = lfsu.finiteElement();
203  const FEV &feV = lfsv.finiteElement();
204 
205  typedef typename FEU::Traits::Basis BasisU;
206  typedef typename FEV::Traits::Basis BasisV;
207 
208  const BasisU &basisU = feU.basis();
209  const BasisV &basisV = feV.basis();
210 
211  typedef typename Params::RangeField RF;
212 
213  typedef typename BasisU::Traits::DomainField DF;
214  static const std::size_t dimDLocal = BasisU::Traits::dimDomainLocal;
215 
216  typedef typename BasisU::Traits::Jacobian JacobianU;
217  typedef typename BasisV::Traits::Jacobian JacobianV;
218 
219  typedef JacobianToCurl<JacobianU> J2CU;
220  typedef JacobianToCurl<JacobianV> J2CV;
221 
222  static const J2CU &j2CU = J2CU();
223  static const J2CV &j2CV = J2CV();
224 
225  typedef typename J2CU::Curl CurlU;
226  typedef typename J2CV::Curl CurlV;
227 
228  static_assert(J2CU::dimCurl == J2CV::dimCurl, "Curl dimension "
229  "of ansatz and test functions must match in "
230  "VectorWave::R0");
231 
232  // select quadrature rule
233  typedef QuadratureRules<DF,dimDLocal> QRs;
234  typedef QuadratureRule<DF,dimDLocal> QR;
235  typedef typename QR::const_iterator QRIterator;
236  GeometryType gt = eg.geometry().type();
237  const QR& rule = QRs::rule(gt,qorder);
238 
239  std::vector<JacobianU> jacobianU(lfsu.size());
240  std::vector<JacobianV> jacobianV(lfsv.size());
241 
242  std::vector<CurlU> curlU(lfsu.size());
243  std::vector<CurlV> curlV(lfsv.size());
244 
245  // loop over quadrature points
246  for(QRIterator it=rule.begin(); it != rule.end(); ++it) {
247  // curls of basefunctions
248  basisU.evaluateJacobian(it->position(), jacobianU);
249  basisV.evaluateJacobian(it->position(), jacobianV);
250 
251  for(std::size_t j = 0; j < lfsu.size(); ++j)
252  j2CU(jacobianU[j], curlU[j]);
253  for(std::size_t i = 0; i < lfsv.size(); ++i)
254  j2CV(jacobianV[i], curlV[i]);
255 
256  RF factor = it->weight()
257  * eg.geometry().integrationElement(it->position())
258  / params.mu(eg.entity(), it->position());
259 
260  for(std::size_t j = 0; j < lfsu.size(); ++j)
261  for(std::size_t i = 0; i < lfsv.size(); ++i)
262  mat(i,j) += factor * (curlU[j] * curlV[i]);
263  }
264  }
265 
267  void setTime(typename Params::Time time) {
268  params.setTime(time);
269  IBase::setTime(time);
270  }
271  };
272 
275 
292  template<class Params>
293  class R1 :
295  public InstationaryLocalOperatorDefaultMethods<typename Params::Time>
296  {
297  public:
299 
306  R1(Params &params_, std::size_t qorder_ = 0) {}
307  };
308 
311 
332  template<class Params>
333  class R2 :
334  public FullVolumePattern,
336  public InstationaryLocalOperatorDefaultMethods<typename Params::Time>,
337  public JacobianBasedAlphaVolume<R2<Params> >
338  {
340  IBase;
341 
342  Params &params;
343  std::size_t qorder;
344 
345  public:
346  // pattern assembly flags
347  enum { doPatternVolume = true };
348  enum { doAlphaVolume = true };
349 
351 
358  R2(Params &params_, std::size_t qorder_ = 2)
359  : params(params_), qorder(qorder_)
360  {}
361 
362  template<typename EG, typename LFSU, typename X, typename LFSV,
363  typename M>
364  void jacobian_volume(const EG& eg, const LFSU& lfsu, const X& x,
365  const LFSV& lfsv, M& mat) const
366  {
367  // domain and range field type
368  typedef typename LFSU::Traits::FiniteElementType FEU;
369  typedef typename LFSV::Traits::FiniteElementType FEV;
370 
371  const FEU &feU = lfsu.finiteElement();
372  const FEV &feV = lfsv.finiteElement();
373 
374  typedef typename FEU::Traits::Basis BasisU;
375  typedef typename FEV::Traits::Basis BasisV;
376 
377  const BasisU &basisU = feU.basis();
378  const BasisV &basisV = feV.basis();
379 
380  typedef typename Params::RangeField RF;
381 
382  typedef typename BasisU::Traits::DomainField DF;
383  static const std::size_t dimDLocal = BasisU::Traits::dimDomainLocal;
384 
385  typedef typename BasisU::Traits::Range RangeU;
386  typedef typename BasisV::Traits::Range RangeV;
387 
388  // select quadrature rule
389  typedef QuadratureRules<DF,dimDLocal> QRs;
390  typedef QuadratureRule<DF,dimDLocal> QR;
391  typedef typename QR::const_iterator QRIterator;
392  GeometryType gt = eg.geometry().type();
393  const QR& rule = QRs::rule(gt,qorder);
394 
395  std::vector<RangeU> phiU(lfsu.size());
396  std::vector<RangeV> phiV(lfsv.size());
397 
398  // loop over quadrature points
399  for(QRIterator it=rule.begin(); it != rule.end(); ++it) {
400  // curls of basefunctions
401  basisU.evaluateFunction(it->position(), phiU);
402  basisV.evaluateFunction(it->position(), phiV);
403 
404  RF factor = it->weight()
405  * eg.geometry().integrationElement(it->position())
406  * params.epsilon(eg.entity(), it->position());
407 
408  for(std::size_t j = 0; j < lfsu.size(); ++j)
409  for(std::size_t i = 0; i < lfsv.size(); ++i)
410  mat(i,j) += factor * (phiU[j] * phiV[i]);
411  }
412  }
413 
415  void setTime(typename Params::Time time) {
416  params.setTime(time);
417  IBase::setTime(time);
418  }
419  };
420 
421  } // namespace VectorWave
422 
424  } // namespace PDELab
425 } // namespace Dune
426 
427 #endif // DUNE_PDELAB_LOCALOPERATOR_VECTORWAVE_HH
Local operator for the vector wave problem, no-temporal-derivatives part.
Definition: vectorwave.hh:164
RangeField mu(const Element &e, const Domain &xl) const
evaluate magnetic permeability
Definition: vectorwave.hh:90
Implement alpha_volume() based on jacobian_volume()
Definition: defaultimp.hh:611
static const std::size_t dimension
export dimension (both domain and range)
Definition: vectorwave.hh:49
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: vectorwave.hh:195
void setTime(R t_)
set time for subsequent evaluation
Definition: idefault.hh:104
bool muChanged(Time t1, Time t2) const
Definition: vectorwave.hh:135
const E & e
Definition: interpolate.hh:172
RangeField epsilon(const Element &e, const Domain &xl) const
evaluate dielectric permittivity
Definition: vectorwave.hh:73
RF muGlobal(const Domain &) const
Definition: vectorwave.hh:134
void setTime(typename Params::Time time)
set time on the function object
Definition: vectorwave.hh:415
void setTime(typename Params::Time time)
set time on the parameter object
Definition: vectorwave.hh:267
GV::ctype DomainField
field type of domain
Definition: vectorwave.hh:52
bool muChanged(Time t1, Time t2) const
Whether mu changes between the given time steps.
Definition: vectorwave.hh:98
Default flags for all local operators.
Definition: flags.hh:18
R2(Params &params_, std::size_t qorder_=2)
Construct a local operator object.
Definition: vectorwave.hh:358
sparsity pattern generator
Definition: pattern.hh:13
FieldVector< RangeField, dimension > Range
vector type of range
Definition: vectorwave.hh:59
R1(Params &params_, std::size_t qorder_=0)
Construct a local operator object.
Definition: vectorwave.hh:306
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: vectorwave.hh:364
Local operator for the vector wave problem, second-temporal-derivatives part.
Definition: vectorwave.hh:333
Time_ Time
export type of temporal values
Definition: vectorwave.hh:46
ConstantParameters(RF epsilon, RF mu)
Definition: vectorwave.hh:128
Definition: vectorwave.hh:348
Definition: vectorwave.hh:179
RF epsilonGlobal(const Domain &) const
Definition: vectorwave.hh:131
void setTime(const Time &time)
set the time for subsequent evaluation
Definition: vectorwave.hh:104
bool epsilonChanged(Time t1, Time t2) const
Definition: vectorwave.hh:132
Definition: adaptivity.hh:27
FieldVector< DomainField, dimension > Domain
vector type of domain
Definition: vectorwave.hh:54
RF RangeField
field type of range
Definition: vectorwave.hh:57
R0(Params &params_, std::size_t qorder_=0)
Construct a local operator object.
Definition: vectorwave.hh:189
Local operator for the vector wave problem, one-temporal-derivative part.
Definition: vectorwave.hh:293
extract the curl of a function from the jacobian of that function
Definition: jacobiantocurl.hh:27
Parameter class interface for the vector wave local operators.
Definition: vectorwave.hh:44
GV::template Codim< 0 >::Entity Element
element type of GridView
Definition: vectorwave.hh:62
bool epsilonChanged(Time t1, Time t2) const
Whether epsilon changes between the given time steps.
Definition: vectorwave.hh:81
const EG & eg
Definition: constraints.hh:280
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
Homogenous parameter class for the vector wave local operators.
Definition: vectorwave.hh:121