Toggle navigation
Documentation
Big picture
The finite element method
The data structure
Not-so-quick guide
Optimisation
Order of action functions
Example codes and tutorials
List of example codes and tutorials
Meshing
Solvers
MPI parallel processing
Post-processing/visualisation
Other
Change log
Creating documentation
Coding conventions
Index
FAQ
Installation
Installation guide
Copyright
About
People
Contact/Get involved
Publications
Acknowledgements
Picture show
Go
src
flux_transport
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-2024 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
30
namespace
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
<>
38
double
EulerEquations<1>::Default_Gamma_Value
= 1.4;
39
template
<>
40
double
EulerEquations<2>::Default_Gamma_Value
= 1.4;
41
template
<>
42
double
EulerEquations<3>::Default_Gamma_Value
= 1.4;
43
44
45
//======================================================
46
/// Calculate the pressure value from the unknowns
47
//=====================================================
48
template
<
unsigned
DIM>
49
double
EulerEquations<DIM>::pressure
(
const
Vector<double>
& u)
const
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>
76
void
EulerEquations<DIM>::flux
(
const
Vector<double>
& u,
77
DenseMatrix<double>
& f)
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>
177
Gauss<1, NNODE_1D>
178
DGSpectralEulerElement<2, NNODE_1D>::Default_face_integration_scheme
;
179
180
template
class
DGSpectralEulerElement<2, 2>
;
181
template
class
DGSpectralEulerElement<2, 3>
;
182
template
class
DGSpectralEulerElement<2, 4>
;
183
184
}
// namespace oomph
s
static char t char * s
Definition:
cfortran.h:568
i
cstr elem_len * i
Definition:
cfortran.h:603
oomph::DGSpectralEulerElement
General DGEulerClass. Establish the template parameters.
Definition:
euler_elements.h:766
oomph::DenseMatrix< double >
oomph::EulerEquations
Base class for Euler equations.
Definition:
euler_elements.h:47
oomph::EulerEquations::pressure
double pressure(const Vector< double > &u) const
Calculate the pressure from the unknowns.
Definition:
euler_elements.cc:49
oomph::EulerEquations::flux
void flux(const Vector< double > &u, DenseMatrix< double > &f)
Return the flux as a function of the unknown.
Definition:
euler_elements.cc:76
oomph::EulerEquations::output
void output(std::ostream &outfile)
Output function: x,y,u or x,y,z,u.
Definition:
euler_elements.h:203
oomph::Gauss< 1, NNODE_1D >
oomph::Vector< double >
euler_elements.h
oomph
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition:
advection_diffusion_elements.cc:30