Commit 08b211b5 authored by Dave Moxey's avatar Dave Moxey

Merge branch 'fix/DG-P0' into 'master'

Fix compressible solver with NUMMODES=1

See merge request nektar/nektar!868
parents 158ebf21 0fd6e75d
...@@ -104,6 +104,7 @@ v5.0.0 ...@@ -104,6 +104,7 @@ v5.0.0
- Allow performing axi-symmetric Euler and NS simulations (!771, !866) - Allow performing axi-symmetric Euler and NS simulations (!771, !866)
- Add ability to use an exponential filtering for stabilization with - Add ability to use an exponential filtering for stabilization with
seg, quad and hex elements (!771, !862) seg, quad and hex elements (!771, !862)
- Fix compressible solver with NUMMODES=1 (!868)
- Introduce equations of state to account for real gas effects (!880) - Introduce equations of state to account for real gas effects (!880)
**APESolver:** **APESolver:**
......
...@@ -196,6 +196,13 @@ namespace Polylib { ...@@ -196,6 +196,13 @@ namespace Polylib {
z[0] = 0.0; z[0] = 0.0;
w[0] = 2.0; w[0] = 2.0;
} }
else if( np == 2 ){
z[0] = -1.0;
z[1] = 1.0;
w[0] = 1.0;
w[1] = 1.0;
}
else{ else{
register int i; register int i;
double fac, one = 1.0, apb = alpha + beta, two = 2.0; double fac, one = 1.0, apb = alpha + beta, two = 2.0;
...@@ -678,8 +685,7 @@ namespace Polylib { ...@@ -678,8 +685,7 @@ namespace Polylib {
void Dglj(double *D, const double *z, const int np, const double alpha, void Dglj(double *D, const double *z, const int np, const double alpha,
const double beta) const double beta)
{ {
if (np <= 1){
if (np <= 0){
D[0] = 0.0; D[0] = 0.0;
} }
else{ else{
......
...@@ -868,7 +868,18 @@ namespace Nektar ...@@ -868,7 +868,18 @@ namespace Nektar
const SpatialDomains::GeomFactorsSharedPtr & geomFactors = const SpatialDomains::GeomFactorsSharedPtr & geomFactors =
GetGeom()->GetMetricInfo(); GetGeom()->GetMetricInfo();
SpatialDomains::GeomType type = geomFactors->GetGtype(); SpatialDomains::GeomType type = geomFactors->GetGtype();
LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys(); LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
for(i = 0; i < ptsKeys.size(); ++i)
{
// Need at least 2 points for computing normals
if (ptsKeys[i].GetNumPoints() == 1)
{
LibUtilities::PointsKey pKey(2, ptsKeys[i].GetPointsType());
ptsKeys[i] = pKey;
}
}
const Array<TwoD, const NekDouble> & df = geomFactors->GetDerivFactors(ptsKeys); const Array<TwoD, const NekDouble> & df = geomFactors->GetDerivFactors(ptsKeys);
const Array<OneD, const NekDouble> & jac = geomFactors->GetJac(ptsKeys); const Array<OneD, const NekDouble> & jac = geomFactors->GetJac(ptsKeys);
...@@ -950,9 +961,9 @@ namespace Nektar ...@@ -950,9 +961,9 @@ namespace Nektar
{ {
int j, k; int j, k;
int nqe0 = m_base[0]->GetNumPoints(); int nqe0 = ptsKeys[0].GetNumPoints();
int nqe1 = m_base[1]->GetNumPoints(); int nqe1 = ptsKeys[0].GetNumPoints();
int nqe2 = m_base[2]->GetNumPoints(); int nqe2 = ptsKeys[0].GetNumPoints();
int nqe01 = nqe0*nqe1; int nqe01 = nqe0*nqe1;
int nqe02 = nqe0*nqe2; int nqe02 = nqe0*nqe2;
int nqe12 = nqe1*nqe2; int nqe12 = nqe1*nqe2;
......
...@@ -703,7 +703,18 @@ namespace Nektar ...@@ -703,7 +703,18 @@ namespace Nektar
{ {
const SpatialDomains::GeomFactorsSharedPtr &geomFactors = const SpatialDomains::GeomFactorsSharedPtr &geomFactors =
GetGeom()->GetMetricInfo(); GetGeom()->GetMetricInfo();
LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys(); LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
for(int i = 0; i < ptsKeys.size(); ++i)
{
// Need at least 2 points for computing normals
if (ptsKeys[i].GetNumPoints() == 1)
{
LibUtilities::PointsKey pKey(2, ptsKeys[i].GetPointsType());
ptsKeys[i] = pKey;
}
}
SpatialDomains::GeomType type = geomFactors->GetGtype(); SpatialDomains::GeomType type = geomFactors->GetGtype();
const Array<TwoD, const NekDouble> &df = geomFactors->GetDerivFactors(ptsKeys); const Array<TwoD, const NekDouble> &df = geomFactors->GetDerivFactors(ptsKeys);
const Array<OneD, const NekDouble> &jac = geomFactors->GetJac(ptsKeys); const Array<OneD, const NekDouble> &jac = geomFactors->GetJac(ptsKeys);
......
...@@ -713,7 +713,18 @@ namespace Nektar ...@@ -713,7 +713,18 @@ namespace Nektar
{ {
const SpatialDomains::GeomFactorsSharedPtr &geomFactors = const SpatialDomains::GeomFactorsSharedPtr &geomFactors =
GetGeom()->GetMetricInfo(); GetGeom()->GetMetricInfo();
LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys(); LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
for(int i = 0; i < ptsKeys.size(); ++i)
{
// Need at least 2 points for computing normals
if (ptsKeys[i].GetNumPoints() == 1)
{
LibUtilities::PointsKey pKey(2, ptsKeys[i].GetPointsType());
ptsKeys[i] = pKey;
}
}
SpatialDomains::GeomType type = geomFactors->GetGtype(); SpatialDomains::GeomType type = geomFactors->GetGtype();
const Array<TwoD, const NekDouble> &df = geomFactors->GetDerivFactors(ptsKeys); const Array<TwoD, const NekDouble> &df = geomFactors->GetDerivFactors(ptsKeys);
const Array<OneD, const NekDouble> &jac = geomFactors->GetJac(ptsKeys); const Array<OneD, const NekDouble> &jac = geomFactors->GetJac(ptsKeys);
......
...@@ -300,6 +300,12 @@ namespace Nektar ...@@ -300,6 +300,12 @@ namespace Nektar
fill(outarray.get(), outarray.get()+m_ncoeffs, 0.0 ); fill(outarray.get(), outarray.get()+m_ncoeffs, 0.0 );
if(nmodes[0] == 1 && nmodes[1] == 1)
{
outarray[0] = inarray[0];
return;
}
Array<OneD, NekDouble> physEdge[4]; Array<OneD, NekDouble> physEdge[4];
Array<OneD, NekDouble> coeffEdge[4]; Array<OneD, NekDouble> coeffEdge[4];
StdRegions::Orientation orient[4]; StdRegions::Orientation orient[4];
...@@ -1194,7 +1200,18 @@ namespace Nektar ...@@ -1194,7 +1200,18 @@ namespace Nektar
const SpatialDomains::GeomFactorsSharedPtr & geomFactors = const SpatialDomains::GeomFactorsSharedPtr & geomFactors =
GetGeom()->GetMetricInfo(); GetGeom()->GetMetricInfo();
SpatialDomains::GeomType type = geomFactors->GetGtype(); SpatialDomains::GeomType type = geomFactors->GetGtype();
LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys(); LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
for(i = 0; i < ptsKeys.size(); ++i)
{
// Need at least 2 points for computing normals
if (ptsKeys[i].GetNumPoints() == 1)
{
LibUtilities::PointsKey pKey(2, ptsKeys[i].GetPointsType());
ptsKeys[i] = pKey;
}
}
const Array<TwoD, const NekDouble> & df = geomFactors->GetDerivFactors(ptsKeys); const Array<TwoD, const NekDouble> & df = geomFactors->GetDerivFactors(ptsKeys);
const Array<OneD, const NekDouble> & jac = geomFactors->GetJac(ptsKeys); const Array<OneD, const NekDouble> & jac = geomFactors->GetJac(ptsKeys);
int nqe; int nqe;
...@@ -1425,7 +1442,7 @@ namespace Nektar ...@@ -1425,7 +1442,7 @@ namespace Nektar
// interpolate Jacobian and invert // interpolate Jacobian and invert
LibUtilities::Interp1D( LibUtilities::Interp1D(
from_key,jac, m_base[0]->GetPointsKey(), work); from_key,jac, m_base[0]->GetPointsKey(), work);
Vmath::Sdiv(nq,1.0,&work[0],1,&work[0],1); Vmath::Sdiv(nqe,1.0,&work[0],1,&work[0],1);
// interpolate // interpolate
for (i = 0; i < GetCoordim(); ++i) for (i = 0; i < GetCoordim(); ++i)
......
...@@ -864,22 +864,28 @@ cout<<"deps/dx ="<<inarray_d0[i]<<" deps/dy="<<inarray_d1[i]<<endl; ...@@ -864,22 +864,28 @@ cout<<"deps/dx ="<<inarray_d0[i]<<" deps/dy="<<inarray_d1[i]<<endl;
{ {
// Assume that input is also Gll_Lagrange // Assume that input is also Gll_Lagrange
// but no way to check; // but no way to check;
LibUtilities::PointsKey p0( LibUtilities::PointsKey f0(
nummodes[mode_offset], nummodes[mode_offset],
LibUtilities::eGaussLobattoLegendre); LibUtilities::eGaussLobattoLegendre);
LibUtilities::PointsKey t0(
m_base[0]->GetNumModes(),
LibUtilities::eGaussLobattoLegendre);
LibUtilities::Interp1D( LibUtilities::Interp1D(
p0,data, m_base[0]->GetPointsKey(), coeffs); f0,data, t0, coeffs);
} }
break; break;
case LibUtilities::eGauss_Lagrange: case LibUtilities::eGauss_Lagrange:
{ {
// Assume that input is also Gauss_Lagrange // Assume that input is also Gauss_Lagrange
// but no way to check; // but no way to check;
LibUtilities::PointsKey p0( LibUtilities::PointsKey f0(
nummodes[mode_offset], nummodes[mode_offset],
LibUtilities::eGaussGaussLegendre); LibUtilities::eGaussGaussLegendre);
LibUtilities::PointsKey t0(
m_base[0]->GetNumModes(),
LibUtilities::eGaussGaussLegendre);
LibUtilities::Interp1D( LibUtilities::Interp1D(
p0,data, m_base[0]->GetPointsKey(), coeffs); f0,data, t0, coeffs);
} }
break; break;
default: default:
......
...@@ -738,7 +738,18 @@ namespace Nektar ...@@ -738,7 +738,18 @@ namespace Nektar
int i; int i;
const SpatialDomains::GeomFactorsSharedPtr &geomFactors = const SpatialDomains::GeomFactorsSharedPtr &geomFactors =
GetGeom()->GetMetricInfo(); GetGeom()->GetMetricInfo();
LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys(); LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
for(int i = 0; i < ptsKeys.size(); ++i)
{
// Need at least 2 points for computing normals
if (ptsKeys[i].GetNumPoints() == 1)
{
LibUtilities::PointsKey pKey(2, ptsKeys[i].GetPointsType());
ptsKeys[i] = pKey;
}
}
SpatialDomains::GeomType type = geomFactors->GetGtype(); SpatialDomains::GeomType type = geomFactors->GetGtype();
const Array<TwoD, const NekDouble> &df = geomFactors->GetDerivFactors(ptsKeys); const Array<TwoD, const NekDouble> &df = geomFactors->GetDerivFactors(ptsKeys);
const Array<OneD, const NekDouble> &jac = geomFactors->GetJac(ptsKeys); const Array<OneD, const NekDouble> &jac = geomFactors->GetJac(ptsKeys);
......
...@@ -270,6 +270,12 @@ namespace Nektar ...@@ -270,6 +270,12 @@ namespace Nektar
fill(outarray.get(), outarray.get()+m_ncoeffs, 0.0 ); fill(outarray.get(), outarray.get()+m_ncoeffs, 0.0 );
if(nmodes[0] == 1 && nmodes[1] == 1)
{
outarray[0] = inarray[0];
return;
}
Array<OneD, NekDouble> physEdge[3]; Array<OneD, NekDouble> physEdge[3];
Array<OneD, NekDouble> coeffEdge[3]; Array<OneD, NekDouble> coeffEdge[3];
for(i = 0; i < 3; i++) for(i = 0; i < 3; i++)
...@@ -828,7 +834,18 @@ namespace Nektar ...@@ -828,7 +834,18 @@ namespace Nektar
{ {
int i; int i;
const SpatialDomains::GeomFactorsSharedPtr & geomFactors = GetGeom()->GetMetricInfo(); const SpatialDomains::GeomFactorsSharedPtr & geomFactors = GetGeom()->GetMetricInfo();
LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys(); LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
for(i = 0; i < ptsKeys.size(); ++i)
{
// Need at least 2 points for computing normals
if (ptsKeys[i].GetNumPoints() == 1)
{
LibUtilities::PointsKey pKey(2, ptsKeys[i].GetPointsType());
ptsKeys[i] = pKey;
}
}
const SpatialDomains::GeomType type = geomFactors->GetGtype(); const SpatialDomains::GeomType type = geomFactors->GetGtype();
const Array<TwoD, const NekDouble> & df = geomFactors->GetDerivFactors(ptsKeys); const Array<TwoD, const NekDouble> & df = geomFactors->GetDerivFactors(ptsKeys);
const Array<OneD, const NekDouble> & jac = geomFactors->GetJac(ptsKeys); const Array<OneD, const NekDouble> & jac = geomFactors->GetJac(ptsKeys);
...@@ -943,7 +960,7 @@ namespace Nektar ...@@ -943,7 +960,7 @@ namespace Nektar
// interpolate Jacobian and invert // interpolate Jacobian and invert
LibUtilities::Interp1D(from_key,jac,m_base[0]->GetPointsKey(),work); LibUtilities::Interp1D(from_key,jac,m_base[0]->GetPointsKey(),work);
Vmath::Sdiv(nq,1.0,&work[0],1,&work[0],1); Vmath::Sdiv(nqe,1.0,&work[0],1,&work[0],1);
// interpolate // interpolate
for(i = 0; i < GetCoordim(); ++i) for(i = 0; i < GetCoordim(); ++i)
......
...@@ -198,9 +198,6 @@ namespace Nektar ...@@ -198,9 +198,6 @@ namespace Nektar
virtual void v_WriteVtkPieceHeader(std::ostream &outfile, int expansion, int istrip); virtual void v_WriteVtkPieceHeader(std::ostream &outfile, int expansion, int istrip);
int m_firstIntEl;
Array<OneD, NekDouble> m_normSign;
}; };
/// Empty ExpList1D object. /// Empty ExpList1D object.
......
...@@ -111,10 +111,11 @@ Array<OneD, NekDouble> AdvectionSystem::GetElmtCFLVals(void) ...@@ -111,10 +111,11 @@ Array<OneD, NekDouble> AdvectionSystem::GetElmtCFLVals(void)
stdVelocity = v_GetMaxStdVelocity(); stdVelocity = v_GetMaxStdVelocity();
Array<OneD, NekDouble> cfl(nelmt, 0.0); Array<OneD, NekDouble> cfl(nelmt, 0.0);
NekDouble order;
for(int el = 0; el < nelmt; ++el) for(int el = 0; el < nelmt; ++el)
{ {
cfl[el] = m_timestep*(stdVelocity[el] * cLambda * order = max(expOrder[el]-1, 1);
(expOrder[el]-1) * (expOrder[el]-1)); cfl[el] = m_timestep*(stdVelocity[el] * cLambda * order * order);
} }
return cfl; return cfl;
......
...@@ -1417,6 +1417,7 @@ void MeshGraph::SetExpansionsToEvenlySpacedPoints(int npoints) ...@@ -1417,6 +1417,7 @@ void MeshGraph::SetExpansionsToEvenlySpacedPoints(int npoints)
{ {
npts = bkeyold.GetNumModes(); npts = bkeyold.GetNumModes();
} }
npts = max(npts, 2);
const LibUtilities::PointsKey pkey( const LibUtilities::PointsKey pkey(
npts, LibUtilities::ePolyEvenlySpaced); npts, LibUtilities::ePolyEvenlySpaced);
......
...@@ -120,7 +120,19 @@ namespace Nektar ...@@ -120,7 +120,19 @@ namespace Nektar
{ {
return false; return false;
} }
if( LibUtilities::ShapeTypeDimMap[lhs.m_shapeType] <
LibUtilities::ShapeTypeDimMap[rhs.m_shapeType])
{
return true;
}
if( LibUtilities::ShapeTypeDimMap[lhs.m_shapeType] >
LibUtilities::ShapeTypeDimMap[rhs.m_shapeType])
{
return false;
}
for(unsigned int i = 0; i < LibUtilities::ShapeTypeDimMap[lhs.m_shapeType]; ++i) for(unsigned int i = 0; i < LibUtilities::ShapeTypeDimMap[lhs.m_shapeType]; ++i)
{ {
if(lhs.m_base[i].get() < rhs.m_base[i].get()) if(lhs.m_base[i].get() < rhs.m_base[i].get())
......
...@@ -530,16 +530,16 @@ namespace Nektar ...@@ -530,16 +530,16 @@ namespace Nektar
Collections::Collection c(CollExp, impTypes); Collections::Collection c(CollExp, impTypes);
const int nq = Exp->GetTotPoints(); const int nq = Exp->GetTotPoints();
Array<OneD, NekDouble> xc(nq), yc(nq); Array<OneD, NekDouble> xc(nq);
Array<OneD, NekDouble> phys(nelmts*nq),tmp; Array<OneD, NekDouble> phys(nelmts*nq),tmp;
Array<OneD, NekDouble> diff1(nelmts*nq); Array<OneD, NekDouble> diff1(nelmts*nq);
Array<OneD, NekDouble> diff2(nelmts*nq); Array<OneD, NekDouble> diff2(nelmts*nq);
Exp->GetCoords(xc, yc); Exp->GetCoords(xc);
for (int i = 0; i < nq; ++i) for (int i = 0; i < nq; ++i)
{ {
phys[i] = sin(xc[i])*cos(yc[i]); phys[i] = sin(xc[i]);
} }
Exp->PhysDeriv(phys, diff1); Exp->PhysDeriv(phys, diff1);
...@@ -586,12 +586,12 @@ namespace Nektar ...@@ -586,12 +586,12 @@ namespace Nektar
Collections::Collection c(CollExp, impTypes); Collections::Collection c(CollExp, impTypes);
const int nq = Exp->GetTotPoints(); const int nq = Exp->GetTotPoints();
Array<OneD, NekDouble> xc(nq),yc(nq); Array<OneD, NekDouble> xc(nq);
Array<OneD, NekDouble> phys(nq),tmp; Array<OneD, NekDouble> phys(nq),tmp;
Array<OneD, NekDouble> diff1(nq); Array<OneD, NekDouble> diff1(nq);
Array<OneD, NekDouble> diff2(nq); Array<OneD, NekDouble> diff2(nq);
Exp->GetCoords(xc,yc); Exp->GetCoords(xc);
for (int i = 0; i < nq; ++i) for (int i = 0; i < nq; ++i)
{ {
...@@ -639,16 +639,16 @@ namespace Nektar ...@@ -639,16 +639,16 @@ namespace Nektar
Collections::Collection c(CollExp, impTypes); Collections::Collection c(CollExp, impTypes);
const int nq = Exp->GetTotPoints(); const int nq = Exp->GetTotPoints();
Array<OneD, NekDouble> xc(nq), yc(nq); Array<OneD, NekDouble> xc(nq);
Array<OneD, NekDouble> phys(nelmts*nq),tmp; Array<OneD, NekDouble> phys(nelmts*nq),tmp;
Array<OneD, NekDouble> diff1(nelmts*nq); Array<OneD, NekDouble> diff1(nelmts*nq);
Array<OneD, NekDouble> diff2(nelmts*nq); Array<OneD, NekDouble> diff2(nelmts*nq);
Exp->GetCoords(xc, yc); Exp->GetCoords(xc);
for (int i = 0; i < nq; ++i) for (int i = 0; i < nq; ++i)
{ {
phys[i] = sin(xc[i])*cos(yc[i]); phys[i] = sin(xc[i]);
} }
Exp->PhysDeriv(phys, diff1); Exp->PhysDeriv(phys, diff1);
...@@ -699,16 +699,16 @@ namespace Nektar ...@@ -699,16 +699,16 @@ namespace Nektar
Collections::Collection c(CollExp, impTypes); Collections::Collection c(CollExp, impTypes);
const int nq = Exp->GetTotPoints(); const int nq = Exp->GetTotPoints();
Array<OneD, NekDouble> xc(nq), yc(nq); Array<OneD, NekDouble> xc(nq);
Array<OneD, NekDouble> phys(nelmts*nq),tmp; Array<OneD, NekDouble> phys(nelmts*nq),tmp;
Array<OneD, NekDouble> diff1(nelmts*nq); Array<OneD, NekDouble> diff1(nelmts*nq);
Array<OneD, NekDouble> diff2(nelmts*nq); Array<OneD, NekDouble> diff2(nelmts*nq);
Exp->GetCoords(xc, yc); Exp->GetCoords(xc);
for (int i = 0; i < nq; ++i) for (int i = 0; i < nq; ++i)
{ {
phys[i] = sin(xc[i])*cos(yc[i]); phys[i] = sin(xc[i]);
} }
Exp->PhysDeriv(phys, diff1); Exp->PhysDeriv(phys, diff1);
...@@ -760,13 +760,13 @@ namespace Nektar ...@@ -760,13 +760,13 @@ namespace Nektar
Array<OneD, NekDouble> coeffs1(nm); Array<OneD, NekDouble> coeffs1(nm);
Array<OneD, NekDouble> coeffs2(nm); Array<OneD, NekDouble> coeffs2(nm);
Array<OneD, NekDouble> xc(nq), yc(nq); Array<OneD, NekDouble> xc(nq);
Exp->GetCoords(xc, yc); Exp->GetCoords(xc);
for (int i = 0; i < nq; ++i) for (int i = 0; i < nq; ++i)
{ {
phys1[i] = sin(xc[i])*cos(yc[i]); phys1[i] = sin(xc[i]);
} }
// Standard routines // Standard routines
...@@ -816,16 +816,16 @@ namespace Nektar ...@@ -816,16 +816,16 @@ namespace Nektar
const int nq = Exp->GetTotPoints(); const int nq = Exp->GetTotPoints();
const int nm = Exp->GetNcoeffs(); const int nm = Exp->GetNcoeffs();
Array<OneD, NekDouble> xc(nq), yc(nq),tmp,tmp1; Array<OneD, NekDouble> xc(nq),tmp,tmp1;
Array<OneD, NekDouble> phys1(nelmts*nq); Array<OneD, NekDouble> phys1(nelmts*nq);
Array<OneD, NekDouble> coeffs1(nelmts*nm); Array<OneD, NekDouble> coeffs1(nelmts*nm);
Array<OneD, NekDouble> coeffs2(nelmts*nm); Array<OneD, NekDouble> coeffs2(nelmts*nm);
Exp->GetCoords(xc, yc); Exp->GetCoords(xc);
for (int i = 0; i < nq; ++i) for (int i = 0; i < nq; ++i)
{ {
phys1[i] = sin(xc[i])*cos(yc[i]); phys1[i] = sin(xc[i]);
} }
Exp->IProductWRTDerivBase(0, phys1, coeffs1); Exp->IProductWRTDerivBase(0, phys1, coeffs1);
...@@ -879,13 +879,13 @@ namespace Nektar ...@@ -879,13 +879,13 @@ namespace Nektar
Array<OneD, NekDouble> coeffs1(nm); Array<OneD, NekDouble> coeffs1(nm);
Array<OneD, NekDouble> coeffs2(nm); Array<OneD, NekDouble> coeffs2(nm);
Array<OneD, NekDouble> xc(nq), yc(nq); Array<OneD, NekDouble> xc(nq);
Exp->GetCoords(xc, yc); Exp->GetCoords(xc);
for (int i = 0; i < nq; ++i) for (int i = 0; i < nq; ++i)
{ {
phys1[i] = sin(xc[i])*cos(yc[i]); phys1[i] = sin(xc[i]);
} }
// Standard routines // Standard routines
...@@ -935,16 +935,16 @@ namespace Nektar ...@@ -935,16 +935,16 @@ namespace Nektar
const int nq = Exp->GetTotPoints(); const int nq = Exp->GetTotPoints();
const int nm = Exp->GetNcoeffs(); const int nm = Exp->GetNcoeffs();
Array<OneD, NekDouble> xc(nq), yc(nq),tmp,tmp1; Array<OneD, NekDouble> xc(nq),tmp,tmp1;
Array<OneD, NekDouble> phys1(nelmts*nq); Array<OneD, NekDouble> phys1(nelmts*nq);
Array<OneD, NekDouble> coeffs1(nelmts*nm); Array<OneD, NekDouble> coeffs1(nelmts*nm);
Array<OneD, NekDouble> coeffs2(nelmts*nm); Array<OneD, NekDouble> coeffs2(nelmts*nm);
Exp->GetCoords(xc, yc); Exp->GetCoords(xc);
for (int i = 0; i < nq; ++i) for (int i = 0; i < nq; ++i)
{ {
phys1[i] = sin(xc[i])*cos(yc[i]); phys1[i] = sin(xc[i]);
} }
Exp->IProductWRTDerivBase(0, phys1, coeffs1); Exp->IProductWRTDerivBase(0, phys1, coeffs1);
...@@ -998,13 +998,13 @@ namespace Nektar ...@@ -998,13 +998,13 @@ namespace Nektar
Array<OneD, NekDouble> coeffs1(nm); Array<OneD, NekDouble> coeffs1(nm);
Array<OneD, NekDouble> coeffs2(nm); Array<OneD, NekDouble> coeffs2(nm);
Array<OneD, NekDouble> xc(nq), yc(nq); Array<OneD, NekDouble> xc(nq);
Exp->GetCoords(xc, yc); Exp->GetCoords(xc);
for (int i = 0; i < nq; ++i) for (int i = 0; i < nq; ++i)
{ {
phys1[i] = sin(xc[i])*cos(yc[i]); phys1[i] = sin(xc[i]);
} }
// Standard routines // Standard routines
...@@ -1053,16 +1053,16 @@ namespace Nektar ...@@ -1053,16 +1053,16 @@ namespace Nektar
const int nq = Exp->GetTotPoints(); const int nq = Exp->GetTotPoints();
const int nm = Exp->GetNcoeffs(); const int nm = Exp->GetNcoeffs();
Array<OneD, NekDouble> xc(nq), yc(nq),tmp,tmp1; Array<OneD, NekDouble> xc(nq),tmp,tmp1;
Array<OneD, NekDouble> phys1(nelmts*nq); Array<OneD, NekDouble> phys1(nelmts*nq);
Array<OneD, NekDouble> coeffs1(nelmts*nm); Array<OneD, NekDouble> coeffs1(nelmts*nm);
Array<OneD, NekDouble> coeffs2(nelmts*nm); Array<OneD, NekDouble> coeffs2(nelmts*nm);
Exp->GetCoords(xc, yc); Exp->GetCoords(xc);
for (int i = 0; i < nq; ++i) for (int i = 0; i < nq; ++i)
{ {
phys1[i] = sin(xc[i])*cos(yc[i]); phys1[i] = sin(xc[i]);
} }
Exp->IProductWRTDerivBase(0, phys1, coeffs1); Exp->IProductWRTDerivBase(0, phys1, coeffs1);
......
...@@ -67,6 +67,7 @@ IF( NEKTAR_SOLVER_COMPRESSIBLE_FLOW ) ...@@ -67,6 +67,7 @@ IF( NEKTAR_SOLVER_COMPRESSIBLE_FLOW )
ADD_NEKTAR_TEST(CylinderSubsonic_P3) ADD_NEKTAR_TEST(CylinderSubsonic_P3)
ADD_NEKTAR_TEST(CylinderSubsonic_P8 LENGTHY) ADD_NEKTAR_TEST(CylinderSubsonic_P8 LENGTHY)
ADD_NEKTAR_TEST(Euler1D) ADD_NEKTAR_TEST(Euler1D)
ADD_NEKTAR_TEST(IsentropicVortex16_P1)