dune-pdelab  2.4-dev
maxwelldg.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil -*-
2 #ifndef DUNE_PDELAB_MAXWELLDG_HH
3 #define DUNE_PDELAB_MAXWELLDG_HH
4 
5 #include<vector>
6 
7 #include<dune/common/exceptions.hh>
8 #include<dune/common/fvector.hh>
9 
10 #include <dune/geometry/referenceelements.hh>
11 
19 
20 #include"maxwellparameter.hh"
21 
22 namespace Dune {
23  namespace PDELab {
24 
25 
26  template<int dim>
28  {};
29 
35  template<>
37  {
38  enum { dim = 3 };
39  public:
40 
50  template<typename T1, typename T2, typename T3>
51  static void eigenvalues (T1 eps, T1 mu, const Dune::FieldVector<T2,2*dim>& e)
52  {
53  T1 s = 1.0/sqrt(mu*eps); //speed of light s = 1/sqrt(\mu \epsilon)
54  e[0] = s;
55  e[1] = s;
56  e[2] = -s;
57  e[3] = -s;
58  e[4] = 0;
59  e[5] = 0;
60  }
61 
73  template<typename T1, typename T2, typename T3>
74  static void eigenvectors (T1 eps, T1 mu, const Dune::FieldVector<T2,dim>& n, Dune::FieldMatrix<T3,2*dim,2*dim>& R)
75  {
76  T1 a=n[0], b=n[1], c=n[2];
77 
78  Dune::FieldVector<T2,dim> re, im;
79  if (std::abs(c)<0.5)
80  {
81  re[0]=a*c; re[1]=b*c; re[2]=c*c-1;
82  im[0]=-b; im[1]=a; im[2]=0;
83  }
84  else
85  {
86  re[0]=a*b; re[1]=b*b-1; re[2]=b*c;
87  im[0]=c; im[1]=0.0; im[2]=-a;
88  }
89 
90  // \lambda_0,1 = s
91  R[0][0] = re[0]; R[0][1] = -im[0];
92  R[1][0] = re[1]; R[1][1] = -im[1];
93  R[2][0] = re[2]; R[2][1] = -im[2];
94  R[3][0] = im[0]; R[3][1] = re[0];
95  R[4][0] = im[1]; R[4][1] = re[1];
96  R[5][0] = im[2]; R[5][1] = re[2];
97 
98  // \lambda_2,3 = -s
99  R[0][2] = im[0]; R[0][3] = re[0];
100  R[1][2] = im[1]; R[1][3] = re[1];
101  R[2][2] = im[2]; R[2][3] = re[2];
102  R[3][2] = re[0]; R[3][3] = -im[0];
103  R[4][2] = re[1]; R[4][3] = -im[1];
104  R[5][2] = re[2]; R[5][3] = -im[2];
105 
106  // \lambda_4,5 = 0
107  R[0][4] = a; R[0][5] = 0;
108  R[1][4] = b; R[1][5] = 0;
109  R[2][4] = c; R[2][5] = 0;
110  R[3][4] = 0; R[3][5] = a;
111  R[4][4] = 0; R[4][5] = b;
112  R[5][4] = 0; R[5][5] = c;
113 
114  // apply scaling
115  T1 weps=sqrt(eps);
116  T1 wmu=sqrt(mu);
117  for (std::size_t i=0; i<3; i++)
118  for (std::size_t j=0; j<6; j++)
119  R[i][j] *= weps;
120  for (std::size_t i=3; i<6; i++)
121  for (std::size_t j=0; j<6; j++)
122  R[i][j] *= wmu;
123 
124  return;
125 
126  // if (std::abs(std::abs(c)-1)<1e-10)
127  // {
128  // if (c>0)
129  // {
130  // // \lambda_0,1 = s
131  // R[0][0] = 0; R[0][1] = 1;
132  // R[1][0] = -1; R[1][1] = 0;
133  // R[2][0] = 0; R[2][1] = 0;
134  // R[3][0] = 1; R[3][1] = 0;
135  // R[4][0] = 0; R[4][1] = 1;
136  // R[5][0] = 0; R[5][1] = 0;
137 
138  // // \lambda_2,3 = -s
139  // R[0][2] = -1; R[0][3] = 0;
140  // R[1][2] = 0; R[1][3] = 1;
141  // R[2][2] = 0; R[2][3] = 0;
142  // R[3][2] = 0; R[3][3] = 1;
143  // R[4][2] = 1; R[4][3] = 0;
144  // R[5][2] = 0; R[5][3] = 0;
145  // }
146  // else
147  // {
148  // // \lambda_0,1 = s
149  // R[0][0] = -1; R[0][1] = 0;
150  // R[1][0] = 0; R[1][1] = 1;
151  // R[2][0] = 0; R[2][1] = 0;
152  // R[3][0] = 0; R[3][1] = 1;
153  // R[4][0] = 1; R[4][1] = 0;
154  // R[5][0] = 0; R[5][1] = 0;
155 
156  // // \lambda_2,3 = -s
157  // R[0][2] = 0; R[0][3] = 1;
158  // R[1][2] = -1; R[1][3] = 0;
159  // R[2][2] = 0; R[2][3] = 0;
160  // R[3][2] = 1; R[3][3] = 0;
161  // R[4][2] = 0; R[4][3] = 1;
162  // R[5][2] = 0; R[5][3] = 0;
163  // }
164  // }
165  // else if (std::abs(std::abs(b)-1)<1e-10)
166  // {
167  // if (b>0)
168  // {
169  // // \lambda_0,1 = s
170  // R[0][0] = -1; R[0][1] = 0;
171  // R[1][0] = 0; R[1][1] = 0;
172  // R[2][0] = 0; R[2][1] = 1;
173  // R[3][0] = 0; R[3][1] = 1;
174  // R[4][0] = 0; R[4][1] = 0;
175  // R[5][0] = 1; R[5][1] = 0;
176 
177  // // \lambda_2,3 = -s
178  // R[0][2] = 0; R[0][3] = 1;
179  // R[1][2] = 0; R[1][3] = 0;
180  // R[2][2] = -1; R[2][3] = 0;
181  // R[3][2] = 1; R[3][3] = 0;
182  // R[4][2] = 0; R[4][3] = 0;
183  // R[5][2] = 0; R[5][3] = 1;
184  // }
185  // else
186  // {
187  // // \lambda_0,1 = s
188  // R[0][0] = 0; R[0][1] = 1;
189  // R[1][0] = 0; R[1][1] = 0;
190  // R[2][0] = -1; R[2][1] = 0;
191  // R[3][0] = 1; R[3][1] = 0;
192  // R[4][0] = 0; R[4][1] = 0;
193  // R[5][0] = 0; R[5][1] = 1;
194 
195  // // \lambda_2,3 = -s
196  // R[0][2] = -1; R[0][3] = 0;
197  // R[1][2] = 0; R[1][3] = 0;
198  // R[2][2] = 0; R[2][3] = 1;
199  // R[3][2] = 0; R[3][3] = 1;
200  // R[4][2] = 0; R[4][3] = 0;
201  // R[5][2] = 1; R[5][3] = 0;
202  // }
203  // }
204  // else if (std::abs(std::abs(a)-1)<1e-10)
205  // {
206  // if (a>0)
207  // {
208  // // \lambda_0,1 = s
209  // R[0][0] = 0; R[0][1] = 0;
210  // R[1][0] = 0; R[1][1] = 1;
211  // R[2][0] = -1; R[2][1] = 0;
212  // R[3][0] = 0; R[3][1] = 0;
213  // R[4][0] = 1; R[4][1] = 0;
214  // R[5][0] = 0; R[5][1] = 1;
215 
216  // // \lambda_2,3 = -s
217  // R[0][2] = 0; R[0][3] = 0;
218  // R[1][2] = -1; R[1][3] = 0;
219  // R[2][2] = 0; R[2][3] = 1;
220  // R[3][2] = 0; R[3][3] = 0;
221  // R[4][2] = 0; R[4][3] = 1;
222  // R[5][2] = 1; R[5][3] = 0;
223  // }
224  // else
225  // {
226  // // \lambda_0,1 = s
227  // R[0][0] = 0; R[0][1] = 0;
228  // R[1][0] = -1; R[1][1] = 0;
229  // R[2][0] = 0; R[2][1] = 1;
230  // R[3][0] = 0; R[3][1] = 0;
231  // R[4][0] = 0; R[4][1] = 1;
232  // R[5][0] = 1; R[5][1] = 0;
233 
234  // // \lambda_2,3 = -s
235  // R[0][2] = 0; R[0][3] = 0;
236  // R[1][2] = 0; R[1][3] = 1;
237  // R[2][2] = -1; R[2][3] = 0;
238  // R[3][2] = 0; R[3][3] = 0;
239  // R[4][2] = 1; R[4][3] = 0;
240  // R[5][2] = 0; R[5][3] = 1;
241  // }
242  // }
243  // else
244  // {
245  // DUNE_THROW(Dune::Exception,"need axiparallel grids for now!");
246 
247  // // \lambda_0,1 = s
248  // R[0][0] = b; R[0][1] = -(b*b+c*c);
249  // R[1][0] = -a; R[1][1] = a*b;
250  // R[2][0] = 0; R[2][1] = a*c;
251  // R[3][0] = a*c; R[3][1] = 0;
252  // R[4][0] = b*c; R[4][1] = -c;
253  // R[5][0] = -(a*a+b*b); R[5][1] = b;
254 
255  // // \lambda_2,3 = -s
256  // R[0][2] = -b; R[0][3] = -(b*b+c*c);
257  // R[1][2] = a; R[1][3] = a*b;
258  // R[2][2] = 0; R[2][3] = a*c;
259  // R[3][2] = a*c; R[3][3] = 0;
260  // R[4][2] = b*c; R[4][3] = c;
261  // R[5][2] = -(a*a+b*b); R[5][3] = -b;
262  // }
263 
264  // // \lambda_4,5 = 0
265  // R[0][4] = 0; R[0][5] = a;
266  // R[1][4] = 0; R[1][5] = b;
267  // R[2][4] = 0; R[2][5] = c;
268  // R[3][4] = a; R[3][5] = 0;
269  // R[4][4] = b; R[4][5] = 0;
270  // R[5][4] = c; R[5][5] = 0;
271 
272  // // apply scaling
273  // T1 weps=sqrt(eps);
274  // T1 wmu=sqrt(mu);
275  // for (std::size_t i=0; i<3; i++)
276  // for (std::size_t j=0; j<6; j++)
277  // R[i][j] *= weps;
278  // for (std::size_t i=3; i<6; i++)
279  // for (std::size_t j=0; j<6; j++)
280  // R[i][j] *= wmu;
281 
282  //std::cout << R << std::endl;
283 
284  }
285  };
286 
301  template<typename T, typename FEM>
303  public NumericalJacobianApplyVolume<DGMaxwellSpatialOperator<T,FEM> >,
304  public NumericalJacobianVolume<DGMaxwellSpatialOperator<T,FEM> >,
305  public NumericalJacobianApplySkeleton<DGMaxwellSpatialOperator<T,FEM> >,
306  public NumericalJacobianSkeleton<DGMaxwellSpatialOperator<T,FEM> >,
307  public NumericalJacobianApplyBoundary<DGMaxwellSpatialOperator<T,FEM> >,
308  public NumericalJacobianBoundary<DGMaxwellSpatialOperator<T,FEM> >,
309  public FullSkeletonPattern,
310  public FullVolumePattern,
312  public InstationaryLocalOperatorDefaultMethods<typename T::Traits::RangeFieldType>
313  {
314  enum { dim = T::Traits::GridViewType::dimension };
315 
316  public:
317  // pattern assembly flags
318  enum { doPatternVolume = true };
319  enum { doPatternSkeleton = true };
320 
321  // residual assembly flags
322  enum { doAlphaVolume = true };
323  enum { doAlphaSkeleton = true };
324  enum { doAlphaBoundary = true };
325  enum { doLambdaVolume = true };
326 
327  // ! constructor
328  DGMaxwellSpatialOperator (T& param_, int overintegration_=0)
329  : param(param_), overintegration(overintegration_), cache(20)
330  {
331  }
332 
333  // volume integral depending on test and ansatz functions
334  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
335  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
336  {
337  // get types
338  using DGSpace = typename LFSV::template Child<0>::Type;
339  using DF = typename DGSpace::Traits::FiniteElementType::Traits
340  ::LocalBasisType::Traits::DomainFieldType;
341  using RF = typename DGSpace::Traits::FiniteElementType::Traits
342  ::LocalBasisType::Traits::RangeFieldType;
343  using size_type = typename DGSpace::Traits::SizeType;
344 
345  // paranoia check number of number of components
346  static_assert(LFSV::CHILDREN == dim*2,
347  "need exactly dim*2 components!");
348 
349  // get local function space that is identical for all components
350  const auto& dgspace = lfsv.template child<0>();
351 
352  // select quadrature rule
353  const int order = dgspace.finiteElement().localBasis().order();
354  const int intorder = overintegration+2*order;
355  auto gt = eg.geometry().type();
356 
357  // intermediate results
358  std::vector<Dune::FieldVector<RF,dim> > gradphi(dgspace.size());
359 
360  // evaluate parameters (assumed constant per element)
361  const auto& localcenter =
362  Dune::ReferenceElements<DF,dim>::general(gt).position(0,0);
363  auto mu = param.mu(eg.entity(),localcenter);
364  auto eps = param.eps(eg.entity(),localcenter);
365  auto sigma = param.sigma(eg.entity(),localcenter);
366  auto muinv = 1.0/mu;
367  auto epsinv = 1.0/eps;
368 
369  //std::cout << "alpha_volume center=" << eg.geometry().center() << std::endl;
370 
371  // loop over quadrature points
372  for (const auto &qp : Dune::QuadratureRules<DF,dim>::rule(gt,intorder))
373  {
374  // evaluate basis functions
375  const auto& phi = cache[order].evaluateFunction
376  (qp.position(), dgspace.finiteElement().localBasis());
377 
378  // evaluate state vector u
379  Dune::FieldVector<RF,dim*2> u(0.0);
380  for (size_type k=0; k<dim*2; k++) // for all components
381  for (size_type j=0; j<dgspace.size(); j++) // for all basis functions
382  u[k] += x(lfsv.child(k),j)*phi[j];
383  //std::cout << " u at " << qp.position() << " : " << u << std::endl;
384 
385  // evaluate gradient of basis functions (we assume Galerkin method lfsu=lfsv)
386  const auto& js = cache[order].evaluateJacobian
387  (qp.position(), dgspace.finiteElement().localBasis());
388 
389  // compute global gradients
390  const auto &jac =
391  eg.geometry().jacobianInverseTransposed(qp.position());
392  for (size_type i=0; i<dgspace.size(); i++)
393  jac.mv(js[i][0],gradphi[i]);
394 
395  // integrate
396  RF factor = qp.weight()
397  * eg.geometry().integrationElement(qp.position());
398 
399  Dune::FieldMatrix<RF,dim*2,dim> F;
400  F[0][0] = 0; F[0][1] = -muinv*u[5]; F[0][2] = muinv*u[4];
401  F[1][0] = muinv*u[5]; F[1][1] = 0; F[1][2] = -muinv*u[3];
402  F[2][0] =-muinv*u[4]; F[2][1] = muinv*u[3]; F[2][2] = 0;
403  F[3][0] = 0; F[3][1] = epsinv*u[2]; F[3][2] = -epsinv*u[1];
404  F[4][0] = -epsinv*u[2]; F[4][1] = 0; F[4][2] = epsinv*u[0];
405  F[5][0] = epsinv*u[1]; F[5][1] = -epsinv*u[0]; F[5][2] = 0;
406 
407  // for all components of the system
408  for (size_type i=0; i<dim*2; i++)
409  // for all test functions of this component
410  for (size_type k=0; k<dgspace.size(); k++)
411  // for all dimensions
412  for (size_type j=0; j<dim; j++)
413  r.accumulate(lfsv.child(i),k,-F[i][j]*gradphi[k][j]*factor);
414 
415  // for the first half of the system
416  for (size_type i=0; i<dim; i++)
417  // for all test functions of this component
418  for (size_type k=0; k<dgspace.size(); k++)
419  r.accumulate(lfsv.child(i),k,(sigma/eps)*u[i]*phi[k]*factor);
420 
421  // std::cout << " residual: ";
422  // for (size_type i=0; i<r.size(); i++) std::cout << r[i] << " ";
423  // std::cout << std::endl;
424  }
425  }
426 
427  // skeleton integral depending on test and ansatz functions
428  // each face is only visited ONCE!
429  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
430  void alpha_skeleton (const IG& ig,
431  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
432  const LFSU& lfsu_n, const X& x_n, const LFSV& lfsv_n,
433  R& r_s, R& r_n) const
434  {
435  using std::max;
436  using std::sqrt;
437 
438  // get types
439  using DGSpace = typename LFSV::template Child<0>::Type;
440  using DF = typename DGSpace::Traits::FiniteElementType::
441  Traits::LocalBasisType::Traits::DomainFieldType;
442  using RF = typename DGSpace::Traits::FiniteElementType::
443  Traits::LocalBasisType::Traits::RangeFieldType;
444  using size_type = typename DGSpace::Traits::SizeType;
445 
446  // get local function space that is identical for all components
447  const auto& dgspace_s = lfsv_s.template child<0>();
448  const auto& dgspace_n = lfsv_n.template child<0>();
449 
450  // normal: assume faces are planar
451  const auto& n_F = ig.centerUnitOuterNormal();
452 
453  // evaluate speed of sound (assumed constant per element)
454  const auto& inside_local =
455  Dune::ReferenceElements<DF,dim>::general(ig.inside().type()).position(0,0);
456  const auto& outside_local =
457  Dune::ReferenceElements<DF,dim>::general(ig.outside().type()).position(0,0);
458  // This is wrong -- this approach (with A- and A+) does not allow
459  // position-dependent eps and mu, so we should not allow the parameter
460  // class to specify them in a position-dependent manner. See my
461  // (Jorrit Fahlke) dissertation on how to do it right.
462  auto mu_s = param.mu(ig.inside(),inside_local);
463  auto mu_n = param.mu(ig.outside(),outside_local);
464  auto eps_s = param.eps(ig.inside(),inside_local);
465  auto eps_n = param.eps(ig.outside(),outside_local);
466  //auto sigma_s = param.sigma(ig.inside(),inside_local);
467  //auto sigma_n = param.sigma(ig.outside(),outside_local);
468 
469  // compute A+ (outgoing waves)
470  Dune::FieldMatrix<DF,dim*2,dim*2> R_s;
471  MaxwellEigenvectors<dim>::eigenvectors(eps_s,mu_s,n_F,R_s);
472  Dune::FieldMatrix<DF,dim*2,dim*2> Dplus_s(0.0);
473  Dplus_s[0][0] = 1.0/sqrt(eps_s*mu_s);
474  Dplus_s[1][1] = 1.0/sqrt(eps_s*mu_s);
475  Dune::FieldMatrix<DF,dim*2,dim*2> Aplus_s(R_s);
476  Aplus_s.rightmultiply(Dplus_s);
477  R_s.invert();
478  Aplus_s.rightmultiply(R_s);
479 
480  // compute A- (incoming waves)
481  Dune::FieldMatrix<DF,dim*2,dim*2> R_n;
482  MaxwellEigenvectors<dim>::eigenvectors(eps_n,mu_n,n_F,R_n);
483  Dune::FieldMatrix<DF,dim*2,dim*2> Dminus_n(0.0);
484  Dminus_n[2][2] = -1.0/sqrt(eps_n*mu_n);
485  Dminus_n[3][3] = -1.0/sqrt(eps_n*mu_n);
486  Dune::FieldMatrix<DF,dim*2,dim*2> Aminus_n(R_n);
487  Aminus_n.rightmultiply(Dminus_n);
488  R_n.invert();
489  Aminus_n.rightmultiply(R_n);
490 
491  // select quadrature rule
492  const int order_s = dgspace_s.finiteElement().localBasis().order();
493  const int order_n = dgspace_n.finiteElement().localBasis().order();
494  const int intorder = overintegration+1+2*max(order_s,order_n);
495  auto gtface = ig.geometry().type();
496 
497  // std::cout << "alpha_skeleton center=" << ig.geometry().center() << std::endl;
498 
499  // loop over quadrature points
500  for (const auto &qp :
501  Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder))
502  {
503  // position of quadrature point in local coordinates of elements
504  const auto& iplocal_s = ig.geometryInInside().global(qp.position());
505  const auto& iplocal_n = ig.geometryInOutside().global(qp.position());
506 
507  // evaluate basis functions
508  const auto& phi_s = cache[order_s].evaluateFunction(iplocal_s,dgspace_s.finiteElement().localBasis());
509  const auto& phi_n = cache[order_n].evaluateFunction(iplocal_n,dgspace_n.finiteElement().localBasis());
510 
511  // evaluate u from inside and outside
512  Dune::FieldVector<RF,dim*2> u_s(0.0);
513  for (size_type i=0; i<dim*2; i++) // for all components
514  for (size_type k=0; k<dgspace_s.size(); k++) // for all basis functions
515  u_s[i] += x_s(lfsv_s.child(i),k)*phi_s[k];
516  Dune::FieldVector<RF,dim*2> u_n(0.0);
517  for (size_type i=0; i<dim*2; i++) // for all components
518  for (size_type k=0; k<dgspace_n.size(); k++) // for all basis functions
519  u_n[i] += x_n(lfsv_n.child(i),k)*phi_n[k];
520 
521  // compute numerical flux at integration point
522  Dune::FieldVector<RF,dim*2> f(0.0);
523  Aplus_s.umv(u_s,f);
524  // std::cout << " after A_plus*u_s " << f << std::endl;
525  Aminus_n.umv(u_n,f);
526  // std::cout << " after A_minus*u_n " << f << std::endl;
527 
528  //std::cout << "f=" << f << " u_s=" << u_s << " u_n=" << u_n << std::endl;
529 
530  // integrate
531  RF factor = qp.weight() * ig.geometry().integrationElement(qp.position());
532  for (size_type k=0; k<dgspace_s.size(); k++) // loop over all vector-valued (!) basis functions (with identical components)
533  for (size_type i=0; i<dim*2; i++) // loop over all components
534  r_s.accumulate(lfsv_s.child(i),k,f[i]*phi_s[k]*factor);
535  for (size_type k=0; k<dgspace_n.size(); k++) // loop over all vector-valued (!) basis functions (with identical components)
536  for (size_type i=0; i<dim*2; i++) // loop over all components
537  r_n.accumulate(lfsv_n.child(i),k,-f[i]*phi_n[k]*factor);
538  }
539 
540  // std::cout << " residual_s: ";
541  // for (auto v : r_s) std::cout << v << " ";
542  // std::cout << std::endl;
543  // std::cout << " residual_n: ";
544  // for (auto v : r_n) std::cout << v << " ";
545  // std::cout << std::endl;
546  }
547 
548  // skeleton integral depending on test and ansatz functions
549  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
550  void alpha_boundary (const IG& ig,
551  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
552  R& r_s) const
553  {
554  // get types
555  using DGSpace = typename LFSV::template Child<0>::Type;
556  using DF = typename DGSpace::Traits::FiniteElementType::
557  Traits::LocalBasisType::Traits::DomainFieldType;
558  using RF = typename DGSpace::Traits::FiniteElementType::
559  Traits::LocalBasisType::Traits::RangeFieldType;
560  using size_type = typename DGSpace::Traits::SizeType;
561 
562  // get local function space that is identical for all components
563  const auto& dgspace_s = lfsv_s.template child<0>();
564 
565  // normal: assume faces are planar
566  const auto& n_F = ig.centerUnitOuterNormal();
567 
568  // evaluate speed of sound (assumed constant per element)
569  const auto& inside_local =
570  Dune::ReferenceElements<DF,dim>::general(ig.inside().type()).position(0,0);
571  auto mu_s = param.mu(ig.inside(),inside_local);
572  auto eps_s = param.eps(ig.inside(),inside_local);
573  //auto sigma_s = param.sigma(ig.inside(),inside_local);
574 
575  // compute A+ (outgoing waves)
576  Dune::FieldMatrix<DF,dim*2,dim*2> R_s;
577  MaxwellEigenvectors<dim>::eigenvectors(eps_s,mu_s,n_F,R_s);
578  Dune::FieldMatrix<DF,dim*2,dim*2> Dplus_s(0.0);
579  Dplus_s[0][0] = 1.0/sqrt(eps_s*mu_s);
580  Dplus_s[1][1] = 1.0/sqrt(eps_s*mu_s);
581  Dune::FieldMatrix<DF,dim*2,dim*2> Aplus_s(R_s);
582  Aplus_s.rightmultiply(Dplus_s);
583  R_s.invert();
584  Aplus_s.rightmultiply(R_s);
585 
586  // compute A- (incoming waves)
587  Dune::FieldMatrix<DF,dim*2,dim*2> R_n;
588  MaxwellEigenvectors<dim>::eigenvectors(eps_s,mu_s,n_F,R_n);
589  Dune::FieldMatrix<DF,dim*2,dim*2> Dminus_n(0.0);
590  Dminus_n[2][2] = -1.0/sqrt(eps_s*mu_s);
591  Dminus_n[3][3] = -1.0/sqrt(eps_s*mu_s);
592  Dune::FieldMatrix<DF,dim*2,dim*2> Aminus_n(R_n);
593  Aminus_n.rightmultiply(Dminus_n);
594  R_n.invert();
595  Aminus_n.rightmultiply(R_n);
596 
597  // select quadrature rule
598  const int order_s = dgspace_s.finiteElement().localBasis().order();
599  const int intorder = overintegration+1+2*order_s;
600  auto gtface = ig.geometry().type();
601 
602  // std::cout << "alpha_boundary center=" << ig.geometry().center() << std::endl;
603 
604  // loop over quadrature points
605  for(const auto &qp :
606  Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder))
607  {
608  // position of quadrature point in local coordinates of elements
609  const auto& iplocal_s = ig.geometryInInside().global(qp.position());
610 
611  // evaluate basis functions
612  const auto& phi_s = cache[order_s].evaluateFunction
613  (iplocal_s,dgspace_s.finiteElement().localBasis());
614 
615  // evaluate u from inside and outside
616  Dune::FieldVector<RF,dim*2> u_s(0.0);
617  for (size_type i=0; i<dim*2; i++) // for all components
618  for (size_type k=0; k<dgspace_s.size(); k++) // for all basis functions
619  u_s[i] += x_s(lfsv_s.child(i),k)*phi_s[k];
620  // std::cout << " u_s " << u_s << std::endl;
621 
622  // evaluate boundary condition
623  Dune::FieldVector<RF,dim*2> u_n(param.g(ig.intersection(),qp.position(),u_s));
624  // std::cout << " u_n " << u_n << " bc: " << param.g(ig.intersection(),qp.position(),u_s) << std::endl;
625 
626  // compute numerical flux at integration point
627  Dune::FieldVector<RF,dim*2> f(0.0);
628  Aplus_s.umv(u_s,f);
629  // std::cout << " after A_plus*u_s " << f << std::endl;
630  Aminus_n.umv(u_n,f);
631  // std::cout << " after A_minus*u_n " << f << std::endl;
632 
633  // integrate
634  RF factor = qp.weight() * ig.geometry().integrationElement(qp.position());
635  for (size_type k=0; k<dgspace_s.size(); k++) // loop over all vector-valued (!) basis functions (with identical components)
636  for (size_type i=0; i<dim*2; i++) // loop over all components
637  r_s.accumulate(lfsv_s.child(i),k,f[i]*phi_s[k]*factor);
638  }
639 
640  // std::cout << " residual_s: ";
641  // for (size_type i=0; i<r_s.size(); i++) std::cout << r_s[i] << " ";
642  // std::cout << std::endl;
643  }
644 
645  // volume integral depending only on test functions
646  template<typename EG, typename LFSV, typename R>
647  void lambda_volume (const EG& eg, const LFSV& lfsv, R& r) const
648  {
649  // get types
650  using DGSpace = typename LFSV::template Child<0>::Type;
651  using DF = typename DGSpace::Traits::FiniteElementType::
652  Traits::LocalBasisType::Traits::DomainFieldType;
653  using size_type = typename DGSpace::Traits::SizeType;
654 
655  // get local function space that is identical for all components
656  const auto& dgspace = lfsv.template child<0>();
657 
658  // select quadrature rule
659  const int order_s = dgspace.finiteElement().localBasis().order();
660  const int intorder = overintegration+2*order_s;
661  Dune::GeometryType gt = eg.geometry().type();
662 
663  // loop over quadrature points
664  for (const auto &qp : Dune::QuadratureRules<DF,dim>::rule(gt,intorder))
665  {
666  // evaluate right hand side
667  auto j = param.j(eg.entity(),qp.position());
668 
669  // evaluate basis functions
670  const auto& phi = cache[order_s].evaluateFunction(qp.position(),dgspace.finiteElement().localBasis());
671 
672  // integrate
673  auto factor = qp.weight() * eg.geometry().integrationElement(qp.position());
674  for (size_type k=0; k<dim*2; k++) // for all components
675  for (size_type i=0; i<dgspace.size(); i++) // for all test functions of this component
676  r.accumulate(lfsv.child(k),i,-j[k]*phi[i]*factor);
677  }
678  }
679 
681  void setTime (typename T::Traits::RangeFieldType t)
682  {
683  }
684 
686  void preStep (typename T::Traits::RangeFieldType time, typename T::Traits::RangeFieldType dt,
687  int stages)
688  {
689  }
690 
692  void preStage (typename T::Traits::RangeFieldType time, int r)
693  {
694  }
695 
697  void postStage ()
698  {
699  }
700 
702  typename T::Traits::RangeFieldType suggestTimestep (typename T::Traits::RangeFieldType dt) const
703  {
704  return dt;
705  }
706 
707  private:
708  T& param;
709  int overintegration;
710  typedef typename FEM::Traits::FiniteElementType::Traits::LocalBasisType LocalBasisType;
712  std::vector<Cache> cache;
713  };
714 
715 
716 
728  template<typename T, typename FEM>
730  public NumericalJacobianApplyVolume<DGMaxwellTemporalOperator<T,FEM> >,
732  public InstationaryLocalOperatorDefaultMethods<typename T::Traits::RangeFieldType>
733  {
734  enum { dim = T::Traits::GridViewType::dimension };
735  public:
736  // pattern assembly flags
737  enum { doPatternVolume = true };
738 
739  // residual assembly flags
740  enum { doAlphaVolume = true };
741 
742  DGMaxwellTemporalOperator (T& param_, int overintegration_=0)
743  : param(param_), overintegration(overintegration_), cache(20)
744  {}
745 
746  // define sparsity pattern of operator representation
747  template<typename LFSU, typename LFSV, typename LocalPattern>
748  void pattern_volume (const LFSU& lfsu, const LFSV& lfsv,
749  LocalPattern& pattern) const
750  {
751  // paranoia check number of number of components
752  static_assert(LFSU::CHILDREN==LFSV::CHILDREN, "need U=V!");
753  static_assert(LFSV::CHILDREN==dim*2, "need exactly dim*2 components!");
754 
755  for (size_t k=0; k<LFSV::CHILDREN; k++)
756  for (size_t i=0; i<lfsv.child(k).size(); ++i)
757  for (size_t j=0; j<lfsu.child(k).size(); ++j)
758  pattern.addLink(lfsv.child(k),i,lfsu.child(k),j);
759  }
760 
761  // volume integral depending on test and ansatz functions
762  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
763  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
764  {
765  // get types
766  using DGSpace = typename LFSV::template Child<0>::Type;
767  using DF = typename DGSpace::Traits::FiniteElementType::
768  Traits::LocalBasisType::Traits::DomainFieldType;
769  using RF = typename DGSpace::Traits::FiniteElementType::
770  Traits::LocalBasisType::Traits::RangeFieldType;
771  using size_type = typename DGSpace::Traits::SizeType;
772 
773  // get local function space that is identical for all components
774  const auto& dgspace = lfsv.template child<0>();
775 
776  // select quadrature rule
777  const int order = dgspace.finiteElement().localBasis().order();
778  const int intorder = overintegration+2*order;
779  auto gt = eg.geometry().type();
780 
781  // loop over quadrature points
782  for (const auto& qp : Dune::QuadratureRules<DF,dim>::rule(gt,intorder))
783  {
784  // evaluate basis functions
785  const auto& phi = cache[order].evaluateFunction(qp.position(),dgspace.finiteElement().localBasis());
786 
787  // evaluate u
788  Dune::FieldVector<RF,dim*2> u(0.0);
789  for (size_type k=0; k<dim*2; k++) // for all components
790  for (size_type j=0; j<dgspace.size(); j++) // for all basis functions
791  u[k] += x(lfsv.child(k),j)*phi[j];
792 
793  // integrate
794  RF factor = qp.weight() * eg.geometry().integrationElement(qp.position());
795  for (size_type k=0; k<dim*2; k++) // for all components
796  for (size_type i=0; i<dgspace.size(); i++) // for all test functions of this component
797  r.accumulate(lfsv.child(k),i,u[k]*phi[i]*factor);
798  }
799  }
800 
801  // jacobian of volume term
802  template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
803  void jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv,
804  M& mat) const
805  {
806  // get types
807  using DGSpace = typename LFSV::template Child<0>::Type;
808  using DF = typename DGSpace::Traits::FiniteElementType::
809  Traits::LocalBasisType::Traits::DomainFieldType;
810  using RF = typename DGSpace::Traits::FiniteElementType::
811  Traits::LocalBasisType::Traits::RangeFieldType;
812  using size_type = typename DGSpace::Traits::SizeType;
813 
814  // get local function space that is identical for all components
815  const auto& dgspace = lfsv.template child<0>();
816 
817  // select quadrature rule
818  const int order = dgspace.finiteElement().localBasis().order();
819  const int intorder = overintegration+2*order;
820  auto gt = eg.geometry().type();
821 
822  // loop over quadrature points
823  for (const auto& qp : Dune::QuadratureRules<DF,dim>::rule(gt,intorder))
824  {
825  // evaluate basis functions
826  const auto& phi = cache[order].evaluateFunction(qp.position(),dgspace.finiteElement().localBasis());
827 
828  // integrate
829  RF factor = qp.weight() * eg.geometry().integrationElement(qp.position());
830  for (size_type k=0; k<dim*2; k++) // for all components
831  for (size_type i=0; i<dgspace.size(); i++) // for all test functions of this component
832  for (size_type j=0; j<dgspace.size(); j++) // for all ansatz functions of this component
833  mat.accumulate(lfsv.child(k),i,lfsu.child(k),j,phi[j]*phi[i]*factor);
834  }
835  }
836 
837  private:
838  T& param;
839  int overintegration;
840  typedef typename FEM::Traits::FiniteElementType::Traits::LocalBasisType LocalBasisType;
842  std::vector<Cache> cache;
843  };
844 
845  }
846 }
847 
848 #endif
static void eigenvalues(T1 eps, T1 mu, const Dune::FieldVector< T2, 2 *dim > &e)
Definition: maxwelldg.hh:51
const E & e
Definition: interpolate.hh:172
Implement jacobian_apply_volume() based on alpha_volume()
Definition: defaultimp.hh:321
Definition: maxwelldg.hh:729
void preStage(typename T::Traits::RangeFieldType time, int r)
to be called once before each stage
Definition: maxwelldg.hh:692
T::Traits::RangeFieldType suggestTimestep(typename T::Traits::RangeFieldType dt) const
to be called once before each stage
Definition: maxwelldg.hh:702
const IG & ig
Definition: constraints.hh:147
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: maxwelldg.hh:763
static void eigenvectors(T1 eps, T1 mu, const Dune::FieldVector< T2, dim > &n, Dune::FieldMatrix< T3, 2 *dim, 2 *dim > &R)
Definition: maxwelldg.hh:74
Definition: maxwelldg.hh:27
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: maxwelldg.hh:335
Default flags for all local operators.
Definition: flags.hh:18
void alpha_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, R &r_s) const
Definition: maxwelldg.hh:550
sparsity pattern generator
Definition: pattern.hh:13
DGMaxwellSpatialOperator(T &param_, int overintegration_=0)
Definition: maxwelldg.hh:328
Definition: maxwelldg.hh:302
Implement jacobian_apply_skeleton() based on alpha_skeleton()
Definition: defaultimp.hh:445
Definition: adaptivity.hh:27
Implement jacobian_volume() based on alpha_volume()
Definition: defaultimp.hh:32
Implement jacobian_boundary() based on alpha_boundary()
Definition: defaultimp.hh:251
static const int dim
Definition: adaptivity.hh:83
void alpha_skeleton(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, const LFSU &lfsu_n, const X &x_n, const LFSV &lfsv_n, R &r_s, R &r_n) const
Definition: maxwelldg.hh:430
DGMaxwellTemporalOperator(T &param_, int overintegration_=0)
Definition: maxwelldg.hh:742
const std::string s
Definition: function.hh:1102
sparsity pattern generator
Definition: pattern.hh:29
void preStep(typename T::Traits::RangeFieldType time, typename T::Traits::RangeFieldType dt, int stages)
to be called once before each time step
Definition: maxwelldg.hh:686
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: maxwelldg.hh:803
void pattern_volume(const LFSU &lfsu, const LFSV &lfsv, LocalPattern &pattern) const
Definition: maxwelldg.hh:748
void postStage()
to be called once at the end of each stage
Definition: maxwelldg.hh:697
store values of basis functions and gradients in a cache
Definition: localbasiscache.hh:15
Implement jacobian_skeleton() based on alpha_skeleton()
Definition: defaultimp.hh:157
void setTime(typename T::Traits::RangeFieldType t)
set time in parameter class
Definition: maxwelldg.hh:681
const EG & eg
Definition: constraints.hh:280
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: maxwelldg.hh:647
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
Implement jacobian_apply_boundary() based on alpha_boundary()
Definition: defaultimp.hh:538