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