2#ifndef DUNE_PDELAB_LOCALOPERATOR_MAXWELLDG_HH
3#define DUNE_PDELAB_LOCALOPERATOR_MAXWELLDG_HH
9#include<dune/common/indices.hh>
11#include <dune/geometry/referenceelements.hh>
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>
22#include"maxwellparameter.hh"
29 class MaxwellEigenvectors
38 class MaxwellEigenvectors<3>
52 template<
typename T1,
typename T2,
typename T3>
56 T1 s = 1.0/sqrt(mu*eps);
76 template<
typename T1,
typename T2,
typename T3>
80 T1 a=n[0], b=n[1], c=n[2];
85 re[0]=a*c; re[1]=b*c; re[2]=c*c-1;
86 im[0]=-b; im[1]=a; im[2]=0;
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;
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];
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];
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;
121 for (std::size_t i=0; i<3; i++)
122 for (std::size_t j=0; j<6; j++)
124 for (std::size_t i=3; i<6; i++)
125 for (std::size_t j=0; j<6; j++)
315 template<
typename T,
typename FEM>
328 enum { dim = T::Traits::GridViewType::dimension };
332 enum { doPatternVolume =
true };
333 enum { doPatternSkeleton =
true };
336 enum { doAlphaVolume =
true };
337 enum { doAlphaSkeleton =
true };
338 enum { doAlphaBoundary =
true };
339 enum { doLambdaVolume =
true };
343 : param(param_), overintegration(overintegration_), cache(20)
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
354 using RF =
typename DGSpace::Traits::FiniteElementType::Traits
355 ::LocalBasisType::Traits::RangeFieldType;
356 using size_type =
typename DGSpace::Traits::SizeType;
360 "need exactly dim*2 components!");
363 const auto& dgspace =
child(lfsv,
_0);
366 const auto& cell = eg.entity();
369 auto geo = eg.geometry();
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);
378 auto epsinv = 1.0/eps;
383 typename EG::Geometry::JacobianInverseTransposed jac;
387 std::vector<Dune::FieldVector<RF,dim> > gradphi(dgspace.size());
390 const int order = dgspace.finiteElement().localBasis().order();
391 const int intorder = overintegration+2*order;
392 for (
const auto &qp : quadratureRule(geo,intorder))
395 const auto& phi = cache[order].evaluateFunction
396 (qp.position(), dgspace.finiteElement().localBasis());
399 for (size_type k=0; k<dim*2; k++){
401 for (size_type j=0; j<dgspace.size(); j++)
402 u[k] += x(lfsv.child(k),j)*phi[j];
407 const auto& js = cache[order].evaluateJacobian
408 (qp.position(), dgspace.finiteElement().localBasis());
411 jac = geo.jacobianInverseTransposed(qp.position());
412 for (size_type i=0; i<dgspace.size(); i++)
413 jac.mv(js[i][0],gradphi[i]);
416 auto factor = qp.weight() * geo.integrationElement(qp.position());
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;
429 for (size_type i=0; i<dim*2; i++)
431 for (size_type k=0; k<dgspace.size(); k++)
433 for (size_type j=0; j<dim; j++)
434 r.accumulate(lfsv.child(i),k,-F[i][j]*gradphi[k][j]*factor);
437 for (size_type i=0; i<dim; i++)
439 for (size_type k=0; k<dgspace.size(); k++)
440 r.accumulate(lfsv.child(i),k,(sigma/eps)*u[i]*phi[k]*factor);
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
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;
469 const auto& dgspace_s =
child(lfsv_s,
_0);
470 const auto& dgspace_n =
child(lfsv_n,
_0);
473 const auto& cell_inside = ig.inside();
474 const auto& cell_outside = ig.outside();
477 auto geo = ig.geometry();
478 auto geo_inside = cell_inside.geometry();
479 auto geo_outside = cell_outside.geometry();
482 auto geo_in_inside = ig.geometryInInside();
483 auto geo_in_outside = ig.geometryInOutside();
488 auto local_inside = ref_el_inside.position(0,0);
489 auto local_outside = ref_el_outside.position(0,0);
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);
502 const auto& n_F = ig.centerUnitOuterNormal();
506 MaxwellEigenvectors<dim>::eigenvectors(eps_s,mu_s,n_F,R_s);
508 Dplus_s[0][0] = 1.0/sqrt(eps_s*mu_s);
509 Dplus_s[1][1] = 1.0/sqrt(eps_s*mu_s);
517 MaxwellEigenvectors<dim>::eigenvectors(eps_n,mu_n,n_F,R_n);
519 Dminus_n[2][2] = -1.0/sqrt(eps_n*mu_n);
520 Dminus_n[3][3] = -1.0/sqrt(eps_n*mu_n);
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))
540 const auto& iplocal_s = geo_in_inside.global(qp.position());
541 const auto& iplocal_n = geo_in_outside.global(qp.position());
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());
548 for (size_type i=0; i<dim*2; i++){
550 for (size_type k=0; k<dgspace_s.size(); k++)
551 u_s[i] += x_s(lfsv_s.child(i),k)*phi_s[k];
553 for (size_type i=0; i<dim*2; i++){
555 for (size_type k=0; k<dgspace_n.size(); k++)
556 u_n[i] += x_n(lfsv_n.child(i),k)*phi_n[k];
569 auto factor = qp.weight() * geo.integrationElement(qp.position());
570 for (size_type k=0; k<dgspace_s.size(); k++)
571 for (size_type i=0; i<dim*2; i++)
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++)
574 for (size_type i=0; i<dim*2; i++)
575 r_n.accumulate(lfsv_n.child(i),k,-f[i]*phi_n[k]*factor);
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,
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;
603 const auto& dgspace_s =
child(lfsv_s,
_0);
606 const auto& cell_inside = ig.inside();
609 auto geo = ig.geometry();
610 auto geo_inside = cell_inside.geometry();
613 auto geo_in_inside = ig.geometryInInside();
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);
623 const auto& n_F = ig.centerUnitOuterNormal();
627 MaxwellEigenvectors<dim>::eigenvectors(eps_s,mu_s,n_F,R_s);
629 Dplus_s[0][0] = 1.0/sqrt(eps_s*mu_s);
630 Dplus_s[1][1] = 1.0/sqrt(eps_s*mu_s);
638 MaxwellEigenvectors<dim>::eigenvectors(eps_s,mu_s,n_F,R_n);
640 Dminus_n[2][2] = -1.0/sqrt(eps_s*mu_s);
641 Dminus_n[3][3] = -1.0/sqrt(eps_s*mu_s);
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))
660 const auto& iplocal_s = geo_in_inside.global(qp.position());
663 const auto& phi_s = cache[order_s].evaluateFunction
664 (iplocal_s,dgspace_s.finiteElement().localBasis());
667 for (size_type i=0; i<dim*2; i++){
669 for (size_type k=0; k<dgspace_s.size(); k++)
670 u_s[i] += x_s(lfsv_s.child(i),k)*phi_s[k];
675 u_n = (param.g(ig.intersection(),qp.position(),u_s));
686 auto factor = qp.weight() * geo.integrationElement(qp.position());
687 for (size_type k=0; k<dgspace_s.size(); k++)
688 for (size_type i=0; i<dim*2; i++)
689 r_s.accumulate(lfsv_s.child(i),k,f[i]*phi_s[k]*factor);
698 template<
typename EG,
typename LFSV,
typename R>
699 void lambda_volume (
const EG& eg,
const LFSV& lfsv, R& r)
const
704 using size_type =
typename DGSpace::Traits::SizeType;
707 const auto& dgspace =
child(lfsv,
_0);
710 const auto& cell = eg.entity();
713 auto geo = eg.geometry();
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))
721 auto j = param.j(cell,qp.position());
724 const auto& phi = cache[order_s].evaluateFunction(qp.position(),dgspace.finiteElement().localBasis());
727 auto factor = qp.weight() * geo.integrationElement(qp.position());
728 for (size_type k=0; k<dim*2; k++)
729 for (size_type i=0; i<dgspace.size(); i++)
730 r.accumulate(lfsv.child(k),i,-j[k]*phi[i]*factor);
735 void setTime (
typename T::Traits::RangeFieldType t)
740 void preStep (
typename T::Traits::RangeFieldType time,
typename T::Traits::RangeFieldType dt,
746 void preStage (
typename T::Traits::RangeFieldType time,
int r)
756 typename T::Traits::RangeFieldType
suggestTimestep (
typename T::Traits::RangeFieldType dt)
const
764 typedef typename FEM::Traits::FiniteElementType::Traits::LocalBasisType LocalBasisType;
766 std::vector<Cache> cache;
782 template<
typename T,
typename FEM>
788 enum { dim = T::Traits::GridViewType::dimension };
791 enum { doPatternVolume =
true };
794 enum { doAlphaVolume =
true };
797 : param(param_), overintegration(overintegration_), cache(20)
801 template<
typename LFSU,
typename LFSV,
typename LocalPattern>
802 void pattern_volume (
const LFSU& lfsu,
const LFSV& lfsv,
803 LocalPattern& pattern)
const
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);
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
822 using RF =
typename DGSpace::Traits::FiniteElementType::
823 Traits::LocalBasisType::Traits::RangeFieldType;
824 using size_type =
typename DGSpace::Traits::SizeType;
827 const auto& dgspace =
child(lfsv,
_0);
830 auto geo = eg.geometry();
836 const int order = dgspace.finiteElement().localBasis().order();
837 const int intorder = overintegration+2*order;
838 for (
const auto& qp : quadratureRule(geo,intorder))
841 const auto& phi = cache[order].evaluateFunction(qp.position(),dgspace.finiteElement().localBasis());
844 for (size_type k=0; k<dim*2; k++){
846 for (size_type j=0; j<dgspace.size(); j++)
847 u[k] += x(lfsv.child(k),j)*phi[j];
851 auto factor = qp.weight() * geo.integrationElement(qp.position());
852 for (size_type k=0; k<dim*2; k++)
853 for (size_type i=0; i<dgspace.size(); i++)
854 r.accumulate(lfsv.child(k),i,u[k]*phi[i]*factor);
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,
866 using size_type =
typename DGSpace::Traits::SizeType;
869 const auto& dgspace =
child(lfsv,
_0);
872 auto geo = eg.geometry();
875 const int order = dgspace.finiteElement().localBasis().order();
876 const int intorder = overintegration+2*order;
877 for (
const auto& qp : quadratureRule(geo,intorder))
880 const auto& phi = cache[order].evaluateFunction(qp.position(),dgspace.finiteElement().localBasis());
883 auto factor = qp.weight() * geo.integrationElement(qp.position());
884 for (size_type k=0; k<dim*2; k++)
885 for (size_type i=0; i<dgspace.size(); i++)
886 for (size_type j=0; j<dgspace.size(); j++)
887 mat.accumulate(lfsv.child(k),i,lfsu.child(k),j,phi[j]*phi[i]*factor);
894 typedef typename FEM::Traits::FiniteElementType::Traits::LocalBasisType LocalBasisType;
896 std::vector<Cache> cache;
void invert(bool doPivoting=true)
Compute inverse.
void umv(const X &x, Y &y) const
y += A x
Definition: densematrix.hh:406
A dense n x m matrix.
Definition: fmatrix.hh:117
FieldMatrix & rightmultiply(const FieldMatrix< K, r, c > &M)
Multiplies M from the right to this matrix.
Definition: fmatrix.hh:333
vector space out of a tensor product of fields.
Definition: fvector.hh:91
Spatial local operator for discontinuous Galerkin method for Maxwells Equations.
Definition: maxwelldg.hh:327
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
T::Traits::RangeFieldType suggestTimestep(typename T::Traits::RangeFieldType dt) const
to be called once before each stage
Definition: maxwelldg.hh:756
void setTime(typename T::Traits::RangeFieldType t)
set time in parameter class
Definition: maxwelldg.hh:735
void preStage(typename T::Traits::RangeFieldType time, int r)
to be called once before each stage
Definition: maxwelldg.hh:746
void postStage()
to be called once at the end of each stage
Definition: maxwelldg.hh:751
Definition: maxwelldg.hh:787
sparsity pattern generator
Definition: pattern.hh:30
sparsity pattern generator
Definition: pattern.hh:14
Default class for additional methods in instationary local operators.
Definition: idefault.hh:90
store values of basis functions and gradients in a cache
Definition: localbasiscache.hh:18
Default flags for all local operators.
Definition: flags.hh:19
static void eigenvalues(T1 eps, T1 mu, const Dune::FieldVector< T2, 2 *dim > &e)
Definition: maxwelldg.hh:53
static void eigenvectors(T1 eps, T1 mu, const Dune::FieldVector< T2, dim > &n, Dune::FieldMatrix< T3, 2 *dim, 2 *dim > &R)
Definition: maxwelldg.hh:77
Implements linear and nonlinear versions of jacobian_apply_boundary() based on alpha_boundary()
Definition: numericaljacobianapply.hh:287
Implements linear and nonlinear versions of jacobian_apply_skeleton() based on alpha_skeleton()
Definition: numericaljacobianapply.hh:182
Implements linear and nonlinear versions of jacobian_apply_volume() based on alpha_volume()
Definition: numericaljacobianapply.hh:34
Implement jacobian_boundary() based on alpha_boundary()
Definition: numericaljacobian.hh:251
Implement jacobian_skeleton() based on alpha_skeleton()
Definition: numericaljacobian.hh:157
Implement jacobian_volume() based on alpha_volume()
Definition: numericaljacobian.hh:32
A few common exception classes.
Implements a vector constructed from a given type representing a field and a compile-time given size.
constexpr index_constant< 0 > _0
Compile time index with value 0.
Definition: indices.hh:52
unspecified value type referenceElement(T &&... t)
Returns a reference element for the objects t....
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
std::size_t degree(const Node &node)
Returns the degree of node as run time information.
Definition: nodeinterface.hh:79
decltype(Node::degree()) StaticDegree
Returns the statically known degree of the given Node type as a std::integral_constant.
Definition: nodeinterface.hh:107
Namespace with predefined compile time indices for the range [0,19].
Definition: indices.hh:50
Dune namespace.
Definition: alignedallocator.hh:13