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