Commit 06c67cce authored by Dave Moxey's avatar Dave Moxey
Browse files

Merge remote-tracking branch 'upstream/master' into feature/hdf5-mesh-4

parents c379e224 796ba082
......@@ -27,6 +27,7 @@ v5.0.0
function definitions for the Absorption Forcing (!769)
- Improve performance of DisContField2D::v_ExtractTracePhys (!824)
- Fix small bug in Jacobian Energy (!857)
- fix variable name overriding in file functions (!870)
- Adds CFI CAD engine back-end (!864)
- Adds CFI Mesh IO support (!864)
- Cleanup of CAD system data structures (!864)
......@@ -62,6 +63,7 @@ v5.0.0
**FieldConvert**:
- Add input module for Semtex field files (!777)
- Fixed interppoints module (!760)
- Fix OutputTecplot in 2DH1D (!818)
- Move StreamFunction utility to a FieldConvert module (!809)
- Extend wss module to compressible flows (!810)
- Allow explicitly setting bool options of FieldConvert modules as false (!811)
......@@ -74,7 +76,7 @@ v5.0.0
**IncNavierStokesSolver**
- Replace steady-state check based on difference of norms by check based on
norm of the difference, to be consistent with the compressible solver (!832)
- Updated SVV to allow for the DGKernel extension (!851)
- Updated SVV to allow for the DGKernel extension (!851)
**CompressibleFlowSolver**
- Add 3D regression tests (!567)
......@@ -83,6 +85,10 @@ v5.0.0
- Add ability to use an exponential filtering for stabilization with
seg, quad and hex elements (!771, !862)
**APESolver:**
- Added two new boundary conditions to the APE system: RiemannInvariantBC
and WhiteNoise (!782)
**Documentation**:
- Added the developer-guide repository as a submodule (!751)
......@@ -92,8 +98,11 @@ v4.4.2
- Fix evaluation of points (e.g. HistoryPoints, Interpolation to pts) close to
the interface of two elements (!836)
- Fix deadlock in Hdf5 with homogeneous expansions (!858)
- Fix a few memory leaks in polylib (!863)
- Fix a crash when Interpolator is called on an empty field (!869)
- Fix petsc compile without MPI (!873)
- Fix calculation of BLPoints (!892)
- Fix deadlock in DiffusionLDG (!885)
- Fix uninitialised coefficients in DirectFull solver (!898)
**NekMesh**
......@@ -139,6 +148,7 @@ v4.4.1
IMEXGear, CNAB, 2nd order IMEX-DIRK, 3rd order IMEX-DIRK (!854)
- Fix bug due to subtractive cancellation in polylib routines (!778)
**FieldConvert:**
- Fix issue with field ordering in the interppointdatatofld module (!754)
- Fix issue with FieldConvert when range flag used (!761)
......
......@@ -51,11 +51,57 @@ Currently, only \inltt{APEUpwind} supported.
\subsection{Functions}
\begin{itemize}
\item \inltt{BaseFlow} Baseflow $(\overline{u}_i, \overline{p}, \overline{\rho})$ defined by the variables \inltt{u0,v0,w0,p0,rho0}
\item \inltt{Source} Source term $\overline{c}^2 q_c$
\item \inltt{InitialConditions}
\end{itemize}
\subsection{Boundary Conditions}
In addition to plain Dirichlet and Neumann boundary conditions, the APESolver features a wall boundray condition, a non-reflecting boundary and a white noise boundary condition.
\begin{itemize}
\item Rigid (Slip-) Wall Boundary Condition
\begin{lstlisting}[style=XmlStyle]
<BOUNDARYCONDITIONS>
<REGION REF="0">
<D VAR="p" USERDEFINEDTYPE="Wall" VALUE="0" />
<D VAR="u" USERDEFINEDTYPE="Wall" VALUE="0" />
<D VAR="v" USERDEFINEDTYPE="Wall" VALUE="0" />
<D VAR="w" USERDEFINEDTYPE="Wall" VALUE="0" />
</REGION>
</BOUNDARYCONDITIONS>
\end{lstlisting}
This BC imposes zero wall-normal perturbation velocity in a way that is more robust than using a Dirichlet boundary condition directly.
\item Non-Reflecting Boundary Condition
\begin{lstlisting}[style=XmlStyle]
<BOUNDARYCONDITIONS>
<REGION REF="0">
<D VAR="p" USERDEFINEDTYPE="RiemannInvariantBC"/>
<D VAR="u" USERDEFINEDTYPE="RiemannInvariantBC"/>
<D VAR="v" USERDEFINEDTYPE="RiemannInvariantBC"/>
<D VAR="w" USERDEFINEDTYPE="RiemannInvariantBC"/>
</REGION>
</BOUNDARYCONDITIONS>
\end{lstlisting}
The Riemann-Invariant BC approximates a non-reflecting (r.g. Farfield) boundary condition by setting incoming invariants to zero.
\item White Noise Boundary Condition
\begin{lstlisting}[style=XmlStyle]
<BOUNDARYCONDITIONS>
<REGION REF="0">
<D VAR="p" USERDEFINEDTYPE="Wall" VALUE="10" />
<D VAR="u" USERDEFINEDTYPE="Wall" VALUE="10" />
<D VAR="v" USERDEFINEDTYPE="Wall" VALUE="10" />
<D VAR="w" USERDEFINEDTYPE="Wall" VALUE="10" />
</REGION>
</BOUNDARYCONDITIONS>
\end{lstlisting}
The white noise BC imposes a stochastic, uniform pressure at the boundary. The implementation uses a Mersenne-Twister pseudo random number generator to generate white Gaussian noise.
The standard deviation $\sigma$ of the pressure is specified by the \inltt{VALUE} attribute.
\end{itemize}
\section{Examples}
\subsection{Aeroacoustic Wave Propagation}
In this section we explain how to set up a simple simulation of aeroacoustics in
......@@ -108,10 +154,6 @@ eventual source terms using the following functions:
<E VAR="rho0" VALUE="Rho0" />
</FUNCTION>
<FUNCTION NAME="Source">
<E VAR="S" VALUE="0" />
</FUNCTION>
<!-- Gaussian pulse located at the origin -->
<FUNCTION NAME="InitialConditions">
<E VAR="p" VALUE="100*exp(-32*(((x)*(x))+((y)*(y))))" />
......
......@@ -239,11 +239,10 @@ void OutputTecplot::OutputFromExp(po::variables_map &vm)
// Calculate coordinate dimension
int nBases = m_f->m_exp[0]->GetExp(0)->GetNumBases();
MultiRegions::ExpansionType HomoExpType = m_f->m_exp[0]->GetExpType();
m_coordim = m_f->m_exp[0]->GetExp(0)->GetCoordim();
if (HomoExpType == MultiRegions::e3DH1D)
if (m_f->m_numHomogeneousDir > 0)
{
int nPlanes = m_f->m_exp[0]->GetZIDs().num_elements();
if (nPlanes == 1) // halfMode case
......@@ -252,17 +251,12 @@ void OutputTecplot::OutputFromExp(po::variables_map &vm)
}
else
{
nBases += 1;
m_coordim += 1;
nBases += m_f->m_numHomogeneousDir;
m_coordim += m_f->m_numHomogeneousDir;
NekDouble tmp = m_numBlocks * (nPlanes - 1);
m_numBlocks = (int)tmp;
}
}
else if (HomoExpType == MultiRegions::e3DH2D)
{
nBases += 2;
m_coordim += 2;
}
m_zoneType = (TecplotZoneType)(2*(nBases-1) + 1);
......@@ -935,18 +929,45 @@ void OutputTecplot::CalculateConnectivity()
if (nbase == 1)
{
int cnt2 = 0;
int np0 = m_f->m_exp[0]->GetExp(i)->GetNumPoints(0);
int cnt2 = 0;
int np0 = m_f->m_exp[0]->GetExp(i)->GetNumPoints(0);
int nPlanes = 1;
Array<OneD, int> conn(2 * (np0 - 1));
for (k = 1; k < np0; ++k)
if (m_f->m_exp[0]->GetExpType() == MultiRegions::e2DH1D)
{
conn[cnt2++] = cnt + k;
conn[cnt2++] = cnt + k - 1;
nPlanes = m_f->m_exp[0]->GetZIDs().num_elements();
if (nPlanes > 1)
{
int totPoints = m_f->m_exp[0]->GetPlane(0)->GetTotPoints();
Array<OneD, int> conn(4 * (np0 - 1) * (nPlanes - 1));
for (int n = 1; n < nPlanes; ++n)
{
for (k = 1; k < np0; ++k)
{
conn[cnt2++] = cnt + (n - 1) * totPoints + k;
conn[cnt2++] = cnt + (n - 1) * totPoints + k - 1;
conn[cnt2++] = cnt + n * totPoints + k - 1;
conn[cnt2++] = cnt + n * totPoints + k;
}
}
m_conn[i] = conn;
}
}
m_conn[i] = conn;
if (nPlanes == 1)
{
Array<OneD, int> conn(2 * (np0 - 1));
for (k = 1; k < np0; ++k)
{
conn[cnt2++] = cnt + k;
conn[cnt2++] = cnt + k - 1;
}
m_conn[i] = conn;
}
}
else if (nbase == 2)
{
......
......@@ -63,12 +63,12 @@ namespace Nektar
}
else
{
NekDouble a = 2.0 * (1.0-r) / (1.0 - pow(r,(double)npts));
NekDouble a = 2.0 * (1.0-r) / (1.0 - pow(r,(double)(npts-1)));
m_points[0][0] = -1.0;
for (unsigned int i = 1; i < npts; ++i)
{
m_points[0][i] = m_points[0][i-1] + a*pow(r,(double)i);
m_points[0][i] = m_points[0][i-1] + a*pow(r,(double)(i-1));
}
m_points[0][npts-1] = 1.0;
......
This diff is collapsed.
......@@ -77,11 +77,12 @@ namespace Nektar
e0D,
e1D,
e2D,
e2DH1D,
e3DH1D,
e3DH2D,
e3D,
eNoType
};
};
MultiRegions::Direction const DirCartesianMap[] =
{
......
......@@ -47,6 +47,7 @@ namespace Nektar
ExpList2DHomogeneous1D::ExpList2DHomogeneous1D():
ExpListHomogeneous1D()
{
SetExpType(e2DH1D);
}
// Constructor for ExpList2DHomogeneous1D to act as a Explist2D field
......@@ -59,6 +60,7 @@ namespace Nektar
const Array<OneD, ExpListSharedPtr> &planes)
: ExpListHomogeneous1D(pSession,HomoBasis,lhom,useFFT,dealiasing)
{
SetExpType(e2DH1D);
int i, n, cnt, nel;
ASSERTL1(m_planes.num_elements() == planes.num_elements(),
......@@ -98,6 +100,7 @@ namespace Nektar
const Collections::ImplementationType ImpType):
ExpListHomogeneous1D(pSession,HomoBasis,lhom,useFFT,dealiasing)
{
SetExpType(e2DH1D);
int n, j, nel;
ExpList1DSharedPtr plane_zero;
......@@ -140,6 +143,7 @@ namespace Nektar
const ExpList2DHomogeneous1D &In):
ExpListHomogeneous1D(In)
{
SetExpType(e2DH1D);
ExpList1DSharedPtr zero_plane =
std::dynamic_pointer_cast<ExpList1D> (In.m_planes[0]);
......
......@@ -499,7 +499,7 @@ void SessionFunction::EvaluatePts(string pFieldName,
vector<string> fieldNames = outPts->GetFieldNames();
for (fieldInd = 0; fieldInd < fieldNames.size(); ++fieldInd)
{
if (outPts->GetFieldName(fieldInd) == pFieldName)
if (outPts->GetFieldName(fieldInd) == fileVar)
{
break;
}
......
......@@ -265,9 +265,11 @@ namespace Nektar
// if Vn >= 0, uflux = uFwd at Neumann, i.e.,
// edge::eForward, if V*n<0 <=> V*n_F<0, pick uflux = uFwd
// edge::eBackward, if V*n<0 <=> V*n_B>=0, pick uflux = uFwd
Array<OneD, NekDouble > uplus(nTracePts);
fields[i]->ExtractTracePhys(ufield[i], uplus);
if(fields[0]->GetBndCondExpansions().num_elements())
{
v_WeakPenaltyforScalar(fields, i, ufield[i], fluxtemp);
v_WeakPenaltyforScalar(fields, i, ufield[i], uplus, fluxtemp);
}
for (j = 0; j < nDim; ++j)
......@@ -296,6 +298,7 @@ namespace Nektar
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const int var,
const Array<OneD, const NekDouble> &ufield,
const Array<OneD, const NekDouble> &uplus,
Array<OneD, NekDouble> &penaltyflux)
{
int i, e, id1, id2;
......@@ -304,10 +307,7 @@ namespace Nektar
int nBndEdgePts, nBndEdges;
int cnt = 0;
int nBndRegions = fields[var]->GetBndCondExpansions().num_elements();
int nTracePts = fields[0]->GetTrace()->GetTotPoints();
Array<OneD, NekDouble > uplus(nTracePts);
fields[var]->ExtractTracePhys(ufield, uplus);
for (i = 0; i < nBndRegions; ++i)
{
// Number of boundary expansion related to that region
......@@ -434,11 +434,14 @@ namespace Nektar
qfluxtemp, 1,
qfluxtemp, 1);
Array<OneD, NekDouble > qtemp(nTracePts);
fields[i]->ExtractTracePhys(qfield[j][i], qtemp);
// Imposing weak boundary condition with flux
if (fields[0]->GetBndCondExpansions().num_elements())
{
v_WeakPenaltyforVector(fields, i, j,
qfield[j][i],
qfield[j][i], qtemp,
qfluxtemp, C11);
}
......@@ -465,15 +468,13 @@ namespace Nektar
const int var,
const int dir,
const Array<OneD, const NekDouble> &qfield,
const Array<OneD, const NekDouble> &qtemp,
Array<OneD, NekDouble> &penaltyflux,
NekDouble C11)
{
int i, e, id1, id2;
int nBndEdges, nBndEdgePts;
int nBndRegions = fields[var]->GetBndCondExpansions().num_elements();
int nTracePts = fields[0]->GetTrace()->GetTotPoints();
Array<OneD, NekDouble > qtemp(nTracePts);
int cnt = 0;
/*
......@@ -486,8 +487,6 @@ namespace Nektar
fields[0]->GetTrace()->GetNormals(m_traceNormals);
*/
fields[var]->ExtractTracePhys(qfield, qtemp);
for (i = 0; i < nBndRegions; ++i)
{
nBndEdges = fields[var]->
......
......@@ -83,6 +83,7 @@ namespace Nektar
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const int var,
const Array<OneD, const NekDouble> &ufield,
const Array<OneD, const NekDouble> &uplus,
Array<OneD, NekDouble> &penaltyflux);
virtual void v_NumFluxforVector(
......@@ -96,6 +97,7 @@ namespace Nektar
const int var,
const int dir,
const Array<OneD, const NekDouble> &qfield,
const Array<OneD, const NekDouble> &qtemp,
Array<OneD, NekDouble> &penaltyflux,
NekDouble C11);
};
......
......@@ -2430,7 +2430,6 @@ namespace Nektar
}
else if(mkey.ConstFactorExists(eFactorSVVDGKerDiffCoeff)) // Rodrigo/Mansoor's DG Kernel
{
NekDouble cutoff = mkey.GetConstFactor(eFactorSVVCutoffRatio);
NekDouble SvvDiffCoeff =
mkey.GetConstFactor(eFactorSVVDGKerDiffCoeff)*
mkey.GetConstFactor(eFactorSVVDiffCoeff);
......
......@@ -2053,7 +2053,6 @@ namespace Nektar
}
else if(mkey.ConstFactorExists(eFactorSVVDGKerDiffCoeff)) // Rodrigo/Mansoor's DG Kernel
{
NekDouble cutoff = mkey.GetConstFactor(eFactorSVVCutoffRatio);
NekDouble SvvDiffCoeff =
mkey.GetConstFactor(eFactorSVVDGKerDiffCoeff)*
mkey.GetConstFactor(eFactorSVVDiffCoeff);
......
......@@ -1962,7 +1962,6 @@ namespace Nektar
}
else if(mkey.ConstFactorExists(eFactorSVVDGKerDiffCoeff)) // Rodrigo/Mansoor's DG Kernel
{
NekDouble cutoff = mkey.GetConstFactor(eFactorSVVCutoffRatio);
NekDouble SvvDiffCoeff =
mkey.GetConstFactor(eFactorSVVDGKerDiffCoeff)*
mkey.GetConstFactor(eFactorSVVDiffCoeff);
......
......@@ -1581,7 +1581,6 @@ namespace Nektar
}
else if(mkey.ConstFactorExists(eFactorSVVDGKerDiffCoeff)) // Rodrigo/mansoor's DG kernel
{
NekDouble cutoff = mkey.GetConstFactor(eFactorSVVCutoffRatio);
NekDouble SvvDiffCoeff =
mkey.GetConstFactor(eFactorSVVDGKerDiffCoeff)*
mkey.GetConstFactor(eFactorSVVDiffCoeff);
......
......@@ -2146,7 +2146,6 @@ namespace Nektar
}
else if(mkey.ConstFactorExists(eFactorSVVDGKerDiffCoeff)) // Rodrigo/Mansoor's DG Kernel
{
NekDouble cutoff = mkey.GetConstFactor(eFactorSVVCutoffRatio);
NekDouble SvvDiffCoeff =
mkey.GetConstFactor(eFactorSVVDGKerDiffCoeff)*
mkey.GetConstFactor(eFactorSVVDiffCoeff);
......
......@@ -1390,7 +1390,6 @@ namespace Nektar
int qb = m_base[1]->GetNumPoints();
int nmodes_a = m_base[0]->GetNumModes();
int nmodes_b = m_base[1]->GetNumModes();
int nmodes = min(nmodes_a,nmodes_b);
// Declare orthogonal basis.
LibUtilities::PointsKey pa(qa,m_base[0]->GetPointsType());
......@@ -1427,7 +1426,6 @@ namespace Nektar
}
else if(mkey.ConstFactorExists(eFactorSVVDGKerDiffCoeff)) // Rodrigo/mansoor's DG kernel
{
NekDouble cutoff = mkey.GetConstFactor(eFactorSVVCutoffRatio);
NekDouble SvvDiffCoeff =
mkey.GetConstFactor(eFactorSVVDGKerDiffCoeff)*
mkey.GetConstFactor(eFactorSVVDiffCoeff);
......
......@@ -36,6 +36,9 @@
#include <iostream>
#include <boost/random/variate_generator.hpp>
#include <boost/random/normal_distribution.hpp>
#include <MultiRegions/AssemblyMap/AssemblyMapDG.h>
#include <APESolver/EquationSystems/APE.h>
#include <LocalRegions/TriExp.h>
......@@ -147,10 +150,10 @@ void APE::v_InitObject()
// Define the normal velocity fields
if (m_fields[0]->GetTrace())
{
m_traceBasefield = Array<OneD, Array<OneD, NekDouble> > (m_spacedim + 2);
for (int i = 0; i < m_spacedim + 2; i++)
m_bfTrace = Array<OneD, Array<OneD, NekDouble> > (m_spacedim + 2);
for (int i = 0; i < m_spacedim + 2; ++i)
{
m_traceBasefield[i] = Array<OneD, NekDouble>(GetTraceNpoints());
m_bfTrace[i] = Array<OneD, NekDouble>(GetTraceNpoints(), 0.0);
}
}
......@@ -168,9 +171,9 @@ void APE::v_InitObject()
m_riemannSolver = SolverUtils::GetRiemannSolverFactory().CreateInstance(
riemName);
m_riemannSolver->SetVector("N", &APE::GetNormals, this);
m_riemannSolver->SetVector("basefield", &APE::GetBasefield, this);
m_riemannSolver->SetVector("basefield", &APE::GetBfTrace, this);
m_riemannSolver->SetAuxVec("vecLocs", &APE::GetVecLocs, this);
m_riemannSolver->SetParam("Gamma", &APE::GetGamma, this);
m_riemannSolver->SetParam("Gamma", &APE::GetGamma, this);
// Set up advection operator
string advName;
......@@ -190,6 +193,9 @@ void APE::v_InitObject()
{
ASSERTL0(false, "Implicit APE not set up.");
}
m_whiteNoiseBC_lastUpdate = -1.0;
m_whiteNoiseBC_p = 0.0;
}
......@@ -364,18 +370,45 @@ void APE::SetBoundaryConditions(Array<OneD, Array<OneD, NekDouble> > &inarray,
Fwd[i] = Array<OneD, NekDouble>(nTracePts);
m_fields[i]->ExtractTracePhys(inarray[i], Fwd[i]);
}
Array<OneD, Array<OneD, NekDouble> > bfFwd = GetBfTrace();
// loop over Boundary Regions
for(int n = 0; n < m_fields[0]->GetBndConditions().num_elements(); ++n)
for (int n = 0; n < m_fields[0]->GetBndConditions().num_elements(); ++n)
{
// Wall Boundary Condition
if (boost::iequals(m_fields[0]->GetBndConditions()[n]->GetUserDefined(),"Wall"))
std::string userDefStr =
m_fields[0]->GetBndConditions()[n]->GetUserDefined();
if (!userDefStr.empty())
{
WallBC(n, cnt, Fwd, inarray);
// Wall Boundary Condition
if (boost::iequals(userDefStr, "Wall"))
{
WallBC(n, cnt, Fwd, inarray);
}
else if (boost::iequals(userDefStr, "WhiteNoise"))
{
WhiteNoiseBC(n, cnt, Fwd, bfFwd, inarray);
}
else if (boost::iequals(userDefStr, "RiemannInvariantBC"))
{
RiemannInvariantBC(n, cnt, Fwd, bfFwd, inarray);
}
else if (boost::iequals(userDefStr, "TimeDependent"))
{
for (int i = 0; i < nvariables; ++i)
{
varName = m_session->GetVariable(i);
m_fields[i]->EvaluateBoundaryConditions(time, varName);
}
}
else
{
string errmsg = "Unrecognised boundary condition: ";
errmsg += userDefStr;
ASSERTL0(false, errmsg.c_str());
}
}
// Time Dependent Boundary Condition (specified in meshfile)
if (m_fields[0]->GetBndConditions()[n]->IsTimeDependent())
else
{
for (int i = 0; i < nvariables; ++i)
{
......@@ -383,7 +416,8 @@ void APE::SetBoundaryConditions(Array<OneD, Array<OneD, NekDouble> > &inarray,
m_fields[i]->EvaluateBoundaryConditions(time, varName);
}
}
cnt +=m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
cnt += m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
}
}
......@@ -448,6 +482,200 @@ void APE::WallBC(int bcRegion, int cnt,
}
/**
* @brief Outflow characteristic boundary conditions for compressible
* flow problems.
*/
void APE::RiemannInvariantBC(int bcRegion,
int cnt,
Array<OneD, Array<OneD, NekDouble> > &Fwd,
Array<OneD, Array<OneD, NekDouble> > &BfFwd,
Array<OneD, Array<OneD, NekDouble> > &physarray)
{
int id1, id2, nBCEdgePts;
int nVariables = physarray.num_elements();
const Array<OneD, const int> &traceBndMap = m_fields[0]->GetTraceBndMap();
int eMax = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
for (int e = 0; e < eMax; ++e)
{
nBCEdgePts = m_fields[0]
->GetBndCondExpansions()[bcRegion]
->GetExp(e)
->GetTotPoints();
id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetPhys_Offset(e);
id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt + e]);
// Calculate (v.n)
Array<OneD, NekDouble> Vn(nBCEdgePts, 0.0);
for (int i = 0; i < m_spacedim; ++i)
{
Vmath::Vvtvp(nBCEdgePts,
&Fwd[1 + i][id2], 1,
&m_traceNormals[i][id2], 1,
&Vn[0], 1,
&Vn[0], 1);
}
// Calculate (v0.n)
Array<OneD, NekDouble> Vn0(nBCEdgePts, 0.0);
for (int i = 0; i < m_spacedim; ++i)
{
Vmath::Vvtvp(nBCEdgePts,
&BfFwd[2 + i][id2], 1,
&m_traceNormals[i][id2], 1,
&Vn0[0], 1,
&Vn0[0], 1);
}
for (int i = 0; i < nBCEdgePts; ++i)
{
NekDouble c = sqrt(m_gamma * BfFwd[0][id2 + i] / BfFwd[1][id2 + i]);
NekDouble l0 = Vn0[i] + c;
NekDouble l1 = Vn0[i] - c;
NekDouble h0, h1;