69{
70 int i, j;
71 int nTracePts =
m_fields[0]->GetTrace()->GetNpoints();
72 int nVariables = physarray.size();
74
75 const Array<OneD, const int> &traceBndMap =
m_fields[0]->GetTraceBndMap();
76
77
78
79 Array<OneD, NekDouble> Vn(nTracePts, 0.0);
80 Array<OneD, NekDouble> Vel(nTracePts, 0.0);
81 for (i = 0; i < nDimensions; ++i)
82 {
83 Vmath::Vdiv(nTracePts, Fwd[i + 1], 1, Fwd[0], 1, Vel, 1);
85 }
86
87
88
89 Array<OneD, NekDouble> absVel(nTracePts, 0.0);
90 m_varConv->GetAbsoluteVelocity(Fwd, absVel);
91
92
93 Array<OneD, NekDouble> soundSpeed(nTracePts);
94 m_varConv->GetSoundSpeed(Fwd, soundSpeed);
95
96
97 Array<OneD, NekDouble> Mach(nTracePts, 0.0);
98 Vmath::Vdiv(nTracePts, Vn, 1, soundSpeed, 1, Mach, 1);
100
101
102 int e, id1, id2, npts, pnt;
104
105
106 for (e = 0;
108 {
111 ->GetExp(e)
112 ->GetTotPoints();
113 id1 =
115 id2 =
117
118
120 Array<OneD, NekDouble> rho(npts, Fwd[0] + id2);
121 Array<OneD, NekDouble> Ei(npts);
123
124
125 for (i = 0; i < npts; i++)
126 {
127 pnt = id2 + i;
128
129
130 if (Mach[pnt] < 0.99)
131 {
132
134 for (j = 1; j < nVariables - 1; ++j)
135 {
136 Ek += 0.5 * (Fwd[j][pnt] * Fwd[j][pnt]) / Fwd[0][pnt];
137 }
138
139 rhoeb = Fwd[0][pnt] * Ei[i] + Ek;
140
141
142 for (j = 0; j < nVariables - 1; ++j)
143 {
146 ->UpdatePhys())[id1 + i] = Fwd[j][pnt];
147 }
148
151 ->UpdatePhys())[id1 + i] =
152 2.0 * rhoeb - Fwd[nVariables - 1][pnt];
153 }
154
155 else
156 {
157 for (j = 0; j < nVariables; ++j)
158 {
159
162 ->UpdatePhys())[id1 + i] = Fwd[j][pnt];
163 }
164 }
165 }
166 }
167}
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.