dune-pdelab  2.0.0
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 <algorithm>
8 #include <cstddef>
9 #include <vector>
10 
11 #include <dune/common/fvector.hh>
12 #include <dune/common/static_assert.hh>
13 
14 #include <dune/geometry/type.hh>
15 #include <dune/geometry/quadraturerules.hh>
16 
23 
24 namespace Dune {
25  namespace PDELab {
29 
30  namespace VectorWave {
31 
33 
46  template<class Imp, class GV, class RF = double, class Time_ = double>
47  struct Parameters {
49  typedef Time_ Time;
50 
52  static const std::size_t dimension = GV::dimension;
53 
55  typedef typename GV::ctype DomainField;
57  typedef FieldVector<DomainField, dimension> Domain;
58 
60  typedef RF RangeField;
62  typedef FieldVector<RangeField, dimension> Range;
63 
65  typedef typename GV::template Codim<0>::Entity Element;
66 
69 
76  RangeField epsilon(const Element& e, const Domain& xl) const
77  { return asImp().epsilonGlobal(e.geometry().global(xl)); }
79 
84  bool epsilonChanged(Time t1, Time t2) const { return true; }
86 
93  RangeField mu(const Element& e, const Domain& xl) const
94  { return asImp().muGlobal(e.geometry().global(xl)); }
96 
101  bool muChanged(Time t1, Time t2) const { return true; }
102 
104 
107  void setTime(const Time &time) { }
108 
109  private:
110  const Imp &asImp() const { return *static_cast<const Imp*>(this); }
111  };
112 
114 
123  template<class GV, class RF = double, class Time = double>
125  public Parameters<ConstantParameters<GV, RF, Time>, GV, RF, Time>
126  {
127  RF epsilon_;
128  RF mu_;
129 
130  public:
131  ConstantParameters(RF epsilon, RF mu) : epsilon_(epsilon), mu_(mu) { }
132 
133  template<typename Domain>
134  RF epsilonGlobal(const Domain &) const { return epsilon_; }
135  bool epsilonChanged(Time t1, Time t2) const { return false; }
136  template<typename Domain>
137  RF muGlobal(const Domain &) const { return mu_; }
138  bool muChanged(Time t1, Time t2) const { return false; }
139  };
140 
143 
166  template<class Params>
167  class R0 :
168  public FullVolumePattern,
170  public InstationaryLocalOperatorDefaultMethods<typename Params::Time>,
171  public JacobianBasedAlphaVolume<R0<Params> >
172  {
174  IBase;
175 
176  Params &params;
177  std::size_t qorder;
178 
179  public:
180  // pattern assembly flags
181  enum { doPatternVolume = true };
182  enum { doAlphaVolume = true };
183 
185 
192  R0(Params &params_, std::size_t qorder_ = 0) :
193  params(params_), qorder(qorder_)
194  {}
195 
196  template<typename EG, typename LFSU, typename X, typename LFSV,
197  typename M>
198  void jacobian_volume(const EG& eg, const LFSU& lfsu, const X& x,
199  const LFSV& lfsv, M& mat) const
200  {
201  // domain and range field type
202  typedef typename LFSU::Traits::FiniteElementType FEU;
203  typedef typename LFSV::Traits::FiniteElementType FEV;
204 
205  const FEU &feU = lfsu.finiteElement();
206  const FEV &feV = lfsv.finiteElement();
207 
208  typedef typename FEU::Traits::Basis BasisU;
209  typedef typename FEV::Traits::Basis BasisV;
210 
211  const BasisU &basisU = feU.basis();
212  const BasisV &basisV = feV.basis();
213 
214  typedef typename Params::RangeField RF;
215 
216  typedef typename BasisU::Traits::DomainField DF;
217  static const std::size_t dimDLocal = BasisU::Traits::dimDomainLocal;
218 
219  typedef typename BasisU::Traits::Jacobian JacobianU;
220  typedef typename BasisV::Traits::Jacobian JacobianV;
221 
222  typedef JacobianToCurl<JacobianU> J2CU;
223  typedef JacobianToCurl<JacobianV> J2CV;
224 
225  static const J2CU &j2CU = J2CU();
226  static const J2CV &j2CV = J2CV();
227 
228  typedef typename J2CU::Curl CurlU;
229  typedef typename J2CV::Curl CurlV;
230 
231  dune_static_assert(J2CU::dimCurl == J2CV::dimCurl, "Curl dimension "
232  "of ansatz and test functions must match in "
233  "VectorWave::R0");
234 
235  // select quadrature rule
236  typedef QuadratureRules<DF,dimDLocal> QRs;
237  typedef QuadratureRule<DF,dimDLocal> QR;
238  typedef typename QR::const_iterator QRIterator;
239  GeometryType gt = eg.geometry().type();
240  const QR& rule = QRs::rule(gt,qorder);
241 
242  std::vector<JacobianU> jacobianU(lfsu.size());
243  std::vector<JacobianV> jacobianV(lfsv.size());
244 
245  std::vector<CurlU> curlU(lfsu.size());
246  std::vector<CurlV> curlV(lfsv.size());
247 
248  // loop over quadrature points
249  for(QRIterator it=rule.begin(); it != rule.end(); ++it) {
250  // curls of basefunctions
251  basisU.evaluateJacobian(it->position(), jacobianU);
252  basisV.evaluateJacobian(it->position(), jacobianV);
253 
254  for(std::size_t j = 0; j < lfsu.size(); ++j)
255  j2CU(jacobianU[j], curlU[j]);
256  for(std::size_t i = 0; i < lfsv.size(); ++i)
257  j2CV(jacobianV[i], curlV[i]);
258 
259  RF factor = it->weight()
260  * eg.geometry().integrationElement(it->position())
261  / params.mu(eg.entity(), it->position());
262 
263  for(std::size_t j = 0; j < lfsu.size(); ++j)
264  for(std::size_t i = 0; i < lfsv.size(); ++i)
265  mat(i,j) += factor * (curlU[j] * curlV[i]);
266  }
267  }
268 
270  void setTime(typename Params::Time time) {
271  params.setTime(time);
272  IBase::setTime(time);
273  }
274  };
275 
278 
295  template<class Params>
296  class R1 :
298  public InstationaryLocalOperatorDefaultMethods<typename Params::Time>
299  {
300  public:
302 
309  R1(Params &params_, std::size_t qorder_ = 0) {}
310  };
311 
314 
335  template<class Params>
336  class R2 :
337  public FullVolumePattern,
339  public InstationaryLocalOperatorDefaultMethods<typename Params::Time>,
340  public JacobianBasedAlphaVolume<R2<Params> >
341  {
343  IBase;
344 
345  Params &params;
346  std::size_t qorder;
347 
348  public:
349  // pattern assembly flags
350  enum { doPatternVolume = true };
351  enum { doAlphaVolume = true };
352 
354 
361  R2(Params &params_, std::size_t qorder_ = 2)
362  : params(params_), qorder(qorder_)
363  {}
364 
365  template<typename EG, typename LFSU, typename X, typename LFSV,
366  typename M>
367  void jacobian_volume(const EG& eg, const LFSU& lfsu, const X& x,
368  const LFSV& lfsv, M& mat) const
369  {
370  // domain and range field type
371  typedef typename LFSU::Traits::FiniteElementType FEU;
372  typedef typename LFSV::Traits::FiniteElementType FEV;
373 
374  const FEU &feU = lfsu.finiteElement();
375  const FEV &feV = lfsv.finiteElement();
376 
377  typedef typename FEU::Traits::Basis BasisU;
378  typedef typename FEV::Traits::Basis BasisV;
379 
380  const BasisU &basisU = feU.basis();
381  const BasisV &basisV = feV.basis();
382 
383  typedef typename Params::RangeField RF;
384 
385  typedef typename BasisU::Traits::DomainField DF;
386  static const std::size_t dimDLocal = BasisU::Traits::dimDomainLocal;
387 
388  typedef typename BasisU::Traits::Range RangeU;
389  typedef typename BasisV::Traits::Range RangeV;
390 
391  // select quadrature rule
392  typedef QuadratureRules<DF,dimDLocal> QRs;
393  typedef QuadratureRule<DF,dimDLocal> QR;
394  typedef typename QR::const_iterator QRIterator;
395  GeometryType gt = eg.geometry().type();
396  const QR& rule = QRs::rule(gt,qorder);
397 
398  std::vector<RangeU> phiU(lfsu.size());
399  std::vector<RangeV> phiV(lfsv.size());
400 
401  // loop over quadrature points
402  for(QRIterator it=rule.begin(); it != rule.end(); ++it) {
403  // curls of basefunctions
404  basisU.evaluateFunction(it->position(), phiU);
405  basisV.evaluateFunction(it->position(), phiV);
406 
407  RF factor = it->weight()
408  * eg.geometry().integrationElement(it->position())
409  * params.epsilon(eg.entity(), it->position());
410 
411  for(std::size_t j = 0; j < lfsu.size(); ++j)
412  for(std::size_t i = 0; i < lfsv.size(); ++i)
413  mat(i,j) += factor * (phiU[j] * phiV[i]);
414  }
415  }
416 
418  void setTime(typename Params::Time time) {
419  params.setTime(time);
420  IBase::setTime(time);
421  }
422  };
423 
425  template<class Params, class Step = int, class Time = double>
426  class CachePolicy :
427  public MultiStepCachePolicy<Step, Time>
428  {
430  const Params &params;
431 
432  protected:
433  using Base::currentStep;
434  using Base::endTime;
435  using Base::dt;
436 
437  public:
439  CachePolicy(const Params &params_) : params(params_) { }
440 
442  virtual bool isAffine(std::size_t order, Step step) const
443  { return true; }
445  virtual bool isComposedAffine(Step step) const
446  { return true; }
448  virtual bool hasPureLinearAlpha(std::size_t order, Step step) const
449  { return true; }
450  virtual bool canReuseJacobian(std::size_t order,
451  Step requested, Step available) const
452  {
453  Time t1 = endTime - dt*(currentStep-requested);
454  Time t2 = endTime - dt*(currentStep-available);
455  if(t1 > t2) std::swap(t1, t2);
456 
457  switch(order) {
458  case 0: return !params.muChanged(t1, t2);
459  case 1: return true;
460  case 2: return !params.epsilonChanged(t1, t2);
461  }
462 
463  DUNE_THROW(InvalidStateException,
464  "VectorWave::CachePolicy::canReuseJacobian(): Invalid "
465  "temporal derivative order=" << order);
466  }
467  virtual bool canReuseZeroResidual(std::size_t order,
468  Step requested, Step available) const
469  { return true; }
470  virtual bool canReuseComposedJacobian(Step requested,
471  Step available) const
472  {
473  Time t1 = endTime - dt*(currentStep-requested);
474  Time t2 = endTime - dt*(currentStep-available);
475  if(t1 > t2) std::swap(t1, t2);
476 
477  return !params.muChanged(t1, t2) && !params.epsilonChanged(t1, t2);
478  }
479  };
480 
481  } // namespace VectorWave
482 
484  } // namespace PDELab
485 } // namespace Dune
486 
487 #endif // DUNE_PDELAB_LOCALOPERATOR_VECTORWAVE_HH
Time_ Time
export type of temporal values
Definition: vectorwave.hh:49
RangeField mu(const Element &e, const Domain &xl) const
evaluate magnetic permeability
Definition: vectorwave.hh:93
bool epsilonChanged(Time t1, Time t2) const
Definition: vectorwave.hh:135
Local operator for the vector wave problem, one-temporal-derivative part.
Definition: vectorwave.hh:296
R1(Params &params_, std::size_t qorder_=0)
Construct a local operator object.
Definition: vectorwave.hh:309
ConstantParameters(RF epsilon, RF mu)
Definition: vectorwave.hh:131
Default flags for all local operators.
Definition: flags.hh:18
virtual bool canReuseJacobian(std::size_t order, Step requested, Step available) const
Whether a jacobian can be reused.
Definition: vectorwave.hh:450
CachePolicy(const Params &params_)
construct a CachePolicy object
Definition: vectorwave.hh:439
void setTime(typename Params::Time time)
set time on the function object
Definition: vectorwave.hh:418
void setTime(typename Params::Time time)
set time on the parameter object
Definition: vectorwave.hh:270
Time dt
Time step size.
Definition: cache.hh:75
Parameter class interface for the vector wave local operators.
Definition: vectorwave.hh:47
Implement alpha_volume() based on jacobian_volume()
Definition: defaultimp.hh:611
Definition: vectorwave.hh:182
RF RangeField
field type of range
Definition: vectorwave.hh:60
Step currentStep
The number of the step set via preStep().
Definition: cache.hh:61
virtual bool isComposedAffine(Step step) const
The composed operator is affine.
Definition: vectorwave.hh:445
bool muChanged(Time t1, Time t2) const
Definition: vectorwave.hh:138
MultiStepCachePolicy for VectorWave operators.
Definition: vectorwave.hh:426
Local operator for the vector wave problem, no-temporal-derivatives part.
Definition: vectorwave.hh:167
RF epsilonGlobal(const Domain &) const
Definition: vectorwave.hh:134
bool muChanged(Time t1, Time t2) const
Whether mu changes between the given time steps.
Definition: vectorwave.hh:101
FieldVector< RangeField, dimension > Range
vector type of range
Definition: vectorwave.hh:62
Policy class for the MultiStepCache.
Definition: cache.hh:46
virtual bool hasPureLinearAlpha(std::size_t order, Step step) const
All local operators have pure linear alpha_*()
Definition: vectorwave.hh:448
extract the curl of a function from the jacobian of that function
Definition: jacobiantocurl.hh:28
R0(Params &params_, std::size_t qorder_=0)
Construct a local operator object.
Definition: vectorwave.hh:192
GV::ctype DomainField
field type of domain
Definition: vectorwave.hh:55
bool epsilonChanged(Time t1, Time t2) const
Whether epsilon changes between the given time steps.
Definition: vectorwave.hh:84
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
Time endTime
The end-time of the step set via preStep().
Definition: cache.hh:68
GV::template Codim< 0 >::Entity Element
element type of GridView
Definition: vectorwave.hh:65
RangeField epsilon(const Element &e, const Domain &xl) const
evaluate dielectric permittivity
Definition: vectorwave.hh:76
Local operator for the vector wave problem, second-temporal-derivatives part.
Definition: vectorwave.hh:336
R2(Params &params_, std::size_t qorder_=2)
Construct a local operator object.
Definition: vectorwave.hh:361
static const std::size_t dimension
export dimension (both domain and range)
Definition: vectorwave.hh:52
sparsity pattern generator
Definition: pattern.hh:14
virtual bool canReuseZeroResidual(std::size_t order, Step requested, Step available) const
Whether a zero-residual can be reused.
Definition: vectorwave.hh:467
void setTime(const Time &time)
set the time for subsequent evaluation
Definition: vectorwave.hh:107
const EG & eg
Definition: common/constraints.hh:277
Homogenous parameter class for the vector wave local operators.
Definition: vectorwave.hh:124
Definition: vectorwave.hh:351
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: vectorwave.hh:198
RF muGlobal(const Domain &) const
Definition: vectorwave.hh:137
const E & e
Definition: interpolate.hh:172
FieldVector< DomainField, dimension > Domain
vector type of domain
Definition: vectorwave.hh:57
virtual bool canReuseComposedJacobian(Step requested, Step available) const
Whether a composed jacobian can be reused.
Definition: vectorwave.hh:470
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: vectorwave.hh:367
virtual bool isAffine(std::size_t order, Step step) const
All component operators are affine.
Definition: vectorwave.hh:442
void setTime(R t_)
set time for subsequent evaluation
Definition: idefault.hh:104