DUNE PDELab (2.7)

maxwelldg.hh
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
13#include<dune/pdelab/common/quadraturerules.hh>
14#include<dune/pdelab/common/function.hh>
15#include<dune/pdelab/common/geometrywrapper.hh>
16#include<dune/pdelab/finiteelement/localbasiscache.hh>
17#include<dune/pdelab/localoperator/defaultimp.hh>
18#include<dune/pdelab/localoperator/flags.hh>
19#include<dune/pdelab/localoperator/idefault.hh>
20#include<dune/pdelab/localoperator/pattern.hh>
21
22#include"maxwellparameter.hh"
23
24namespace Dune {
25 namespace PDELab {
26
27
28 template<int dim>
29 class MaxwellEigenvectors
30 {};
31
37 template<>
38 class MaxwellEigenvectors<3>
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>
314 class DGMaxwellSpatialOperator :
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,
323 public LocalOperatorDefaultFlags,
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;
762 typedef Dune::PDELab::LocalBasisCache<LocalBasisType> Cache;
763 std::vector<Cache> cache;
764 };
765
766
767
779 template<typename T, typename FEM>
780 class DGMaxwellTemporalOperator :
781 public NumericalJacobianApplyVolume<DGMaxwellTemporalOperator<T,FEM> >,
782 public LocalOperatorDefaultFlags,
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
803 static_assert(TypeTree::StaticDegree<LFSU>::value==TypeTree::StaticDegree<LFSV>::value, "need U=V!");
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;
892 typedef Dune::PDELab::LocalBasisCache<LocalBasisType> Cache;
893 std::vector<Cache> cache;
894 };
895
896 }
897}
898
899#endif // DUNE_PDELAB_LOCALOPERATOR_MAXWELLDG_HH
Dune::DenseMatrix::invert
void invert(bool doPivoting=true)
Compute inverse.
Dune::DenseMatrix::umv
void umv(const X &x, Y &y) const
y += A x
Definition: densematrix.hh:433
Dune::FieldMatrix
A dense n x m matrix.
Definition: fmatrix.hh:69
Dune::FieldMatrix::rightmultiply
FieldMatrix & rightmultiply(const FieldMatrix< K, r, c > &M)
Multiplies M from the right to this matrix.
Definition: fmatrix.hh:235
Dune::FieldVector
vector space out of a tensor product of fields.
Definition: fvector.hh:96
Dune::PDELab::DGMaxwellSpatialOperator
Spatial local operator for discontinuous Galerkin method for Maxwells Equations.
Definition: maxwelldg.hh:325
Dune::PDELab::DGMaxwellSpatialOperator::preStep
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
Dune::PDELab::DGMaxwellSpatialOperator::suggestTimestep
T::Traits::RangeFieldType suggestTimestep(typename T::Traits::RangeFieldType dt) const
to be called once before each stage
Definition: maxwelldg.hh:753
Dune::PDELab::DGMaxwellSpatialOperator::setTime
void setTime(typename T::Traits::RangeFieldType t)
set time in parameter class
Definition: maxwelldg.hh:732
Dune::PDELab::DGMaxwellSpatialOperator::preStage
void preStage(typename T::Traits::RangeFieldType time, int r)
to be called once before each stage
Definition: maxwelldg.hh:743
Dune::PDELab::DGMaxwellSpatialOperator::postStage
void postStage()
to be called once at the end of each stage
Definition: maxwelldg.hh:748
Dune::PDELab::DGMaxwellTemporalOperator
Definition: maxwelldg.hh:784
Dune::PDELab::FullSkeletonPattern
sparsity pattern generator
Definition: pattern.hh:30
Dune::PDELab::FullVolumePattern
sparsity pattern generator
Definition: pattern.hh:14
Dune::PDELab::InstationaryLocalOperatorDefaultMethods
Default class for additional methods in instationary local operators.
Definition: idefault.hh:90
Dune::PDELab::LocalBasisCache
store values of basis functions and gradients in a cache
Definition: localbasiscache.hh:18
Dune::PDELab::LocalOperatorDefaultFlags
Default flags for all local operators.
Definition: flags.hh:19
Dune::PDELab::MaxwellEigenvectors< 3 >::eigenvalues
static void eigenvalues(T1 eps, T1 mu, const Dune::FieldVector< T2, 2 *dim > &e)
Definition: maxwelldg.hh:53
Dune::PDELab::MaxwellEigenvectors< 3 >::eigenvectors
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
Dune::PDELab::NumericalJacobianApplyBoundary
Implements linear and nonlinear versions of jacobian_apply_boundary() based on alpha_boundary()
Definition: numericaljacobianapply.hh:287
Dune::PDELab::NumericalJacobianApplySkeleton
Implements linear and nonlinear versions of jacobian_apply_skeleton() based on alpha_skeleton()
Definition: numericaljacobianapply.hh:182
Dune::PDELab::NumericalJacobianApplyVolume
Implements linear and nonlinear versions of jacobian_apply_volume() based on alpha_volume()
Definition: numericaljacobianapply.hh:34
Dune::PDELab::NumericalJacobianBoundary
Implement jacobian_boundary() based on alpha_boundary()
Definition: numericaljacobian.hh:251
Dune::PDELab::NumericalJacobianSkeleton
Implement jacobian_skeleton() based on alpha_skeleton()
Definition: numericaljacobian.hh:157
Dune::PDELab::NumericalJacobianVolume
Implement jacobian_volume() based on alpha_volume()
Definition: numericaljacobian.hh:32
exceptions.hh
A few common exception classes.
fvector.hh
Implements a vector constructed from a given type representing a field and a compile-time given size.
Dune::Indices::_0
constexpr index_constant< 0 > _0
Compile time index with value 0.
Definition: indices.hh:51
Dune::referenceElement
unspecified value type referenceElement(T &&... t)
Returns a reference element for the objects t....
Dune::TypeTree::degree
std::size_t degree(const Node &node)
Returns the degree of node as run time information.
Definition: nodeinterface.hh:71
Dune::TypeTree::StaticDegree
std::integral_constant< std::size_t, degree(static_cast< std::decay_t< Node > * >(nullptr), NodeTag< std::decay_t< Node > >()) > StaticDegree
Returns the statically known degree of the given Node type as a std::integral_constant.
Definition: nodeinterface.hh:105
Dune::Simd::Overloads::max
auto max(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::max()
Definition: defaults.hh:79
Dune::TypeTree::Child
typename impl::_Child< Node, indices... >::type Child
Template alias for the type of a child node given by a list of child indices.
Definition: childextraction.hh:276
Dune::TypeTree::child
ImplementationDefined child(Node &&node, Indices... indices)
Extracts the child of a node given by a sequence of compile-time and run-time indices.
Definition: childextraction.hh:179
Dune::Indices
Namespace with predefined compile time indices for the range [0,19].
Definition: indices.hh:49
Dune
Dune namespace.
Definition: alignedallocator.hh:14
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden & Uni Heidelberg  |  generated with Hugo v0.111.3 (Mar 3, 23:33, 2026)