dune-pdelab  2.4-dev
linearacousticsdg.hh
Go to the documentation of this file.
1 // -*- tab-width: 8; indent-tabs-mode: nil -*-
2 #ifndef DUNE_PDELAB_LINEARACOUSTICSDG_HH
3 #define DUNE_PDELAB_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>
17 
19 
20 namespace Dune {
21  namespace PDELab {
22 
23 
24  template<int dim>
26  {};
27 
31  template<>
33  {
34  enum { dim = 1 };
35  public:
41  template<typename T1, typename T2, typename T3>
42  static void eigenvectors_transposed (T1 c, const Dune::FieldVector<T2,dim>& n, Dune::FieldMatrix<T3,dim+1,dim+1>& RT)
43  {
44  RT[0][0] = 1; RT[0][1] = c*n[0];
45  RT[1][0] = -1; RT[1][1] = c*n[0];
46  }
47 
48  template<typename T1, typename T2, typename T3>
49  static void eigenvectors (T1 c, const Dune::FieldVector<T2,dim>& n, Dune::FieldMatrix<T3,dim+1,dim+1>& RT)
50  {
51  RT[0][0] = 1; RT[1][0] = c*n[0];
52  RT[0][1] = -1; RT[1][1] = c*n[0];
53  }
54 
55  };
56 
60  template<>
62  {
63  enum { dim = 2 };
64  public:
70  template<typename T1, typename T2, typename T3>
71  static void eigenvectors_transposed (T1 c, const Dune::FieldVector<T2,dim>& n, Dune::FieldMatrix<T3,dim+1,dim+1>& RT)
72  {
73  RT[0][0] = 0; RT[0][1] = -n[1]; RT[0][2] = n[0];
74  RT[1][0] = 1; RT[1][1] = c*n[0]; RT[1][2] = c*n[1];
75  RT[2][0] = -1; RT[2][1] = c*n[0]; RT[2][2] = c*n[1];
76  }
77 
78  template<typename T1, typename T2, typename T3>
79  static void eigenvectors (T1 c, const Dune::FieldVector<T2,dim>& n, Dune::FieldMatrix<T3,dim+1,dim+1>& RT)
80  {
81  RT[0][0] = 0; RT[1][0] = -n[1]; RT[2][0] = n[0];
82  RT[0][1] = 1; RT[1][1] = c*n[0]; RT[2][1] = c*n[1];
83  RT[0][2] = -1; RT[1][2] = c*n[0]; RT[2][2] = c*n[1];
84  }
85  };
86 
90  template<>
92  {
93  enum { dim = 3 };
94  public:
100  template<typename T1, typename T2, typename T3>
101  static void eigenvectors_transposed (T1 c, const Dune::FieldVector<T2,dim>& n, Dune::FieldMatrix<T3,dim+1,dim+1>& RT)
102  {
103  Dune::FieldVector<T2,dim> s;
104  s[0] = n[1]-n[2]; s[1] = n[2]-n[0]; s[2] = n[0]-n[1];
105  if (s.two_norm()<1e-5)
106  {
107  s[0] = n[1]+n[2]; s[1] = n[2]-n[0]; s[2] = -(n[0]+n[1]);
108  }
109 
110  Dune::FieldVector<T2,dim> t; // compute cross product s * n
111  t[0] = n[1]*s[2] - n[2]*s[1];
112  t[1] = n[2]*s[0] - n[0]*s[2];
113  t[2] = n[0]*s[1] - n[1]*s[0];
114 
115  RT[0][0] = 0; RT[0][1] = s[0]; RT[0][2] = s[1]; RT[0][3] = s[2];
116  RT[1][0] = 0; RT[1][1] = t[0]; RT[1][2] = t[1]; RT[1][3] = t[2];
117  RT[2][0] = 1; RT[2][1] = c*n[0]; RT[2][2] = c*n[1]; RT[2][3] = c*n[2];
118  RT[3][0] = -1; RT[3][1] = c*n[0]; RT[3][2] = c*n[1]; RT[3][3] = c*n[2];
119  }
120 
121  template<typename T1, typename T2, typename T3>
122  static void eigenvectors (T1 c, const Dune::FieldVector<T2,dim>& n, Dune::FieldMatrix<T3,dim+1,dim+1>& RT)
123  {
124  Dune::FieldVector<T2,dim> s;
125  s[0] = n[1]-n[2]; s[1] = n[2]-n[0]; s[2] = n[0]-n[1];
126  if (s.two_norm()<1e-5)
127  {
128  s[0] = n[1]+n[2]; s[1] = n[2]-n[0]; s[2] = -(n[0]+n[1]);
129  }
130 
131  Dune::FieldVector<T2,dim> t; // compute cross product s * n
132  t[0] = n[1]*s[2] - n[2]*s[1];
133  t[1] = n[2]*s[0] - n[0]*s[2];
134  t[2] = n[0]*s[1] - n[1]*s[0];
135 
136  RT[0][0] = 0; RT[1][0] = s[0]; RT[2][0] = s[1]; RT[3][0] = s[2];
137  RT[0][1] = 0; RT[1][1] = t[0]; RT[2][1] = t[1]; RT[3][1] = t[2];
138  RT[0][2] = 1; RT[1][2] = c*n[0]; RT[2][2] = c*n[1]; RT[3][2] = c*n[2];
139  RT[0][3] = -1; RT[1][3] = c*n[0]; RT[2][3] = c*n[1]; RT[3][3] = c*n[2];
140  }
141  };
142 
159  template<typename T, typename FEM>
161  public NumericalJacobianApplyVolume<DGLinearAcousticsSpatialOperator<T,FEM> >,
162  public NumericalJacobianVolume<DGLinearAcousticsSpatialOperator<T,FEM> >,
163  public NumericalJacobianApplySkeleton<DGLinearAcousticsSpatialOperator<T,FEM> >,
164  public NumericalJacobianSkeleton<DGLinearAcousticsSpatialOperator<T,FEM> >,
165  public NumericalJacobianApplyBoundary<DGLinearAcousticsSpatialOperator<T,FEM> >,
166  public NumericalJacobianBoundary<DGLinearAcousticsSpatialOperator<T,FEM> >,
167  public FullSkeletonPattern,
168  public FullVolumePattern,
170  public InstationaryLocalOperatorDefaultMethods<typename T::Traits::RangeFieldType>
171  {
172  enum { dim = T::Traits::GridViewType::dimension };
173 
174  public:
175  // pattern assembly flags
176  enum { doPatternVolume = true };
177  enum { doPatternSkeleton = true };
178 
179  // residual assembly flags
180  enum { doAlphaVolume = true };
181  enum { doAlphaSkeleton = true };
182  enum { doAlphaBoundary = true };
183  enum { doLambdaVolume = true };
184 
185  // ! constructor
186  DGLinearAcousticsSpatialOperator (T& param_, int overintegration_=0)
187  : param(param_), overintegration(overintegration_), cache(20)
188  {
189  }
190 
191  // volume integral depending on test and ansatz functions
192  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
193  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
194  {
195  // get types
196  typedef typename LFSV::template Child<0>::Type DGSpace;
197  typedef typename DGSpace::Traits::FiniteElementType::
198  Traits::LocalBasisType::Traits::DomainFieldType DF;
199  typedef typename DGSpace::Traits::FiniteElementType::
200  Traits::LocalBasisType::Traits::RangeFieldType RF;
201  typedef typename DGSpace::Traits::FiniteElementType::
202  Traits::LocalBasisType::Traits::RangeType RangeType;
203  typedef typename DGSpace::Traits::FiniteElementType::
204  Traits::LocalBasisType::Traits::JacobianType JacobianType;
205  typedef typename DGSpace::Traits::SizeType size_type;
206 
207  // paranoia check number of number of components
208  if (LFSV::CHILDREN!=dim+1)
209  DUNE_THROW(Dune::Exception,"need exactly dim+1 components!");
210 
211  // get local function space that is identical for all components
212  const DGSpace& dgspace = lfsv.template child<0>();
213 
214  auto geometry = eg.geometry();
215 
216  // select quadrature rule
217  const int order = dgspace.finiteElement().localBasis().order();
218  const int intorder = overintegration+2*order;
219  Dune::GeometryType gt = geometry.type();
220 
221  // transformation
222  typename EG::Geometry::JacobianInverseTransposed jac;
223 
224  // evaluate speed of sound (assumed constant per element)
225  Dune::FieldVector<DF,dim> localcenter = Dune::ReferenceElements<DF,dim>::general(gt).position(0,0);
226  RF c2 = param.c(eg.entity(),localcenter);
227  c2 = c2*c2; // square it
228 
229  // std::cout << "alpha_volume center=" << geometry.center() << std::endl;
230 
231  // loop over quadrature points
232  for (const auto& ip : Dune::QuadratureRules<DF,dim>::rule(gt,intorder))
233  {
234  // evaluate basis functions
235  const std::vector<RangeType>& phi = cache[order].evaluateFunction(ip.position(),dgspace.finiteElement().localBasis());
236 
237  // evaluate u
238  Dune::FieldVector<RF,dim+1> u(0.0);
239  for (size_type k=0; k<=dim; k++) // for all components
240  for (size_type j=0; j<dgspace.size(); j++) // for all basis functions
241  u[k] += x(lfsv.child(k),j)*phi[j];
242  // std::cout << " u at " << ip.position() << " : " << u << std::endl;
243 
244  // evaluate gradient of basis functions (we assume Galerkin method lfsu=lfsv)
245  const std::vector<JacobianType>& js = cache[order].evaluateJacobian(ip.position(),dgspace.finiteElement().localBasis());
246 
247  // compute global gradients
248  jac = geometry.jacobianInverseTransposed(ip.position());
249  std::vector<Dune::FieldVector<RF,dim> > gradphi(dgspace.size());
250  for (size_type i=0; i<dgspace.size(); i++)
251  jac.mv(js[i][0],gradphi[i]);
252 
253  // integrate
254  RF factor = ip.weight() * geometry.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  typedef typename LFSV::template Child<0>::Type DGSpace;
280  typedef typename DGSpace::Traits::FiniteElementType::
281  Traits::LocalBasisType::Traits::DomainFieldType DF;
282  typedef typename DGSpace::Traits::FiniteElementType::
283  Traits::LocalBasisType::Traits::RangeFieldType RF;
284  typedef typename DGSpace::Traits::FiniteElementType::
285  Traits::LocalBasisType::Traits::RangeType RangeType;
286  typedef typename DGSpace::Traits::SizeType size_type;
287 
288  // get local function space that is identical for all components
289  const DGSpace& dgspace_s = lfsv_s.template child<0>();
290  const DGSpace& dgspace_n = lfsv_n.template child<0>();
291 
292  // normal: assume faces are planar
293  const Dune::FieldVector<DF,dim> n_F = ig.centerUnitOuterNormal();
294 
295  // store inside and outside cells
296  auto inside_cell = ig.inside();
297  auto outside_cell = ig.outside();
298 
299  // evaluate speed of sound (assumed constant per element)
300  const Dune::FieldVector<DF,dim>&
301  inside_local = Dune::ReferenceElements<DF,dim>::general(inside_cell.type()).position(0,0);
302  const Dune::FieldVector<DF,dim>&
303  outside_local = Dune::ReferenceElements<DF,dim>::general(outside_cell.type()).position(0,0);
304  RF c_s = param.c(inside_cell,inside_local);
305  RF c_n = param.c(outside_cell,outside_local);
306 
307  // compute A+ (outgoing waves)
308  Dune::FieldMatrix<DF,dim+1,dim+1> RT;
310  Dune::FieldVector<DF,dim+1> alpha;
311  for (int i=0; i<=dim; i++) alpha[i] = RT[dim-1][i]; // row dim-1 corresponds to eigenvalue +c
312  Dune::FieldVector<DF,dim+1> unit(0.0);
313  unit[dim-1] = 1.0;
314  Dune::FieldVector<DF,dim+1> beta;
315  RT.solve(beta,unit);
316  Dune::FieldMatrix<DF,dim+1,dim+1> A_plus_s;
317  for (int i=0; i<=dim; i++)
318  for (int j=0; j<=dim; j++)
319  A_plus_s[i][j] = c_s*alpha[i]*beta[j];
320 
321  // compute A- (incoming waves)
323  for (int i=0; i<=dim; i++) alpha[i] = RT[dim][i]; // row dim corresponds to eigenvalue -c
324  unit = 0.0;
325  unit[dim] = 1.0;
326  RT.solve(beta,unit);
327  Dune::FieldMatrix<DF,dim+1,dim+1> A_minus_n;
328  for (int i=0; i<=dim; i++)
329  for (int j=0; j<=dim; j++)
330  A_minus_n[i][j] = -c_n*alpha[i]*beta[j];
331 
332  auto geometry = ig.geometry();
333 
334  // select quadrature rule
335  const int order_s = dgspace_s.finiteElement().localBasis().order();
336  const int order_n = dgspace_n.finiteElement().localBasis().order();
337  const int intorder = overintegration+1+2*std::max(order_s,order_n);
338  Dune::GeometryType gtface = geometry.type();
339 
340  // std::cout << "alpha_skeleton center=" << geometry.center() << std::endl;
341 
342  // loop over quadrature points
343  for (const auto& ip : Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder))
344  {
345  // position of quadrature point in local coordinates of elements
346  Dune::FieldVector<DF,dim> iplocal_s = ig.geometryInInside().global(ip.position());
347  Dune::FieldVector<DF,dim> iplocal_n = ig.geometryInOutside().global(ip.position());
348 
349  // evaluate basis functions
350  const std::vector<RangeType>& phi_s = cache[order_s].evaluateFunction(iplocal_s,dgspace_s.finiteElement().localBasis());
351  const std::vector<RangeType>& phi_n = cache[order_n].evaluateFunction(iplocal_n,dgspace_n.finiteElement().localBasis());
352 
353  // evaluate u from inside and outside
354  Dune::FieldVector<RF,dim+1> u_s(0.0);
355  for (size_type i=0; i<=dim; i++) // for all components
356  for (size_type k=0; k<dgspace_s.size(); k++) // for all basis functions
357  u_s[i] += x_s(lfsv_s.child(i),k)*phi_s[k];
358  Dune::FieldVector<RF,dim+1> u_n(0.0);
359  for (size_type i=0; i<=dim; i++) // for all components
360  for (size_type k=0; k<dgspace_n.size(); k++) // for all basis functions
361  u_n[i] += x_n(lfsv_n.child(i),k)*phi_n[k];
362 
363  // compute numerical flux at integration point
364  Dune::FieldVector<RF,dim+1> f(0.0);
365  A_plus_s.umv(u_s,f);
366  // std::cout << " after A_plus*u_s " << f << std::endl;
367  A_minus_n.umv(u_n,f);
368  // std::cout << " after A_minus*u_n " << f << std::endl;
369 
370  // integrate
371  RF factor = ip.weight() * geometry.integrationElement(ip.position());
372  for (size_type k=0; k<dgspace_s.size(); k++) // loop over all vector-valued (!) basis functions (with identical components)
373  for (size_type i=0; i<=dim; i++) // loop over all components
374  r_s.accumulate(lfsv_s.child(i),k, f[i]*phi_s[k]*factor);
375  for (size_type k=0; k<dgspace_n.size(); k++) // loop over all vector-valued (!) basis functions (with identical components)
376  for (size_type i=0; i<=dim; i++) // loop over all components
377  r_n.accumulate(lfsv_n.child(i),k, - f[i]*phi_n[k]*factor);
378  }
379 
380  // std::cout << " residual_s: ";
381  // for (size_type i=0; i<r_s.size(); i++) std::cout << r_s[i] << " ";
382  // std::cout << std::endl;
383  // std::cout << " residual_n: ";
384  // for (size_type i=0; i<r_n.size(); i++) std::cout << r_n[i] << " ";
385  // std::cout << std::endl;
386 
387  }
388 
389  // skeleton integral depending on test and ansatz functions
390  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
391  void alpha_boundary (const IG& ig,
392  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
393  R& r_s) const
394  {
395  // get types
396  typedef typename LFSV::template Child<0>::Type DGSpace;
397  typedef typename DGSpace::Traits::FiniteElementType::
398  Traits::LocalBasisType::Traits::DomainFieldType DF;
399  typedef typename DGSpace::Traits::FiniteElementType::
400  Traits::LocalBasisType::Traits::RangeFieldType RF;
401  typedef typename DGSpace::Traits::FiniteElementType::
402  Traits::LocalBasisType::Traits::RangeType RangeType;
403  typedef typename DGSpace::Traits::SizeType size_type;
404 
405  // get local function space that is identical for all components
406  const DGSpace& dgspace_s = lfsv_s.template child<0>();
407 
408  // normal: assume faces are planar
409  const Dune::FieldVector<DF,dim> n_F = ig.centerUnitOuterNormal();
410 
411  auto cell = ig.inside();
412 
413  // evaluate speed of sound (assumed constant per element)
414  const Dune::FieldVector<DF,dim>&
415  inside_local = Dune::ReferenceElements<DF,dim>::general(cell.type()).position(0,0);
416  RF c_s = param.c(cell,inside_local);
417 
418  // compute A+ (outgoing waves)
419  Dune::FieldMatrix<DF,dim+1,dim+1> RT;
421  Dune::FieldVector<DF,dim+1> alpha;
422  for (int i=0; i<=dim; i++) alpha[i] = RT[dim-1][i]; // row dim-1 corresponds to eigenvalue +c
423  Dune::FieldVector<DF,dim+1> unit(0.0);
424  unit[dim-1] = 1.0;
425  Dune::FieldVector<DF,dim+1> beta;
426  RT.solve(beta,unit);
427  Dune::FieldMatrix<DF,dim+1,dim+1> A_plus_s;
428  for (int i=0; i<=dim; i++)
429  for (int j=0; j<=dim; j++)
430  A_plus_s[i][j] = c_s*alpha[i]*beta[j];
431 
432  // compute A- (incoming waves)
434  for (int i=0; i<=dim; i++) alpha[i] = RT[dim][i]; // row dim corresponds to eigenvalue -c
435  unit = 0.0;
436  unit[dim] = 1.0;
437  RT.solve(beta,unit);
438  Dune::FieldMatrix<DF,dim+1,dim+1> A_minus_n;
439  for (int i=0; i<=dim; i++)
440  for (int j=0; j<=dim; j++)
441  A_minus_n[i][j] = -c_s*alpha[i]*beta[j];
442 
443  auto geometry = ig.geometry();
444 
445  // select quadrature rule
446  const int order_s = dgspace_s.finiteElement().localBasis().order();
447  const int intorder = overintegration+1+2*order_s;
448  Dune::GeometryType gtface = geometry.type();
449 
450  // std::cout << "alpha_boundary center=" << geometry.center() << std::endl;
451 
452  // loop over quadrature points
453  for (const auto& ip : Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder))
454  {
455  // position of quadrature point in local coordinates of elements
456  Dune::FieldVector<DF,dim> iplocal_s = ig.geometryInInside().global(ip.position());
457 
458  // evaluate basis functions
459  const std::vector<RangeType>& phi_s = cache[order_s].evaluateFunction(iplocal_s,dgspace_s.finiteElement().localBasis());
460 
461  // evaluate u from inside and outside
462  Dune::FieldVector<RF,dim+1> u_s(0.0);
463  for (size_type i=0; i<=dim; i++) // for all components
464  for (size_type k=0; k<dgspace_s.size(); k++) // for all basis functions
465  u_s[i] += x_s(lfsv_s.child(i),k)*phi_s[k];
466  // std::cout << " u_s " << u_s << std::endl;
467 
468  // evaluate boundary condition
469  Dune::FieldVector<RF,dim+1> u_n(param.g(ig.intersection(),ip.position(),u_s));
470  // std::cout << " u_n " << u_n << " bc: " << param.g(ig.intersection(),ip.position(),u_s) << std::endl;
471 
472  // compute numerical flux at integration point
473  Dune::FieldVector<RF,dim+1> f(0.0);
474  A_plus_s.umv(u_s,f);
475  // std::cout << " after A_plus*u_s " << f << std::endl;
476  A_minus_n.umv(u_n,f);
477  // std::cout << " after A_minus*u_n " << f << std::endl;
478 
479  // integrate
480  RF factor = ip.weight() * geometry.integrationElement(ip.position());
481  for (size_type k=0; k<dgspace_s.size(); k++) // loop over all vector-valued (!) basis functions (with identical components)
482  for (size_type i=0; i<=dim; i++) // loop over all components
483  r_s.accumulate(lfsv_s.child(i),k, f[i]*phi_s[k]*factor);
484  }
485 
486  // std::cout << " residual_s: ";
487  // for (size_type i=0; i<r_s.size(); i++) std::cout << r_s[i] << " ";
488  // std::cout << std::endl;
489  }
490 
491  // volume integral depending only on test functions
492  template<typename EG, typename LFSV, typename R>
493  void lambda_volume (const EG& eg, const LFSV& lfsv, R& r) const
494  {
495  // get types
496  typedef typename LFSV::template Child<0>::Type DGSpace;
497  typedef typename DGSpace::Traits::FiniteElementType::
498  Traits::LocalBasisType::Traits::DomainFieldType DF;
499  typedef typename DGSpace::Traits::FiniteElementType::
500  Traits::LocalBasisType::Traits::RangeFieldType RF;
501  typedef typename DGSpace::Traits::FiniteElementType::
502  Traits::LocalBasisType::Traits::RangeType RangeType;
503  typedef typename DGSpace::Traits::SizeType size_type;
504 
505  // get local function space that is identical for all components
506  const DGSpace& dgspace = lfsv.template child<0>();
507 
508  auto geometry = eg.geometry();
509 
510  // select quadrature rule
511  const int order_s = dgspace.finiteElement().localBasis().order();
512  const int intorder = overintegration+2*order_s;
513  Dune::GeometryType gt = geometry.type();
514 
515  // loop over quadrature points
516  for (const auto& ip : Dune::QuadratureRules<DF,dim>::rule(gt,intorder))
517  {
518  // evaluate right hand side
519  Dune::FieldVector<RF,dim+1> q(param.q(eg.entity(),ip.position()));
520 
521  // evaluate basis functions
522  const std::vector<RangeType>& phi = cache[order_s].evaluateFunction(ip.position(),dgspace.finiteElement().localBasis());
523 
524  // integrate
525  RF factor = ip.weight() * geometry.integrationElement(ip.position());
526  for (size_type k=0; k<=dim; k++) // for all components
527  for (size_type i=0; i<dgspace.size(); i++) // for all test functions of this component
528  r.accumulate(lfsv.child(k),i, - q[k]*phi[i]*factor);
529  }
530  }
531 
533  void setTime (typename T::Traits::RangeFieldType t)
534  {
535  }
536 
538  void preStep (typename T::Traits::RangeFieldType time, typename T::Traits::RangeFieldType dt,
539  int stages)
540  {
541  }
542 
544  void preStage (typename T::Traits::RangeFieldType time, int r)
545  {
546  }
547 
549  void postStage ()
550  {
551  }
552 
554  typename T::Traits::RangeFieldType suggestTimestep (typename T::Traits::RangeFieldType dt) const
555  {
556  return dt;
557  }
558 
559  private:
560  T& param;
561  int overintegration;
562  typedef typename FEM::Traits::FiniteElementType::Traits::LocalBasisType LocalBasisType;
564  std::vector<Cache> cache;
565  };
566 
567 
568 
575  template<typename T, typename FEM>
577  public NumericalJacobianApplyVolume<DGLinearAcousticsTemporalOperator<T,FEM> >,
579  public InstationaryLocalOperatorDefaultMethods<typename T::Traits::RangeFieldType>
580  {
581  enum { dim = T::Traits::GridViewType::dimension };
582  public:
583  // pattern assembly flags
584  enum { doPatternVolume = true };
585 
586  // residual assembly flags
587  enum { doAlphaVolume = true };
588 
589  DGLinearAcousticsTemporalOperator (T& param_, int overintegration_=0)
590  : param(param_), overintegration(overintegration_), cache(20)
591  {}
592 
593  // define sparsity pattern of operator representation
594  template<typename LFSU, typename LFSV, typename LocalPattern>
595  void pattern_volume (const LFSU& lfsu, const LFSV& lfsv,
596  LocalPattern& pattern) const
597  {
598  // paranoia check number of number of components
599  if (LFSU::CHILDREN!=LFSV::CHILDREN)
600  DUNE_THROW(Dune::Exception,"need U=V!");
601  if (LFSV::CHILDREN!=dim+1)
602  DUNE_THROW(Dune::Exception,"need exactly dim+1 components!");
603 
604  for (size_t k=0; k<LFSV::CHILDREN; k++)
605  for (size_t i=0; i<lfsv.child(k).size(); ++i)
606  for (size_t j=0; j<lfsu.child(k).size(); ++j)
607  pattern.addLink(lfsv.child(k),i,lfsu.child(k),j);
608  }
609 
610  // volume integral depending on test and ansatz functions
611  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
612  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
613  {
614  // get types
615  typedef typename LFSV::template Child<0>::Type DGSpace;
616  typedef typename DGSpace::Traits::FiniteElementType::
617  Traits::LocalBasisType::Traits::DomainFieldType DF;
618  typedef typename DGSpace::Traits::FiniteElementType::
619  Traits::LocalBasisType::Traits::RangeFieldType RF;
620  typedef typename DGSpace::Traits::FiniteElementType::
621  Traits::LocalBasisType::Traits::RangeType RangeType;
622  typedef typename DGSpace::Traits::SizeType size_type;
623 
624  // get local function space that is identical for all components
625  const DGSpace& dgspace = lfsv.template child<0>();
626 
627  auto geometry = eg.geometry();
628 
629  // select quadrature rule
630  const int order = dgspace.finiteElement().localBasis().order();
631  const int intorder = overintegration+2*order;
632  Dune::GeometryType gt = geometry.type();
633 
634  // loop over quadrature points
635  for (const auto& ip : Dune::QuadratureRules<DF,dim>::rule(gt,intorder))
636  {
637  // evaluate basis functions
638  const std::vector<RangeType>& phi = cache[order].evaluateFunction(ip.position(),dgspace.finiteElement().localBasis());
639 
640  // evaluate u
641  Dune::FieldVector<RF,dim+1> u(0.0);
642  for (size_type k=0; k<=dim; k++) // for all components
643  for (size_type j=0; j<dgspace.size(); j++) // for all basis functions
644  u[k] += x(lfsv.child(k),j)*phi[j];
645 
646  // integrate
647  RF factor = ip.weight() * geometry.integrationElement(ip.position());
648  for (size_type k=0; k<=dim; k++) // for all components
649  for (size_type i=0; i<dgspace.size(); i++) // for all test functions of this component
650  r.accumulate(lfsv.child(k),i, u[k]*phi[i]*factor);
651  }
652  }
653 
654  // jacobian of volume term
655  template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
656  void jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv,
657  M & mat) const
658  {
659  // get types
660  typedef typename LFSV::template Child<0>::Type DGSpace;
661  typedef typename DGSpace::Traits::FiniteElementType::
662  Traits::LocalBasisType::Traits::DomainFieldType DF;
663  typedef typename DGSpace::Traits::FiniteElementType::
664  Traits::LocalBasisType::Traits::RangeFieldType RF;
665  typedef typename DGSpace::Traits::FiniteElementType::
666  Traits::LocalBasisType::Traits::RangeType RangeType;
667  typedef typename DGSpace::Traits::SizeType size_type;
668 
669  // get local function space that is identical for all components
670  const DGSpace& dgspace = lfsv.template child<0>();
671 
672  auto geometry = eg.geometry();
673 
674  // select quadrature rule
675  const int order = dgspace.finiteElement().localBasis().order();
676  const int intorder = overintegration+2*order;
677  Dune::GeometryType gt = geometry.type();
678 
679  // loop over quadrature points
680  for (const auto& ip : Dune::QuadratureRules<DF,dim>::rule(gt,intorder))
681  {
682  // evaluate basis functions
683  const std::vector<RangeType>& phi = cache[order].evaluateFunction(ip.position(),dgspace.finiteElement().localBasis());
684 
685  // integrate
686  RF factor = ip.weight() * geometry.integrationElement(ip.position());
687  for (size_type k=0; k<=dim; k++) // for all components
688  for (size_type i=0; i<dgspace.size(); i++) // for all test functions of this component
689  for (size_type j=0; j<dgspace.size(); j++) // for all ansatz functions of this component
690  mat.accumulate(lfsv.child(k),i,lfsu.child(k),j, phi[j]*phi[i]*factor);
691  }
692  }
693 
694  private:
695  T& param;
696  int overintegration;
697  typedef typename FEM::Traits::FiniteElementType::Traits::LocalBasisType LocalBasisType;
699  std::vector<Cache> cache;
700  };
701 
702  }
703 }
704 
705 #endif
DGLinearAcousticsSpatialOperator(T &param_, int overintegration_=0)
Definition: linearacousticsdg.hh:186
static void eigenvectors_transposed(T1 c, const Dune::FieldVector< T2, dim > &n, Dune::FieldMatrix< T3, dim+1, dim+1 > &RT)
Definition: linearacousticsdg.hh:42
void setTime(typename T::Traits::RangeFieldType t)
set time in parameter class
Definition: linearacousticsdg.hh:533
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: linearacousticsdg.hh:493
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: linearacousticsdg.hh:193
void postStage()
to be called once at the end of each stage
Definition: linearacousticsdg.hh:549
Definition: linearacousticsdg.hh:160
static void eigenvectors_transposed(T1 c, const Dune::FieldVector< T2, dim > &n, Dune::FieldMatrix< T3, dim+1, dim+1 > &RT)
Definition: linearacousticsdg.hh:101
const E & e
Definition: interpolate.hh:172
Implement jacobian_apply_volume() based on alpha_volume()
Definition: defaultimp.hh:321
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
const IG & ig
Definition: constraints.hh:147
static void eigenvectors(T1 c, const Dune::FieldVector< T2, dim > &n, Dune::FieldMatrix< T3, dim+1, dim+1 > &RT)
Definition: linearacousticsdg.hh:79
Default flags for all local operators.
Definition: flags.hh:18
sparsity pattern generator
Definition: pattern.hh:13
static void eigenvectors(T1 c, const Dune::FieldVector< T2, dim > &n, Dune::FieldMatrix< T3, dim+1, dim+1 > &RT)
Definition: linearacousticsdg.hh:122
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: linearacousticsdg.hh:656
static void eigenvectors_transposed(T1 c, const Dune::FieldVector< T2, dim > &n, Dune::FieldMatrix< T3, dim+1, dim+1 > &RT)
Definition: linearacousticsdg.hh:71
Implement jacobian_apply_skeleton() based on alpha_skeleton()
Definition: defaultimp.hh:445
Definition: adaptivity.hh:27
Implement jacobian_volume() based on alpha_volume()
Definition: defaultimp.hh:32
Implement jacobian_boundary() based on alpha_boundary()
Definition: defaultimp.hh:251
static const int dim
Definition: adaptivity.hh:83
T::Traits::RangeFieldType suggestTimestep(typename T::Traits::RangeFieldType dt) const
to be called once before each stage
Definition: linearacousticsdg.hh:554
const std::string s
Definition: function.hh:1102
Definition: linearacousticsdg.hh:25
sparsity pattern generator
Definition: pattern.hh:29
void pattern_volume(const LFSU &lfsu, const LFSV &lfsv, LocalPattern &pattern) const
Definition: linearacousticsdg.hh:595
Definition: linearacousticsdg.hh:576
void preStage(typename T::Traits::RangeFieldType time, int r)
to be called once before each stage
Definition: linearacousticsdg.hh:544
DGLinearAcousticsTemporalOperator(T &param_, int overintegration_=0)
Definition: linearacousticsdg.hh:589
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:538
store values of basis functions and gradients in a cache
Definition: localbasiscache.hh:15
Implement jacobian_skeleton() based on alpha_skeleton()
Definition: defaultimp.hh:157
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:391
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: linearacousticsdg.hh:612
const EG & eg
Definition: constraints.hh:280
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
Implement jacobian_apply_boundary() based on alpha_boundary()
Definition: defaultimp.hh:538
static void eigenvectors(T1 c, const Dune::FieldVector< T2, dim > &n, Dune::FieldMatrix< T3, dim+1, dim+1 > &RT)
Definition: linearacousticsdg.hh:49