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