2 #ifndef DUNE_PDELAB_MAXWELLDG_HH 3 #define DUNE_PDELAB_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 TypeTree::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;
347 static_assert(LFSV::CHILDREN ==
dim*2,
348 "need exactly dim*2 components!");
351 using namespace TypeTree::Indices;
352 const auto& dgspace = child(lfsv,_0);
355 const auto& cell = eg.entity();
358 auto geo = eg.geometry();
362 auto localcenter = ref_el.position(0,0);
363 auto mu = param.mu(cell,localcenter);
364 auto eps = param.eps(cell,localcenter);
365 auto sigma = param.sigma(cell,localcenter);
367 auto epsinv = 1.0/eps;
372 typename EG::Geometry::JacobianInverseTransposed jac;
375 Dune::FieldVector<RF,dim*2> u(0.0);
376 std::vector<Dune::FieldVector<RF,dim> > gradphi(dgspace.size());
379 const int order = dgspace.finiteElement().localBasis().order();
380 const int intorder = overintegration+2*order;
384 const auto& phi = cache[order].evaluateFunction
385 (qp.position(), dgspace.finiteElement().localBasis());
388 for (size_type k=0; k<
dim*2; k++){
390 for (size_type j=0; j<dgspace.size(); j++)
391 u[k] += x(lfsv.child(k),j)*phi[j];
396 const auto& js = cache[order].evaluateJacobian
397 (qp.position(), dgspace.finiteElement().localBasis());
400 jac = geo.jacobianInverseTransposed(qp.position());
401 for (size_type i=0; i<dgspace.size(); i++)
402 jac.mv(js[i][0],gradphi[i]);
405 auto factor = qp.weight() * geo.integrationElement(qp.position());
407 Dune::FieldMatrix<RF,dim*2,dim> F;
408 F[0][0] = 0; F[0][1] = -muinv*u[5]; F[0][2] = muinv*u[4];
409 F[1][0] = muinv*u[5]; F[1][1] = 0; F[1][2] = -muinv*u[3];
410 F[2][0] =-muinv*u[4]; F[2][1] = muinv*u[3]; F[2][2] = 0;
411 F[3][0] = 0; F[3][1] = epsinv*u[2]; F[3][2] = -epsinv*u[1];
412 F[4][0] = -epsinv*u[2]; F[4][1] = 0; F[4][2] = epsinv*u[0];
413 F[5][0] = epsinv*u[1]; F[5][1] = -epsinv*u[0]; F[5][2] = 0;
416 for (size_type i=0; i<dim*2; i++)
418 for (size_type k=0; k<dgspace.size(); k++)
420 for (size_type j=0; j<
dim; j++)
421 r.accumulate(lfsv.child(i),k,-F[i][j]*gradphi[k][j]*factor);
424 for (size_type i=0; i<
dim; i++)
426 for (size_type k=0; k<dgspace.size(); k++)
427 r.accumulate(lfsv.child(i),k,(sigma/eps)*u[i]*phi[k]*factor);
437 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
439 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
440 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
441 R& r_s, R& r_n)
const 447 using namespace TypeTree::Indices;
448 using DGSpace = TypeTree::Child<LFSV,_0>;
449 using DF =
typename DGSpace::Traits::FiniteElementType::
450 Traits::LocalBasisType::Traits::DomainFieldType;
451 using RF =
typename DGSpace::Traits::FiniteElementType::
452 Traits::LocalBasisType::Traits::RangeFieldType;
453 using size_type =
typename DGSpace::Traits::SizeType;
456 using namespace TypeTree::Indices;
457 const auto& dgspace_s = child(lfsv_s,_0);
458 const auto& dgspace_n = child(lfsv_n,_0);
461 const auto& cell_inside = ig.inside();
462 const auto& cell_outside = ig.outside();
465 auto geo = ig.geometry();
466 auto geo_inside = cell_inside.geometry();
467 auto geo_outside = cell_outside.geometry();
470 auto geo_in_inside = ig.geometryInInside();
471 auto geo_in_outside = ig.geometryInOutside();
476 auto local_inside = ref_el_inside.position(0,0);
477 auto local_outside = ref_el_outside.position(0,0);
482 auto mu_s = param.mu(cell_inside,local_inside);
483 auto mu_n = param.mu(cell_outside,local_outside);
484 auto eps_s = param.eps(cell_inside,local_inside);
485 auto eps_n = param.eps(cell_outside,local_outside);
490 const auto& n_F = ig.centerUnitOuterNormal();
493 Dune::FieldMatrix<DF,dim*2,dim*2> R_s;
495 Dune::FieldMatrix<DF,dim*2,dim*2> Dplus_s(0.0);
496 Dplus_s[0][0] = 1.0/sqrt(eps_s*mu_s);
497 Dplus_s[1][1] = 1.0/sqrt(eps_s*mu_s);
498 Dune::FieldMatrix<DF,dim*2,dim*2> Aplus_s(R_s);
499 Aplus_s.rightmultiply(Dplus_s);
501 Aplus_s.rightmultiply(R_s);
504 Dune::FieldMatrix<DF,dim*2,dim*2> R_n;
506 Dune::FieldMatrix<DF,dim*2,dim*2> Dminus_n(0.0);
507 Dminus_n[2][2] = -1.0/sqrt(eps_n*mu_n);
508 Dminus_n[3][3] = -1.0/sqrt(eps_n*mu_n);
509 Dune::FieldMatrix<DF,dim*2,dim*2> Aminus_n(R_n);
510 Aminus_n.rightmultiply(Dminus_n);
512 Aminus_n.rightmultiply(R_n);
515 Dune::FieldVector<RF,dim*2> u_s(0.0);
516 Dune::FieldVector<RF,dim*2> u_n(0.0);
517 Dune::FieldVector<RF,dim*2> f(0.0);
522 const int order_s = dgspace_s.finiteElement().localBasis().order();
523 const int order_n = dgspace_n.finiteElement().localBasis().order();
524 const int intorder = overintegration+1+2*max(order_s,order_n);
528 const auto& iplocal_s = geo_in_inside.global(qp.position());
529 const auto& iplocal_n = geo_in_outside.global(qp.position());
532 const auto& phi_s = cache[order_s].evaluateFunction(iplocal_s,dgspace_s.finiteElement().localBasis());
533 const auto& phi_n = cache[order_n].evaluateFunction(iplocal_n,dgspace_n.finiteElement().localBasis());
536 for (size_type i=0; i<
dim*2; i++){
538 for (size_type k=0; k<dgspace_s.size(); k++)
539 u_s[i] += x_s(lfsv_s.child(i),k)*phi_s[k];
541 for (size_type i=0; i<dim*2; i++){
543 for (size_type k=0; k<dgspace_n.size(); k++)
544 u_n[i] += x_n(lfsv_n.child(i),k)*phi_n[k];
557 auto factor = qp.weight() * geo.integrationElement(qp.position());
558 for (size_type k=0; k<dgspace_s.size(); k++)
559 for (size_type i=0; i<dim*2; i++)
560 r_s.accumulate(lfsv_s.child(i),k,f[i]*phi_s[k]*factor);
561 for (size_type k=0; k<dgspace_n.size(); k++)
562 for (size_type i=0; i<dim*2; i++)
563 r_n.accumulate(lfsv_n.child(i),k,-f[i]*phi_n[k]*factor);
575 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
577 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
581 using namespace TypeTree::Indices;
582 using DGSpace = TypeTree::Child<LFSV,_0>;
583 using DF =
typename DGSpace::Traits::FiniteElementType::
584 Traits::LocalBasisType::Traits::DomainFieldType;
585 using RF =
typename DGSpace::Traits::FiniteElementType::
586 Traits::LocalBasisType::Traits::RangeFieldType;
587 using size_type =
typename DGSpace::Traits::SizeType;
590 using namespace TypeTree::Indices;
591 const auto& dgspace_s = child(lfsv_s,_0);
594 const auto& cell_inside = ig.inside();
597 auto geo = ig.geometry();
598 auto geo_inside = cell_inside.geometry();
601 auto geo_in_inside = ig.geometryInInside();
605 auto local_inside = ref_el_inside.position(0,0);
606 auto mu_s = param.mu(cell_inside,local_inside);
607 auto eps_s = param.eps(cell_inside,local_inside);
611 const auto& n_F = ig.centerUnitOuterNormal();
614 Dune::FieldMatrix<DF,dim*2,dim*2> R_s;
616 Dune::FieldMatrix<DF,dim*2,dim*2> Dplus_s(0.0);
617 Dplus_s[0][0] = 1.0/sqrt(eps_s*mu_s);
618 Dplus_s[1][1] = 1.0/sqrt(eps_s*mu_s);
619 Dune::FieldMatrix<DF,dim*2,dim*2> Aplus_s(R_s);
620 Aplus_s.rightmultiply(Dplus_s);
622 Aplus_s.rightmultiply(R_s);
625 Dune::FieldMatrix<DF,dim*2,dim*2> R_n;
627 Dune::FieldMatrix<DF,dim*2,dim*2> Dminus_n(0.0);
628 Dminus_n[2][2] = -1.0/sqrt(eps_s*mu_s);
629 Dminus_n[3][3] = -1.0/sqrt(eps_s*mu_s);
630 Dune::FieldMatrix<DF,dim*2,dim*2> Aminus_n(R_n);
631 Aminus_n.rightmultiply(Dminus_n);
633 Aminus_n.rightmultiply(R_n);
636 Dune::FieldVector<RF,dim*2> u_s(0.0);
637 Dune::FieldVector<RF,dim*2> u_n;
638 Dune::FieldVector<RF,dim*2> f(0.0);
643 const int order_s = dgspace_s.finiteElement().localBasis().order();
644 const int intorder = overintegration+1+2*order_s;
648 const auto& iplocal_s = geo_in_inside.global(qp.position());
651 const auto& phi_s = cache[order_s].evaluateFunction
652 (iplocal_s,dgspace_s.finiteElement().localBasis());
655 for (size_type i=0; i<
dim*2; i++){
657 for (size_type k=0; k<dgspace_s.size(); k++)
658 u_s[i] += x_s(lfsv_s.child(i),k)*phi_s[k];
663 u_n = (param.g(ig.intersection(),qp.position(),u_s));
674 auto factor = qp.weight() * geo.integrationElement(qp.position());
675 for (size_type k=0; k<dgspace_s.size(); k++)
676 for (size_type i=0; i<dim*2; i++)
677 r_s.accumulate(lfsv_s.child(i),k,f[i]*phi_s[k]*factor);
686 template<
typename EG,
typename LFSV,
typename R>
690 using namespace TypeTree::Indices;
691 using DGSpace = TypeTree::Child<LFSV,_0>;
692 using size_type =
typename DGSpace::Traits::SizeType;
695 using namespace TypeTree::Indices;
696 const auto& dgspace = child(lfsv,_0);
699 const auto& cell = eg.entity();
702 auto geo = eg.geometry();
705 const int order_s = dgspace.finiteElement().localBasis().order();
706 const int intorder = overintegration+2*order_s;
710 auto j = param.j(cell,qp.position());
713 const auto& phi = cache[order_s].evaluateFunction(qp.position(),dgspace.finiteElement().localBasis());
716 auto factor = qp.weight() * geo.integrationElement(qp.position());
717 for (size_type k=0; k<
dim*2; k++)
718 for (size_type i=0; i<dgspace.size(); i++)
719 r.accumulate(lfsv.child(k),i,-j[k]*phi[i]*factor);
724 void setTime (
typename T::Traits::RangeFieldType t)
729 void preStep (
typename T::Traits::RangeFieldType time,
typename T::Traits::RangeFieldType dt,
735 void preStage (
typename T::Traits::RangeFieldType time,
int r)
745 typename T::Traits::RangeFieldType
suggestTimestep (
typename T::Traits::RangeFieldType dt)
const 753 typedef typename FEM::Traits::FiniteElementType::Traits::LocalBasisType LocalBasisType;
755 std::vector<Cache> cache;
771 template<
typename T,
typename FEM>
777 enum {
dim = T::Traits::GridViewType::dimension };
780 enum { doPatternVolume =
true };
783 enum { doAlphaVolume =
true };
786 : param(param_), overintegration(overintegration_), cache(20)
790 template<
typename LFSU,
typename LFSV,
typename LocalPattern>
792 LocalPattern& pattern)
const 795 static_assert(LFSU::CHILDREN==LFSV::CHILDREN,
"need U=V!");
796 static_assert(LFSV::CHILDREN==
dim*2,
"need exactly dim*2 components!");
798 for (
size_t k=0; k<LFSV::CHILDREN; k++)
799 for (
size_t i=0; i<lfsv.child(k).size(); ++i)
800 for (
size_t j=0; j<lfsu.child(k).size(); ++j)
801 pattern.addLink(lfsv.child(k),i,lfsu.child(k),j);
805 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
806 void alpha_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const 809 using namespace TypeTree::Indices;
810 using DGSpace = TypeTree::Child<LFSV,_0>;
811 using RF =
typename DGSpace::Traits::FiniteElementType::
812 Traits::LocalBasisType::Traits::RangeFieldType;
813 using size_type =
typename DGSpace::Traits::SizeType;
816 using namespace TypeTree::Indices;
817 const auto& dgspace = child(lfsv,_0);
820 auto geo = eg.geometry();
823 Dune::FieldVector<RF,dim*2> u(0.0);
826 const int order = dgspace.finiteElement().localBasis().order();
827 const int intorder = overintegration+2*order;
831 const auto& phi = cache[order].evaluateFunction(qp.position(),dgspace.finiteElement().localBasis());
834 for (size_type k=0; k<
dim*2; k++){
836 for (size_type j=0; j<dgspace.size(); j++)
837 u[k] += x(lfsv.child(k),j)*phi[j];
841 auto factor = qp.weight() * geo.integrationElement(qp.position());
842 for (size_type k=0; k<dim*2; k++)
843 for (size_type i=0; i<dgspace.size(); i++)
844 r.accumulate(lfsv.child(k),i,u[k]*phi[i]*factor);
849 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
854 using namespace TypeTree::Indices;
855 using DGSpace = TypeTree::Child<LFSV,_0>;
856 using size_type =
typename DGSpace::Traits::SizeType;
859 using namespace TypeTree::Indices;
860 const auto& dgspace = child(lfsv,_0);
863 auto geo = eg.geometry();
866 const int order = dgspace.finiteElement().localBasis().order();
867 const int intorder = overintegration+2*order;
871 const auto& phi = cache[order].evaluateFunction(qp.position(),dgspace.finiteElement().localBasis());
874 auto factor = qp.weight() * geo.integrationElement(qp.position());
875 for (size_type k=0; k<
dim*2; k++)
876 for (size_type i=0; i<dgspace.size(); i++)
877 for (size_type j=0; j<dgspace.size(); j++)
878 mat.accumulate(lfsv.child(k),i,lfsu.child(k),j,phi[j]*phi[i]*factor);
885 typedef typename FEM::Traits::FiniteElementType::Traits::LocalBasisType LocalBasisType;
887 std::vector<Cache> cache;
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:729
Definition: maxwelldg.hh:772
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 jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: maxwelldg.hh:850
static void eigenvalues(T1 eps, T1 mu, const Dune::FieldVector< T2, 2 *dim > &e)
Definition: maxwelldg.hh:53
sparsity pattern generator
Definition: pattern.hh:13
Implements linear and nonlinear versions of jacobian_apply_skeleton() based on alpha_skeleton() ...
Definition: numericaljacobianapply.hh:180
const IG & ig
Definition: constraints.hh:148
static const int dim
Definition: adaptivity.hh:83
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:576
Implement jacobian_boundary() based on alpha_boundary()
Definition: numericaljacobian.hh:250
ReferenceElementWrapper< ReferenceElement< typename Geometry::ctype, Geometry::mydimension > > referenceElement(const Geometry &geo)
Returns the reference element for the given geometry.
Definition: referenceelements.hh:144
Definition: adaptivity.hh:27
sparsity pattern generator
Definition: pattern.hh:29
T::Traits::RangeFieldType suggestTimestep(typename T::Traits::RangeFieldType dt) const
to be called once before each stage
Definition: maxwelldg.hh:745
Implements linear and nonlinear versions of jacobian_apply_boundary() based on alpha_boundary() ...
Definition: numericaljacobianapply.hh:285
Default flags for all local operators.
Definition: flags.hh:18
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: maxwelldg.hh:806
const Entity & e
Definition: localfunctionspace.hh:111
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
void preStage(typename T::Traits::RangeFieldType time, int r)
to be called once before each stage
Definition: maxwelldg.hh:735
Definition: maxwelldg.hh:29
Implement jacobian_volume() based on alpha_volume()
Definition: numericaljacobian.hh:31
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
void setTime(typename T::Traits::RangeFieldType t)
set time in parameter class
Definition: maxwelldg.hh:724
DGMaxwellTemporalOperator(T ¶m_, int overintegration_=0)
Definition: maxwelldg.hh:785
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:438
void pattern_volume(const LFSU &lfsu, const LFSV &lfsv, LocalPattern &pattern) const
Definition: maxwelldg.hh:791
const std::string s
Definition: function.hh:1102
store values of basis functions and gradients in a cache
Definition: localbasiscache.hh:17
Implements linear and nonlinear versions of jacobian_apply_volume() based on alpha_volume() ...
Definition: numericaljacobianapply.hh:32
DGMaxwellSpatialOperator(T ¶m_, int overintegration_=0)
Definition: maxwelldg.hh:330
void postStage()
to be called once at the end of each stage
Definition: maxwelldg.hh:740
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: maxwelldg.hh:337
Definition: maxwelldg.hh:304
Implement jacobian_skeleton() based on alpha_skeleton()
Definition: numericaljacobian.hh:156
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: maxwelldg.hh:687