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