dune-pdelab  2.5-dev
linearacousticsdg.hh
Go to the documentation of this file.
1 // -*- tab-width: 8; indent-tabs-mode: nil -*-
2 #ifndef DUNE_PDELAB_LOCALOPERATOR_LINEARACOUSTICSDG_HH
3 #define DUNE_PDELAB_LOCALOPERATOR_LINEARACOUSTICSDG_HH
4 
5 #include<vector>
6 
7 #include<dune/common/exceptions.hh>
8 #include<dune/common/fvector.hh>
9 #include<dune/geometry/referenceelements.hh>
10 
20 
22 
23 namespace Dune {
24  namespace PDELab {
25 
26 
27  template<int dim>
29  {};
30 
34  template<>
36  {
37  enum { dim = 1 };
38  public:
44  template<typename T1, typename T2, typename T3>
45  static void eigenvectors_transposed (T1 c, const Dune::FieldVector<T2,dim>& n, Dune::FieldMatrix<T3,dim+1,dim+1>& RT)
46  {
47  RT[0][0] = 1; RT[0][1] = c*n[0];
48  RT[1][0] = -1; RT[1][1] = c*n[0];
49  }
50 
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)
53  {
54  RT[0][0] = 1; RT[1][0] = c*n[0];
55  RT[0][1] = -1; RT[1][1] = c*n[0];
56  }
57 
58  };
59 
63  template<>
65  {
66  enum { dim = 2 };
67  public:
73  template<typename T1, typename T2, typename T3>
74  static void eigenvectors_transposed (T1 c, const Dune::FieldVector<T2,dim>& n, Dune::FieldMatrix<T3,dim+1,dim+1>& RT)
75  {
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];
79  }
80 
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)
83  {
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];
87  }
88  };
89 
93  template<>
95  {
96  enum { dim = 3 };
97  public:
103  template<typename T1, typename T2, typename T3>
104  static void eigenvectors_transposed (T1 c, const Dune::FieldVector<T2,dim>& n, Dune::FieldMatrix<T3,dim+1,dim+1>& RT)
105  {
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()<1e-5)
109  {
110  s[0] = n[1]+n[2]; s[1] = n[2]-n[0]; s[2] = -(n[0]+n[1]);
111  }
112 
113  Dune::FieldVector<T2,dim> t; // compute cross product s * n
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];
117 
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];
122  }
123 
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)
126  {
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()<1e-5)
130  {
131  s[0] = n[1]+n[2]; s[1] = n[2]-n[0]; s[2] = -(n[0]+n[1]);
132  }
133 
134  Dune::FieldVector<T2,dim> t; // compute cross product s * n
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];
138 
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];
143  }
144  };
145 
162  template<typename T, typename FEM>
164  public NumericalJacobianApplyVolume<DGLinearAcousticsSpatialOperator<T,FEM> >,
165  public NumericalJacobianVolume<DGLinearAcousticsSpatialOperator<T,FEM> >,
166  public NumericalJacobianApplySkeleton<DGLinearAcousticsSpatialOperator<T,FEM> >,
167  public NumericalJacobianSkeleton<DGLinearAcousticsSpatialOperator<T,FEM> >,
168  public NumericalJacobianApplyBoundary<DGLinearAcousticsSpatialOperator<T,FEM> >,
169  public NumericalJacobianBoundary<DGLinearAcousticsSpatialOperator<T,FEM> >,
170  public FullSkeletonPattern,
171  public FullVolumePattern,
173  public InstationaryLocalOperatorDefaultMethods<typename T::Traits::RangeFieldType>
174  {
175  enum { dim = T::Traits::GridViewType::dimension };
176 
177  public:
178  // pattern assembly flags
179  enum { doPatternVolume = true };
180  enum { doPatternSkeleton = true };
181 
182  // residual assembly flags
183  enum { doAlphaVolume = true };
184  enum { doAlphaSkeleton = true };
185  enum { doAlphaBoundary = true };
186  enum { doLambdaVolume = true };
187 
188  // ! constructor
189  DGLinearAcousticsSpatialOperator (T& param_, int overintegration_=0)
190  : param(param_), overintegration(overintegration_), cache(20)
191  {
192  }
193 
194  // volume integral depending on test and ansatz functions
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
197  {
198  // Get types
199  using namespace 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;
204 
205  // paranoia check number of number of components
206  if (TypeTree::degree(lfsv)!=dim+1)
207  DUNE_THROW(Dune::Exception,"need exactly dim+1 components!");
208 
209  // get local function space that is identical for all components
210  const auto& dgspace = child(lfsv,_0);
211 
212  // Reference to cell
213  const auto& cell = eg.entity();
214 
215  // Get geometry
216  auto geo = eg.geometry();
217 
218  // evaluate speed of sound (assumed constant per element)
219  auto ref_el = referenceElement(geo);
220  auto localcenter = ref_el.position(0,0);
221  auto c2 = param.c(cell,localcenter);
222  c2 = c2*c2; // square it
223 
224  // Transformation
225  typename EG::Geometry::JacobianInverseTransposed jac;
226 
227  // Initialize vectors outside for loop
228  Dune::FieldVector<RF,dim+1> u(0.0);
229  std::vector<Dune::FieldVector<RF,dim> > gradphi(dgspace.size());
230 
231  // loop over quadrature points
232  const int order = dgspace.finiteElement().localBasis().order();
233  const int intorder = overintegration+2*order;
234  for (const auto& ip : quadratureRule(geo,intorder))
235  {
236  // evaluate basis functions
237  auto& phi = cache[order].evaluateFunction(ip.position(),dgspace.finiteElement().localBasis());
238 
239  // evaluate u
240  u = 0.0;
241  for (size_type k=0; k<=dim; k++) // for all components
242  for (size_type j=0; j<dgspace.size(); j++) // for all basis functions
243  u[k] += x(lfsv.child(k),j)*phi[j];
244  // std::cout << " u at " << ip.position() << " : " << u << std::endl;
245 
246  // evaluate gradient of basis functions (we assume Galerkin method lfsu=lfsv)
247  auto& js = cache[order].evaluateJacobian(ip.position(),dgspace.finiteElement().localBasis());
248 
249  // compute global gradients
250  jac = geo.jacobianInverseTransposed(ip.position());
251  for (size_type i=0; i<dgspace.size(); i++)
252  jac.mv(js[i][0],gradphi[i]);
253 
254  // integrate
255  auto factor = ip.weight() * geo.integrationElement(ip.position());
256  for (size_type k=0; k<dgspace.size(); k++) // loop over all vector-valued (!) basis functions (with identical components)
257  {
258  // component i=0
259  for (size_type j=0; j<dim; j++)
260  r.accumulate(lfsv.child(0),k, - u[j+1]*gradphi[k][j]*factor);
261  // components i=1...d
262  for (size_type i=1; i<=dim; i++)
263  r.accumulate(lfsv.child(i),k, - c2*u[0]*gradphi[k][i-1]*factor);
264  }
265  // std::cout << " residual: ";
266  // for (size_type i=0; i<r.size(); i++) std::cout << r[i] << " ";
267  // std::cout << std::endl;
268  }
269  }
270 
271  // skeleton integral depending on test and ansatz functions
272  // each face is only visited ONCE!
273  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
274  void alpha_skeleton (const IG& ig,
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
278  {
279  // Get types
280  using namespace 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;
287 
288  // Get local function space that is identical for all components
289  const auto& dgspace_s = child(lfsv_s,_0);
290  const auto& dgspace_n = child(lfsv_n,_0);
291 
292  // Normal: assume faces are planar
293  const Dune::FieldVector<DF,dim> n_F = ig.centerUnitOuterNormal();
294 
295  // References to inside and outside cells
296  const auto& cell_inside = ig.inside();
297  const auto& cell_outside = ig.outside();
298 
299  // Get geometries
300  auto geo = ig.geometry();
301  auto geo_inside = cell_inside.geometry();
302  auto geo_outside = cell_outside.geometry();
303 
304  // Get geometry of intersection in local coordinates of cell_inside and cell_outside
305  auto geo_in_inside = ig.geometryInInside();
306  auto geo_in_outside = ig.geometryInOutside();
307 
308  // Evaluate speed of sound (assumed constant per element)
309  auto ref_el_inside = referenceElement(geo_inside);
310  auto ref_el_outside = referenceElement(geo_outside);
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);
315 
316  // Compute A+ (outgoing waves)
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]; // row dim-1 corresponds to eigenvalue +c
321  Dune::FieldVector<DF,dim+1> unit(0.0);
322  unit[dim-1] = 1.0;
323  Dune::FieldVector<DF,dim+1> beta;
324  RT.solve(beta,unit);
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];
329 
330  // Compute A- (incoming waves)
332  for (int i=0; i<=dim; i++) alpha[i] = RT[dim][i]; // row dim corresponds to eigenvalue -c
333  unit = 0.0;
334  unit[dim] = 1.0;
335  RT.solve(beta,unit);
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];
340 
341  // Initialize vectors outside for loop
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);
345 
346  // Loop over quadrature points
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);
350  for (const auto& ip : quadratureRule(geo,intorder))
351  {
352  // Position of quadrature point in local coordinates of elements
353  auto iplocal_s = geo_in_inside.global(ip.position());
354  auto iplocal_n = geo_in_outside.global(ip.position());
355 
356  // Evaluate basis functions
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());
359 
360  // Evaluate u from inside and outside
361  u_s = 0.0;
362  for (size_type i=0; i<=dim; i++) // for all components
363  for (size_type k=0; k<dgspace_s.size(); k++) // for all basis functions
364  u_s[i] += x_s(lfsv_s.child(i),k)*phi_s[k];
365  u_n = 0.0;
366  for (size_type i=0; i<=dim; i++) // for all components
367  for (size_type k=0; k<dgspace_n.size(); k++) // for all basis functions
368  u_n[i] += x_n(lfsv_n.child(i),k)*phi_n[k];
369 
370  // Compute numerical flux at integration point
371  f = 0.0;
372  A_plus_s.umv(u_s,f);
373  // std::cout << " after A_plus*u_s " << f << std::endl;
374  A_minus_n.umv(u_n,f);
375  // std::cout << " after A_minus*u_n " << f << std::endl;
376 
377  // Integrate
378  auto factor = ip.weight() * geo.integrationElement(ip.position());
379  for (size_type k=0; k<dgspace_s.size(); k++) // loop over all vector-valued (!) basis functions (with identical components)
380  for (size_type i=0; i<=dim; i++) // loop over all components
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++) // loop over all vector-valued (!) basis functions (with identical components)
383  for (size_type i=0; i<=dim; i++) // loop over all components
384  r_n.accumulate(lfsv_n.child(i),k, - f[i]*phi_n[k]*factor);
385  }
386 
387  // std::cout << " residual_s: ";
388  // for (size_type i=0; i<r_s.size(); i++) std::cout << r_s[i] << " ";
389  // std::cout << std::endl;
390  // std::cout << " residual_n: ";
391  // for (size_type i=0; i<r_n.size(); i++) std::cout << r_n[i] << " ";
392  // std::cout << std::endl;
393 
394  }
395 
396  // Skeleton integral depending on test and ansatz functions
397  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
398  void alpha_boundary (const IG& ig,
399  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
400  R& r_s) const
401  {
402  // Get types
403  using namespace 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;
410 
411  // Get local function space that is identical for all components
412  const auto& dgspace_s = child(lfsv_s,_0);
413 
414  // Normal: assume faces are planar
415  const Dune::FieldVector<DF,dim> n_F = ig.centerUnitOuterNormal();
416 
417  // Reference to inside cell
418  const auto& cell_inside = ig.inside();
419 
420  // Get geometries
421  auto geo = ig.geometry();
422  auto geo_inside = cell_inside.geometry();
423 
424  // Get geometry of intersection in local coordinates of cell_inside
425  auto geo_in_inside = ig.geometryInInside();
426 
427  // Evaluate speed of sound (assumed constant per element)
428  auto ref_el_inside = referenceElement(geo_inside);
429  auto local_inside = ref_el_inside.position(0,0);
430  auto c_s = param.c(cell_inside,local_inside);
431 
432  // Compute A+ (outgoing waves)
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]; // row dim-1 corresponds to eigenvalue +c
437  Dune::FieldVector<DF,dim+1> unit(0.0);
438  unit[dim-1] = 1.0;
439  Dune::FieldVector<DF,dim+1> beta;
440  RT.solve(beta,unit);
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];
445 
446  // Compute A- (incoming waves)
448  for (int i=0; i<=dim; i++) alpha[i] = RT[dim][i]; // row dim corresponds to eigenvalue -c
449  unit = 0.0;
450  unit[dim] = 1.0;
451  RT.solve(beta,unit);
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];
456 
457  // Initialize vectors outside for loop
458  Dune::FieldVector<RF,dim+1> u_s(0.0);
459  Dune::FieldVector<RF,dim+1> f(0.0);
460 
461  // Loop over quadrature points
462  const int order_s = dgspace_s.finiteElement().localBasis().order();
463  const int intorder = overintegration+1+2*order_s;
464  for (const auto& ip : quadratureRule(geo,intorder))
465  {
466  // Position of quadrature point in local coordinates of elements
467  auto iplocal_s = geo_in_inside.global(ip.position());
468 
469  // Evaluate basis functions
470  auto& phi_s = cache[order_s].evaluateFunction(iplocal_s,dgspace_s.finiteElement().localBasis());
471 
472  // Evaluate u from inside and outside
473  u_s = 0.0;
474  for (size_type i=0; i<=dim; i++) // for all components
475  for (size_type k=0; k<dgspace_s.size(); k++) // for all basis functions
476  u_s[i] += x_s(lfsv_s.child(i),k)*phi_s[k];
477  // std::cout << " u_s " << u_s << std::endl;
478 
479  // Evaluate boundary condition
480  Dune::FieldVector<RF,dim+1> u_n(param.g(ig.intersection(),ip.position(),u_s));
481  // std::cout << " u_n " << u_n << " bc: " << param.g(ig.intersection(),ip.position(),u_s) << std::endl;
482 
483  // Compute numerical flux at integration point
484  f = 0.0;
485  A_plus_s.umv(u_s,f);
486  // std::cout << " after A_plus*u_s " << f << std::endl;
487  A_minus_n.umv(u_n,f);
488  // std::cout << " after A_minus*u_n " << f << std::endl;
489 
490  // Integrate
491  auto factor = ip.weight() * geo.integrationElement(ip.position());
492  for (size_type k=0; k<dgspace_s.size(); k++) // loop over all vector-valued (!) basis functions (with identical components)
493  for (size_type i=0; i<=dim; i++) // loop over all components
494  r_s.accumulate(lfsv_s.child(i),k, f[i]*phi_s[k]*factor);
495  }
496 
497  // std::cout << " residual_s: ";
498  // for (size_type i=0; i<r_s.size(); i++) std::cout << r_s[i] << " ";
499  // std::cout << std::endl;
500  }
501 
502  // Volume integral depending only on test functions
503  template<typename EG, typename LFSV, typename R>
504  void lambda_volume (const EG& eg, const LFSV& lfsv, R& r) const
505  {
506  // Get types
507  using namespace Indices;
508  using DGSpace = TypeTree::Child<LFSV,_0>;
509  using size_type = typename DGSpace::Traits::SizeType;
510 
511  // Get local function space that is identical for all components
512  const DGSpace& dgspace = child(lfsv,_0);
513 
514  // Reference to cell
515  const auto& cell = eg.entity();
516 
517  // Get geometries
518  auto geo = eg.geometry();
519 
520  // Loop over quadrature points
521  const int order_s = dgspace.finiteElement().localBasis().order();
522  const int intorder = overintegration+2*order_s;
523  for (const auto& ip : quadratureRule(geo,intorder))
524  {
525  // Evaluate right hand side
526  auto q(param.q(cell,ip.position()));
527 
528  // Evaluate basis functions
529  auto& phi = cache[order_s].evaluateFunction(ip.position(),dgspace.finiteElement().localBasis());
530 
531  // Integrate
532  auto factor = ip.weight() * geo.integrationElement(ip.position());
533  for (size_type k=0; k<=dim; k++) // for all components
534  for (size_type i=0; i<dgspace.size(); i++) // for all test functions of this component
535  r.accumulate(lfsv.child(k),i, - q[k]*phi[i]*factor);
536  }
537  }
538 
540  void setTime (typename T::Traits::RangeFieldType t)
541  {
542  }
543 
545  void preStep (typename T::Traits::RangeFieldType time, typename T::Traits::RangeFieldType dt,
546  int stages)
547  {
548  }
549 
551  void preStage (typename T::Traits::RangeFieldType time, int r)
552  {
553  }
554 
556  void postStage ()
557  {
558  }
559 
561  typename T::Traits::RangeFieldType suggestTimestep (typename T::Traits::RangeFieldType dt) const
562  {
563  return dt;
564  }
565 
566  private:
567  T& param;
568  int overintegration;
569  using LocalBasisType = typename FEM::Traits::FiniteElementType::Traits::LocalBasisType;
571  std::vector<Cache> cache;
572  };
573 
574 
575 
582  template<typename T, typename FEM>
584  public NumericalJacobianApplyVolume<DGLinearAcousticsTemporalOperator<T,FEM> >,
586  public InstationaryLocalOperatorDefaultMethods<typename T::Traits::RangeFieldType>
587  {
588  enum { dim = T::Traits::GridViewType::dimension };
589  public:
590  // pattern assembly flags
591  enum { doPatternVolume = true };
592 
593  // residual assembly flags
594  enum { doAlphaVolume = true };
595 
596  DGLinearAcousticsTemporalOperator (T& param_, int overintegration_=0)
597  : param(param_), overintegration(overintegration_), cache(20)
598  {}
599 
600  // define sparsity pattern of operator representation
601  template<typename LFSU, typename LFSV, typename LocalPattern>
602  void pattern_volume (const LFSU& lfsu, const LFSV& lfsv,
603  LocalPattern& pattern) const
604  {
605  // paranoia check number of number of components
606  if (TypeTree::degree(lfsu)!=TypeTree::degree(lfsv))
607  DUNE_THROW(Dune::Exception,"need U=V!");
608  if (TypeTree::degree(lfsv)!=dim+1)
609  DUNE_THROW(Dune::Exception,"need exactly dim+1 components!");
610 
611  for (size_t k=0; k<TypeTree::degree(lfsv); 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);
615  }
616 
617  // volume integral depending on test and ansatz functions
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
620  {
621  // get types
622  using namespace 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;
627 
628  // get local function space that is identical for all components
629  const auto& dgspace = child(lfsv,_0);
630 
631  // get geometry
632  auto geo = eg.geometry();
633 
634  // Initialize vectors outside for loop
635  Dune::FieldVector<RF,dim+1> u(0.0);
636 
637  // loop over quadrature points
638  const int order = dgspace.finiteElement().localBasis().order();
639  const int intorder = overintegration+2*order;
640  for (const auto& ip : quadratureRule(geo,intorder))
641  {
642  // evaluate basis functions
643  auto& phi = cache[order].evaluateFunction(ip.position(),dgspace.finiteElement().localBasis());
644 
645  // evaluate u
646  u = 0.0;
647  for (size_type k=0; k<=dim; k++) // for all components
648  for (size_type j=0; j<dgspace.size(); j++) // for all basis functions
649  u[k] += x(lfsv.child(k),j)*phi[j];
650 
651  // integrate
652  auto factor = ip.weight() * geo.integrationElement(ip.position());
653  for (size_type k=0; k<=dim; k++) // for all components
654  for (size_type i=0; i<dgspace.size(); i++) // for all test functions of this component
655  r.accumulate(lfsv.child(k),i, u[k]*phi[i]*factor);
656  }
657  }
658 
659  // jacobian of volume term
660  template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
661  void jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv,
662  M & mat) const
663  {
664  // get types
665  using namespace Indices;
666  using DGSpace = TypeTree::Child<LFSV,_0>;
667  using size_type = typename DGSpace::Traits::SizeType;
668 
669  // get local function space that is identical for all components
670  const auto& dgspace = child(lfsv,_0);
671 
672  // get geometry
673  auto geo = eg.geometry();
674 
675  // loop over quadrature points
676  const int order = dgspace.finiteElement().localBasis().order();
677  const int intorder = overintegration+2*order;
678  for (const auto& ip : quadratureRule(geo,intorder))
679  {
680  // evaluate basis functions
681  auto& phi = cache[order].evaluateFunction(ip.position(),dgspace.finiteElement().localBasis());
682 
683  // integrate
684  auto factor = ip.weight() * geo.integrationElement(ip.position());
685  for (size_type k=0; k<=dim; k++) // for all components
686  for (size_type i=0; i<dgspace.size(); i++) // for all test functions of this component
687  for (size_type j=0; j<dgspace.size(); j++) // for all ansatz functions of this component
688  mat.accumulate(lfsv.child(k),i,lfsu.child(k),j, phi[j]*phi[i]*factor);
689  }
690  }
691 
692  private:
693  T& param;
694  int overintegration;
695  using LocalBasisType = typename FEM::Traits::FiniteElementType::Traits::LocalBasisType;
697  std::vector<Cache> cache;
698  };
699 
700  }
701 }
702 
703 #endif // DUNE_PDELAB_LOCALOPERATOR_LINEARACOUSTICSDG_HH
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
const IG & ig
Definition: constraints.hh:148
static const int dim
Definition: adaptivity.hh:83
DGLinearAcousticsSpatialOperator(T &param_, int overintegration_=0)
Definition: linearacousticsdg.hh:189
void pattern_volume(const LFSU &lfsu, const LFSV &lfsv, LocalPattern &pattern) const
Definition: linearacousticsdg.hh:602
static void eigenvectors_transposed(T1 c, const Dune::FieldVector< T2, dim > &n, Dune::FieldMatrix< T3, dim+1, dim+1 > &RT)
Definition: linearacousticsdg.hh:45
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
static void eigenvectors(T1 c, const Dune::FieldVector< T2, dim > &n, Dune::FieldMatrix< T3, dim+1, dim+1 > &RT)
Definition: linearacousticsdg.hh:82
sparsity pattern generator
Definition: pattern.hh:29
const std::string s
Definition: function.hh:1101
Implement jacobian_skeleton() based on alpha_skeleton()
Definition: numericaljacobian.hh:156
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: linearacousticsdg.hh:619
Definition: linearacousticsdg.hh:28
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 preStage(typename T::Traits::RangeFieldType time, int r)
to be called once before each stage
Definition: linearacousticsdg.hh:551
Implements linear and nonlinear versions of jacobian_apply_boundary() based on alpha_boundary() ...
Definition: numericaljacobianapply.hh:285
Implements linear and nonlinear versions of jacobian_apply_volume() based on alpha_volume() ...
Definition: numericaljacobianapply.hh:32
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:52
Implements linear and nonlinear versions of jacobian_apply_skeleton() based on alpha_skeleton() ...
Definition: numericaljacobianapply.hh:180
store values of basis functions and gradients in a cache
Definition: localbasiscache.hh:17
void postStage()
to be called once at the end of each stage
Definition: linearacousticsdg.hh:556
static void eigenvectors_transposed(T1 c, const Dune::FieldVector< T2, dim > &n, Dune::FieldMatrix< T3, dim+1, dim+1 > &RT)
Definition: linearacousticsdg.hh:104
static void eigenvectors_transposed(T1 c, const Dune::FieldVector< T2, dim > &n, Dune::FieldMatrix< T3, dim+1, dim+1 > &RT)
Definition: linearacousticsdg.hh:74
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: linearacousticsdg.hh:196
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: linearacousticsdg.hh:504
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
sparsity pattern generator
Definition: pattern.hh:13
const Entity & e
Definition: localfunctionspace.hh:111
static void eigenvectors(T1 c, const Dune::FieldVector< T2, dim > &n, Dune::FieldMatrix< T3, dim+1, dim+1 > &RT)
Definition: linearacousticsdg.hh:125
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
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
DGLinearAcousticsTemporalOperator(T &param_, int overintegration_=0)
Definition: linearacousticsdg.hh:596
Implement jacobian_boundary() based on alpha_boundary()
Definition: numericaljacobian.hh:250
Definition: linearacousticsdg.hh:163
Implement jacobian_volume() based on alpha_volume()
Definition: numericaljacobian.hh:31
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: linearacousticsdg.hh:661
T::Traits::RangeFieldType suggestTimestep(typename T::Traits::RangeFieldType dt) const
to be called once before each stage
Definition: linearacousticsdg.hh:561
void setTime(typename T::Traits::RangeFieldType t)
set time in parameter class
Definition: linearacousticsdg.hh:540
ReferenceElementWrapper< ReferenceElement< typename Geometry::ctype, Geometry::mydimension > > referenceElement(const Geometry &geo)
Returns the reference element for the given geometry.
Definition: referenceelements.hh:144
Default flags for all local operators.
Definition: flags.hh:18