DUNE PDELab (git)

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