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> 10 #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++)
303 template<
typename T,
typename FEM>
316 enum {
dim = T::Traits::GridViewType::dimension };
320 enum { doPatternVolume =
true };
321 enum { doPatternSkeleton =
true };
324 enum { doAlphaVolume =
true };
325 enum { doAlphaSkeleton =
true };
326 enum { doAlphaBoundary =
true };
327 enum { doLambdaVolume =
true };
331 : param(param_), overintegration(overintegration_), cache(20)
336 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
337 void alpha_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const 340 using namespace Indices;
341 using DGSpace = TypeTree::Child<LFSV,_0>;
342 using RF =
typename DGSpace::Traits::FiniteElementType::Traits
343 ::LocalBasisType::Traits::RangeFieldType;
344 using size_type =
typename DGSpace::Traits::SizeType;
348 "need exactly dim*2 components!");
351 const auto& dgspace = child(lfsv,_0);
354 const auto& cell = eg.entity();
357 auto geo = eg.geometry();
361 auto localcenter = ref_el.position(0,0);
362 auto mu = param.mu(cell,localcenter);
363 auto eps = param.eps(cell,localcenter);
364 auto sigma = param.sigma(cell,localcenter);
366 auto epsinv = 1.0/eps;
371 typename EG::Geometry::JacobianInverseTransposed jac;
374 Dune::FieldVector<RF,dim*2> u(0.0);
375 std::vector<Dune::FieldVector<RF,dim> > gradphi(dgspace.size());
378 const int order = dgspace.finiteElement().localBasis().order();
379 const int intorder = overintegration+2*order;
383 const auto& phi = cache[order].evaluateFunction
384 (qp.position(), dgspace.finiteElement().localBasis());
387 for (size_type k=0; k<
dim*2; k++){
389 for (size_type j=0; j<dgspace.size(); j++)
390 u[k] += x(lfsv.child(k),j)*phi[j];
395 const auto& js = cache[order].evaluateJacobian
396 (qp.position(), dgspace.finiteElement().localBasis());
399 jac = geo.jacobianInverseTransposed(qp.position());
400 for (size_type i=0; i<dgspace.size(); i++)
401 jac.mv(js[i][0],gradphi[i]);
404 auto factor = qp.weight() * geo.integrationElement(qp.position());
406 Dune::FieldMatrix<RF,dim*2,dim> F;
407 F[0][0] = 0; F[0][1] = -muinv*u[5]; F[0][2] = muinv*u[4];
408 F[1][0] = muinv*u[5]; F[1][1] = 0; F[1][2] = -muinv*u[3];
409 F[2][0] =-muinv*u[4]; F[2][1] = muinv*u[3]; F[2][2] = 0;
410 F[3][0] = 0; F[3][1] = epsinv*u[2]; F[3][2] = -epsinv*u[1];
411 F[4][0] = -epsinv*u[2]; F[4][1] = 0; F[4][2] = epsinv*u[0];
412 F[5][0] = epsinv*u[1]; F[5][1] = -epsinv*u[0]; F[5][2] = 0;
415 for (size_type i=0; i<dim*2; i++)
417 for (size_type k=0; k<dgspace.size(); k++)
419 for (size_type j=0; j<
dim; j++)
420 r.accumulate(lfsv.child(i),k,-F[i][j]*gradphi[k][j]*factor);
423 for (size_type i=0; i<
dim; i++)
425 for (size_type k=0; k<dgspace.size(); k++)
426 r.accumulate(lfsv.child(i),k,(sigma/eps)*u[i]*phi[k]*factor);
436 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
438 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
439 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
440 R& r_s, R& r_n)
const 446 using namespace Indices;
447 using DGSpace = TypeTree::Child<LFSV,_0>;
448 using DF =
typename DGSpace::Traits::FiniteElementType::
449 Traits::LocalBasisType::Traits::DomainFieldType;
450 using RF =
typename DGSpace::Traits::FiniteElementType::
451 Traits::LocalBasisType::Traits::RangeFieldType;
452 using size_type =
typename DGSpace::Traits::SizeType;
455 const auto& dgspace_s = child(lfsv_s,_0);
456 const auto& dgspace_n = child(lfsv_n,_0);
459 const auto& cell_inside = ig.inside();
460 const auto& cell_outside = ig.outside();
463 auto geo = ig.geometry();
464 auto geo_inside = cell_inside.geometry();
465 auto geo_outside = cell_outside.geometry();
468 auto geo_in_inside = ig.geometryInInside();
469 auto geo_in_outside = ig.geometryInOutside();
474 auto local_inside = ref_el_inside.position(0,0);
475 auto local_outside = ref_el_outside.position(0,0);
480 auto mu_s = param.mu(cell_inside,local_inside);
481 auto mu_n = param.mu(cell_outside,local_outside);
482 auto eps_s = param.eps(cell_inside,local_inside);
483 auto eps_n = param.eps(cell_outside,local_outside);
488 const auto& n_F = ig.centerUnitOuterNormal();
491 Dune::FieldMatrix<DF,dim*2,dim*2> R_s;
493 Dune::FieldMatrix<DF,dim*2,dim*2> Dplus_s(0.0);
494 Dplus_s[0][0] = 1.0/sqrt(eps_s*mu_s);
495 Dplus_s[1][1] = 1.0/sqrt(eps_s*mu_s);
496 Dune::FieldMatrix<DF,dim*2,dim*2> Aplus_s(R_s);
497 Aplus_s.rightmultiply(Dplus_s);
499 Aplus_s.rightmultiply(R_s);
502 Dune::FieldMatrix<DF,dim*2,dim*2> R_n;
504 Dune::FieldMatrix<DF,dim*2,dim*2> Dminus_n(0.0);
505 Dminus_n[2][2] = -1.0/sqrt(eps_n*mu_n);
506 Dminus_n[3][3] = -1.0/sqrt(eps_n*mu_n);
507 Dune::FieldMatrix<DF,dim*2,dim*2> Aminus_n(R_n);
508 Aminus_n.rightmultiply(Dminus_n);
510 Aminus_n.rightmultiply(R_n);
513 Dune::FieldVector<RF,dim*2> u_s(0.0);
514 Dune::FieldVector<RF,dim*2> u_n(0.0);
515 Dune::FieldVector<RF,dim*2> f(0.0);
520 const int order_s = dgspace_s.finiteElement().localBasis().order();
521 const int order_n = dgspace_n.finiteElement().localBasis().order();
522 const int intorder = overintegration+1+2*max(order_s,order_n);
526 const auto& iplocal_s = geo_in_inside.global(qp.position());
527 const auto& iplocal_n = geo_in_outside.global(qp.position());
530 const auto& phi_s = cache[order_s].evaluateFunction(iplocal_s,dgspace_s.finiteElement().localBasis());
531 const auto& phi_n = cache[order_n].evaluateFunction(iplocal_n,dgspace_n.finiteElement().localBasis());
534 for (size_type i=0; i<
dim*2; i++){
536 for (size_type k=0; k<dgspace_s.size(); k++)
537 u_s[i] += x_s(lfsv_s.child(i),k)*phi_s[k];
539 for (size_type i=0; i<dim*2; i++){
541 for (size_type k=0; k<dgspace_n.size(); k++)
542 u_n[i] += x_n(lfsv_n.child(i),k)*phi_n[k];
555 auto factor = qp.weight() * geo.integrationElement(qp.position());
556 for (size_type k=0; k<dgspace_s.size(); k++)
557 for (size_type i=0; i<dim*2; i++)
558 r_s.accumulate(lfsv_s.child(i),k,f[i]*phi_s[k]*factor);
559 for (size_type k=0; k<dgspace_n.size(); k++)
560 for (size_type i=0; i<dim*2; i++)
561 r_n.accumulate(lfsv_n.child(i),k,-f[i]*phi_n[k]*factor);
573 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
575 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
579 using namespace Indices;
580 using DGSpace = TypeTree::Child<LFSV,_0>;
581 using DF =
typename DGSpace::Traits::FiniteElementType::
582 Traits::LocalBasisType::Traits::DomainFieldType;
583 using RF =
typename DGSpace::Traits::FiniteElementType::
584 Traits::LocalBasisType::Traits::RangeFieldType;
585 using size_type =
typename DGSpace::Traits::SizeType;
588 const auto& dgspace_s = child(lfsv_s,_0);
591 const auto& cell_inside = ig.inside();
594 auto geo = ig.geometry();
595 auto geo_inside = cell_inside.geometry();
598 auto geo_in_inside = ig.geometryInInside();
602 auto local_inside = ref_el_inside.position(0,0);
603 auto mu_s = param.mu(cell_inside,local_inside);
604 auto eps_s = param.eps(cell_inside,local_inside);
608 const auto& n_F = ig.centerUnitOuterNormal();
611 Dune::FieldMatrix<DF,dim*2,dim*2> R_s;
613 Dune::FieldMatrix<DF,dim*2,dim*2> Dplus_s(0.0);
614 Dplus_s[0][0] = 1.0/sqrt(eps_s*mu_s);
615 Dplus_s[1][1] = 1.0/sqrt(eps_s*mu_s);
616 Dune::FieldMatrix<DF,dim*2,dim*2> Aplus_s(R_s);
617 Aplus_s.rightmultiply(Dplus_s);
619 Aplus_s.rightmultiply(R_s);
622 Dune::FieldMatrix<DF,dim*2,dim*2> R_n;
624 Dune::FieldMatrix<DF,dim*2,dim*2> Dminus_n(0.0);
625 Dminus_n[2][2] = -1.0/sqrt(eps_s*mu_s);
626 Dminus_n[3][3] = -1.0/sqrt(eps_s*mu_s);
627 Dune::FieldMatrix<DF,dim*2,dim*2> Aminus_n(R_n);
628 Aminus_n.rightmultiply(Dminus_n);
630 Aminus_n.rightmultiply(R_n);
633 Dune::FieldVector<RF,dim*2> u_s(0.0);
634 Dune::FieldVector<RF,dim*2> u_n;
635 Dune::FieldVector<RF,dim*2> f(0.0);
640 const int order_s = dgspace_s.finiteElement().localBasis().order();
641 const int intorder = overintegration+1+2*order_s;
645 const auto& iplocal_s = geo_in_inside.global(qp.position());
648 const auto& phi_s = cache[order_s].evaluateFunction
649 (iplocal_s,dgspace_s.finiteElement().localBasis());
652 for (size_type i=0; i<
dim*2; i++){
654 for (size_type k=0; k<dgspace_s.size(); k++)
655 u_s[i] += x_s(lfsv_s.child(i),k)*phi_s[k];
660 u_n = (param.g(ig.intersection(),qp.position(),u_s));
671 auto factor = qp.weight() * geo.integrationElement(qp.position());
672 for (size_type k=0; k<dgspace_s.size(); k++)
673 for (size_type i=0; i<dim*2; i++)
674 r_s.accumulate(lfsv_s.child(i),k,f[i]*phi_s[k]*factor);
683 template<
typename EG,
typename LFSV,
typename R>
687 using namespace Indices;
688 using DGSpace = TypeTree::Child<LFSV,_0>;
689 using size_type =
typename DGSpace::Traits::SizeType;
692 const auto& dgspace = child(lfsv,_0);
695 const auto& cell = eg.entity();
698 auto geo = eg.geometry();
701 const int order_s = dgspace.finiteElement().localBasis().order();
702 const int intorder = overintegration+2*order_s;
706 auto j = param.j(cell,qp.position());
709 const auto& phi = cache[order_s].evaluateFunction(qp.position(),dgspace.finiteElement().localBasis());
712 auto factor = qp.weight() * geo.integrationElement(qp.position());
713 for (size_type k=0; k<
dim*2; k++)
714 for (size_type i=0; i<dgspace.size(); i++)
715 r.accumulate(lfsv.child(k),i,-j[k]*phi[i]*factor);
720 void setTime (
typename T::Traits::RangeFieldType t)
725 void preStep (
typename T::Traits::RangeFieldType time,
typename T::Traits::RangeFieldType dt,
731 void preStage (
typename T::Traits::RangeFieldType time,
int r)
741 typename T::Traits::RangeFieldType
suggestTimestep (
typename T::Traits::RangeFieldType dt)
const 749 typedef typename FEM::Traits::FiniteElementType::Traits::LocalBasisType LocalBasisType;
751 std::vector<Cache> cache;
767 template<
typename T,
typename FEM>
773 enum {
dim = T::Traits::GridViewType::dimension };
776 enum { doPatternVolume =
true };
779 enum { doAlphaVolume =
true };
782 : param(param_), overintegration(overintegration_), cache(20)
786 template<
typename LFSU,
typename LFSV,
typename LocalPattern>
788 LocalPattern& pattern)
const 794 for (
size_t k=0; k<TypeTree::degree(lfsv); k++)
795 for (
size_t i=0; i<lfsv.child(k).size(); ++i)
796 for (
size_t j=0; j<lfsu.child(k).size(); ++j)
797 pattern.addLink(lfsv.child(k),i,lfsu.child(k),j);
801 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
802 void alpha_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const 805 using namespace Indices;
806 using DGSpace = TypeTree::Child<LFSV,_0>;
807 using RF =
typename DGSpace::Traits::FiniteElementType::
808 Traits::LocalBasisType::Traits::RangeFieldType;
809 using size_type =
typename DGSpace::Traits::SizeType;
812 const auto& dgspace = child(lfsv,_0);
815 auto geo = eg.geometry();
818 Dune::FieldVector<RF,dim*2> u(0.0);
821 const int order = dgspace.finiteElement().localBasis().order();
822 const int intorder = overintegration+2*order;
826 const auto& phi = cache[order].evaluateFunction(qp.position(),dgspace.finiteElement().localBasis());
829 for (size_type k=0; k<
dim*2; k++){
831 for (size_type j=0; j<dgspace.size(); j++)
832 u[k] += x(lfsv.child(k),j)*phi[j];
836 auto factor = qp.weight() * geo.integrationElement(qp.position());
837 for (size_type k=0; k<dim*2; k++)
838 for (size_type i=0; i<dgspace.size(); i++)
839 r.accumulate(lfsv.child(k),i,u[k]*phi[i]*factor);
844 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
849 using namespace Indices;
850 using DGSpace = TypeTree::Child<LFSV,_0>;
851 using size_type =
typename DGSpace::Traits::SizeType;
854 const auto& dgspace = child(lfsv,_0);
857 auto geo = eg.geometry();
860 const int order = dgspace.finiteElement().localBasis().order();
861 const int intorder = overintegration+2*order;
865 const auto& phi = cache[order].evaluateFunction(qp.position(),dgspace.finiteElement().localBasis());
868 auto factor = qp.weight() * geo.integrationElement(qp.position());
869 for (size_type k=0; k<
dim*2; k++)
870 for (size_type i=0; i<dgspace.size(); i++)
871 for (size_type j=0; j<dgspace.size(); j++)
872 mat.accumulate(lfsv.child(k),i,lfsu.child(k),j,phi[j]*phi[i]*factor);
879 typedef typename FEM::Traits::FiniteElementType::Traits::LocalBasisType LocalBasisType;
881 std::vector<Cache> cache;
887 #endif // DUNE_PDELAB_LOCALOPERATOR_MAXWELLDG_HH 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
T::Traits::RangeFieldType suggestTimestep(typename T::Traits::RangeFieldType dt) const
to be called once before each stage
Definition: maxwelldg.hh:741
DGMaxwellTemporalOperator(T ¶m_, int overintegration_=0)
Definition: maxwelldg.hh:781
const IG & ig
Definition: constraints.hh:148
static const int dim
Definition: adaptivity.hh:83
void setTime(typename T::Traits::RangeFieldType t)
set time in parameter class
Definition: maxwelldg.hh:720
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
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:802
sparsity pattern generator
Definition: pattern.hh:29
const std::string s
Definition: function.hh:1101
Implement jacobian_skeleton() based on alpha_skeleton()
Definition: numericaljacobian.hh:156
Implements linear and nonlinear versions of jacobian_apply_boundary() based on alpha_boundary() ...
Definition: numericaljacobianapply.hh:285
Implements linear and nonlinear versions of jacobian_apply_volume() based on alpha_volume() ...
Definition: numericaljacobianapply.hh:32
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: maxwelldg.hh:684
void pattern_volume(const LFSU &lfsu, const LFSV &lfsv, LocalPattern &pattern) const
Definition: maxwelldg.hh:787
void postStage()
to be called once at the end of each stage
Definition: maxwelldg.hh:736
DGMaxwellSpatialOperator(T ¶m_, int overintegration_=0)
Definition: maxwelldg.hh:330
Implements linear and nonlinear versions of jacobian_apply_skeleton() based on alpha_skeleton() ...
Definition: numericaljacobianapply.hh:180
store values of basis functions and gradients in a cache
Definition: localbasiscache.hh:17
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: maxwelldg.hh:845
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
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:725
Definition: maxwelldg.hh:768
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: maxwelldg.hh:337
sparsity pattern generator
Definition: pattern.hh:13
void preStage(typename T::Traits::RangeFieldType time, int r)
to be called once before each stage
Definition: maxwelldg.hh:731
const Entity & e
Definition: localfunctionspace.hh:111
static void eigenvalues(T1 eps, T1 mu, const Dune::FieldVector< T2, 2 *dim > &e)
Definition: maxwelldg.hh:53
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:574
Implement jacobian_boundary() based on alpha_boundary()
Definition: numericaljacobian.hh:250
Definition: maxwelldg.hh:29
Implement jacobian_volume() based on alpha_volume()
Definition: numericaljacobian.hh:31
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:437
ReferenceElementWrapper< ReferenceElement< typename Geometry::ctype, Geometry::mydimension > > referenceElement(const Geometry &geo)
Returns the reference element for the given geometry.
Definition: referenceelements.hh:144
Default flags for all local operators.
Definition: flags.hh:18
Definition: maxwelldg.hh:304