euler_elements.h
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Version 1.0; svn revision $LastChangedRevision: 1182 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-05-20 16:50:20 +0100 (Fri, 20 May 2016) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 //Header file for scalar advection elements
31 
32 #ifndef OOMPH_EULER_ELEMENTS_HEADER
33 #define OOMPH_EULER_ELEMENTS_HEADER
34 
35 // Config header generated by autoconfig
36 #ifdef HAVE_CONFIG_H
37 #include <oomph-lib-config.h>
38 #endif
39 
41 #include "../generic/Qspectral_elements.h"
42 #include "../generic/dg_elements.h"
43 
44 namespace oomph
45 {
46 
47 //==============================================================
48 ///Base class for Euler equations
49 //=============================================================
50 template<unsigned DIM>
52 {
53  //\short Default value for gamma (initialised to 1.4)
54  static double Default_Gamma_Value;
55 
56  //Pointer to the specific gas constant
57  double *Gamma_pt;
58 
59  /// \short Storage for the average primitive values
61 
62  /// \short Storage for the approximated gradients
64 
65 protected:
66 
67  /// \short DIM momentum-components, a density and an energy are transported
68  inline unsigned nflux() const {return DIM+2;}
69 
70  /// \short Return the flux as a function of the unknown
71  void flux(const Vector<double> &u, DenseMatrix<double> &f);
72 
73  /// \short Return the flux derivatives as a function of the unknowns
74  //void dflux_du(const Vector<double> &u, RankThreeTensor<double> &df_du);
75 
76 public:
77 
78  ///Constructor
81  {
82  //Set the default value of gamma
84  }
85 
86  ///Destructor
87  virtual ~EulerEquations()
88  {
89  if(this->Average_prim_value!=0)
91  if(this->Average_gradient!=0)
92  {delete[] Average_gradient; Average_gradient=0;}
93  }
94 
95  ///Calculate the pressure from the unknowns
96  double pressure(const Vector<double> &u) const;
97 
98  /// \short Return the value of gamma
99  double gamma() const {return *Gamma_pt;}
100 
101  /// \short Access function for the pointer to gamma
102  double* &gamma_pt() {return Gamma_pt;}
103 
104  /// \short Access function for the pointer to gamma (const version)
105  double* const &gamma_pt() const {return Gamma_pt;}
106 
107  ///\short The number of unknowns at each node is the number of flux components
108  inline unsigned required_nvalue(const unsigned &n) const {return DIM+2;}
109 
110  ///Compute the error and norm of solution integrated over the element
111  ///Does not plot the error in the outfile
112  void compute_error(std::ostream &outfile,
114  exact_solution_pt, const double &t,
115  Vector<double> &error, Vector<double> &norm)
116  {
117  //Find the number of fluxes
118  const unsigned n_flux = this->nflux();
119  //Find the number of nodes
120  const unsigned n_node = this->nnode();
121  //Storage for the shape function and derivatives of shape function
122  Shape psi(n_node);
123  DShape dpsidx(n_node,DIM);
124 
125  //Find the number of integration points
126  unsigned n_intpt = this->integral_pt()->nweight();
127 
128  error.resize(n_flux);
129  norm.resize(n_flux);
130  for(unsigned i=0;i<n_flux;i++)
131  {error[i] = 0.0; norm[i] = 0.0;}
132 
133  Vector<double> s(DIM);
134 
135  //Loop over the integration points
136  for(unsigned ipt=0;ipt<n_intpt;ipt++)
137  {
138  //Get the shape functions at the knot
139  double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
140 
141  //Get the integral weight
142  double W = this->integral_pt()->weight(ipt)*J;
143 
144  //Storage for the local functions
146  Vector<double> interpolated_u(n_flux,0.0);
147 
148  //Loop over the shape functions
149  for(unsigned l=0;l<n_node;l++)
150  {
151  //Locally cache the shape function
152  const double psi_ = psi(l);
153  for(unsigned i=0;i<DIM;i++)
154  {
155  interpolated_x[i] += this->nodal_position(l,i)*psi_;
156  }
157 
158  for(unsigned i=0;i<n_flux;i++)
159  {
160  //Calculate the velocity and tangent vector
161  interpolated_u[i] += this->nodal_value(l,i)*psi_;
162  }
163  }
164 
165  //Now get the exact solution at this value of x and t
166  Vector<double> exact_u(n_flux,0.0);
167  (*exact_solution_pt)(t,interpolated_x,exact_u);
168 
169  //Loop over the unknowns
170  for(unsigned i=0;i<n_flux;i++)
171  {
172  error[i] +=
173  pow((interpolated_u[i] - exact_u[i]),2.0)*W;
174  norm[i] += interpolated_u[i]*interpolated_u[i]*W;
175  }
176  }
177  }
178 
179 
180  /// \short Output function:
181  /// x,y,u or x,y,z,u
182  void output(std::ostream &outfile)
183  {
184  unsigned nplot=5;
185  output(outfile,nplot);
186  }
187 
188  /// \short Output function:
189  /// x,y,u or x,y,z,u at n_plot^DIM plot points
190  void output(std::ostream &outfile, const unsigned &n_plot);
191 
193  {
194  //Find the number of fluxes
195  const unsigned n_flux = this->nflux();
196  //Resize the memory if necessary
197  if(this->Average_prim_value==0)
198  {this->Average_prim_value = new double[n_flux];}
199  //Will move this one day
200  if(this->Average_gradient==0)
201  {this->Average_gradient = new double[n_flux*DIM];}
202  }
203 
204  /// \short Return access to the average gradient
205  double &average_gradient(const unsigned &i, const unsigned &j)
206  {
207  if(Average_gradient==0)
208  {
209  throw OomphLibError("Averages not calculated yet",
210  OOMPH_CURRENT_FUNCTION,
211  OOMPH_EXCEPTION_LOCATION);
212  }
213  return Average_gradient[DIM*i+j];
214  }
215 
216  /// \short Return the average values
217  double &average_prim_value(const unsigned &i)
218  {
219  if(Average_prim_value==0)
220  {
221  throw OomphLibError("Averages not calculated yet",
222  OOMPH_CURRENT_FUNCTION,
223  OOMPH_EXCEPTION_LOCATION);
224  }
225  return Average_prim_value[i];
226  }
227 
228 };
229 
230 
231 template <unsigned DIM, unsigned NNODE_1D>
233  public virtual QSpectralElement<DIM,NNODE_1D>,
234  public virtual EulerEquations<DIM>
235  {
236  public:
237 
238 
239  ///\short Constructor: Call constructors for QElement and
240  /// Advection Diffusion equations
242  EulerEquations<DIM>()
243  { }
244 
245  /// Broken copy constructor
248  {
249  BrokenCopy::broken_copy("QSpectralEulerElement");
250  }
251 
252  /// Broken assignment operator
253 //Commented out broken assignment operator because this can lead to a conflict warning
254 //when used in the virtual inheritence hierarchy. Essentially the compiler doesn't
255 //realise that two separate implementations of the broken function are the same and so,
256 //quite rightly, it shouts.
257  /*void operator=(
258  const QSpectralEulerElement<DIM,NNODE_1D>&)
259  {
260  BrokenCopy::broken_assign("QSpectralEulerElement");
261  }*/
262 
263 
264  /// \short Output function:
265  /// x,y,u or x,y,z,u
266  void output(std::ostream &outfile)
267  {EulerEquations<DIM>::output(outfile);}
268 
269  /// \short Output function:
270  /// x,y,u or x,y,z,u at n_plot^DIM plot points
271  void output(std::ostream &outfile, const unsigned &n_plot)
272  {EulerEquations<DIM>::output(outfile,n_plot);}
273 
274 
275  /*/// \short C-style output function:
276  /// x,y,u or x,y,z,u
277  void output(FILE* file_pt)
278  {
279  EulerEquations<NFLUX,DIM>::output(file_pt);
280  }
281 
282  /// \short C-style output function:
283  /// x,y,u or x,y,z,u at n_plot^DIM plot points
284  void output(FILE* file_pt, const unsigned &n_plot)
285  {
286  EulerEquations<NFLUX,DIM>::output(file_pt,n_plot);
287  }
288 
289  /// \short Output function for an exact solution:
290  /// x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points
291  void output_fct(std::ostream &outfile, const unsigned &n_plot,
292  FiniteElement::SteadyExactSolutionFctPt
293  exact_soln_pt)
294  {
295  EulerEquations<NFLUX,DIM>::
296  output_fct(outfile,n_plot,exact_soln_pt);}
297 
298 
299  /// \short Output function for a time-dependent exact solution.
300  /// x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points
301  /// (Calls the steady version)
302  void output_fct(std::ostream &outfile, const unsigned &n_plot,
303  const double& time,
304  FiniteElement::UnsteadyExactSolutionFctPt
305  exact_soln_pt)
306  {
307  EulerEquations<NFLUX,DIM>::
308  output_fct(outfile,n_plot,time,exact_soln_pt);
309  }*/
310 
311 
312 protected:
313 
314  /// Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
316  const Vector<double> &s,
317  Shape &psi,
318  DShape &dpsidx,
319  Shape &test,
320  DShape &dtestdx) const;
321 
322  /// \short Shape, test functions & derivs. w.r.t. to global coords. at
323  /// integration point ipt. Return Jacobian.
325  const unsigned& ipt,
326  Shape &psi,
327  DShape &dpsidx,
328  Shape &test,
329  DShape &dtestdx)
330  const;
331 
332 };
333 
334 //Inline functions:
335 
336 
337 //======================================================================
338 /// \short Define the shape functions and test functions and derivatives
339 /// w.r.t. global coordinates and return Jacobian of mapping.
340 ///
341 /// Galerkin: Test functions = shape functions
342 //======================================================================
343 template<unsigned DIM, unsigned NNODE_1D>
346  Shape &psi,
347  DShape &dpsidx,
348  Shape &test,
349  DShape &dtestdx) const
350 {
351  //Call the geometrical shape functions and derivatives
352  double J = this->dshape_eulerian(s,psi,dpsidx);
353 
354  //Loop over the test functions and derivatives and set them equal to the
355  //shape functions
356  for(unsigned i=0;i<NNODE_1D;i++)
357  {
358  test[i] = psi[i];
359  for(unsigned j=0;j<DIM;j++)
360  {
361  dtestdx(i,j) = dpsidx(i,j);
362  }
363  }
364 
365  //Return the jacobian
366  return J;
367 }
368 
369 
370 
371 //======================================================================
372 /// Define the shape functions and test functions and derivatives
373 /// w.r.t. global coordinates and return Jacobian of mapping.
374 ///
375 /// Galerkin: Test functions = shape functions
376 //======================================================================
377 template<unsigned DIM, unsigned NNODE_1D>
380  const unsigned &ipt,
381  Shape &psi,
382  DShape &dpsidx,
383  Shape &test,
384  DShape &dtestdx) const
385 {
386  //Call the geometrical shape functions and derivatives
387  double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
388 
389  //Set the test functions equal to the shape functions (pointer copy)
390  test = psi;
391  dtestdx = dpsidx;
392 
393  //Return the jacobian
394  return J;
395 }
396 
397 
398 ////////////////////////////////////////////////////////////////////////
399 ////////////////////////////////////////////////////////////////////////
400 ////////////////////////////////////////////////////////////////////////
401 
402 
403 
404 //=======================================================================
405 /// \short Face geometry for the QEulerElement elements:
406 /// The spatial dimension of the face elements is one lower than that
407 /// of the bulk element but they have the same number of points along
408 /// their 1D edges.
409 //=======================================================================
410 template<unsigned DIM, unsigned NNODE_1D>
412  QSpectralEulerElement<DIM,NNODE_1D> >:
413 public virtual QSpectralElement<DIM-1,NNODE_1D>
414 {
415 
416  public:
417 
418  /// \short Constructor: Call the constructor for the
419  /// appropriate lower-dimensional QElement
420  FaceGeometry() : QElement<DIM-1,NNODE_1D>() {}
421 
422 };
423 
424 
425 
426 
427  //======================================================================
428  /// FaceElement for Discontinuous Galerkin Problems
429  //======================================================================
430  template<class ELEMENT>
431  class DGEulerFaceElement : public virtual FaceGeometry<ELEMENT>,
432  public virtual DGFaceElement
433  {
434  unsigned Nflux;
435 
437 
438  public:
439 
440  /// Constructor
441  DGEulerFaceElement(FiniteElement* const &element_pt,
442  const int &face_index) :
443  FaceGeometry<ELEMENT>(), DGFaceElement()
444  {
445  //Attach geometric information to the element
446  //N.B. This function also assigns nbulk_value from required_nvalue
447  //of the bulk element.
448  element_pt->build_face_element(face_index,this);
449  //Set the value of the flux
450  Nflux = 2 + element_pt->dim();
451  //Use the Face integration scheme of the element for 2D elements
452  if(Nflux > 3)
453  {
455  dynamic_cast<ELEMENT*>(element_pt)->face_integration_pt());
456  }
457  }
458 
459  /// Specify the value of nodal zeta from the face geometry
460  /// \short The "global" intrinsic coordinate of the element when
461  /// viewed as part of a geometric object should be given by
462  /// the FaceElement representation, by default (needed to break
463  /// indeterminacy if bulk element is SolidElement)
464  double zeta_nodal(const unsigned &n, const unsigned &k,
465  const unsigned &i) const
466  {return FaceElement::zeta_nodal(n,k,i);}
467 
468 
469  //There is a single required n_flux
470  unsigned required_nflux() {return Nflux;}
471 
472  ///Calculate the normal flux, so this is the dot product of the
473  ///numerical flux with n_out
474  void numerical_flux(const Vector<double> &n_out,
475  const Vector<double> &u_int,
476  const Vector<double> &u_ext,
477  Vector<double> &flux)
478  {
479  //Let's follow the yellow book and use local Lax-Friedrichs
480  //This is almost certainly not the best flux to use ---
481  //further investigation is required here.
482  //Cache the bulk element
483  ELEMENT* cast_bulk_element_pt =
484  dynamic_cast<ELEMENT*>(this->bulk_element_pt());
485 
486  //Find the dimension of the problem
487  const unsigned dim = cast_bulk_element_pt->dim();
488  //Find the number of variables
489  const unsigned n_flux = this->Nflux;
490 
491  const double gamma = cast_bulk_element_pt->gamma();
492 
493  //Now let's find the fluxes
494  DenseMatrix<double> flux_int(n_flux,dim);
495  DenseMatrix<double> flux_ext(n_flux,dim);
496 
497  cast_bulk_element_pt->flux(u_int,flux_int);
498  cast_bulk_element_pt->flux(u_ext,flux_ext);
499 
500  DenseMatrix<double> flux_av(n_flux,dim);
501  for(unsigned i=0;i<n_flux;i++)
502  {
503  for(unsigned j=0;j<dim;j++)
504  {
505  flux_av(i,j) = 0.5*(flux_int(i,j) + flux_ext(i,j));
506  }
507  }
508 
509  flux.initialise(0.0);
510  //Now set the first part of the numerical flux
511  for(unsigned i=0;i<n_flux;i++)
512  {
513  for(unsigned j=0;j<dim;j++)
514  {
515  flux[i] += flux_av(i,j)*n_out[j];
516  }
517  }
518 
519  //Now let's find the normal jumps in the fluxes
520  Vector<double> jump(n_flux);
521  //Now we take the normal jump in the quantities
522  //The first two are scalars
523  for(unsigned i=0;i<2;i++)
524  {
525  jump[i] = u_int[i] - u_ext[i];
526  }
527 
528  //The next are vectors so treat them accordingly ??
529  //By taking the component dotted with the normal component
530  double velocity_jump = 0.0;
531  for(unsigned j=0;j<dim;j++)
532  {
533  velocity_jump += (u_int[2+j] - u_ext[2+j])*n_out[j];
534  }
535  //Now we multiply by the outer normal to get the vector flux
536  for(unsigned j=0;j<dim;j++)
537  {
538  jump[2+j] = velocity_jump*n_out[j];
539  }
540 
541 
542  //Let's find the Roe average
543  /* Vector<double> u_average(n_flux);
544  double sum = sqrt(u_int[0]) + sqrt(u_ext[0]);
545  for(unsigned i=0;i<n_flux;i++)
546  {
547  u_average[i] = (sqrt(u_int[0])*u_int[i] + sqrt(u_ext[0])*u_ext[i])/sum;
548  }
549 
550  //Find the average pressure
551  double p = cast_bulk_element_pt->pressure(u_average);
552 
553  //Now do the normal velocity
554  double vel = 0.0;
555  for(unsigned j=0;j<dim;j++)
556  {
557  //vel += u_average[2+j]*u_average[2+j];
558  vel += u_average[2+j]*n_out[j];
559  }
560  //vel = sqrt(vel)/u_average[0];
561  vel /= u_average[0];
562 
563  double arg = std::fabs(gamma*p/u_average[0]);*/
564 
565  //Let's calculate the internal and external pressures and then enthalpies
566  double p_int = cast_bulk_element_pt->pressure(u_int);
567  double p_ext = cast_bulk_element_pt->pressure(u_ext);
568 
569  //Limit the pressures to zero if necessary, but keep the energy the
570  //same
571  if(p_int < 0)
572  {oomph_info << "Negative int pressure" << p_int << "\n"; p_int=0.0;}
573  if(p_ext < 0)
574  {oomph_info << "Negative ext pressure " << p_ext << "\n"; p_ext=0.0;}
575 
576  double H_int = (u_int[1] + p_int)/u_int[0];
577  double H_ext = (u_ext[1] + p_ext)/u_ext[0];
578 
579  //Also calculate the velocities
580  Vector<double> vel_int(2), vel_ext(2);
581  for(unsigned j=0;j<dim;j++)
582  {
583  vel_int[j] = u_int[2+j]/u_int[0];
584  vel_ext[j] = u_ext[2+j]/u_ext[0];
585  }
586 
587  //Now we calculate the Roe averaged values
588  Vector<double> vel_average(dim);
589  double s_int = sqrt(u_int[0]);
590  double s_ext = sqrt(u_ext[0]);
591  double sum = s_int + s_ext;
592 
593  //Velocities
594  for(unsigned j=0;j<dim;j++)
595  {
596  vel_average[j] = (s_int*vel_int[j] + s_ext*vel_ext[j])/sum;
597  }
598 
599  //The enthalpy
600  double H_average = (s_int*H_int + s_ext*H_ext)/sum;
601 
602  //Now calculate the square of the sound speed
603  double arg = H_average;
604  for(unsigned j=0;j<dim;j++) {arg -= 0.5*vel_average[j];}
605  arg *= (gamma - 1.0);
606  //Get the local sound speed
607  if(arg < 0.0)
608  {oomph_info << "Square of sound speed is negative!\n"; arg = 0.0;}
609  double a = sqrt(arg);
610 
611  //Calculate the normal average velocity
612  //Now do the normal velocity
613  double vel = 0.0;
614  for(unsigned j=0;j<dim;j++) {vel += vel_average[j]*n_out[j];}
615 
616  Vector<double> eigA(3,0.0);
617 
618  eigA[0] = vel - a;
619  eigA[1] = vel;
620  eigA[2] = vel + a;
621 
622  double lambda = std::fabs(eigA[0]);
623  for(unsigned i=1;i<3;i++)
624  {
625  if(std::fabs(eigA[i]) > lambda) {lambda = std::fabs(eigA[i]);}
626  }
627 
628 
629  //Find the magnitude of the external velocity
630  /*double vel_mag = 0.0;
631  for(unsigned j=0;j<dim;j++) {vel_mag += u_ext[2+j]*u_ext[2+j];}
632  vel_mag = sqrt(vel_mag)/u_ext[0];
633 
634  //Get the pressure
635  double p = cast_bulk_element_pt->pressure(u_ext);
636  double lambda_ext = vel_mag + sqrt(std::fabs(gamma*p/u_ext[0]));
637 
638  //Let's do the same for the internal one
639  vel_mag = 0.0;
640  for(unsigned j=0;j<dim;j++) {vel_mag += u_int[2+j]*u_int[2+j];}
641  vel_mag = sqrt(vel_mag)/u_int[0];
642 
643  //Get the pressure
644  p = cast_bulk_element_pt->pressure(u_int);
645  double lambda_int = vel_mag + sqrt(std::fabs(gamma*p/u_int[0]));
646 
647  //Now take the largest
648  double lambda = lambda_int;
649  if(lambda_ext > lambda) {lambda = lambda_ext;}*/
650 
651  for(unsigned i=0;i<n_flux;i++) {flux[i] += 0.5*lambda*jump[i];}
652  }
653 
654 
655  };
656 
657 
658  //======================================================================
659  /// \short FaceElement for Discontinuous Galerkin Problems with reflection
660  /// boundary conditions
661  //======================================================================
662  template<class ELEMENT>
663  class DGEulerFaceReflectionElement : public virtual FaceGeometry<ELEMENT>,
664  public virtual DGFaceElement
665  {
666  unsigned Nflux;
667 
668  public:
669 
670  /// Constructor
672  const int &face_index) :
673  FaceGeometry<ELEMENT>(), DGFaceElement()
674  {
675  //Attach geometric information to the element
676  //N.B. This function also assigns nbulk_value from required_nvalue
677  //of the bulk element.
678  element_pt->build_face_element(face_index,this);
679  //Set the value of the flux
680  Nflux = 2 + element_pt->dim();
681  }
682 
683 
684  /// Specify the value of nodal zeta from the face geometry
685  /// \short The "global" intrinsic coordinate of the element when
686  /// viewed as part of a geometric object should be given by
687  /// the FaceElement representation, by default (needed to break
688  /// indeterminacy if bulk element is SolidElement)
689  double zeta_nodal(const unsigned &n, const unsigned &k,
690  const unsigned &i) const
691  {return FaceElement::zeta_nodal(n,k,i);}
692 
693  //There is a single required n_flux
694  unsigned required_nflux() {return Nflux;}
695 
696  ///We overload interpolated_u to reflect
698  Vector<double> &u)
699  {
700  //Get the standard interpolated_u
702 
703  //Now do the reflection condition for the velocities
704  //Find dot product of normal and velocities
705  const unsigned nodal_dim = this->nodal_dimension();
706  Vector<double> n(nodal_dim);
707  this->outer_unit_normal(s,n);
708 
709  double dot = 0.0;
710  for(unsigned j=0;j<nodal_dim;j++)
711  {
712  dot += n[j]*u[2+j];
713  }
714 
715  //Now subtract
716  for(unsigned j=0;j<nodal_dim;j++)
717  {
718  u[2+j] -= 2.0*dot*n[j];
719  }
720  }
721 };
722 
723 
724 
725 //=================================================================
726 ///General DGEulerClass. Establish the template parameters
727 //===================================================================
728 template<unsigned DIM, unsigned NNODE_1D>
730 {
731 
732 };
733 
734 
735 //==================================================================
736 //Specialization for 1D DG Advection element
737 //==================================================================
738 template<unsigned NNODE_1D>
739 class DGSpectralEulerElement<1,NNODE_1D> :
740  public QSpectralEulerElement<1,NNODE_1D>,
741  public DGElement
742 {
743  friend class
745 
746  Vector<double> Inverse_mass_diagonal;
747 
748 public:
749 
750  ///Overload the required number of fluxes for the DGElement
751  unsigned required_nflux() {return this->nflux();}
752 
753  //Calculate averages
754  void calculate_element_averages(double* &average_value)
756 
757 
758  //Constructor
760  DGElement()
761  {
762  //Need to up the order of integration for the accurate resolution
763  //of the quadratic non-linearities
764  this->set_integration_scheme(new Gauss<1,NNODE_1D>);
765  }
766 
767  //Dummy
768  Integral* face_integration_pt() const {return 0;}
769 
770 
772  {
773  }
774 
775  void build_all_faces()
776  {
777  //Make the two faces
778  Face_element_pt.resize(2);
779  //Make the face on the left
780  Face_element_pt[0] =
781  new DGEulerFaceElement<
783  (this,-1);
784  //Make the face on the right
785  Face_element_pt[1] =
786  new DGEulerFaceElement<
787  DGSpectralEulerElement<1,NNODE_1D> >
788  (this,+1);
789  }
790 
791 
792  ///\short Compute the residuals for the Navier--Stokes equations;
793  /// flag=1(or 0): do (or don't) compute the Jacobian as well.
794  void fill_in_generic_residual_contribution_flux_transport(
795  Vector<double> &residuals, DenseMatrix<double> &jacobian,
796  DenseMatrix<double> &mass_matrix, unsigned flag)
797  {
800  jacobian,
801  mass_matrix,
802  flag);
803 
804  this->add_flux_contributions_to_residuals(residuals,jacobian,flag);
805  }
806 
807 
808 //============================================================================
809 ///Function that returns the current value of the residuals
810 ///multiplied by the inverse mass matrix (virtual so that it can be overloaded
811 ///specific elements in which time and memory saving tricks can be applied)
812 //============================================================================
813 /* void get_inverse_mass_matrix_times_residuals(Vector<double> &minv_res)
814 {
815  //Now let's assemble stuff
816  const unsigned n_dof = this->ndof();
817 
818  //Resize and initialise the vector that will holds the residuals
819  minv_res.resize(n_dof);
820  for(unsigned n=0;n<n_dof;n++) {minv_res[n] = 0.0;}
821 
822  //If we are recycling the mass matrix
823  if(Mass_matrix_reuse_is_enabled && Mass_matrix_has_been_computed)
824  {
825  //Get the residuals
826  this->fill_in_contribution_to_residuals(minv_res);
827  }
828  //Otherwise
829  else
830  {
831  //Temporary mass matrix
832  DenseDoubleMatrix M(n_dof,n_dof,0.0);
833 
834  //Get the local mass matrix and residuals
835  this->fill_in_contribution_to_mass_matrix(minv_res,M);
836 
837  //Store the diagonal entries
838  Inverse_mass_diagonal.clear();
839  for(unsigned n=0;n<n_dof;n++) {Inverse_mass_diagonal.push_back(1.0/M(n,n));}
840 
841  //The mass matrix has been computed
842  Mass_matrix_has_been_computed=true;
843  }
844 
845  for(unsigned n=0;n<n_dof;n++) {minv_res[n] *= Inverse_mass_diagonal[n];}
846 
847 }*/
848 
849 
850 
851 };
852 
853 
854 //=======================================================================
855 /// Face geometry of the 1D DG elements
856 //=======================================================================
857 template<unsigned NNODE_1D>
859 public virtual PointElement
860 {
861  public:
863 };
864 
865 
866 
867 //==================================================================
868 ///Specialisation for 2D DG Elements
869 //==================================================================
870 template<unsigned NNODE_1D>
871 class DGSpectralEulerElement<2,NNODE_1D> :
872  public QSpectralEulerElement<2,NNODE_1D>,
873  public DGElement
874 {
875  friend class
877 
878  static Gauss<1,NNODE_1D> Default_face_integration_scheme;
879 
880 public:
881 
882  ///Overload the required number of fluxes for the DGElement
883  unsigned required_nflux() {return this->nflux();}
884 
885  //Calculate averages
886  void calculate_element_averages(double* &average_value)
888 
889  //Constructor
891  QSpectralEulerElement<2,NNODE_1D>(),
892  DGElement()
893  {
894  //Need to up the order of integration for the accurate resolution
895  //of the quadratic non-linearities
896  this->set_integration_scheme(new GaussLobattoLegendre<2,3*NNODE_1D/2>);
897  }
898 
901  Integral* face_integration_pt() const
902  {return &Default_face_integration_scheme;}
903 
904  void build_all_faces()
905  {
906  Face_element_pt.resize(4);
907  Face_element_pt[0] =
908  new DGEulerFaceElement<
910  (this,2);
911  Face_element_pt[1] =
912  new DGEulerFaceElement<
913  DGSpectralEulerElement<2,NNODE_1D> >
914  (this,1);
915  Face_element_pt[2] =
916  new DGEulerFaceElement<
917  DGSpectralEulerElement<2,NNODE_1D> >
918  (this,-2);
919  Face_element_pt[3] =
920  new DGEulerFaceElement<
921  DGSpectralEulerElement<2,NNODE_1D> >
922  (this,-1);
923  }
924 
925 
926  ///\short Compute the residuals for the Navier--Stokes equations;
927  /// flag=1(or 0): do (or don't) compute the Jacobian as well.
928  void fill_in_generic_residual_contribution_flux_transport(
929  Vector<double> &residuals, DenseMatrix<double> &jacobian,
930  DenseMatrix<double> &mass_matrix, unsigned flag)
931  {
934  jacobian,
935  mass_matrix,
936  flag);
937 
938  this->add_flux_contributions_to_residuals(residuals,jacobian,flag);
939  }
940 
941 
942 };
943 
944 
945 //=======================================================================
946 /// Face geometry of the DG elements
947 //=======================================================================
948 template<unsigned NNODE_1D>
950  public virtual QSpectralElement<1,NNODE_1D>
951 {
952  public:
953  FaceGeometry() : QSpectralElement<1,NNODE_1D>() {}
954 };
955 
956 
957 }
958 
959 #endif
double dshape_and_dtest_eulerian_flux_transport(const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
FaceElement for Discontinuous Galerkin Problems.
double gamma() const
Return the value of gamma.
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement. ...
void flux(const Vector< double > &u, DenseMatrix< double > &f)
Return the flux as a function of the unknown.
double & average_prim_value(const unsigned &i)
Return the average values.
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
The "global" intrinsic coordinate of the element when viewed as part of a geometric object should be ...
cstr elem_len * i
Definition: cfortran.h:607
double dshape_and_dtest_eulerian_at_knot_flux_transport(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Shape, test functions & derivs. w.r.t. to global coords. at integration point ipt. Return Jacobian.
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
Definition: elements.h:4367
virtual void set_integration_scheme(Integral *const &integral_pt)
Set the spatial integration scheme.
Definition: elements.cc:3148
static double Default_Gamma_Value
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
Definition: elements.h:2458
QSpectralEulerElement()
Constructor: Call constructors for QElement and Advection Diffusion equations.
A general Finite Element class.
Definition: elements.h:1271
char t
Definition: cfortran.h:572
A Base class for DGElements.
Definition: dg_elements.h:160
EulerEquations()
Return the flux derivatives as a function of the unknowns.
OomphInfo oomph_info
unsigned nodal_dimension() const
Return the required Eulerian dimension of the nodes in this element.
Definition: elements.h:2358
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
The "global" intrinsic coordinate of the element when viewed as part of a geometric object should be ...
void outer_unit_normal(const Vector< double > &s, Vector< double > &unit_normal) const
Compute outer unit normal at the specified local coordinate.
Definition: elements.cc:5695
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
unsigned required_nflux()
Set the number of flux components.
FaceElement for Discontinuous Galerkin Problems with reflection boundary conditions.
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3841
unsigned required_nflux()
Set the number of flux components.
virtual void interpolated_u(const Vector< double > &s, Vector< double > &f)
Return the interpolated values of the unknown fluxes.
Definition: dg_elements.cc:179
void output(std::ostream &outfile, const unsigned &n_plot)
Output function: x,y,u or x,y,z,u at n_plot^DIM plot points.
unsigned required_nvalue(const unsigned &n) const
The number of unknowns at each node is the number of flux components.
double pressure(const Vector< double > &u) const
Calculate the pressure from the unknowns.
virtual ~EulerEquations()
Destructor.
static char t char * s
Definition: cfortran.h:572
void output(std::ostream &outfile)
Broken assignment operator.
DGEulerFaceReflectionElement(FiniteElement *const &element_pt, const int &face_index)
Constructor.
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1900
QSpectralEulerElement(const QSpectralEulerElement< DIM, NNODE_1D > &dummy)
Broken copy constructor.
General DGEulerClass. Establish the template parameters.
void initialise(const _Tp &__value)
Iterate over all values and set to the desired value.
Definition: Vector.h:171
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
In a FaceElement, the "global" intrinsic coordinate of the element along the boundary, when viewed as part of a compound geometric object is specified using the boundary coordinate defined by the mesh. Note: Boundary coordinates will have been set up when creating the underlying mesh, and their values will have been stored at the nodes.
Definition: elements.h:4251
double *const & gamma_pt() const
Access function for the pointer to gamma (const version)
unsigned nflux() const
DIM momentum-components, a density and an energy are transported.
double nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. If the node is hanging, the appropriate interpolation is ...
Definition: elements.h:2215
double dot(const Vector< double > &a, const Vector< double > &b)
Probably not always best/fastest because not optimised for dimension but useful...
Definition: Vector.h:289
Specialisation for 2D DG Elements.
virtual void fill_in_generic_residual_contribution_flux_transport(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
Compute the residuals for the Navier–Stokes equations; flag=1(or 0): do (or don't) compute the Jacobi...
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
Base class for Discontinuous Galerkin Faces. These are responsible for calculating the normal fluxes ...
Definition: dg_elements.h:53
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2470
FiniteElement *& bulk_element_pt()
Pointer to higher-dimensional "bulk" element.
Definition: elements.h:4470
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
void(* UnsteadyExactSolutionFctPt)(const double &, const Vector< double > &, Vector< double > &)
Function pointer for function that computes Vector-valued time-dependent function as ...
Definition: elements.h:1726
void output(std::ostream &outfile)
Output function: x,y,u or x,y,z,u.
virtual void build_face_element(const int &face_index, FaceElement *face_element_pt)
Function for building a lower dimensional FaceElement on the specified face of the FiniteElement...
Definition: elements.cc:4922
void calculate_element_averages(double *&average_values)
Compute the average values of the fluxes.
Base class for Euler equations.
double * Average_gradient
Storage for the approximated gradients.
double * Average_prim_value
Storage for the average primitive values.
DGEulerFaceElement(FiniteElement *const &element_pt, const int &face_index)
Constructor.
void numerical_flux(const Vector< double > &n_out, const Vector< double > &u_int, const Vector< double > &u_ext, Vector< double > &flux)
double & average_gradient(const unsigned &i, const unsigned &j)
Return access to the average gradient.
double *& gamma_pt()
Access function for the pointer to gamma.
void interpolated_u(const Vector< double > &s, Vector< double > &u)
We overload interpolated_u to reflect.
void compute_error(std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_solution_pt, const double &t, Vector< double > &error, Vector< double > &norm)
virtual double dshape_eulerian_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsidx) const
Return the geometric shape functions and also first derivatives w.r.t. global coordinates at the ipt-...
Definition: elements.cc:3252