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