74{
75 boost::ignore_unused(time);
76
77 int i, j;
78 int nTracePts =
m_fields[0]->GetTrace()->GetNpoints();
79 int nVariables = physarray.size();
81
82 const Array<OneD, const int> &traceBndMap =
m_fields[0]->GetTraceBndMap();
83
84
85
86 Array<OneD, NekDouble> Vn(nTracePts, 0.0);
87 Array<OneD, NekDouble> Vel(nTracePts, 0.0);
88 for (i = 0; i < nDimensions; ++i)
89 {
90 Vmath::Vdiv(nTracePts, Fwd[i + 1], 1, Fwd[0], 1, Vel, 1);
92 }
93
94
95
96 Array<OneD, NekDouble> absVel(nTracePts, 0.0);
97 m_varConv->GetAbsoluteVelocity(Fwd, absVel);
98
99
100 Array<OneD, NekDouble>
pressure(nTracePts);
101 Array<OneD, NekDouble> soundSpeed(nTracePts);
102
104 m_varConv->GetSoundSpeed(Fwd, soundSpeed);
105
106
107 Array<OneD, NekDouble> Mach(nTracePts, 0.0);
108 Vmath::Vdiv(nTracePts, Vn, 1, soundSpeed, 1, Mach, 1);
110
111
112 int e, id1, id2, npts, pnt;
114
115
116 for (e = 0;
118 {
121 ->GetExp(e)
122 ->GetTotPoints();
123 id1 =
125 id2 =
127
128
129 Array<OneD, NekDouble> tmpPressure(npts,
pressure + id2);
130 Array<OneD, NekDouble> rho(npts, Fwd[0] + id2);
131 Array<OneD, NekDouble> Ei(npts);
132 m_varConv->GetEFromRhoP(rho, tmpPressure, Ei);
133
134
135 for (i = 0; i < npts; i++)
136 {
137 pnt = id2 + i;
138
139
140 if (Mach[pnt] < 0.99)
141 {
142
143 for (j = 0; j < nVariables - 1; ++j)
144 {
148 }
149
150
152 for (j = 1; j < nVariables - 1; ++j)
153 {
154 Ek += 0.5 *
158 }
159
160 rhoeb = Fwd[0][pnt] * Ei[i] + Ek;
161
164 ->UpdatePhys())[id1 + i] = rhoeb;
165 }
166
167 else
168 {
169 for (j = 0; j < nVariables; ++j)
170 {
171
174 ->UpdatePhys())[id1 + i] = Fwd[j][pnt];
175 }
176 }
177 }
178 }
179}
int m_spacedim
Space dimension.
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Trace normals.
VariableConverterSharedPtr m_varConv
Auxiliary object to convert variables.
void Vabs(int n, const T *x, const int incx, T *y, const int incy)
vabs: y = |x|
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.