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// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26// Non-inline member function of the flux transport elements class
27
28#include "euler_elements.h"
29
30namespace oomph
31{
32 //===========================================================
33 /// Set the default value of Gamma to be 1.4 in all three
34 /// dimension specialisations. This form seems to be required for
35 /// gcc 4.4.4, rather than a more general templated version.
36 //===========================================================
37 template<>
39 template<>
41 template<>
43
44
45 //======================================================
46 /// Calculate the pressure value from the unknowns
47 //=====================================================
48 template<unsigned DIM>
50 {
51 // Initialise the pressure to zero
52 double p = 0.0;
53 // Subtract off the momentum components
54 for (unsigned j = 0; j < DIM; j++)
55 {
56 p -= u[2 + j] * u[2 + j];
57 }
58 // Multiply by half and divide by the extra density component
59 p *= 0.5 / u[0];
60 // Now add on the energy
61 p += u[1];
62 // Finaly multiply by gamma minus 1
63 p *= (this->gamma() - 1);
64
65 // return the pressure
66 return p;
67 }
68
69
70 //=========================================================
71 /// Return the flux as a function of the unknowns
72 /// The unknowns are stored as density, energy and then
73 /// the velocity components
74 //=========================================================
75 template<unsigned DIM>
78 {
79 // The density flux is the momentum
80 for (unsigned j = 0; j < DIM; j++)
81 {
82 f(0, j) = u[2 + j];
83 }
84 // The energy flux is given by the velocity component multiplied by
85 // E + p
86 // Find the pressure
87 double p = pressure(u);
88
89 // The we can do the energy fluxes
90 for (unsigned j = 0; j < DIM; j++)
91 {
92 f(1, j) = u[2 + j] * (u[1] + p) / u[0];
93 }
94
95 // Now the momentum fluxes
96 for (unsigned i = 0; i < DIM; i++)
97 {
98 for (unsigned j = 0; j < DIM; j++)
99 {
100 f(2 + i, j) = u[2 + j] * u[2 + i] / u[0];
101 }
102 // Add the additional diagonal terms
103 f(2 + i, i) += p;
104 }
105 }
106
107
108 //====================================================================
109 /// Output function, print the values of all unknowns
110 //==================================================================
111 template<unsigned DIM>
112 void EulerEquations<DIM>::output(std::ostream& outfile, const unsigned& nplot)
113 {
114 // Find the number of fluxes
115 const unsigned n_flux = this->nflux();
116
117 // Vector of local coordinates
118 Vector<double> s(DIM);
119 // Vector of values
120 Vector<double> u(n_flux, 0.0);
121
122 // Tecplot header info
123 outfile << this->tecplot_zone_string(nplot);
124
125 // Loop over plot points
126 unsigned num_plot_points = this->nplot_points(nplot);
127 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
128 {
129 // Get local coordinates of plot point
130 this->get_s_plot(iplot, nplot, s);
131
132 // Coordinates
133 for (unsigned i = 0; i < DIM; i++)
134 {
135 outfile << this->interpolated_x(s, i) << " ";
136 }
137
138 // Values
139 for (unsigned i = 0; i < n_flux; i++)
140 {
141 u[i] = this->interpolated_u_flux_transport(s, i);
142 outfile << u[i] << " ";
143 }
144
145 // Now output the velocity
146 for (unsigned j = 0; j < DIM; j++)
147 {
148 outfile << u[2 + j] / u[0] << " ";
149 }
150
151 // Also the pressure
152 outfile << pressure(u);
153
154 outfile << std::endl;
155 }
156 outfile << std::endl;
157
158 // Write tecplot footer (e.g. FE connectivity lists)
159 this->write_tecplot_zone_footer(outfile, nplot);
160 }
161
162
163 //======================================================================
164 /// Return the flux derivatives as a function of the unknowns
165 //=====================================================================
166 /*template<unsigned DIM>
167 void EulerEquations<DIM>::
168 dflux_du(const Vector<double> &u, RankThreeTensor<double> &df_du)
169 {
170 }*/
171
172 template class EulerEquations<1>;
173 template class EulerEquations<2>;
174 template class EulerEquations<3>;
175
176 template<unsigned NNODE_1D>
179
180 template class DGSpectralEulerElement<2, 2>;
181 template class DGSpectralEulerElement<2, 3>;
182 template class DGSpectralEulerElement<2, 4>;
183
184} // namespace oomph
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
General DGEulerClass. Establish the template parameters.
Base class for Euler equations.
double pressure(const Vector< double > &u) const
Calculate the pressure from the unknowns.
void flux(const Vector< double > &u, DenseMatrix< double > &f)
Return the flux as a function of the unknown.
void output(std::ostream &outfile)
Output function: x,y,u or x,y,z,u.
Class for multidimensional Gaussian integration rules.
Definition: integral.h:145
//////////////////////////////////////////////////////////////////// ////////////////////////////////...