2#ifndef DUNE_PDELAB_LOCALOPERATOR_LINEARACOUSTICSDG_HH
3#define DUNE_PDELAB_LOCALOPERATOR_LINEARACOUSTICSDG_HH
9#include<dune/geometry/referenceelements.hh>
11#include<dune/pdelab/common/quadraturerules.hh>
12#include<dune/pdelab/common/geometrywrapper.hh>
13#include<dune/pdelab/common/function.hh>
14#include<dune/pdelab/localoperator/pattern.hh>
15#include<dune/pdelab/localoperator/flags.hh>
16#include<dune/pdelab/localoperator/idefault.hh>
17#include<dune/pdelab/localoperator/defaultimp.hh>
18#include<dune/pdelab/finiteelement/localbasiscache.hh>
20#include"linearacousticsparameter.hh"
27 class LinearAcousticsEigenvectors
34 class LinearAcousticsEigenvectors<1>
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>
53 RT[0][0] = 1; RT[1][0] = c*n[0];
54 RT[0][1] = -1; RT[1][1] = c*n[0];
63 class LinearAcousticsEigenvectors<2>
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>
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];
93 class LinearAcousticsEigenvectors<3>
102 template<
typename T1,
typename T2,
typename T3>
106 s[0] = n[1]-n[2]; s[1] = n[2]-n[0]; s[2] = n[0]-n[1];
107 if (s.two_norm()<1e-5)
109 s[0] = n[1]+n[2]; s[1] = n[2]-n[0]; s[2] = -(n[0]+n[1]);
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>
127 s[0] = n[1]-n[2]; s[1] = n[2]-n[0]; s[2] = n[0]-n[1];
128 if (s.two_norm()<1e-5)
130 s[0] = n[1]+n[2]; s[1] = n[2]-n[0]; s[2] = -(n[0]+n[1]);
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 };
178 enum { doPatternVolume =
true };
179 enum { doPatternSkeleton =
true };
182 enum { doAlphaVolume =
true };
183 enum { doAlphaSkeleton =
true };
184 enum { doAlphaBoundary =
true };
185 enum { doLambdaVolume =
true };
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
200 using RF =
typename DGSpace::Traits::FiniteElementType::
201 Traits::LocalBasisType::Traits::RangeFieldType;
202 using size_type =
typename DGSpace::Traits::SizeType;
209 const auto& dgspace =
child(lfsv,
_0);
212 const auto& cell = eg.entity();
215 auto geo = eg.geometry();
219 auto localcenter = ref_el.position(0,0);
220 auto c2 = param.c(cell,localcenter);
224 typename EG::Geometry::JacobianInverseTransposed jac;
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;
233 for (
const auto& ip : quadratureRule(geo,intorder))
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>
273 void alpha_skeleton (
const IG& ig,
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
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);
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();
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);
317 LinearAcousticsEigenvectors<dim>::eigenvectors_transposed(c_s,n_F,RT);
319 for (
int i=0; i<=dim; i++) alpha[i] = RT[dim-1][i];
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];
330 LinearAcousticsEigenvectors<dim>::eigenvectors_transposed(c_n,n_F,RT);
331 for (
int i=0; i<=dim; i++) alpha[i] = RT[dim][i];
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];
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);
349 for (
const auto& ip : quadratureRule(geo,intorder))
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>
397 void alpha_boundary (
const IG& ig,
398 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
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);
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();
428 auto local_inside = ref_el_inside.position(0,0);
429 auto c_s = param.c(cell_inside,local_inside);
433 LinearAcousticsEigenvectors<dim>::eigenvectors_transposed(c_s,n_F,RT);
435 for (
int i=0; i<=dim; i++) alpha[i] = RT[dim-1][i];
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];
446 LinearAcousticsEigenvectors<dim>::eigenvectors_transposed(c_s,n_F,RT);
447 for (
int i=0; i<=dim; i++) alpha[i] = RT[dim][i];
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];
461 const int order_s = dgspace_s.finiteElement().localBasis().order();
462 const int intorder = overintegration+1+2*order_s;
463 for (
const auto& ip : quadratureRule(geo,intorder))
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];
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>
503 void lambda_volume (
const EG& eg,
const LFSV& lfsv, R& r)
const
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;
522 for (
const auto& ip : quadratureRule(geo,intorder))
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 };
590 enum { doPatternVolume =
true };
593 enum { doAlphaVolume =
true };
596 : param(param_), overintegration(overintegration_), cache(20)
600 template<
typename LFSU,
typename LFSV,
typename LocalPattern>
601 void pattern_volume (
const LFSU& lfsu,
const LFSV& lfsv,
602 LocalPattern& pattern)
const
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
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();
637 const int order = dgspace.finiteElement().localBasis().order();
638 const int intorder = overintegration+2*order;
639 for (
const auto& ip : quadratureRule(geo,intorder))
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>
660 void jacobian_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
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;
677 for (
const auto& ip : quadratureRule(geo,intorder))
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;
void solve(V1 &x, const V2 &b, bool doPivoting=true) const
Solve system A x = b.
void umv(const X &x, Y &y) const
y += A x
Definition: densematrix.hh:434
Base class for Dune-Exceptions.
Definition: exceptions.hh:94
A dense n x m matrix.
Definition: fmatrix.hh:69
vector space out of a tensor product of fields.
Definition: fvector.hh:95
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
void preStage(typename T::Traits::RangeFieldType time, int r)
to be called once before each stage
Definition: linearacousticsdg.hh:550
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
void postStage()
to be called once at the end of each stage
Definition: linearacousticsdg.hh:555
Definition: linearacousticsdg.hh:586
sparsity pattern generator
Definition: pattern.hh:30
sparsity pattern generator
Definition: pattern.hh:14
Default class for additional methods in instationary local operators.
Definition: idefault.hh:90
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_transposed(T1 c, const Dune::FieldVector< T2, dim > &n, Dune::FieldMatrix< T3, dim+1, dim+1 > &RT)
Definition: linearacousticsdg.hh:73
static void eigenvectors_transposed(T1 c, const Dune::FieldVector< T2, dim > &n, Dune::FieldMatrix< T3, dim+1, dim+1 > &RT)
Definition: linearacousticsdg.hh:103
store values of basis functions and gradients in a cache
Definition: localbasiscache.hh:18
Default flags for all local operators.
Definition: flags.hh:19
Implements linear and nonlinear versions of jacobian_apply_boundary() based on alpha_boundary()
Definition: numericaljacobianapply.hh:287
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_volume() based on alpha_volume()
Definition: numericaljacobianapply.hh:34
Implement jacobian_boundary() based on alpha_boundary()
Definition: numericaljacobian.hh:251
Implement jacobian_skeleton() based on alpha_skeleton()
Definition: numericaljacobian.hh:157
Implement jacobian_volume() based on alpha_volume()
Definition: numericaljacobian.hh:32
A few common exception classes.
Implements a vector constructed from a given type representing a field and a compile-time given size.
constexpr index_constant< 0 > _0
Compile time index with value 0.
Definition: indices.hh:51
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
unspecified value type referenceElement(T &&... t)
Returns a reference element for the objects t....
std::size_t degree(const Node &node)
Returns the degree of node as run time information.
Definition: nodeinterface.hh:76
auto max(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::max()
Definition: defaults.hh:79
Namespace with predefined compile time indices for the range [0,19].
Definition: indices.hh:49
Dune namespace.
Definition: alignedallocator.hh:11