2 #ifndef DUNE_PDELAB_LOCALOPERATOR_MAXWELLDG_HH
3 #define DUNE_PDELAB_LOCALOPERATOR_MAXWELLDG_HH
7 #include<dune/common/exceptions.hh>
8 #include<dune/common/fvector.hh>
9 #include<dune/common/indices.hh>
11 #include <dune/geometry/referenceelements.hh>
52 template<
typename T1,
typename T2,
typename T3>
53 static void eigenvalues (T1 eps, T1 mu,
const Dune::FieldVector<T2,2*dim>&
e)
55 T1
s = 1.0/sqrt(mu*eps);
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)
78 T1 a=n[0], b=n[1], c=n[2];
80 Dune::FieldVector<T2,dim> re, im;
83 re[0]=a*c; re[1]=b*c; re[2]=c*c-1;
84 im[0]=-b; im[1]=a; im[2]=0;
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;
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];
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];
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;
119 for (std::size_t i=0; i<3; i++)
120 for (std::size_t j=0; j<6; j++)
122 for (std::size_t i=3; i<6; i++)
123 for (std::size_t j=0; j<6; j++)
313 template<
typename T,
typename FEM>
326 enum { dim = T::Traits::GridViewType::dimension };
341 : param(param_), overintegration(overintegration_), cache(20)
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
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;
358 "need exactly dim*2 components!");
361 const auto& dgspace = child(lfsv,_0);
364 const auto& cell = eg.entity();
367 auto geo = eg.geometry();
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);
376 auto epsinv = 1.0/eps;
381 typename EG::Geometry::JacobianInverseTransposed jac;
384 Dune::FieldVector<RF,dim*2> u(0.0);
385 std::vector<Dune::FieldVector<RF,dim> > gradphi(dgspace.size());
388 const int order = dgspace.finiteElement().localBasis().order();
389 const int intorder = overintegration+2*order;
393 const auto& phi = cache[order].evaluateFunction
394 (qp.position(), dgspace.finiteElement().localBasis());
397 for (size_type k=0; k<dim*2; k++){
399 for (size_type j=0; j<dgspace.size(); j++)
400 u[k] += x(lfsv.child(k),j)*phi[j];
405 const auto& js = cache[order].evaluateJacobian
406 (qp.position(), dgspace.finiteElement().localBasis());
409 jac = geo.jacobianInverseTransposed(qp.position());
410 for (size_type i=0; i<dgspace.size(); i++)
411 jac.mv(js[i][0],gradphi[i]);
414 auto factor = qp.weight() * geo.integrationElement(qp.position());
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;
427 for (size_type i=0; i<dim*2; i++)
429 for (size_type k=0; k<dgspace.size(); k++)
431 for (size_type j=0; j<dim; j++)
432 r.accumulate(lfsv.child(i),k,-F[i][j]*gradphi[k][j]*factor);
435 for (size_type i=0; i<dim; i++)
437 for (size_type k=0; k<dgspace.size(); k++)
438 r.accumulate(lfsv.child(i),k,(sigma/eps)*u[i]*phi[k]*factor);
448 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
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
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;
467 const auto& dgspace_s = child(lfsv_s,_0);
468 const auto& dgspace_n = child(lfsv_n,_0);
471 const auto& cell_inside =
ig.inside();
472 const auto& cell_outside =
ig.outside();
475 auto geo =
ig.geometry();
476 auto geo_inside = cell_inside.geometry();
477 auto geo_outside = cell_outside.geometry();
480 auto geo_in_inside =
ig.geometryInInside();
481 auto geo_in_outside =
ig.geometryInOutside();
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);
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);
500 const auto& n_F =
ig.centerUnitOuterNormal();
503 Dune::FieldMatrix<DF,dim*2,dim*2> 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);
511 Aplus_s.rightmultiply(R_s);
514 Dune::FieldMatrix<DF,dim*2,dim*2> 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);
522 Aminus_n.rightmultiply(R_n);
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);
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);
538 const auto& iplocal_s = geo_in_inside.global(qp.position());
539 const auto& iplocal_n = geo_in_outside.global(qp.position());
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());
546 for (size_type i=0; i<dim*2; i++){
548 for (size_type k=0; k<dgspace_s.size(); k++)
549 u_s[i] += x_s(lfsv_s.child(i),k)*phi_s[k];
551 for (size_type i=0; i<dim*2; i++){
553 for (size_type k=0; k<dgspace_n.size(); k++)
554 u_n[i] += x_n(lfsv_n.child(i),k)*phi_n[k];
567 auto factor = qp.weight() * geo.integrationElement(qp.position());
568 for (size_type k=0; k<dgspace_s.size(); k++)
569 for (size_type i=0; i<dim*2; i++)
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++)
572 for (size_type i=0; i<dim*2; i++)
573 r_n.accumulate(lfsv_n.child(i),k,-f[i]*phi_n[k]*factor);
585 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
587 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
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;
600 const auto& dgspace_s = child(lfsv_s,_0);
603 const auto& cell_inside =
ig.inside();
606 auto geo =
ig.geometry();
607 auto geo_inside = cell_inside.geometry();
610 auto geo_in_inside =
ig.geometryInInside();
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);
620 const auto& n_F =
ig.centerUnitOuterNormal();
623 Dune::FieldMatrix<DF,dim*2,dim*2> 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);
631 Aplus_s.rightmultiply(R_s);
634 Dune::FieldMatrix<DF,dim*2,dim*2> 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);
642 Aminus_n.rightmultiply(R_n);
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);
652 const int order_s = dgspace_s.finiteElement().localBasis().order();
653 const int intorder = overintegration+1+2*order_s;
657 const auto& iplocal_s = geo_in_inside.global(qp.position());
660 const auto& phi_s = cache[order_s].evaluateFunction
661 (iplocal_s,dgspace_s.finiteElement().localBasis());
664 for (size_type i=0; i<dim*2; i++){
666 for (size_type k=0; k<dgspace_s.size(); k++)
667 u_s[i] += x_s(lfsv_s.child(i),k)*phi_s[k];
672 u_n = (param.g(
ig.intersection(),qp.position(),u_s));
683 auto factor = qp.weight() * geo.integrationElement(qp.position());
684 for (size_type k=0; k<dgspace_s.size(); k++)
685 for (size_type i=0; i<dim*2; i++)
686 r_s.accumulate(lfsv_s.child(i),k,f[i]*phi_s[k]*factor);
695 template<
typename EG,
typename LFSV,
typename R>
699 using namespace Indices;
700 using DGSpace = TypeTree::Child<LFSV,0>;
701 using size_type =
typename DGSpace::Traits::SizeType;
704 const auto& dgspace = child(lfsv,_0);
707 const auto& cell = eg.entity();
710 auto geo = eg.geometry();
713 const int order_s = dgspace.finiteElement().localBasis().order();
714 const int intorder = overintegration+2*order_s;
718 auto j = param.j(cell,qp.position());
721 const auto& phi = cache[order_s].evaluateFunction(qp.position(),dgspace.finiteElement().localBasis());
724 auto factor = qp.weight() * geo.integrationElement(qp.position());
725 for (size_type k=0; k<dim*2; k++)
726 for (size_type i=0; i<dgspace.size(); i++)
727 r.accumulate(lfsv.child(k),i,-j[k]*phi[i]*factor);
732 void setTime (
typename T::Traits::RangeFieldType t)
737 void preStep (
typename T::Traits::RangeFieldType time,
typename T::Traits::RangeFieldType dt,
743 void preStage (
typename T::Traits::RangeFieldType time,
int r)
753 typename T::Traits::RangeFieldType
suggestTimestep (
typename T::Traits::RangeFieldType dt)
const
761 typedef typename FEM::Traits::FiniteElementType::Traits::LocalBasisType LocalBasisType;
763 std::vector<Cache> cache;
779 template<
typename T,
typename FEM>
785 enum { dim = T::Traits::GridViewType::dimension };
794 : param(param_), overintegration(overintegration_), cache(20)
798 template<
typename LFSU,
typename LFSV,
typename LocalPattern>
800 LocalPattern& pattern)
const
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);
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
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;
824 const auto& dgspace = child(lfsv,_0);
827 auto geo = eg.geometry();
830 Dune::FieldVector<RF,dim*2> u(0.0);
833 const int order = dgspace.finiteElement().localBasis().order();
834 const int intorder = overintegration+2*order;
838 const auto& phi = cache[order].evaluateFunction(qp.position(),dgspace.finiteElement().localBasis());
841 for (size_type k=0; k<dim*2; k++){
843 for (size_type j=0; j<dgspace.size(); j++)
844 u[k] += x(lfsv.child(k),j)*phi[j];
848 auto factor = qp.weight() * geo.integrationElement(qp.position());
849 for (size_type k=0; k<dim*2; k++)
850 for (size_type i=0; i<dgspace.size(); i++)
851 r.accumulate(lfsv.child(k),i,u[k]*phi[i]*factor);
856 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
861 using namespace Indices;
862 using DGSpace = TypeTree::Child<LFSV,0>;
863 using size_type =
typename DGSpace::Traits::SizeType;
866 const auto& dgspace = child(lfsv,_0);
869 auto geo = eg.geometry();
872 const int order = dgspace.finiteElement().localBasis().order();
873 const int intorder = overintegration+2*order;
877 const auto& phi = cache[order].evaluateFunction(qp.position(),dgspace.finiteElement().localBasis());
880 auto factor = qp.weight() * geo.integrationElement(qp.position());
881 for (size_type k=0; k<dim*2; k++)
882 for (size_type i=0; i<dgspace.size(); i++)
883 for (size_type j=0; j<dgspace.size(); j++)
884 mat.accumulate(lfsv.child(k),i,lfsu.child(k),j,phi[j]*phi[i]*factor);
891 typedef typename FEM::Traits::FiniteElementType::Traits::LocalBasisType LocalBasisType;
893 std::vector<Cache> cache;
static const int dim
Definition: adaptivity.hh:84
const std::string s
Definition: function.hh:843
const IG & ig
Definition: constraints.hh:149
const Entity & e
Definition: localfunctionspace.hh:121
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
QuadratureRuleWrapper< QuadratureRule< typename Geometry::ctype, Geometry::mydimension > > quadratureRule(const Geometry &geo, std::size_t order, QuadratureType::Enum quadrature_type=QuadratureType::GaussLegendre)
Returns a quadrature rule for the given geometry.
Definition: quadraturerules.hh:117
store values of basis functions and gradients in a cache
Definition: localbasiscache.hh:18
Default flags for all local operators.
Definition: flags.hh:19
Default class for additional methods in instationary local operators.
Definition: idefault.hh:90
Definition: maxwelldg.hh:30
static void eigenvalues(T1 eps, T1 mu, const Dune::FieldVector< T2, 2 *dim > &e)
Definition: maxwelldg.hh:53
static void eigenvectors(T1 eps, T1 mu, const Dune::FieldVector< T2, dim > &n, Dune::FieldMatrix< T3, 2 *dim, 2 *dim > &R)
Definition: maxwelldg.hh:76
Spatial local operator for discontinuous Galerkin method for Maxwells Equations.
Definition: maxwelldg.hh:325
@ doPatternVolume
Definition: maxwelldg.hh:330
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
DGMaxwellSpatialOperator(T ¶m_, int overintegration_=0)
Definition: maxwelldg.hh:340
@ doAlphaBoundary
Definition: maxwelldg.hh:336
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: maxwelldg.hh:696
void alpha_skeleton(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, const LFSU &lfsu_n, const X &x_n, const LFSV &lfsv_n, R &r_s, R &r_n) const
Definition: maxwelldg.hh:449
@ doAlphaVolume
Definition: maxwelldg.hh:334
void alpha_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, R &r_s) const
Definition: maxwelldg.hh:586
T::Traits::RangeFieldType suggestTimestep(typename T::Traits::RangeFieldType dt) const
to be called once before each stage
Definition: maxwelldg.hh:753
@ doPatternSkeleton
Definition: maxwelldg.hh:331
void setTime(typename T::Traits::RangeFieldType t)
set time in parameter class
Definition: maxwelldg.hh:732
@ doAlphaSkeleton
Definition: maxwelldg.hh:335
@ doLambdaVolume
Definition: maxwelldg.hh:337
void preStage(typename T::Traits::RangeFieldType time, int r)
to be called once before each stage
Definition: maxwelldg.hh:743
void postStage()
to be called once at the end of each stage
Definition: maxwelldg.hh:748
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: maxwelldg.hh:347
Definition: maxwelldg.hh:784
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: maxwelldg.hh:857
void pattern_volume(const LFSU &lfsu, const LFSV &lfsv, LocalPattern &pattern) const
Definition: maxwelldg.hh:799
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: maxwelldg.hh:814
DGMaxwellTemporalOperator(T ¶m_, int overintegration_=0)
Definition: maxwelldg.hh:793
Implement jacobian_volume() based on alpha_volume()
Definition: numericaljacobian.hh:32
Implement jacobian_skeleton() based on alpha_skeleton()
Definition: numericaljacobian.hh:157
Implement jacobian_boundary() based on alpha_boundary()
Definition: numericaljacobian.hh:251
Implements linear and nonlinear versions of jacobian_apply_volume() based on alpha_volume()
Definition: numericaljacobianapply.hh:34
Implements linear and nonlinear versions of jacobian_apply_skeleton() based on alpha_skeleton()
Definition: numericaljacobianapply.hh:182
Implements linear and nonlinear versions of jacobian_apply_boundary() based on alpha_boundary()
Definition: numericaljacobianapply.hh:287
sparsity pattern generator
Definition: pattern.hh:14
sparsity pattern generator
Definition: pattern.hh:30
static const unsigned int value
Definition: gridfunctionspace/tags.hh:139