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