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>
50 template<
typename T1,
typename T2,
typename T3>
51 static void eigenvalues (T1 eps, T1 mu,
const Dune::FieldVector<T2,2*dim>&
e)
53 T1
s = 1.0/sqrt(mu*eps);
73 template<
typename T1,
typename T2,
typename T3>
74 static void eigenvectors (T1 eps, T1 mu,
const Dune::FieldVector<T2,dim>& n, Dune::FieldMatrix<T3,2*dim,2*dim>& R)
76 T1 a=n[0], b=n[1], c=n[2];
78 Dune::FieldVector<T2,dim> re, im;
81 re[0]=a*c; re[1]=b*c; re[2]=c*c-1;
82 im[0]=-b; im[1]=a; im[2]=0;
86 re[0]=a*b; re[1]=b*b-1; re[2]=b*c;
87 im[0]=c; im[1]=0.0; im[2]=-a;
91 R[0][0] = re[0]; R[0][1] = -im[0];
92 R[1][0] = re[1]; R[1][1] = -im[1];
93 R[2][0] = re[2]; R[2][1] = -im[2];
94 R[3][0] = im[0]; R[3][1] = re[0];
95 R[4][0] = im[1]; R[4][1] = re[1];
96 R[5][0] = im[2]; R[5][1] = re[2];
99 R[0][2] = im[0]; R[0][3] = re[0];
100 R[1][2] = im[1]; R[1][3] = re[1];
101 R[2][2] = im[2]; R[2][3] = re[2];
102 R[3][2] = re[0]; R[3][3] = -im[0];
103 R[4][2] = re[1]; R[4][3] = -im[1];
104 R[5][2] = re[2]; R[5][3] = -im[2];
107 R[0][4] = a; R[0][5] = 0;
108 R[1][4] = b; R[1][5] = 0;
109 R[2][4] = c; R[2][5] = 0;
110 R[3][4] = 0; R[3][5] = a;
111 R[4][4] = 0; R[4][5] = b;
112 R[5][4] = 0; R[5][5] = c;
117 for (std::size_t i=0; i<3; i++)
118 for (std::size_t j=0; j<6; j++)
120 for (std::size_t i=3; i<6; i++)
121 for (std::size_t j=0; j<6; j++)
301 template<
typename T,
typename FEM>
314 enum { dim = T::Traits::GridViewType::dimension };
329 : param(param_), overintegration(overintegration_), cache(20)
334 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
335 void alpha_volume (
const EG&
eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
338 using DGSpace =
typename LFSV::template Child<0>::Type;
339 using DF =
typename DGSpace::Traits::FiniteElementType::Traits
340 ::LocalBasisType::Traits::DomainFieldType;
341 using RF =
typename DGSpace::Traits::FiniteElementType::Traits
342 ::LocalBasisType::Traits::RangeFieldType;
343 using size_type =
typename DGSpace::Traits::SizeType;
346 static_assert(LFSV::CHILDREN == dim*2,
347 "need exactly dim*2 components!");
350 const auto& dgspace = lfsv.template child<0>();
353 const int order = dgspace.finiteElement().localBasis().order();
354 const int intorder = overintegration+2*order;
355 auto gt = eg.geometry().type();
358 std::vector<Dune::FieldVector<RF,dim> > gradphi(dgspace.size());
361 const auto& localcenter =
362 Dune::ReferenceElements<DF,dim>::general(gt).position(0,0);
363 auto mu = param.mu(eg.entity(),localcenter);
364 auto eps = param.eps(eg.entity(),localcenter);
365 auto sigma = param.sigma(eg.entity(),localcenter);
367 auto epsinv = 1.0/eps;
372 for (
const auto &qp : Dune::QuadratureRules<DF,dim>::rule(gt,intorder))
375 const auto& phi = cache[order].evaluateFunction
376 (qp.position(), dgspace.finiteElement().localBasis());
379 Dune::FieldVector<RF,dim*2> u(0.0);
380 for (size_type k=0; k<dim*2; k++)
381 for (size_type j=0; j<dgspace.size(); j++)
382 u[k] += x(lfsv.child(k),j)*phi[j];
386 const auto& js = cache[order].evaluateJacobian
387 (qp.position(), dgspace.finiteElement().localBasis());
391 eg.geometry().jacobianInverseTransposed(qp.position());
392 for (size_type i=0; i<dgspace.size(); i++)
393 jac.mv(js[i][0],gradphi[i]);
396 RF factor = qp.weight()
397 * eg.geometry().integrationElement(qp.position());
399 Dune::FieldMatrix<RF,dim*2,dim> F;
400 F[0][0] = 0; F[0][1] = -muinv*u[5]; F[0][2] = muinv*u[4];
401 F[1][0] = muinv*u[5]; F[1][1] = 0; F[1][2] = -muinv*u[3];
402 F[2][0] =-muinv*u[4]; F[2][1] = muinv*u[3]; F[2][2] = 0;
403 F[3][0] = 0; F[3][1] = epsinv*u[2]; F[3][2] = -epsinv*u[1];
404 F[4][0] = -epsinv*u[2]; F[4][1] = 0; F[4][2] = epsinv*u[0];
405 F[5][0] = epsinv*u[1]; F[5][1] = -epsinv*u[0]; F[5][2] = 0;
408 for (size_type i=0; i<dim*2; i++)
410 for (size_type k=0; k<dgspace.size(); k++)
412 for (size_type j=0; j<dim; j++)
413 r.accumulate(lfsv.child(i),k,-F[i][j]*gradphi[k][j]*factor);
416 for (size_type i=0; i<dim; i++)
418 for (size_type k=0; k<dgspace.size(); k++)
419 r.accumulate(lfsv.child(i),k,(sigma/eps)*u[i]*phi[k]*factor);
429 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
431 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
432 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
433 R& r_s, R& r_n)
const
439 using DGSpace =
typename LFSV::template Child<0>::Type;
440 using DF =
typename DGSpace::Traits::FiniteElementType::
441 Traits::LocalBasisType::Traits::DomainFieldType;
442 using RF =
typename DGSpace::Traits::FiniteElementType::
443 Traits::LocalBasisType::Traits::RangeFieldType;
444 using size_type =
typename DGSpace::Traits::SizeType;
447 const auto& dgspace_s = lfsv_s.template child<0>();
448 const auto& dgspace_n = lfsv_n.template child<0>();
451 const auto& n_F = ig.centerUnitOuterNormal();
454 const auto& inside_local =
455 Dune::ReferenceElements<DF,dim>::general(ig.inside().type()).position(0,0);
456 const auto& outside_local =
457 Dune::ReferenceElements<DF,dim>::general(ig.outside().type()).position(0,0);
462 auto mu_s = param.mu(ig.inside(),inside_local);
463 auto mu_n = param.mu(ig.outside(),outside_local);
464 auto eps_s = param.eps(ig.inside(),inside_local);
465 auto eps_n = param.eps(ig.outside(),outside_local);
470 Dune::FieldMatrix<DF,dim*2,dim*2> R_s;
472 Dune::FieldMatrix<DF,dim*2,dim*2> Dplus_s(0.0);
473 Dplus_s[0][0] = 1.0/sqrt(eps_s*mu_s);
474 Dplus_s[1][1] = 1.0/sqrt(eps_s*mu_s);
475 Dune::FieldMatrix<DF,dim*2,dim*2> Aplus_s(R_s);
476 Aplus_s.rightmultiply(Dplus_s);
478 Aplus_s.rightmultiply(R_s);
481 Dune::FieldMatrix<DF,dim*2,dim*2> R_n;
483 Dune::FieldMatrix<DF,dim*2,dim*2> Dminus_n(0.0);
484 Dminus_n[2][2] = -1.0/sqrt(eps_n*mu_n);
485 Dminus_n[3][3] = -1.0/sqrt(eps_n*mu_n);
486 Dune::FieldMatrix<DF,dim*2,dim*2> Aminus_n(R_n);
487 Aminus_n.rightmultiply(Dminus_n);
489 Aminus_n.rightmultiply(R_n);
492 const int order_s = dgspace_s.finiteElement().localBasis().order();
493 const int order_n = dgspace_n.finiteElement().localBasis().order();
494 const int intorder = overintegration+1+2*max(order_s,order_n);
495 auto gtface = ig.geometry().type();
500 for (
const auto &qp :
501 Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder))
504 const auto& iplocal_s = ig.geometryInInside().global(qp.position());
505 const auto& iplocal_n = ig.geometryInOutside().global(qp.position());
508 const auto& phi_s = cache[order_s].evaluateFunction(iplocal_s,dgspace_s.finiteElement().localBasis());
509 const auto& phi_n = cache[order_n].evaluateFunction(iplocal_n,dgspace_n.finiteElement().localBasis());
512 Dune::FieldVector<RF,dim*2> u_s(0.0);
513 for (size_type i=0; i<dim*2; i++)
514 for (size_type k=0; k<dgspace_s.size(); k++)
515 u_s[i] += x_s(lfsv_s.child(i),k)*phi_s[k];
516 Dune::FieldVector<RF,dim*2> u_n(0.0);
517 for (size_type i=0; i<dim*2; i++)
518 for (size_type k=0; k<dgspace_n.size(); k++)
519 u_n[i] += x_n(lfsv_n.child(i),k)*phi_n[k];
522 Dune::FieldVector<RF,dim*2> f(0.0);
531 RF factor = qp.weight() * ig.geometry().integrationElement(qp.position());
532 for (size_type k=0; k<dgspace_s.size(); k++)
533 for (size_type i=0; i<dim*2; i++)
534 r_s.accumulate(lfsv_s.child(i),k,f[i]*phi_s[k]*factor);
535 for (size_type k=0; k<dgspace_n.size(); k++)
536 for (size_type i=0; i<dim*2; i++)
537 r_n.accumulate(lfsv_n.child(i),k,-f[i]*phi_n[k]*factor);
549 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
551 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
555 using DGSpace =
typename LFSV::template Child<0>::Type;
556 using DF =
typename DGSpace::Traits::FiniteElementType::
557 Traits::LocalBasisType::Traits::DomainFieldType;
558 using RF =
typename DGSpace::Traits::FiniteElementType::
559 Traits::LocalBasisType::Traits::RangeFieldType;
560 using size_type =
typename DGSpace::Traits::SizeType;
563 const auto& dgspace_s = lfsv_s.template child<0>();
566 const auto& n_F = ig.centerUnitOuterNormal();
569 const auto& inside_local =
570 Dune::ReferenceElements<DF,dim>::general(ig.inside().type()).position(0,0);
571 auto mu_s = param.mu(ig.inside(),inside_local);
572 auto eps_s = param.eps(ig.inside(),inside_local);
576 Dune::FieldMatrix<DF,dim*2,dim*2> R_s;
578 Dune::FieldMatrix<DF,dim*2,dim*2> Dplus_s(0.0);
579 Dplus_s[0][0] = 1.0/sqrt(eps_s*mu_s);
580 Dplus_s[1][1] = 1.0/sqrt(eps_s*mu_s);
581 Dune::FieldMatrix<DF,dim*2,dim*2> Aplus_s(R_s);
582 Aplus_s.rightmultiply(Dplus_s);
584 Aplus_s.rightmultiply(R_s);
587 Dune::FieldMatrix<DF,dim*2,dim*2> R_n;
589 Dune::FieldMatrix<DF,dim*2,dim*2> Dminus_n(0.0);
590 Dminus_n[2][2] = -1.0/sqrt(eps_s*mu_s);
591 Dminus_n[3][3] = -1.0/sqrt(eps_s*mu_s);
592 Dune::FieldMatrix<DF,dim*2,dim*2> Aminus_n(R_n);
593 Aminus_n.rightmultiply(Dminus_n);
595 Aminus_n.rightmultiply(R_n);
598 const int order_s = dgspace_s.finiteElement().localBasis().order();
599 const int intorder = overintegration+1+2*order_s;
600 auto gtface = ig.geometry().type();
606 Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder))
609 const auto& iplocal_s = ig.geometryInInside().global(qp.position());
612 const auto& phi_s = cache[order_s].evaluateFunction
613 (iplocal_s,dgspace_s.finiteElement().localBasis());
616 Dune::FieldVector<RF,dim*2> u_s(0.0);
617 for (size_type i=0; i<dim*2; i++)
618 for (size_type k=0; k<dgspace_s.size(); k++)
619 u_s[i] += x_s(lfsv_s.child(i),k)*phi_s[k];
623 Dune::FieldVector<RF,dim*2> u_n(param.g(ig.intersection(),qp.position(),u_s));
627 Dune::FieldVector<RF,dim*2> f(0.0);
634 RF factor = qp.weight() * ig.geometry().integrationElement(qp.position());
635 for (size_type k=0; k<dgspace_s.size(); k++)
636 for (size_type i=0; i<dim*2; i++)
637 r_s.accumulate(lfsv_s.child(i),k,f[i]*phi_s[k]*factor);
646 template<
typename EG,
typename LFSV,
typename R>
650 using DGSpace =
typename LFSV::template Child<0>::Type;
651 using DF =
typename DGSpace::Traits::FiniteElementType::
652 Traits::LocalBasisType::Traits::DomainFieldType;
653 using size_type =
typename DGSpace::Traits::SizeType;
656 const auto& dgspace = lfsv.template child<0>();
659 const int order_s = dgspace.finiteElement().localBasis().order();
660 const int intorder = overintegration+2*order_s;
661 Dune::GeometryType gt = eg.geometry().type();
664 for (
const auto &qp : Dune::QuadratureRules<DF,dim>::rule(gt,intorder))
667 auto j = param.j(eg.entity(),qp.position());
670 const auto& phi = cache[order_s].evaluateFunction(qp.position(),dgspace.finiteElement().localBasis());
673 auto factor = qp.weight() * eg.geometry().integrationElement(qp.position());
674 for (size_type k=0; k<dim*2; k++)
675 for (size_type i=0; i<dgspace.size(); i++)
676 r.accumulate(lfsv.child(k),i,-j[k]*phi[i]*factor);
681 void setTime (
typename T::Traits::RangeFieldType t)
686 void preStep (
typename T::Traits::RangeFieldType time,
typename T::Traits::RangeFieldType dt,
692 void preStage (
typename T::Traits::RangeFieldType time,
int r)
702 typename T::Traits::RangeFieldType
suggestTimestep (
typename T::Traits::RangeFieldType dt)
const
710 typedef typename FEM::Traits::FiniteElementType::Traits::LocalBasisType LocalBasisType;
712 std::vector<Cache> cache;
728 template<
typename T,
typename FEM>
734 enum { dim = T::Traits::GridViewType::dimension };
743 : param(param_), overintegration(overintegration_), cache(20)
747 template<
typename LFSU,
typename LFSV,
typename LocalPattern>
749 LocalPattern& pattern)
const
752 static_assert(LFSU::CHILDREN==LFSV::CHILDREN,
"need U=V!");
753 static_assert(LFSV::CHILDREN==dim*2,
"need exactly dim*2 components!");
755 for (
size_t k=0; k<LFSV::CHILDREN; k++)
756 for (
size_t i=0; i<lfsv.child(k).size(); ++i)
757 for (
size_t j=0; j<lfsu.child(k).size(); ++j)
758 pattern.addLink(lfsv.child(k),i,lfsu.child(k),j);
762 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
763 void alpha_volume (
const EG&
eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
766 using DGSpace =
typename LFSV::template Child<0>::Type;
767 using DF =
typename DGSpace::Traits::FiniteElementType::
768 Traits::LocalBasisType::Traits::DomainFieldType;
769 using RF =
typename DGSpace::Traits::FiniteElementType::
770 Traits::LocalBasisType::Traits::RangeFieldType;
771 using size_type =
typename DGSpace::Traits::SizeType;
774 const auto& dgspace = lfsv.template child<0>();
777 const int order = dgspace.finiteElement().localBasis().order();
778 const int intorder = overintegration+2*order;
779 auto gt = eg.geometry().type();
782 for (
const auto& qp : Dune::QuadratureRules<DF,dim>::rule(gt,intorder))
785 const auto& phi = cache[order].evaluateFunction(qp.position(),dgspace.finiteElement().localBasis());
788 Dune::FieldVector<RF,dim*2> u(0.0);
789 for (size_type k=0; k<dim*2; k++)
790 for (size_type j=0; j<dgspace.size(); j++)
791 u[k] += x(lfsv.child(k),j)*phi[j];
794 RF factor = qp.weight() * eg.geometry().integrationElement(qp.position());
795 for (size_type k=0; k<dim*2; k++)
796 for (size_type i=0; i<dgspace.size(); i++)
797 r.accumulate(lfsv.child(k),i,u[k]*phi[i]*factor);
802 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
807 using DGSpace =
typename LFSV::template Child<0>::Type;
808 using DF =
typename DGSpace::Traits::FiniteElementType::
809 Traits::LocalBasisType::Traits::DomainFieldType;
810 using RF =
typename DGSpace::Traits::FiniteElementType::
811 Traits::LocalBasisType::Traits::RangeFieldType;
812 using size_type =
typename DGSpace::Traits::SizeType;
815 const auto& dgspace = lfsv.template child<0>();
818 const int order = dgspace.finiteElement().localBasis().order();
819 const int intorder = overintegration+2*order;
820 auto gt = eg.geometry().type();
823 for (
const auto& qp : Dune::QuadratureRules<DF,dim>::rule(gt,intorder))
826 const auto& phi = cache[order].evaluateFunction(qp.position(),dgspace.finiteElement().localBasis());
829 RF factor = qp.weight() * eg.geometry().integrationElement(qp.position());
830 for (size_type k=0; k<dim*2; k++)
831 for (size_type i=0; i<dgspace.size(); i++)
832 for (size_type j=0; j<dgspace.size(); j++)
833 mat.accumulate(lfsv.child(k),i,lfsu.child(k),j,phi[j]*phi[i]*factor);
840 typedef typename FEM::Traits::FiniteElementType::Traits::LocalBasisType LocalBasisType;
842 std::vector<Cache> cache;
static void eigenvalues(T1 eps, T1 mu, const Dune::FieldVector< T2, 2 *dim > &e)
Definition: maxwelldg.hh:51
Definition: maxwelldg.hh:323
Definition: maxwelldg.hh:324
Definition: maxwelldg.hh:737
const E & e
Definition: interpolate.hh:172
Implement jacobian_apply_volume() based on alpha_volume()
Definition: defaultimp.hh:321
Definition: maxwelldg.hh:729
void preStage(typename T::Traits::RangeFieldType time, int r)
to be called once before each stage
Definition: maxwelldg.hh:692
T::Traits::RangeFieldType suggestTimestep(typename T::Traits::RangeFieldType dt) const
to be called once before each stage
Definition: maxwelldg.hh:702
const IG & ig
Definition: constraints.hh:147
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: maxwelldg.hh:763
static void eigenvectors(T1 eps, T1 mu, const Dune::FieldVector< T2, dim > &n, Dune::FieldMatrix< T3, 2 *dim, 2 *dim > &R)
Definition: maxwelldg.hh:74
Definition: maxwelldg.hh:27
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: maxwelldg.hh:335
Default flags for all local operators.
Definition: flags.hh:18
Definition: maxwelldg.hh:322
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:550
sparsity pattern generator
Definition: pattern.hh:13
DGMaxwellSpatialOperator(T ¶m_, int overintegration_=0)
Definition: maxwelldg.hh:328
Definition: maxwelldg.hh:302
Implement jacobian_apply_skeleton() based on alpha_skeleton()
Definition: defaultimp.hh:445
Definition: maxwelldg.hh:325
Definition: adaptivity.hh:27
Implement jacobian_volume() based on alpha_volume()
Definition: defaultimp.hh:32
Definition: maxwelldg.hh:318
Implement jacobian_boundary() based on alpha_boundary()
Definition: defaultimp.hh:251
static const int dim
Definition: adaptivity.hh:83
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:430
DGMaxwellTemporalOperator(T ¶m_, int overintegration_=0)
Definition: maxwelldg.hh:742
const std::string s
Definition: function.hh:1102
sparsity pattern generator
Definition: pattern.hh:29
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:686
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: maxwelldg.hh:803
void pattern_volume(const LFSU &lfsu, const LFSV &lfsv, LocalPattern &pattern) const
Definition: maxwelldg.hh:748
void postStage()
to be called once at the end of each stage
Definition: maxwelldg.hh:697
store values of basis functions and gradients in a cache
Definition: localbasiscache.hh:15
Implement jacobian_skeleton() based on alpha_skeleton()
Definition: defaultimp.hh:157
Definition: maxwelldg.hh:319
void setTime(typename T::Traits::RangeFieldType t)
set time in parameter class
Definition: maxwelldg.hh:681
const EG & eg
Definition: constraints.hh:280
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: maxwelldg.hh:647
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
Definition: maxwelldg.hh:740
Implement jacobian_apply_boundary() based on alpha_boundary()
Definition: defaultimp.hh:538