euler_elements.cc
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: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
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 //Non-inline member function of the flux transport elements class
31 
32 #include "euler_elements.h"
33 
34 namespace oomph
35 {
36  //===========================================================
37  /// Set the default value of Gamma to be 1.4 in all three
38  /// dimension specialisations. This form seems to be required for
39  /// gcc 4.4.4, rather than a more general templated version.
40  //===========================================================
41  template<>
43  template<>
45  template<>
47 
48 
49 
50  //======================================================
51  ///Calculate the pressure value from the unknowns
52  //=====================================================
53  template<unsigned DIM>
55  {
56  //Initialise the pressure to zero
57  double p = 0.0;
58  //Subtract off the momentum components
59  for(unsigned j=0;j<DIM;j++) {p -= u[2+j]*u[2+j];}
60  //Multiply by half and divide by the extra density component
61  p *= 0.5/u[0];
62  //Now add on the energy
63  p += u[1];
64  //Finaly multiply by gamma minus 1
65  p *= (this->gamma() - 1);
66 
67  //return the pressure
68  return p;
69  }
70 
71 
72  //=========================================================
73  /// \short Return the flux as a function of the unknowns
74  /// The unknowns are stored as density, energy and then
75  /// the velocity components
76  //=========================================================
77  template<unsigned DIM>
80  {
81  //The density flux is the momentum
82  for(unsigned j=0;j<DIM;j++) {f(0,j) = u[2+j];}
83  //The energy flux is given by the velocity component multiplied by
84  //E + p
85  //Find the pressure
86  double p = pressure(u);
87 
88  //The we can do the energy fluxes
89  for(unsigned j=0;j<DIM;j++) {f(1,j) = u[2+j]*(u[1] + p)/u[0];}
90 
91  //Now the momentum fluxes
92  for(unsigned i=0;i<DIM;i++)
93  {
94  for(unsigned j=0;j<DIM;j++)
95  {
96  f(2+i,j) = u[2+j]*u[2+i]/u[0];
97  }
98  //Add the additional diagonal terms
99  f(2+i,i) += p;
100  }
101  }
102 
103 
104  //====================================================================
105  ///Output function, print the values of all unknowns
106  //==================================================================
107  template<unsigned DIM>
109  output(std::ostream &outfile, const unsigned &nplot)
110  {
111  //Find the number of fluxes
112  const unsigned n_flux = this->nflux();
113 
114  //Vector of local coordinates
115  Vector<double> s(DIM);
116  //Vector of values
117  Vector<double> u(n_flux,0.0);
118 
119  // Tecplot header info
120  outfile << this->tecplot_zone_string(nplot);
121 
122  // Loop over plot points
123  unsigned num_plot_points = this->nplot_points(nplot);
124  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
125  {
126 
127  // Get local coordinates of plot point
128  this->get_s_plot(iplot,nplot,s);
129 
130  // Coordinates
131  for(unsigned i=0;i<DIM;i++)
132  {
133  outfile << this->interpolated_x(s,i) << " ";
134  }
135 
136  // Values
137  for(unsigned i=0;i<n_flux;i++)
138  {
139  u[i] = this->interpolated_u_flux_transport(s,i);
140  outfile << u[i] << " ";
141  }
142 
143  //Now output the velocity
144  for(unsigned j=0;j<DIM;j++)
145  {
146  outfile << u[2+j]/u[0] << " ";
147  }
148 
149  //Also the pressure
150  outfile << pressure(u);
151 
152  outfile << std::endl;
153  }
154  outfile << std::endl;
155 
156  // Write tecplot footer (e.g. FE connectivity lists)
157  this->write_tecplot_zone_footer(outfile,nplot);
158 
159  }
160 
161 
162  //======================================================================
163  /// \short Return the flux derivatives as a function of the unknowns
164  //=====================================================================
165  /*template<unsigned DIM>
166  void EulerEquations<DIM>::
167  dflux_du(const Vector<double> &u, RankThreeTensor<double> &df_du)
168  {
169  }*/
170 
171  template class EulerEquations<1>;
172  template class EulerEquations<2>;
173  template class EulerEquations<3>;
174 
175 template<unsigned NNODE_1D>
178 
179  template class DGSpectralEulerElement<2,2>;
180  template class DGSpectralEulerElement<2,3>;
181  template class DGSpectralEulerElement<2,4>;
182 
183 }
void flux(const Vector< double > &u, DenseMatrix< double > &f)
Return the flux as a function of the unknown.
cstr elem_len * i
Definition: cfortran.h:607
double pressure(const Vector< double > &u) const
Calculate the pressure from the unknowns.
static char t char * s
Definition: cfortran.h:572
General DGEulerClass. Establish the template parameters.
void output(std::ostream &outfile)
Output function: x,y,u or x,y,z,u.
Base class for Euler equations.