DUNE PDELab (git)

linearacousticsdg.hh
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
11#include<dune/pdelab/common/quadraturerules.hh>
12#include<dune/pdelab/common/geometrywrapper.hh>
13#include<dune/pdelab/common/function.hh>
14#include<dune/pdelab/localoperator/pattern.hh>
15#include<dune/pdelab/localoperator/flags.hh>
16#include<dune/pdelab/localoperator/idefault.hh>
17#include<dune/pdelab/localoperator/defaultimp.hh>
18#include<dune/pdelab/finiteelement/localbasiscache.hh>
19
20#include"linearacousticsparameter.hh"
21
22namespace Dune {
23 namespace PDELab {
24
25
26 template<int dim>
27 class LinearAcousticsEigenvectors
28 {};
29
33 template<>
34 class LinearAcousticsEigenvectors<1>
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<>
63 class LinearAcousticsEigenvectors<2>
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<>
93 class LinearAcousticsEigenvectors<3>
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>
162 class DGLinearAcousticsSpatialOperator :
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,
171 public LocalOperatorDefaultFlags,
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;
317 LinearAcousticsEigenvectors<dim>::eigenvectors_transposed(c_s,n_F,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)
330 LinearAcousticsEigenvectors<dim>::eigenvectors_transposed(c_n,n_F,RT);
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;
433 LinearAcousticsEigenvectors<dim>::eigenvectors_transposed(c_s,n_F,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)
446 LinearAcousticsEigenvectors<dim>::eigenvectors_transposed(c_s,n_F,RT);
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;
569 using Cache = Dune::PDELab::LocalBasisCache<LocalBasisType>;
570 std::vector<Cache> cache;
571 };
572
573
574
581 template<typename T, typename FEM>
582 class DGLinearAcousticsTemporalOperator :
583 public NumericalJacobianApplyVolume<DGLinearAcousticsTemporalOperator<T,FEM> >,
584 public LocalOperatorDefaultFlags,
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;
695 using Cache = Dune::PDELab::LocalBasisCache<LocalBasisType>;
696 std::vector<Cache> cache;
697 };
698
699 }
700}
701
702#endif // DUNE_PDELAB_LOCALOPERATOR_LINEARACOUSTICSDG_HH
Dune::DenseMatrix::solve
void solve(V1 &x, const V2 &b, bool doPivoting=true) const
Solve system A x = b.
Dune::DenseMatrix::umv
void umv(const X &x, Y &y) const
y += A x
Definition: densematrix.hh:406
Dune::Exception
Base class for Dune-Exceptions.
Definition: exceptions.hh:96
Dune::FieldMatrix
A dense n x m matrix.
Definition: fmatrix.hh:117
Dune::FieldVector
vector space out of a tensor product of fields.
Definition: fvector.hh:91
Dune::PDELab::DGLinearAcousticsSpatialOperator
Definition: linearacousticsdg.hh:173
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::DGLinearAcousticsSpatialOperator::preStage
void preStage(typename T::Traits::RangeFieldType time, int r)
to be called once before each stage
Definition: linearacousticsdg.hh:550
Dune::PDELab::DGLinearAcousticsSpatialOperator::setTime
void setTime(typename T::Traits::RangeFieldType t)
set time in parameter class
Definition: linearacousticsdg.hh:539
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::DGLinearAcousticsSpatialOperator::postStage
void postStage()
to be called once at the end of each stage
Definition: linearacousticsdg.hh:555
Dune::PDELab::DGLinearAcousticsTemporalOperator
Definition: linearacousticsdg.hh:586
Dune::PDELab::FullSkeletonPattern
sparsity pattern generator
Definition: pattern.hh:30
Dune::PDELab::FullVolumePattern
sparsity pattern generator
Definition: pattern.hh:14
Dune::PDELab::InstationaryLocalOperatorDefaultMethods
Default class for additional methods in instationary local operators.
Definition: idefault.hh:90
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
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::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
Dune::PDELab::LocalBasisCache
store values of basis functions and gradients in a cache
Definition: localbasiscache.hh:18
Dune::PDELab::LocalOperatorDefaultFlags
Default flags for all local operators.
Definition: flags.hh:19
Dune::PDELab::NumericalJacobianApplyBoundary
Implements linear and nonlinear versions of jacobian_apply_boundary() based on alpha_boundary()
Definition: numericaljacobianapply.hh:287
Dune::PDELab::NumericalJacobianApplySkeleton
Implements linear and nonlinear versions of jacobian_apply_skeleton() based on alpha_skeleton()
Definition: numericaljacobianapply.hh:182
Dune::PDELab::NumericalJacobianApplyVolume
Implements linear and nonlinear versions of jacobian_apply_volume() based on alpha_volume()
Definition: numericaljacobianapply.hh:34
Dune::PDELab::NumericalJacobianBoundary
Implement jacobian_boundary() based on alpha_boundary()
Definition: numericaljacobian.hh:251
Dune::PDELab::NumericalJacobianSkeleton
Implement jacobian_skeleton() based on alpha_skeleton()
Definition: numericaljacobian.hh:157
Dune::PDELab::NumericalJacobianVolume
Implement jacobian_volume() based on alpha_volume()
Definition: numericaljacobian.hh:32
exceptions.hh
A few common exception classes.
fvector.hh
Implements a vector constructed from a given type representing a field and a compile-time given size.
Dune::Indices::_0
constexpr index_constant< 0 > _0
Compile time index with value 0.
Definition: indices.hh:52
DUNE_THROW
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
Dune::referenceElement
unspecified value type referenceElement(T &&... t)
Returns a reference element for the objects t....
Dune::Hybrid::max
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
Dune::TypeTree::degree
std::size_t degree(const Node &node)
Returns the degree of node as run time information.
Definition: nodeinterface.hh:79
Dune::TypeTree::Child
typename impl::_Child< Node, indices... >::type Child
Template alias for the type of a child node given by a list of child indices.
Definition: childextraction.hh:225
Dune::TypeTree::child
ImplementationDefined child(Node &&node, Indices... indices)
Extracts the child of a node given by a sequence of compile-time and run-time indices.
Definition: childextraction.hh:128
Dune::Indices
Namespace with predefined compile time indices for the range [0,19].
Definition: indices.hh:50
Dune
Dune namespace.
Definition: alignedallocator.hh:13
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 23, 23:29, 2024)