dune-pdelab  2.0.0
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_ELECTRODYNAMIC_HH
4 #define DUNE_PDELAB_ELECTRODYNAMIC_HH
5 
6 #include<vector>
7 
8 #include<dune/common/exceptions.hh>
9 #include<dune/common/fmatrix.hh>
10 #include<dune/common/fvector.hh>
11 #include<dune/common/static_assert.hh>
12 #include<dune/common/typetraits.hh>
13 
14 #include<dune/geometry/type.hh>
15 #include<dune/geometry/quadraturerules.hh>
16 
17 #include"defaultimp.hh"
18 #include"pattern.hh"
19 #include"flags.hh"
20 
21 namespace Dune {
22  namespace PDELab {
26 
28 
37  template<typename Eps>
39  : public FullVolumePattern
41  , public JacobianBasedAlphaVolume<Electrodynamic_T<Eps> >
42  {
43  public:
44 
45  // pattern assembly flags
46  enum { doPatternVolume = true };
47 
48  enum { doAlphaVolume = true };
49 
51 
58  Electrodynamic_T(const Eps &eps_, int qorder_ = 2)
59  : eps(eps_)
60  , qorder(qorder_)
61  {}
62 
66  template<typename EG, typename LFS, typename X, typename M>
67  void jacobian_volume (const EG& eg, const LFS& lfsu, const X& x,
68  const LFS& lfsv, M& mat) const
69  {
70  // domain and range field type
71  typedef typename LFS::Traits::FiniteElementType::
72  Traits::Basis::Traits BasisTraits;
73 
74  typedef typename BasisTraits::DomainField DF;
75  static const unsigned dimDLocal = BasisTraits::dimDomainLocal;
76 
77  typedef typename BasisTraits::RangeField RF;
78  typedef typename BasisTraits::Range Range;
79  static const unsigned dimR = BasisTraits::dimRange;
80 
81  // static checks
82  dune_static_assert(dimR == 3 || dimR == 2,
83  "Works only in 2D or 3D");
84  dune_static_assert
86  "Grids ctype and Finite Elements DomainFieldType must match");
87 
88  // select quadrature rule
89  typedef Dune::QuadratureRule<DF,dimDLocal> QR;
90  typedef Dune::QuadratureRules<DF,dimDLocal> QRs;
91  Dune::GeometryType gt = eg.geometry().type();
92  const QR& rule = QRs::rule(gt,qorder);
93 
94  // loop over quadrature points
95  for(typename QR::const_iterator it=rule.begin();
96  it!=rule.end(); ++it) {
97  // values of basefunctions
98  std::vector<Range> phi(lfsu.size());
99  lfsu.finiteElement().basis().evaluateFunction(it->position(),phi);
100 
101  // calculate T
102  typename Eps::Traits::RangeType epsval;
103  eps.evaluate(eg.entity(), it->position(), epsval);
104 
105  RF factor = it->weight()
106  * eg.geometry().integrationElement(it->position()) * epsval;
107 
108  for(unsigned i = 0; i < lfsu.size(); ++i)
109  for(unsigned j = 0; j < lfsu.size(); ++j)
110  mat.accumulate(lfsv,i,lfsu,j,factor * (phi[i] * phi[j]));
111  }
112  }
113 
114  private:
115  const Eps &eps;
116  const int qorder;
117  };
118 
120 
130  template<typename Mu>
132  : public FullVolumePattern
134  , public JacobianBasedAlphaVolume<Electrodynamic_S<Mu> >
135  {
137  template <unsigned d>
138  struct CurlTraits {
139  static const unsigned dim =
140  d == 1 ? 2 :
141  d == 2 ? 1 :
142  /*else*/ 3;
143  };
144 
145  template<typename RF>
146  static void
147  jacobianToCurl(FieldVector<RF, 1> &curl,
148  const FieldMatrix<RF, 2, 2> &jacobian)
149  {
150  curl[0] = jacobian[1][0] - jacobian[0][1];
151  }
152  template<typename RF>
153  static void
154  jacobianToCurl(FieldVector<RF, 3> &curl,
155  const FieldMatrix<RF, 3, 3> &jacobian)
156  {
157  for(unsigned i = 0; i < 3; ++i)
158  curl[i] = jacobian[(i+2)%3][(i+1)%3] - jacobian[(i+1)%3][(i+2)%3];
159  }
160 
161  public:
162 
163  // pattern assembly flags
164  enum { doPatternVolume = true };
165 
166  enum { doAlphaVolume = true };
167 
169 
176  Electrodynamic_S(const Mu &mu_, int qorder_ = 2)
177  : mu(mu_)
178  , qorder(qorder_)
179  {}
180 
184  template<typename EG, typename LFS, typename X, typename M>
185  void jacobian_volume (const EG& eg, const LFS& lfsu, const X& x,
186  const LFS& lfsv, M& mat) const
187  {
188  // domain and range field type
189  typedef typename LFS::Traits::FiniteElementType::Traits::Basis::Traits
190  BasisTraits;
191 
192  typedef typename BasisTraits::DomainField DF;
193  static const unsigned dimDLocal = BasisTraits::dimDomainLocal;
194 
195  typedef typename BasisTraits::RangeField RF;
196  static const unsigned dimR = BasisTraits::dimRange;
197 
198  typedef typename BasisTraits::Jacobian Jacobian;
200 
201  // static checks
202  dune_static_assert(dimR == 3 || dimR == 2,
203  "Works only in 2D or 3D");
204  dune_static_assert
206  "Grids ctype and Finite Elements DomainFieldType must match");
207 
208  // select quadrature rule
209  typedef Dune::QuadratureRule<DF,dimDLocal> QR;
210  typedef Dune::QuadratureRules<DF,dimDLocal> QRs;
211  Dune::GeometryType gt = eg.geometry().type();
212  const QR& rule = QRs::rule(gt,qorder);
213 
214  // loop over quadrature points
215  for(typename QR::const_iterator it=rule.begin();
216  it!=rule.end(); ++it) {
217  // curl of the basefunctions
218  std::vector<Jacobian> J(lfsu.size());
219  lfsu.finiteElement().basis().evaluateJacobian(it->position(),J);
220 
221  std::vector<Curl> rotphi(lfsu.size());
222  for(unsigned i = 0; i < lfsu.size(); ++i)
223  jacobianToCurl(rotphi[i], J[i]);
224 
225  // calculate S
226  typename Mu::Traits::RangeType muval;
227  mu.evaluate(eg.entity(), it->position(), muval);
228 
229  RF factor = it->weight()
230  * eg.geometry().integrationElement(it->position()) / muval;
231 
232  for(unsigned i = 0; i < lfsu.size(); ++i)
233  for(unsigned j = 0; j < lfsu.size(); ++j)
234  mat.accumulate(lfsv,i,lfsu,j,factor * (rotphi[i] * rotphi[j]));
235 
236  }
237  }
238 
239  private:
240  const Mu &mu;
241  const int qorder;
242  };
243 
245  } // namespace PDELab
246 } // namespace Dune
247 
248 #endif // DUNE_PDELAB_ELECTRODYNAMIC_HH
void jacobian_volume(const EG &eg, const LFS &lfsu, const X &x, const LFS &lfsv, M &mat) const
Definition: electrodynamic.hh:185
Default flags for all local operators.
Definition: flags.hh:18
void jacobian_volume(const EG &eg, const LFS &lfsu, const X &x, const LFS &lfsv, M &mat) const
Definition: electrodynamic.hh:67
Definition: electrodynamic.hh:48
Contruct matrix T for the Electrodynamic operator.
Definition: electrodynamic.hh:38
static const int dim
Definition: adaptivity.hh:82
Implement alpha_volume() based on jacobian_volume()
Definition: defaultimp.hh:611
Electrodynamic_T(const Eps &eps_, int qorder_=2)
Construct an Electrodynamic_T localoperator.
Definition: electrodynamic.hh:58
Electrodynamic_S(const Mu &mu_, int qorder_=2)
Construct an Electrodynamic_S localoperator.
Definition: electrodynamic.hh:176
Contruct matrix S for the Electrodynamic operator.
Definition: electrodynamic.hh:131
Definition: electrodynamic.hh:166
static const unsigned int value
Definition: gridfunctionspace/tags.hh:175
sparsity pattern generator
Definition: pattern.hh:14
const EG & eg
Definition: common/constraints.hh:277
Definition: electrodynamic.hh:46
Definition: electrodynamic.hh:164