2 #ifndef DUNE_PDELAB_LOCALOPERATOR_LINEARACOUSTICSDG_HH
3 #define DUNE_PDELAB_LOCALOPERATOR_LINEARACOUSTICSDG_HH
7 #include<dune/common/exceptions.hh>
8 #include<dune/common/fvector.hh>
9 #include<dune/geometry/referenceelements.hh>
43 template<
typename T1,
typename T2,
typename T3>
46 RT[0][0] = 1; RT[0][1] = c*n[0];
47 RT[1][0] = -1; RT[1][1] = c*n[0];
50 template<
typename T1,
typename T2,
typename T3>
51 static void eigenvectors (T1 c,
const Dune::FieldVector<T2,dim>& n, Dune::FieldMatrix<T3,dim+1,dim+1>& RT)
53 RT[0][0] = 1; RT[1][0] = c*n[0];
54 RT[0][1] = -1; RT[1][1] = c*n[0];
72 template<
typename T1,
typename T2,
typename T3>
75 RT[0][0] = 0; RT[0][1] = -n[1]; RT[0][2] = n[0];
76 RT[1][0] = 1; RT[1][1] = c*n[0]; RT[1][2] = c*n[1];
77 RT[2][0] = -1; RT[2][1] = c*n[0]; RT[2][2] = c*n[1];
80 template<
typename T1,
typename T2,
typename T3>
81 static void eigenvectors (T1 c,
const Dune::FieldVector<T2,dim>& n, Dune::FieldMatrix<T3,dim+1,dim+1>& RT)
83 RT[0][0] = 0; RT[1][0] = -n[1]; RT[2][0] = n[0];
84 RT[0][1] = 1; RT[1][1] = c*n[0]; RT[2][1] = c*n[1];
85 RT[0][2] = -1; RT[1][2] = c*n[0]; RT[2][2] = c*n[1];
102 template<
typename T1,
typename T2,
typename T3>
105 Dune::FieldVector<T2,dim>
s;
106 s[0] = n[1]-n[2];
s[1] = n[2]-n[0];
s[2] = n[0]-n[1];
107 if (
s.two_norm()<1
e-5)
109 s[0] = n[1]+n[2];
s[1] = n[2]-n[0];
s[2] = -(n[0]+n[1]);
112 Dune::FieldVector<T2,dim> t;
113 t[0] = n[1]*
s[2] - n[2]*
s[1];
114 t[1] = n[2]*
s[0] - n[0]*
s[2];
115 t[2] = n[0]*
s[1] - n[1]*
s[0];
117 RT[0][0] = 0; RT[0][1] =
s[0]; RT[0][2] =
s[1]; RT[0][3] =
s[2];
118 RT[1][0] = 0; RT[1][1] = t[0]; RT[1][2] = t[1]; RT[1][3] = t[2];
119 RT[2][0] = 1; RT[2][1] = c*n[0]; RT[2][2] = c*n[1]; RT[2][3] = c*n[2];
120 RT[3][0] = -1; RT[3][1] = c*n[0]; RT[3][2] = c*n[1]; RT[3][3] = c*n[2];
123 template<
typename T1,
typename T2,
typename T3>
124 static void eigenvectors (T1 c,
const Dune::FieldVector<T2,dim>& n, Dune::FieldMatrix<T3,dim+1,dim+1>& RT)
126 Dune::FieldVector<T2,dim>
s;
127 s[0] = n[1]-n[2];
s[1] = n[2]-n[0];
s[2] = n[0]-n[1];
128 if (
s.two_norm()<1
e-5)
130 s[0] = n[1]+n[2];
s[1] = n[2]-n[0];
s[2] = -(n[0]+n[1]);
133 Dune::FieldVector<T2,dim> t;
134 t[0] = n[1]*
s[2] - n[2]*
s[1];
135 t[1] = n[2]*
s[0] - n[0]*
s[2];
136 t[2] = n[0]*
s[1] - n[1]*
s[0];
138 RT[0][0] = 0; RT[1][0] =
s[0]; RT[2][0] =
s[1]; RT[3][0] =
s[2];
139 RT[0][1] = 0; RT[1][1] = t[0]; RT[2][1] = t[1]; RT[3][1] = t[2];
140 RT[0][2] = 1; RT[1][2] = c*n[0]; RT[2][2] = c*n[1]; RT[3][2] = c*n[2];
141 RT[0][3] = -1; RT[1][3] = c*n[0]; RT[2][3] = c*n[1]; RT[3][3] = c*n[2];
161 template<
typename T,
typename FEM>
174 enum { dim = T::Traits::GridViewType::dimension };
189 : param(param_), overintegration(overintegration_), cache(20)
194 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
195 void alpha_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
198 using namespace Indices;
199 using DGSpace = TypeTree::Child<LFSV,_0>;
200 using RF =
typename DGSpace::Traits::FiniteElementType::
201 Traits::LocalBasisType::Traits::RangeFieldType;
202 using size_type =
typename DGSpace::Traits::SizeType;
205 if (TypeTree::degree(lfsv)!=dim+1)
206 DUNE_THROW(Dune::Exception,
"need exactly dim+1 components!");
209 const auto& dgspace = child(lfsv,_0);
212 const auto& cell = eg.entity();
215 auto geo = eg.geometry();
218 auto ref_el = referenceElement(geo);
219 auto localcenter = ref_el.position(0,0);
220 auto c2 = param.c(cell,localcenter);
224 typename EG::Geometry::JacobianInverseTransposed jac;
227 Dune::FieldVector<RF,dim+1> u(0.0);
228 std::vector<Dune::FieldVector<RF,dim> > gradphi(dgspace.size());
231 const int order = dgspace.finiteElement().localBasis().order();
232 const int intorder = overintegration+2*order;
236 auto& phi = cache[order].evaluateFunction(ip.position(),dgspace.finiteElement().localBasis());
240 for (size_type k=0; k<=dim; k++)
241 for (size_type j=0; j<dgspace.size(); j++)
242 u[k] += x(lfsv.child(k),j)*phi[j];
246 auto& js = cache[order].evaluateJacobian(ip.position(),dgspace.finiteElement().localBasis());
249 jac = geo.jacobianInverseTransposed(ip.position());
250 for (size_type i=0; i<dgspace.size(); i++)
251 jac.mv(js[i][0],gradphi[i]);
254 auto factor = ip.weight() * geo.integrationElement(ip.position());
255 for (size_type k=0; k<dgspace.size(); k++)
258 for (size_type j=0; j<dim; j++)
259 r.accumulate(lfsv.child(0),k, - u[j+1]*gradphi[k][j]*factor);
261 for (size_type i=1; i<=dim; i++)
262 r.accumulate(lfsv.child(i),k, - c2*u[0]*gradphi[k][i-1]*factor);
272 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
274 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
275 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
276 R& r_s, R& r_n)
const
279 using namespace Indices;
280 using DGSpace = TypeTree::Child<LFSV,_0>;
281 using DF =
typename DGSpace::Traits::FiniteElementType::
282 Traits::LocalBasisType::Traits::DomainFieldType;
283 using RF =
typename DGSpace::Traits::FiniteElementType::
284 Traits::LocalBasisType::Traits::RangeFieldType;
285 using size_type =
typename DGSpace::Traits::SizeType;
288 const auto& dgspace_s = child(lfsv_s,_0);
289 const auto& dgspace_n = child(lfsv_n,_0);
292 const Dune::FieldVector<DF,dim> n_F =
ig.centerUnitOuterNormal();
295 const auto& cell_inside =
ig.inside();
296 const auto& cell_outside =
ig.outside();
299 auto geo =
ig.geometry();
300 auto geo_inside = cell_inside.geometry();
301 auto geo_outside = cell_outside.geometry();
304 auto geo_in_inside =
ig.geometryInInside();
305 auto geo_in_outside =
ig.geometryInOutside();
308 auto ref_el_inside = referenceElement(geo_inside);
309 auto ref_el_outside = referenceElement(geo_outside);
310 auto local_inside = ref_el_inside.position(0,0);
311 auto local_outside = ref_el_outside.position(0,0);
312 auto c_s = param.c(cell_inside,local_inside);
313 auto c_n = param.c(cell_outside,local_outside);
316 Dune::FieldMatrix<DF,dim+1,dim+1> RT;
318 Dune::FieldVector<DF,dim+1> alpha;
319 for (
int i=0; i<=dim; i++) alpha[i] = RT[dim-1][i];
320 Dune::FieldVector<DF,dim+1> unit(0.0);
322 Dune::FieldVector<DF,dim+1> beta;
324 Dune::FieldMatrix<DF,dim+1,dim+1> A_plus_s;
325 for (
int i=0; i<=dim; i++)
326 for (
int j=0; j<=dim; j++)
327 A_plus_s[i][j] = c_s*alpha[i]*beta[j];
331 for (
int i=0; i<=dim; i++) alpha[i] = RT[dim][i];
335 Dune::FieldMatrix<DF,dim+1,dim+1> A_minus_n;
336 for (
int i=0; i<=dim; i++)
337 for (
int j=0; j<=dim; j++)
338 A_minus_n[i][j] = -c_n*alpha[i]*beta[j];
341 Dune::FieldVector<RF,dim+1> u_s(0.0);
342 Dune::FieldVector<RF,dim+1> u_n(0.0);
343 Dune::FieldVector<RF,dim+1> f(0.0);
346 const int order_s = dgspace_s.finiteElement().localBasis().order();
347 const int order_n = dgspace_n.finiteElement().localBasis().order();
348 const int intorder = overintegration+1+2*std::max(order_s,order_n);
352 auto iplocal_s = geo_in_inside.global(ip.position());
353 auto iplocal_n = geo_in_outside.global(ip.position());
356 auto& phi_s = cache[order_s].evaluateFunction(iplocal_s,dgspace_s.finiteElement().localBasis());
357 auto& phi_n = cache[order_n].evaluateFunction(iplocal_n,dgspace_n.finiteElement().localBasis());
361 for (size_type i=0; i<=dim; i++)
362 for (size_type k=0; k<dgspace_s.size(); k++)
363 u_s[i] += x_s(lfsv_s.child(i),k)*phi_s[k];
365 for (size_type i=0; i<=dim; i++)
366 for (size_type k=0; k<dgspace_n.size(); k++)
367 u_n[i] += x_n(lfsv_n.child(i),k)*phi_n[k];
373 A_minus_n.umv(u_n,f);
377 auto factor = ip.weight() * geo.integrationElement(ip.position());
378 for (size_type k=0; k<dgspace_s.size(); k++)
379 for (size_type i=0; i<=dim; i++)
380 r_s.accumulate(lfsv_s.child(i),k, f[i]*phi_s[k]*factor);
381 for (size_type k=0; k<dgspace_n.size(); k++)
382 for (size_type i=0; i<=dim; i++)
383 r_n.accumulate(lfsv_n.child(i),k, - f[i]*phi_n[k]*factor);
396 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
398 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
402 using namespace Indices;
403 using DGSpace = TypeTree::Child<LFSV,_0>;
404 using DF =
typename DGSpace::Traits::FiniteElementType::
405 Traits::LocalBasisType::Traits::DomainFieldType;
406 using RF =
typename DGSpace::Traits::FiniteElementType::
407 Traits::LocalBasisType::Traits::RangeFieldType;
408 using size_type =
typename DGSpace::Traits::SizeType;
411 const auto& dgspace_s = child(lfsv_s,_0);
414 const Dune::FieldVector<DF,dim> n_F =
ig.centerUnitOuterNormal();
417 const auto& cell_inside =
ig.inside();
420 auto geo =
ig.geometry();
421 auto geo_inside = cell_inside.geometry();
424 auto geo_in_inside =
ig.geometryInInside();
427 auto ref_el_inside = referenceElement(geo_inside);
428 auto local_inside = ref_el_inside.position(0,0);
429 auto c_s = param.c(cell_inside,local_inside);
432 Dune::FieldMatrix<DF,dim+1,dim+1> RT;
434 Dune::FieldVector<DF,dim+1> alpha;
435 for (
int i=0; i<=dim; i++) alpha[i] = RT[dim-1][i];
436 Dune::FieldVector<DF,dim+1> unit(0.0);
438 Dune::FieldVector<DF,dim+1> beta;
440 Dune::FieldMatrix<DF,dim+1,dim+1> A_plus_s;
441 for (
int i=0; i<=dim; i++)
442 for (
int j=0; j<=dim; j++)
443 A_plus_s[i][j] = c_s*alpha[i]*beta[j];
447 for (
int i=0; i<=dim; i++) alpha[i] = RT[dim][i];
451 Dune::FieldMatrix<DF,dim+1,dim+1> A_minus_n;
452 for (
int i=0; i<=dim; i++)
453 for (
int j=0; j<=dim; j++)
454 A_minus_n[i][j] = -c_s*alpha[i]*beta[j];
457 Dune::FieldVector<RF,dim+1> u_s(0.0);
458 Dune::FieldVector<RF,dim+1> f(0.0);
461 const int order_s = dgspace_s.finiteElement().localBasis().order();
462 const int intorder = overintegration+1+2*order_s;
466 auto iplocal_s = geo_in_inside.global(ip.position());
469 auto& phi_s = cache[order_s].evaluateFunction(iplocal_s,dgspace_s.finiteElement().localBasis());
473 for (size_type i=0; i<=dim; i++)
474 for (size_type k=0; k<dgspace_s.size(); k++)
475 u_s[i] += x_s(lfsv_s.child(i),k)*phi_s[k];
479 Dune::FieldVector<RF,dim+1> u_n(param.g(
ig.intersection(),ip.position(),u_s));
486 A_minus_n.umv(u_n,f);
490 auto factor = ip.weight() * geo.integrationElement(ip.position());
491 for (size_type k=0; k<dgspace_s.size(); k++)
492 for (size_type i=0; i<=dim; i++)
493 r_s.accumulate(lfsv_s.child(i),k, f[i]*phi_s[k]*factor);
502 template<
typename EG,
typename LFSV,
typename R>
506 using namespace Indices;
507 using DGSpace = TypeTree::Child<LFSV,_0>;
508 using size_type =
typename DGSpace::Traits::SizeType;
511 const DGSpace& dgspace = child(lfsv,_0);
514 const auto& cell = eg.entity();
517 auto geo = eg.geometry();
520 const int order_s = dgspace.finiteElement().localBasis().order();
521 const int intorder = overintegration+2*order_s;
525 auto q(param.q(cell,ip.position()));
528 auto& phi = cache[order_s].evaluateFunction(ip.position(),dgspace.finiteElement().localBasis());
531 auto factor = ip.weight() * geo.integrationElement(ip.position());
532 for (size_type k=0; k<=dim; k++)
533 for (size_type i=0; i<dgspace.size(); i++)
534 r.accumulate(lfsv.child(k),i, - q[k]*phi[i]*factor);
539 void setTime (
typename T::Traits::RangeFieldType t)
544 void preStep (
typename T::Traits::RangeFieldType time,
typename T::Traits::RangeFieldType dt,
550 void preStage (
typename T::Traits::RangeFieldType time,
int r)
560 typename T::Traits::RangeFieldType
suggestTimestep (
typename T::Traits::RangeFieldType dt)
const
568 using LocalBasisType =
typename FEM::Traits::FiniteElementType::Traits::LocalBasisType;
570 std::vector<Cache> cache;
581 template<
typename T,
typename FEM>
587 enum { dim = T::Traits::GridViewType::dimension };
596 : param(param_), overintegration(overintegration_), cache(20)
600 template<
typename LFSU,
typename LFSV,
typename LocalPattern>
602 LocalPattern& pattern)
const
605 if (TypeTree::degree(lfsu)!=TypeTree::degree(lfsv))
606 DUNE_THROW(Dune::Exception,
"need U=V!");
607 if (TypeTree::degree(lfsv)!=dim+1)
608 DUNE_THROW(Dune::Exception,
"need exactly dim+1 components!");
610 for (
size_t k=0; k<TypeTree::degree(lfsv); k++)
611 for (
size_t i=0; i<lfsv.child(k).size(); ++i)
612 for (
size_t j=0; j<lfsu.child(k).size(); ++j)
613 pattern.addLink(lfsv.child(k),i,lfsu.child(k),j);
617 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
618 void alpha_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
621 using namespace Indices;
622 using DGSpace = TypeTree::Child<LFSV,_0>;
623 using RF =
typename DGSpace::Traits::FiniteElementType::
624 Traits::LocalBasisType::Traits::RangeFieldType;
625 using size_type =
typename DGSpace::Traits::SizeType;
628 const auto& dgspace = child(lfsv,_0);
631 auto geo = eg.geometry();
634 Dune::FieldVector<RF,dim+1> u(0.0);
637 const int order = dgspace.finiteElement().localBasis().order();
638 const int intorder = overintegration+2*order;
642 auto& phi = cache[order].evaluateFunction(ip.position(),dgspace.finiteElement().localBasis());
646 for (size_type k=0; k<=dim; k++)
647 for (size_type j=0; j<dgspace.size(); j++)
648 u[k] += x(lfsv.child(k),j)*phi[j];
651 auto factor = ip.weight() * geo.integrationElement(ip.position());
652 for (size_type k=0; k<=dim; k++)
653 for (size_type i=0; i<dgspace.size(); i++)
654 r.accumulate(lfsv.child(k),i, u[k]*phi[i]*factor);
659 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
664 using namespace Indices;
665 using DGSpace = TypeTree::Child<LFSV,_0>;
666 using size_type =
typename DGSpace::Traits::SizeType;
669 const auto& dgspace = child(lfsv,_0);
672 auto geo = eg.geometry();
675 const int order = dgspace.finiteElement().localBasis().order();
676 const int intorder = overintegration+2*order;
680 auto& phi = cache[order].evaluateFunction(ip.position(),dgspace.finiteElement().localBasis());
683 auto factor = ip.weight() * geo.integrationElement(ip.position());
684 for (size_type k=0; k<=dim; k++)
685 for (size_type i=0; i<dgspace.size(); i++)
686 for (size_type j=0; j<dgspace.size(); j++)
687 mat.accumulate(lfsv.child(k),i,lfsu.child(k),j, phi[j]*phi[i]*factor);
694 using LocalBasisType =
typename FEM::Traits::FiniteElementType::Traits::LocalBasisType;
696 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:123
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: linearacousticsdg.hh:28
static void eigenvectors_transposed(T1 c, const Dune::FieldVector< T2, dim > &n, Dune::FieldMatrix< T3, dim+1, dim+1 > &RT)
Definition: linearacousticsdg.hh:44
static void eigenvectors(T1 c, const Dune::FieldVector< T2, dim > &n, Dune::FieldMatrix< T3, dim+1, dim+1 > &RT)
Definition: linearacousticsdg.hh:51
static void eigenvectors_transposed(T1 c, const Dune::FieldVector< T2, dim > &n, Dune::FieldMatrix< T3, dim+1, dim+1 > &RT)
Definition: linearacousticsdg.hh:73
static void eigenvectors(T1 c, const Dune::FieldVector< T2, dim > &n, Dune::FieldMatrix< T3, dim+1, dim+1 > &RT)
Definition: linearacousticsdg.hh:81
static void eigenvectors_transposed(T1 c, const Dune::FieldVector< T2, dim > &n, Dune::FieldMatrix< T3, dim+1, dim+1 > &RT)
Definition: linearacousticsdg.hh:103
static void eigenvectors(T1 c, const Dune::FieldVector< T2, dim > &n, Dune::FieldMatrix< T3, dim+1, dim+1 > &RT)
Definition: linearacousticsdg.hh:124
Definition: linearacousticsdg.hh:173
void preStep(typename T::Traits::RangeFieldType time, typename T::Traits::RangeFieldType dt, int stages)
to be called once before each time step
Definition: linearacousticsdg.hh:544
@ doPatternVolume
Definition: linearacousticsdg.hh:178
void preStage(typename T::Traits::RangeFieldType time, int r)
to be called once before each stage
Definition: linearacousticsdg.hh:550
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: linearacousticsdg.hh:195
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: linearacousticsdg.hh:503
@ doPatternSkeleton
Definition: linearacousticsdg.hh:179
void setTime(typename T::Traits::RangeFieldType t)
set time in parameter class
Definition: linearacousticsdg.hh:539
T::Traits::RangeFieldType suggestTimestep(typename T::Traits::RangeFieldType dt) const
to be called once before each stage
Definition: linearacousticsdg.hh:560
@ doAlphaVolume
Definition: linearacousticsdg.hh:182
void alpha_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, R &r_s) const
Definition: linearacousticsdg.hh:397
@ doLambdaVolume
Definition: linearacousticsdg.hh:185
void postStage()
to be called once at the end of each stage
Definition: linearacousticsdg.hh:555
DGLinearAcousticsSpatialOperator(T ¶m_, int overintegration_=0)
Definition: linearacousticsdg.hh:188
@ doAlphaBoundary
Definition: linearacousticsdg.hh:184
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: linearacousticsdg.hh:273
@ doAlphaSkeleton
Definition: linearacousticsdg.hh:183
Definition: linearacousticsdg.hh:586
DGLinearAcousticsTemporalOperator(T ¶m_, int overintegration_=0)
Definition: linearacousticsdg.hh:595
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: linearacousticsdg.hh:618
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: linearacousticsdg.hh:660
void pattern_volume(const LFSU &lfsu, const LFSV &lfsv, LocalPattern &pattern) const
Definition: linearacousticsdg.hh:601
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