2 #ifndef DUNE_PDELAB_LINEARACOUSTICSDG_HH 3 #define DUNE_PDELAB_LINEARACOUSTICSDG_HH 7 #include<dune/common/exceptions.hh> 8 #include<dune/common/fvector.hh> 9 #include<dune/geometry/referenceelements.hh> 44 template<
typename T1,
typename T2,
typename T3>
47 RT[0][0] = 1; RT[0][1] = c*n[0];
48 RT[1][0] = -1; RT[1][1] = c*n[0];
51 template<
typename T1,
typename T2,
typename T3>
52 static void eigenvectors (T1 c,
const Dune::FieldVector<T2,dim>& n, Dune::FieldMatrix<T3,dim+1,dim+1>& RT)
54 RT[0][0] = 1; RT[1][0] = c*n[0];
55 RT[0][1] = -1; RT[1][1] = c*n[0];
73 template<
typename T1,
typename T2,
typename T3>
76 RT[0][0] = 0; RT[0][1] = -n[1]; RT[0][2] = n[0];
77 RT[1][0] = 1; RT[1][1] = c*n[0]; RT[1][2] = c*n[1];
78 RT[2][0] = -1; RT[2][1] = c*n[0]; RT[2][2] = c*n[1];
81 template<
typename T1,
typename T2,
typename T3>
82 static void eigenvectors (T1 c,
const Dune::FieldVector<T2,dim>& n, Dune::FieldMatrix<T3,dim+1,dim+1>& RT)
84 RT[0][0] = 0; RT[1][0] = -n[1]; RT[2][0] = n[0];
85 RT[0][1] = 1; RT[1][1] = c*n[0]; RT[2][1] = c*n[1];
86 RT[0][2] = -1; RT[1][2] = c*n[0]; RT[2][2] = c*n[1];
103 template<
typename T1,
typename T2,
typename T3>
106 Dune::FieldVector<T2,dim>
s;
107 s[0] = n[1]-n[2]; s[1] = n[2]-n[0]; s[2] = n[0]-n[1];
108 if (s.two_norm()<1
e-5)
110 s[0] = n[1]+n[2]; s[1] = n[2]-n[0]; s[2] = -(n[0]+n[1]);
113 Dune::FieldVector<T2,dim> t;
114 t[0] = n[1]*s[2] - n[2]*s[1];
115 t[1] = n[2]*s[0] - n[0]*s[2];
116 t[2] = n[0]*s[1] - n[1]*s[0];
118 RT[0][0] = 0; RT[0][1] = s[0]; RT[0][2] = s[1]; RT[0][3] = s[2];
119 RT[1][0] = 0; RT[1][1] = t[0]; RT[1][2] = t[1]; RT[1][3] = t[2];
120 RT[2][0] = 1; RT[2][1] = c*n[0]; RT[2][2] = c*n[1]; RT[2][3] = c*n[2];
121 RT[3][0] = -1; RT[3][1] = c*n[0]; RT[3][2] = c*n[1]; RT[3][3] = c*n[2];
124 template<
typename T1,
typename T2,
typename T3>
125 static void eigenvectors (T1 c,
const Dune::FieldVector<T2,dim>& n, Dune::FieldMatrix<T3,dim+1,dim+1>& RT)
127 Dune::FieldVector<T2,dim>
s;
128 s[0] = n[1]-n[2]; s[1] = n[2]-n[0]; s[2] = n[0]-n[1];
129 if (s.two_norm()<1
e-5)
131 s[0] = n[1]+n[2]; s[1] = n[2]-n[0]; s[2] = -(n[0]+n[1]);
134 Dune::FieldVector<T2,dim> t;
135 t[0] = n[1]*s[2] - n[2]*s[1];
136 t[1] = n[2]*s[0] - n[0]*s[2];
137 t[2] = n[0]*s[1] - n[1]*s[0];
139 RT[0][0] = 0; RT[1][0] = s[0]; RT[2][0] = s[1]; RT[3][0] = s[2];
140 RT[0][1] = 0; RT[1][1] = t[0]; RT[2][1] = t[1]; RT[3][1] = t[2];
141 RT[0][2] = 1; RT[1][2] = c*n[0]; RT[2][2] = c*n[1]; RT[3][2] = c*n[2];
142 RT[0][3] = -1; RT[1][3] = c*n[0]; RT[2][3] = c*n[1]; RT[3][3] = c*n[2];
162 template<
typename T,
typename FEM>
175 enum {
dim = T::Traits::GridViewType::dimension };
179 enum { doPatternVolume =
true };
180 enum { doPatternSkeleton =
true };
183 enum { doAlphaVolume =
true };
184 enum { doAlphaSkeleton =
true };
185 enum { doAlphaBoundary =
true };
186 enum { doLambdaVolume =
true };
190 : param(param_), overintegration(overintegration_), cache(20)
195 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
196 void alpha_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const 199 using namespace TypeTree::Indices;
200 using DGSpace = TypeTree::Child<LFSV,_0>;
201 using RF =
typename DGSpace::Traits::FiniteElementType::
202 Traits::LocalBasisType::Traits::RangeFieldType;
203 using size_type =
typename DGSpace::Traits::SizeType;
206 if (LFSV::CHILDREN!=
dim+1)
207 DUNE_THROW(Dune::Exception,
"need exactly dim+1 components!");
210 const auto& dgspace = child(lfsv,_0);
213 const auto& cell = eg.entity();
216 auto geo = eg.geometry();
220 auto localcenter = ref_el.position(0,0);
221 auto c2 = param.c(cell,localcenter);
225 typename EG::Geometry::JacobianInverseTransposed jac;
228 Dune::FieldVector<RF,dim+1> u(0.0);
229 std::vector<Dune::FieldVector<RF,dim> > gradphi(dgspace.size());
232 const int order = dgspace.finiteElement().localBasis().order();
233 const int intorder = overintegration+2*order;
237 auto& phi = cache[order].evaluateFunction(ip.position(),dgspace.finiteElement().localBasis());
241 for (size_type k=0; k<=
dim; k++)
242 for (size_type j=0; j<dgspace.size(); j++)
243 u[k] += x(lfsv.child(k),j)*phi[j];
247 auto& js = cache[order].evaluateJacobian(ip.position(),dgspace.finiteElement().localBasis());
250 jac = geo.jacobianInverseTransposed(ip.position());
251 for (size_type i=0; i<dgspace.size(); i++)
252 jac.mv(js[i][0],gradphi[i]);
255 auto factor = ip.weight() * geo.integrationElement(ip.position());
256 for (size_type k=0; k<dgspace.size(); k++)
259 for (size_type j=0; j<
dim; j++)
260 r.accumulate(lfsv.child(0),k, - u[j+1]*gradphi[k][j]*factor);
262 for (size_type i=1; i<=
dim; i++)
263 r.accumulate(lfsv.child(i),k, - c2*u[0]*gradphi[k][i-1]*factor);
273 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
275 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
276 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
277 R& r_s, R& r_n)
const 280 using namespace TypeTree::Indices;
281 using DGSpace = TypeTree::Child<LFSV,_0>;
282 using DF =
typename DGSpace::Traits::FiniteElementType::
283 Traits::LocalBasisType::Traits::DomainFieldType;
284 using RF =
typename DGSpace::Traits::FiniteElementType::
285 Traits::LocalBasisType::Traits::RangeFieldType;
286 using size_type =
typename DGSpace::Traits::SizeType;
289 const auto& dgspace_s = child(lfsv_s,_0);
290 const auto& dgspace_n = child(lfsv_n,_0);
293 const Dune::FieldVector<DF,dim> n_F = ig.centerUnitOuterNormal();
296 const auto& cell_inside = ig.inside();
297 const auto& cell_outside = ig.outside();
300 auto geo = ig.geometry();
301 auto geo_inside = cell_inside.geometry();
302 auto geo_outside = cell_outside.geometry();
305 auto geo_in_inside = ig.geometryInInside();
306 auto geo_in_outside = ig.geometryInOutside();
311 auto local_inside = ref_el_inside.position(0,0);
312 auto local_outside = ref_el_outside.position(0,0);
313 auto c_s = param.c(cell_inside,local_inside);
314 auto c_n = param.c(cell_outside,local_outside);
317 Dune::FieldMatrix<DF,dim+1,dim+1> RT;
319 Dune::FieldVector<DF,dim+1> alpha;
320 for (
int i=0; i<=
dim; i++) alpha[i] = RT[
dim-1][i];
321 Dune::FieldVector<DF,dim+1> unit(0.0);
323 Dune::FieldVector<DF,dim+1> beta;
325 Dune::FieldMatrix<DF,dim+1,dim+1> A_plus_s;
326 for (
int i=0; i<=
dim; i++)
327 for (
int j=0; j<=
dim; j++)
328 A_plus_s[i][j] = c_s*alpha[i]*beta[j];
332 for (
int i=0; i<=
dim; i++) alpha[i] = RT[
dim][i];
336 Dune::FieldMatrix<DF,dim+1,dim+1> A_minus_n;
337 for (
int i=0; i<=
dim; i++)
338 for (
int j=0; j<=
dim; j++)
339 A_minus_n[i][j] = -c_n*alpha[i]*beta[j];
342 Dune::FieldVector<RF,dim+1> u_s(0.0);
343 Dune::FieldVector<RF,dim+1> u_n(0.0);
344 Dune::FieldVector<RF,dim+1> f(0.0);
347 const int order_s = dgspace_s.finiteElement().localBasis().order();
348 const int order_n = dgspace_n.finiteElement().localBasis().order();
349 const int intorder = overintegration+1+2*std::max(order_s,order_n);
353 auto iplocal_s = geo_in_inside.global(ip.position());
354 auto iplocal_n = geo_in_outside.global(ip.position());
357 auto& phi_s = cache[order_s].evaluateFunction(iplocal_s,dgspace_s.finiteElement().localBasis());
358 auto& phi_n = cache[order_n].evaluateFunction(iplocal_n,dgspace_n.finiteElement().localBasis());
362 for (size_type i=0; i<=
dim; i++)
363 for (size_type k=0; k<dgspace_s.size(); k++)
364 u_s[i] += x_s(lfsv_s.child(i),k)*phi_s[k];
366 for (size_type i=0; i<=
dim; i++)
367 for (size_type k=0; k<dgspace_n.size(); k++)
368 u_n[i] += x_n(lfsv_n.child(i),k)*phi_n[k];
374 A_minus_n.umv(u_n,f);
378 auto factor = ip.weight() * geo.integrationElement(ip.position());
379 for (size_type k=0; k<dgspace_s.size(); k++)
380 for (size_type i=0; i<=
dim; i++)
381 r_s.accumulate(lfsv_s.child(i),k, f[i]*phi_s[k]*factor);
382 for (size_type k=0; k<dgspace_n.size(); k++)
383 for (size_type i=0; i<=
dim; i++)
384 r_n.accumulate(lfsv_n.child(i),k, - f[i]*phi_n[k]*factor);
397 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
399 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
403 using namespace TypeTree::Indices;
404 using DGSpace = TypeTree::Child<LFSV,_0>;
405 using DF =
typename DGSpace::Traits::FiniteElementType::
406 Traits::LocalBasisType::Traits::DomainFieldType;
407 using RF =
typename DGSpace::Traits::FiniteElementType::
408 Traits::LocalBasisType::Traits::RangeFieldType;
409 using size_type =
typename DGSpace::Traits::SizeType;
412 const auto& dgspace_s = child(lfsv_s,_0);
415 const Dune::FieldVector<DF,dim> n_F = ig.centerUnitOuterNormal();
418 const auto& cell_inside = ig.inside();
421 auto geo = ig.geometry();
422 auto geo_inside = cell_inside.geometry();
425 auto geo_in_inside = ig.geometryInInside();
429 auto local_inside = ref_el_inside.position(0,0);
430 auto c_s = param.c(cell_inside,local_inside);
433 Dune::FieldMatrix<DF,dim+1,dim+1> RT;
435 Dune::FieldVector<DF,dim+1> alpha;
436 for (
int i=0; i<=
dim; i++) alpha[i] = RT[
dim-1][i];
437 Dune::FieldVector<DF,dim+1> unit(0.0);
439 Dune::FieldVector<DF,dim+1> beta;
441 Dune::FieldMatrix<DF,dim+1,dim+1> A_plus_s;
442 for (
int i=0; i<=
dim; i++)
443 for (
int j=0; j<=
dim; j++)
444 A_plus_s[i][j] = c_s*alpha[i]*beta[j];
448 for (
int i=0; i<=
dim; i++) alpha[i] = RT[
dim][i];
452 Dune::FieldMatrix<DF,dim+1,dim+1> A_minus_n;
453 for (
int i=0; i<=
dim; i++)
454 for (
int j=0; j<=
dim; j++)
455 A_minus_n[i][j] = -c_s*alpha[i]*beta[j];
458 Dune::FieldVector<RF,dim+1> u_s(0.0);
459 Dune::FieldVector<RF,dim+1> f(0.0);
462 const int order_s = dgspace_s.finiteElement().localBasis().order();
463 const int intorder = overintegration+1+2*order_s;
467 auto iplocal_s = geo_in_inside.global(ip.position());
470 auto& phi_s = cache[order_s].evaluateFunction(iplocal_s,dgspace_s.finiteElement().localBasis());
474 for (size_type i=0; i<=
dim; i++)
475 for (size_type k=0; k<dgspace_s.size(); k++)
476 u_s[i] += x_s(lfsv_s.child(i),k)*phi_s[k];
480 Dune::FieldVector<RF,dim+1> u_n(param.g(ig.intersection(),ip.position(),u_s));
487 A_minus_n.umv(u_n,f);
491 auto factor = ip.weight() * geo.integrationElement(ip.position());
492 for (size_type k=0; k<dgspace_s.size(); k++)
493 for (size_type i=0; i<=
dim; i++)
494 r_s.accumulate(lfsv_s.child(i),k, f[i]*phi_s[k]*factor);
503 template<
typename EG,
typename LFSV,
typename R>
507 using namespace TypeTree::Indices;
508 using DGSpace = TypeTree::Child<LFSV,_0>;
509 using size_type =
typename DGSpace::Traits::SizeType;
512 const DGSpace& dgspace = child(lfsv,_0);
515 const auto& cell = eg.entity();
518 auto geo = eg.geometry();
521 const int order_s = dgspace.finiteElement().localBasis().order();
522 const int intorder = overintegration+2*order_s;
526 auto q(param.q(cell,ip.position()));
529 auto& phi = cache[order_s].evaluateFunction(ip.position(),dgspace.finiteElement().localBasis());
532 auto factor = ip.weight() * geo.integrationElement(ip.position());
533 for (size_type k=0; k<=
dim; k++)
534 for (size_type i=0; i<dgspace.size(); i++)
535 r.accumulate(lfsv.child(k),i, - q[k]*phi[i]*factor);
540 void setTime (
typename T::Traits::RangeFieldType t)
545 void preStep (
typename T::Traits::RangeFieldType time,
typename T::Traits::RangeFieldType dt,
551 void preStage (
typename T::Traits::RangeFieldType time,
int r)
561 typename T::Traits::RangeFieldType
suggestTimestep (
typename T::Traits::RangeFieldType dt)
const 569 using LocalBasisType =
typename FEM::Traits::FiniteElementType::Traits::LocalBasisType;
571 std::vector<Cache> cache;
582 template<
typename T,
typename FEM>
588 enum {
dim = T::Traits::GridViewType::dimension };
591 enum { doPatternVolume =
true };
594 enum { doAlphaVolume =
true };
597 : param(param_), overintegration(overintegration_), cache(20)
601 template<
typename LFSU,
typename LFSV,
typename LocalPattern>
603 LocalPattern& pattern)
const 606 if (LFSU::CHILDREN!=LFSV::CHILDREN)
607 DUNE_THROW(Dune::Exception,
"need U=V!");
608 if (LFSV::CHILDREN!=
dim+1)
609 DUNE_THROW(Dune::Exception,
"need exactly dim+1 components!");
611 for (
size_t k=0; k<LFSV::CHILDREN; k++)
612 for (
size_t i=0; i<lfsv.child(k).size(); ++i)
613 for (
size_t j=0; j<lfsu.child(k).size(); ++j)
614 pattern.addLink(lfsv.child(k),i,lfsu.child(k),j);
618 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
619 void alpha_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const 622 using namespace TypeTree::Indices;
623 using DGSpace = TypeTree::Child<LFSV,_0>;
624 using RF =
typename DGSpace::Traits::FiniteElementType::
625 Traits::LocalBasisType::Traits::RangeFieldType;
626 using size_type =
typename DGSpace::Traits::SizeType;
629 const auto& dgspace = child(lfsv,_0);
632 auto geo = eg.geometry();
635 Dune::FieldVector<RF,dim+1> u(0.0);
638 const int order = dgspace.finiteElement().localBasis().order();
639 const int intorder = overintegration+2*order;
643 auto& phi = cache[order].evaluateFunction(ip.position(),dgspace.finiteElement().localBasis());
647 for (size_type k=0; k<=
dim; k++)
648 for (size_type j=0; j<dgspace.size(); j++)
649 u[k] += x(lfsv.child(k),j)*phi[j];
652 auto factor = ip.weight() * geo.integrationElement(ip.position());
653 for (size_type k=0; k<=
dim; k++)
654 for (size_type i=0; i<dgspace.size(); i++)
655 r.accumulate(lfsv.child(k),i, u[k]*phi[i]*factor);
660 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
665 using namespace TypeTree::Indices;
666 using DGSpace = TypeTree::Child<LFSV,_0>;
667 using size_type =
typename DGSpace::Traits::SizeType;
670 const auto& dgspace = child(lfsv,_0);
673 auto geo = eg.geometry();
676 const int order = dgspace.finiteElement().localBasis().order();
677 const int intorder = overintegration+2*order;
681 auto& phi = cache[order].evaluateFunction(ip.position(),dgspace.finiteElement().localBasis());
684 auto factor = ip.weight() * geo.integrationElement(ip.position());
685 for (size_type k=0; k<=
dim; k++)
686 for (size_type i=0; i<dgspace.size(); i++)
687 for (size_type j=0; j<dgspace.size(); j++)
688 mat.accumulate(lfsv.child(k),i,lfsu.child(k),j, phi[j]*phi[i]*factor);
695 using LocalBasisType =
typename FEM::Traits::FiniteElementType::Traits::LocalBasisType;
697 std::vector<Cache> cache;
Definition: linearacousticsdg.hh:583
static void eigenvectors(T1 c, const Dune::FieldVector< T2, dim > &n, Dune::FieldMatrix< T3, dim+1, dim+1 > &RT)
Definition: linearacousticsdg.hh:82
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: linearacousticsdg.hh:619
sparsity pattern generator
Definition: pattern.hh:13
Implements linear and nonlinear versions of jacobian_apply_skeleton() based on alpha_skeleton() ...
Definition: numericaljacobianapply.hh:180
void preStage(typename T::Traits::RangeFieldType time, int r)
to be called once before each stage
Definition: linearacousticsdg.hh:551
const IG & ig
Definition: constraints.hh:148
static const int dim
Definition: adaptivity.hh:83
static void eigenvectors_transposed(T1 c, const Dune::FieldVector< T2, dim > &n, Dune::FieldMatrix< T3, dim+1, dim+1 > &RT)
Definition: linearacousticsdg.hh:104
Implement jacobian_boundary() based on alpha_boundary()
Definition: numericaljacobian.hh:250
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:398
T::Traits::RangeFieldType suggestTimestep(typename T::Traits::RangeFieldType dt) const
to be called once before each stage
Definition: linearacousticsdg.hh:561
ReferenceElementWrapper< ReferenceElement< typename Geometry::ctype, Geometry::mydimension > > referenceElement(const Geometry &geo)
Returns the reference element for the given geometry.
Definition: referenceelements.hh:144
static void eigenvectors_transposed(T1 c, const Dune::FieldVector< T2, dim > &n, Dune::FieldMatrix< T3, dim+1, dim+1 > &RT)
Definition: linearacousticsdg.hh:74
Definition: adaptivity.hh:27
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:545
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:274
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: linearacousticsdg.hh:196
void postStage()
to be called once at the end of each stage
Definition: linearacousticsdg.hh:556
sparsity pattern generator
Definition: pattern.hh:29
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
Definition: linearacousticsdg.hh:163
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 lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: linearacousticsdg.hh:504
DGLinearAcousticsTemporalOperator(T ¶m_, int overintegration_=0)
Definition: linearacousticsdg.hh:596
Implement jacobian_volume() based on alpha_volume()
Definition: numericaljacobian.hh:31
void setTime(typename T::Traits::RangeFieldType t)
set time in parameter class
Definition: linearacousticsdg.hh:540
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
void pattern_volume(const LFSU &lfsu, const LFSV &lfsv, LocalPattern &pattern) const
Definition: linearacousticsdg.hh:602
DGLinearAcousticsSpatialOperator(T ¶m_, int overintegration_=0)
Definition: linearacousticsdg.hh:189
static void eigenvectors(T1 c, const Dune::FieldVector< T2, dim > &n, Dune::FieldMatrix< T3, dim+1, dim+1 > &RT)
Definition: linearacousticsdg.hh:52
Definition: linearacousticsdg.hh:28
const std::string s
Definition: function.hh:1102
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: linearacousticsdg.hh:661
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
Implement jacobian_skeleton() based on alpha_skeleton()
Definition: numericaljacobian.hh:156
static void eigenvectors(T1 c, const Dune::FieldVector< T2, dim > &n, Dune::FieldMatrix< T3, dim+1, dim+1 > &RT)
Definition: linearacousticsdg.hh:125
static void eigenvectors_transposed(T1 c, const Dune::FieldVector< T2, dim > &n, Dune::FieldMatrix< T3, dim+1, dim+1 > &RT)
Definition: linearacousticsdg.hh:45