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> 53 template<
typename T1,
typename T2,
typename T3>
54 static void eigenvalues (T1 eps, T1 mu,
const Dune::FieldVector<T2,2*dim>&
e)
56 T1
s = 1.0/sqrt(mu*eps);
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)
79 T1 a=n[0], b=n[1], c=n[2];
81 Dune::FieldVector<T2,dim> re, im;
84 re[0]=a*c; re[1]=b*c; re[2]=c*c-1;
85 im[0]=-b; im[1]=a; im[2]=0;
89 re[0]=a*b; re[1]=b*b-1; re[2]=b*c;
90 im[0]=c; im[1]=0.0; im[2]=-a;
94 R[0][0] = re[0]; R[0][1] = -im[0];
95 R[1][0] = re[1]; R[1][1] = -im[1];
96 R[2][0] = re[2]; R[2][1] = -im[2];
97 R[3][0] = im[0]; R[3][1] = re[0];
98 R[4][0] = im[1]; R[4][1] = re[1];
99 R[5][0] = im[2]; R[5][1] = re[2];
102 R[0][2] = im[0]; R[0][3] = re[0];
103 R[1][2] = im[1]; R[1][3] = re[1];
104 R[2][2] = im[2]; R[2][3] = re[2];
105 R[3][2] = re[0]; R[3][3] = -im[0];
106 R[4][2] = re[1]; R[4][3] = -im[1];
107 R[5][2] = re[2]; R[5][3] = -im[2];
110 R[0][4] = a; R[0][5] = 0;
111 R[1][4] = b; R[1][5] = 0;
112 R[2][4] = c; R[2][5] = 0;
113 R[3][4] = 0; R[3][5] = a;
114 R[4][4] = 0; R[4][5] = b;
115 R[5][4] = 0; R[5][5] = c;
120 for (std::size_t i=0; i<3; i++)
121 for (std::size_t j=0; j<6; j++)
123 for (std::size_t i=3; i<6; i++)
124 for (std::size_t j=0; j<6; j++)
314 template<
typename T,
typename FEM>
327 enum {
dim = T::Traits::GridViewType::dimension };
331 enum { doPatternVolume =
true };
332 enum { doPatternSkeleton =
true };
335 enum { doAlphaVolume =
true };
336 enum { doAlphaSkeleton =
true };
337 enum { doAlphaBoundary =
true };
338 enum { doLambdaVolume =
true };
342 : param(param_), overintegration(overintegration_), cache(20)
347 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
348 void alpha_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const 351 using namespace Indices;
352 using DGSpace = TypeTree::Child<LFSV,0>;
353 using RF =
typename DGSpace::Traits::FiniteElementType::Traits
354 ::LocalBasisType::Traits::RangeFieldType;
355 using size_type =
typename DGSpace::Traits::SizeType;
359 "need exactly dim*2 components!");
362 const auto& dgspace = child(lfsv,_0);
365 const auto& cell = eg.entity();
368 auto geo = eg.geometry();
372 auto localcenter = ref_el.position(0,0);
373 auto mu = param.mu(cell,localcenter);
374 auto eps = param.eps(cell,localcenter);
375 auto sigma = param.sigma(cell,localcenter);
377 auto epsinv = 1.0/eps;
382 typename EG::Geometry::JacobianInverseTransposed jac;
385 Dune::FieldVector<RF,dim*2> u(0.0);
386 std::vector<Dune::FieldVector<RF,dim> > gradphi(dgspace.size());
389 const int order = dgspace.finiteElement().localBasis().order();
390 const int intorder = overintegration+2*order;
394 const auto& phi = cache[order].evaluateFunction
395 (qp.position(), dgspace.finiteElement().localBasis());
398 for (size_type k=0; k<
dim*2; k++){
400 for (size_type j=0; j<dgspace.size(); j++)
401 u[k] += x(lfsv.child(k),j)*phi[j];
406 const auto& js = cache[order].evaluateJacobian
407 (qp.position(), dgspace.finiteElement().localBasis());
410 jac = geo.jacobianInverseTransposed(qp.position());
411 for (size_type i=0; i<dgspace.size(); i++)
412 jac.mv(js[i][0],gradphi[i]);
415 auto factor = qp.weight() * geo.integrationElement(qp.position());
417 Dune::FieldMatrix<RF,dim*2,dim> F;
418 static_assert(dim == 3,
"Sorry, the following flux implementation " 419 "can only work for dim == 3");
420 F[0][0] = 0; F[0][1] = -muinv*u[5]; F[0][2] = muinv*u[4];
421 F[1][0] = muinv*u[5]; F[1][1] = 0; F[1][2] = -muinv*u[3];
422 F[2][0] =-muinv*u[4]; F[2][1] = muinv*u[3]; F[2][2] = 0;
423 F[3][0] = 0; F[3][1] = epsinv*u[2]; F[3][2] = -epsinv*u[1];
424 F[4][0] = -epsinv*u[2]; F[4][1] = 0; F[4][2] = epsinv*u[0];
425 F[5][0] = epsinv*u[1]; F[5][1] = -epsinv*u[0]; F[5][2] = 0;
428 for (size_type i=0; i<dim*2; i++)
430 for (size_type k=0; k<dgspace.size(); k++)
432 for (size_type j=0; j<
dim; j++)
433 r.accumulate(lfsv.child(i),k,-F[i][j]*gradphi[k][j]*factor);
436 for (size_type i=0; i<
dim; i++)
438 for (size_type k=0; k<dgspace.size(); k++)
439 r.accumulate(lfsv.child(i),k,(sigma/eps)*u[i]*phi[k]*factor);
449 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
451 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
452 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
453 R& r_s, R& r_n)
const 459 using namespace Indices;
460 using DGSpace = TypeTree::Child<LFSV,0>;
461 using DF =
typename DGSpace::Traits::FiniteElementType::
462 Traits::LocalBasisType::Traits::DomainFieldType;
463 using RF =
typename DGSpace::Traits::FiniteElementType::
464 Traits::LocalBasisType::Traits::RangeFieldType;
465 using size_type =
typename DGSpace::Traits::SizeType;
468 const auto& dgspace_s = child(lfsv_s,_0);
469 const auto& dgspace_n = child(lfsv_n,_0);
472 const auto& cell_inside = ig.inside();
473 const auto& cell_outside = ig.outside();
476 auto geo = ig.geometry();
477 auto geo_inside = cell_inside.geometry();
478 auto geo_outside = cell_outside.geometry();
481 auto geo_in_inside = ig.geometryInInside();
482 auto geo_in_outside = ig.geometryInOutside();
487 auto local_inside = ref_el_inside.position(0,0);
488 auto local_outside = ref_el_outside.position(0,0);
493 auto mu_s = param.mu(cell_inside,local_inside);
494 auto mu_n = param.mu(cell_outside,local_outside);
495 auto eps_s = param.eps(cell_inside,local_inside);
496 auto eps_n = param.eps(cell_outside,local_outside);
501 const auto& n_F = ig.centerUnitOuterNormal();
504 Dune::FieldMatrix<DF,dim*2,dim*2> R_s;
506 Dune::FieldMatrix<DF,dim*2,dim*2> Dplus_s(0.0);
507 Dplus_s[0][0] = 1.0/sqrt(eps_s*mu_s);
508 Dplus_s[1][1] = 1.0/sqrt(eps_s*mu_s);
509 Dune::FieldMatrix<DF,dim*2,dim*2> Aplus_s(R_s);
510 Aplus_s.rightmultiply(Dplus_s);
512 Aplus_s.rightmultiply(R_s);
515 Dune::FieldMatrix<DF,dim*2,dim*2> R_n;
517 Dune::FieldMatrix<DF,dim*2,dim*2> Dminus_n(0.0);
518 Dminus_n[2][2] = -1.0/sqrt(eps_n*mu_n);
519 Dminus_n[3][3] = -1.0/sqrt(eps_n*mu_n);
520 Dune::FieldMatrix<DF,dim*2,dim*2> Aminus_n(R_n);
521 Aminus_n.rightmultiply(Dminus_n);
523 Aminus_n.rightmultiply(R_n);
526 Dune::FieldVector<RF,dim*2> u_s(0.0);
527 Dune::FieldVector<RF,dim*2> u_n(0.0);
528 Dune::FieldVector<RF,dim*2> f(0.0);
533 const int order_s = dgspace_s.finiteElement().localBasis().order();
534 const int order_n = dgspace_n.finiteElement().localBasis().order();
535 const int intorder = overintegration+1+2*max(order_s,order_n);
539 const auto& iplocal_s = geo_in_inside.global(qp.position());
540 const auto& iplocal_n = geo_in_outside.global(qp.position());
543 const auto& phi_s = cache[order_s].evaluateFunction(iplocal_s,dgspace_s.finiteElement().localBasis());
544 const auto& phi_n = cache[order_n].evaluateFunction(iplocal_n,dgspace_n.finiteElement().localBasis());
547 for (size_type i=0; i<
dim*2; i++){
549 for (size_type k=0; k<dgspace_s.size(); k++)
550 u_s[i] += x_s(lfsv_s.child(i),k)*phi_s[k];
552 for (size_type i=0; i<dim*2; i++){
554 for (size_type k=0; k<dgspace_n.size(); k++)
555 u_n[i] += x_n(lfsv_n.child(i),k)*phi_n[k];
568 auto factor = qp.weight() * geo.integrationElement(qp.position());
569 for (size_type k=0; k<dgspace_s.size(); k++)
570 for (size_type i=0; i<dim*2; i++)
571 r_s.accumulate(lfsv_s.child(i),k,f[i]*phi_s[k]*factor);
572 for (size_type k=0; k<dgspace_n.size(); k++)
573 for (size_type i=0; i<dim*2; i++)
574 r_n.accumulate(lfsv_n.child(i),k,-f[i]*phi_n[k]*factor);
586 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
588 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
592 using namespace Indices;
593 using DGSpace = TypeTree::Child<LFSV,0>;
594 using DF =
typename DGSpace::Traits::FiniteElementType::
595 Traits::LocalBasisType::Traits::DomainFieldType;
596 using RF =
typename DGSpace::Traits::FiniteElementType::
597 Traits::LocalBasisType::Traits::RangeFieldType;
598 using size_type =
typename DGSpace::Traits::SizeType;
601 const auto& dgspace_s = child(lfsv_s,_0);
604 const auto& cell_inside = ig.inside();
607 auto geo = ig.geometry();
608 auto geo_inside = cell_inside.geometry();
611 auto geo_in_inside = ig.geometryInInside();
615 auto local_inside = ref_el_inside.position(0,0);
616 auto mu_s = param.mu(cell_inside,local_inside);
617 auto eps_s = param.eps(cell_inside,local_inside);
621 const auto& n_F = ig.centerUnitOuterNormal();
624 Dune::FieldMatrix<DF,dim*2,dim*2> R_s;
626 Dune::FieldMatrix<DF,dim*2,dim*2> Dplus_s(0.0);
627 Dplus_s[0][0] = 1.0/sqrt(eps_s*mu_s);
628 Dplus_s[1][1] = 1.0/sqrt(eps_s*mu_s);
629 Dune::FieldMatrix<DF,dim*2,dim*2> Aplus_s(R_s);
630 Aplus_s.rightmultiply(Dplus_s);
632 Aplus_s.rightmultiply(R_s);
635 Dune::FieldMatrix<DF,dim*2,dim*2> R_n;
637 Dune::FieldMatrix<DF,dim*2,dim*2> Dminus_n(0.0);
638 Dminus_n[2][2] = -1.0/sqrt(eps_s*mu_s);
639 Dminus_n[3][3] = -1.0/sqrt(eps_s*mu_s);
640 Dune::FieldMatrix<DF,dim*2,dim*2> Aminus_n(R_n);
641 Aminus_n.rightmultiply(Dminus_n);
643 Aminus_n.rightmultiply(R_n);
646 Dune::FieldVector<RF,dim*2> u_s(0.0);
647 Dune::FieldVector<RF,dim*2> u_n;
648 Dune::FieldVector<RF,dim*2> f(0.0);
653 const int order_s = dgspace_s.finiteElement().localBasis().order();
654 const int intorder = overintegration+1+2*order_s;
658 const auto& iplocal_s = geo_in_inside.global(qp.position());
661 const auto& phi_s = cache[order_s].evaluateFunction
662 (iplocal_s,dgspace_s.finiteElement().localBasis());
665 for (size_type i=0; i<
dim*2; i++){
667 for (size_type k=0; k<dgspace_s.size(); k++)
668 u_s[i] += x_s(lfsv_s.child(i),k)*phi_s[k];
673 u_n = (param.g(ig.intersection(),qp.position(),u_s));
684 auto factor = qp.weight() * geo.integrationElement(qp.position());
685 for (size_type k=0; k<dgspace_s.size(); k++)
686 for (size_type i=0; i<dim*2; i++)
687 r_s.accumulate(lfsv_s.child(i),k,f[i]*phi_s[k]*factor);
696 template<
typename EG,
typename LFSV,
typename R>
700 using namespace Indices;
701 using DGSpace = TypeTree::Child<LFSV,0>;
702 using size_type =
typename DGSpace::Traits::SizeType;
705 const auto& dgspace = child(lfsv,_0);
708 const auto& cell = eg.entity();
711 auto geo = eg.geometry();
714 const int order_s = dgspace.finiteElement().localBasis().order();
715 const int intorder = overintegration+2*order_s;
719 auto j = param.j(cell,qp.position());
722 const auto& phi = cache[order_s].evaluateFunction(qp.position(),dgspace.finiteElement().localBasis());
725 auto factor = qp.weight() * geo.integrationElement(qp.position());
726 for (size_type k=0; k<
dim*2; k++)
727 for (size_type i=0; i<dgspace.size(); i++)
728 r.accumulate(lfsv.child(k),i,-j[k]*phi[i]*factor);
733 void setTime (
typename T::Traits::RangeFieldType t)
738 void preStep (
typename T::Traits::RangeFieldType time,
typename T::Traits::RangeFieldType dt,
744 void preStage (
typename T::Traits::RangeFieldType time,
int r)
754 typename T::Traits::RangeFieldType
suggestTimestep (
typename T::Traits::RangeFieldType dt)
const 762 typedef typename FEM::Traits::FiniteElementType::Traits::LocalBasisType LocalBasisType;
764 std::vector<Cache> cache;
780 template<
typename T,
typename FEM>
786 enum {
dim = T::Traits::GridViewType::dimension };
789 enum { doPatternVolume =
true };
792 enum { doAlphaVolume =
true };
795 : param(param_), overintegration(overintegration_), cache(20)
799 template<
typename LFSU,
typename LFSV,
typename LocalPattern>
801 LocalPattern& pattern)
const 807 for (
size_t k=0; k<TypeTree::degree(lfsv); k++)
808 for (
size_t i=0; i<lfsv.child(k).size(); ++i)
809 for (
size_t j=0; j<lfsu.child(k).size(); ++j)
810 pattern.addLink(lfsv.child(k),i,lfsu.child(k),j);
814 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
815 void alpha_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const 818 using namespace Indices;
819 using DGSpace = TypeTree::Child<LFSV,0>;
820 using RF =
typename DGSpace::Traits::FiniteElementType::
821 Traits::LocalBasisType::Traits::RangeFieldType;
822 using size_type =
typename DGSpace::Traits::SizeType;
825 const auto& dgspace = child(lfsv,_0);
828 auto geo = eg.geometry();
831 Dune::FieldVector<RF,dim*2> u(0.0);
834 const int order = dgspace.finiteElement().localBasis().order();
835 const int intorder = overintegration+2*order;
839 const auto& phi = cache[order].evaluateFunction(qp.position(),dgspace.finiteElement().localBasis());
842 for (size_type k=0; k<
dim*2; k++){
844 for (size_type j=0; j<dgspace.size(); j++)
845 u[k] += x(lfsv.child(k),j)*phi[j];
849 auto factor = qp.weight() * geo.integrationElement(qp.position());
850 for (size_type k=0; k<dim*2; k++)
851 for (size_type i=0; i<dgspace.size(); i++)
852 r.accumulate(lfsv.child(k),i,u[k]*phi[i]*factor);
857 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
862 using namespace Indices;
863 using DGSpace = TypeTree::Child<LFSV,0>;
864 using size_type =
typename DGSpace::Traits::SizeType;
867 const auto& dgspace = child(lfsv,_0);
870 auto geo = eg.geometry();
873 const int order = dgspace.finiteElement().localBasis().order();
874 const int intorder = overintegration+2*order;
878 const auto& phi = cache[order].evaluateFunction(qp.position(),dgspace.finiteElement().localBasis());
881 auto factor = qp.weight() * geo.integrationElement(qp.position());
882 for (size_type k=0; k<
dim*2; k++)
883 for (size_type i=0; i<dgspace.size(); i++)
884 for (size_type j=0; j<dgspace.size(); j++)
885 mat.accumulate(lfsv.child(k),i,lfsu.child(k),j,phi[j]*phi[i]*factor);
892 typedef typename FEM::Traits::FiniteElementType::Traits::LocalBasisType LocalBasisType;
894 std::vector<Cache> cache;
900 #endif // DUNE_PDELAB_LOCALOPERATOR_MAXWELLDG_HH T::Traits::RangeFieldType suggestTimestep(typename T::Traits::RangeFieldType dt) const
to be called once before each stage
Definition: maxwelldg.hh:754
const IG & ig
Definition: constraints.hh:148
static const int dim
Definition: adaptivity.hh:83
const Entity & e
Definition: localfunctionspace.hh:111
Implement jacobian_skeleton() based on alpha_skeleton()
Definition: numericaljacobian.hh:156
DGMaxwellSpatialOperator(T ¶m_, int overintegration_=0)
Definition: maxwelldg.hh:341
const std::string s
Definition: function.hh:830
Implements linear and nonlinear versions of jacobian_apply_boundary() based on alpha_boundary() ...
Definition: numericaljacobianapply.hh:285
void alpha_skeleton(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, const LFSU &lfsu_n, const X &x_n, const LFSV &lfsv_n, R &r_s, R &r_n) const
Definition: maxwelldg.hh:450
void preStep(typename T::Traits::RangeFieldType time, typename T::Traits::RangeFieldType dt, int stages)
to be called once before each time step
Definition: maxwelldg.hh:738
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: maxwelldg.hh:697
Spatial local operator for discontinuous Galerkin method for Maxwells Equations.
Definition: maxwelldg.hh:315
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
void alpha_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, R &r_s) const
Definition: maxwelldg.hh:587
static void eigenvalues(T1 eps, T1 mu, const Dune::FieldVector< T2, 2 *dim > &e)
Definition: maxwelldg.hh:54
sparsity pattern generator
Definition: pattern.hh:29
Implements linear and nonlinear versions of jacobian_apply_volume() based on alpha_volume() ...
Definition: numericaljacobianapply.hh:32
void preStage(typename T::Traits::RangeFieldType time, int r)
to be called once before each stage
Definition: maxwelldg.hh:744
void pattern_volume(const LFSU &lfsu, const LFSV &lfsv, LocalPattern &pattern) const
Definition: maxwelldg.hh:800
static void eigenvectors(T1 eps, T1 mu, const Dune::FieldVector< T2, dim > &n, Dune::FieldMatrix< T3, 2 *dim, 2 *dim > &R)
Definition: maxwelldg.hh:77
void postStage()
to be called once at the end of each stage
Definition: maxwelldg.hh:749
DGMaxwellTemporalOperator(T ¶m_, int overintegration_=0)
Definition: maxwelldg.hh:794
sparsity pattern generator
Definition: pattern.hh:13
Definition: maxwelldg.hh:781
QuadratureRuleWrapper< QuadratureRule< typename Geometry::ctype, Geometry::mydimension > > quadratureRule(const Geometry &geo, std::size_t order, QuadratureType::Enum quadrature_type=QuadratureType::GaussLegendre)
Returns a quadrature rule for the given geometry.
Definition: quadraturerules.hh:117
Definition: maxwelldg.hh:30
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: maxwelldg.hh:858
ReferenceElementWrapper< ReferenceElement< typename Geometry::ctype, Geometry::mydimension > > referenceElement(const Geometry &geo)
Returns the reference element for the given geometry.
Definition: referenceelements.hh:144
store values of basis functions and gradients in a cache
Definition: localbasiscache.hh:17
Default flags for all local operators.
Definition: flags.hh:18
static const unsigned int value
Definition: gridfunctionspace/tags.hh:139
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: maxwelldg.hh:815
Implements linear and nonlinear versions of jacobian_apply_skeleton() based on alpha_skeleton() ...
Definition: numericaljacobianapply.hh:180
void setTime(typename T::Traits::RangeFieldType t)
set time in parameter class
Definition: maxwelldg.hh:733
Implement jacobian_boundary() based on alpha_boundary()
Definition: numericaljacobian.hh:250
Implement jacobian_volume() based on alpha_volume()
Definition: numericaljacobian.hh:31
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: maxwelldg.hh:348