Commit e036d9f3 authored by Gianmarco Mengaldo's avatar Gianmarco Mengaldo

GAUSS_LAGRANGE and GAUSS_LAGRANGE_SEM basis added.

parent ba1b721f
......@@ -692,14 +692,17 @@ namespace Nektar
}
/** \brief Determine if basis has collocation property,
* i.e. GLL_Lagrange with Lobatto integration of appropriate order.
*/
* i.e. GLL_Lagrange with Lobatto integration of appropriate order,
* Gauss_Lagrange with Gauss integration of appropriate order.
*/
bool BasisKey::Collocation() const
{
return (m_basistype == eGLL_Lagrange &&
GetPointsType() == eGaussLobattoLegendre &&
GetNumModes() == GetNumPoints()) ||
m_basistype == eGauss_Lagrange;
return ((m_basistype == eGLL_Lagrange &&
GetPointsType() == eGaussLobattoLegendre &&
GetNumModes() == GetNumPoints()) ||
(m_basistype == eGauss_Lagrange &&
GetPointsType() == eGaussGaussLegendre &&
GetNumModes() == GetNumPoints()));
}
// BasisKey compared to BasisKey
......
......@@ -51,13 +51,13 @@ namespace Nektar
eModified_C, //!< Principle Modified Functions \f$ \phi^c_{pqr}(z_i) \f$
eFourier, //!< Fourier Expansion \f$ \exp(i p\pi z_i)\f$
eGLL_Lagrange, //!< Lagrange for SEM basis \f$ h_p(z_i) \f$
eGauss_Lagrange, //!< Lagrange Polynomials using the Gauss points \f$ h_p(z_i) \f$
eLegendre, //!< Legendre Polynomials \f$ L_p(z_i) = P^{0,0}_p(z_i)\f$. Same as Ortho_A
eChebyshev, //!< Chebyshev Polynomials \f$ T_p(z_i) = P^{-1/2,-1/2}_p(z_i)\f$
eMonomial, //!< Monomial polynomials \f$ L_p(z_i) = z_i^{p}\f$
eFourierSingleMode, //!< Fourier ModifiedExpansion with just the first mode \f$ \exp(i \pi z_i)\f$
eFourierHalfModeRe, //!< Fourier Modified expansions with just the real part of the first mode \f$ Re[\exp(i \pi z_i)]\f$
eFourierHalfModeIm, //!< Fourier Modified expansions with just the imaginary part of the first mode \f$ Im[\exp(i \pi z_i)]\f$
eGauss_Lagrange, //!< Lagrange Polynomials using the Gauss points \f$ h_p(z_i) \f$
SIZE_BasisType //!< Length of enum list
};
}
......
......@@ -55,13 +55,13 @@ namespace Nektar
"Modified_C",
"Fourier",
"GLL_Lagrange",
"Gauss_Lagrange",
"Legendre",
"Chebyshev",
"Monomial",
"FourierSingleMode",
"FourierHalfModeRe",
"FourierHalfModeIm",
"Gauss_Lagrange"
"FourierHalfModeIm"
};
......
......@@ -1082,6 +1082,9 @@ namespace Nektar
{
switch(edgeExp->GetBasis(0)->GetBasisType())
{
case LibUtilities::eGauss_Lagrange:
reverse( map.get() , map.get()+order_e);
break;
case LibUtilities::eGLL_Lagrange:
reverse( map.get() , map.get()+order_e);
break;
......
......@@ -982,6 +982,9 @@ namespace Nektar
{
switch(edgeExp->GetBasis(0)->GetBasisType())
{
case LibUtilities::eGauss_Lagrange:
reverse( map.get() , map.get()+order_e);
break;
case LibUtilities::eGLL_Lagrange:
reverse( map.get() , map.get()+order_e);
break;
......
......@@ -413,6 +413,7 @@ cout<<"deps/dx ="<<inarray_d0[i]<<" deps/dy="<<inarray_d1[i]<<endl;
switch(m_base[0]->GetBasisType())
{
case LibUtilities::eGauss_Lagrange:
case LibUtilities::eGLL_Lagrange:
{
offset = 1;
......
......@@ -515,33 +515,45 @@ namespace Nektar
{
switch(locSegExp->GetBasisType(0))
{
case LibUtilities::eModified_A:
// reverse vertex order
m_localToGlobalBndMap[cnt] = gid + 1;
m_localToGlobalBndMap[cnt+1] = gid;
for(k = 2; k < order_e; ++k)
case LibUtilities::eModified_A:
{
m_localToGlobalBndMap[k+cnt] = gid + k;
// reverse vertex order
m_localToGlobalBndMap[cnt] = gid + 1;
m_localToGlobalBndMap[cnt+1] = gid;
for (k = 2; k < order_e; ++k)
{
m_localToGlobalBndMap[k+cnt] = gid + k;
}
// negate odd modes
for(k = 3; k < order_e; k+=2)
{
m_localToGlobalBndSign[cnt+k] = -1.0;
}
break;
}
// negate odd modes
for(k = 3; k < order_e; k+=2)
case LibUtilities::eGLL_Lagrange:
{
// reverse order
for(k = 0; k < order_e; ++k)
{
m_localToGlobalBndMap[cnt+order_e-k-1] = gid + k;
}
break;
}
case LibUtilities::eGauss_Lagrange:
{
m_localToGlobalBndSign[cnt+k] = -1.0;
// reverse order
for(k = 0; k < order_e; ++k)
{
m_localToGlobalBndMap[cnt+order_e-k-1] = gid + k;
}
break;
}
break;
case LibUtilities::eGLL_Lagrange:
// reverse order
for(k = 0; k < order_e; ++k)
default:
{
m_localToGlobalBndMap[cnt+order_e-k-1] = gid + k;
ASSERTL0(false,"Boundary type not permitted");
}
break;
default:
ASSERTL0(false,"Boundary type not permitted");
}
}
......
......@@ -773,7 +773,9 @@ namespace Nektar
Array<OneD, NekDouble> tmp(nLocalSolutionPts, 0.0);
// Partition the field vector in local vectors
tmp = field + (n * nLocalSolutionPts);
Vmath::Vcopy(nLocalSolutionPts,
(&field[phys_offset]), 1,
(&tmp[0]), 1);
// Basis definition on each element
Basis = (*m_exp)[n]->GetBasis(0);
......@@ -794,8 +796,8 @@ namespace Nektar
}
// Set the x-coordinate of the standard interface point
interface_coord[0] = vertex_coord;
interface_coord[0] = -1.0;
// Implementation for every points except Gauss points
if (Basis->GetPointsType() != LibUtilities::eGaussGaussLegendre)
{
......@@ -925,35 +927,33 @@ namespace Nektar
}
}
/// Note this routine changes m_trace->m_coeffs space; ; is the same
/// for point expansion (croth)
void DisContField1D::v_AddTraceIntegral(
const Array<OneD, const NekDouble> &Fn,
Array<OneD, NekDouble> &outarray)
{
int p,n,offset, t_offset;
double vertnorm =0.0;
int p, n, offset, t_offset;
double vertnorm = 0.0;
for(n = 0; n < GetExpSize(); ++n)
for (n = 0; n < GetExpSize(); ++n)
{
offset = GetCoeff_Offset(n);
for(p = 0; p < 2; ++p)
{
vertnorm = 0.0;
for (int i=0; i<((*m_exp)[n]->GetVertexNormal(p)).num_elements(); i++)
for (int i = 0; i<((*m_exp)[n]->GetVertexNormal(p)).num_elements(); i++)
{
vertnorm += ((*m_exp)[n]->GetVertexNormal(p))[i][0];
}
t_offset = GetTrace()->GetPhys_Offset(n+p);
if(vertnorm >= 0.0)
if (vertnorm >= 0.0)
{
outarray[offset+(*m_exp)[n]->GetVertexMap(1)] += Fn[t_offset];
}
if(vertnorm < 0.0)
if (vertnorm < 0.0)
{
outarray[offset] -= Fn[t_offset];
}
......
......@@ -679,12 +679,15 @@ namespace Nektar
// Maybe a different techique for the smoothing require
// implementation for modal basis.
ASSERTL0((*m_exp)[0]->GetBasisType(0)
== LibUtilities::eGLL_Lagrange,
ASSERTL0((*m_exp)[0]->GetBasisType(0)
== LibUtilities::eGLL_Lagrange ||
(*m_exp)[0]->GetBasisType(0)
== LibUtilities::eGauss_Lagrange,
"Smoothing is currently not allowed unless you are using "
"a nodal base for efficiency reasons. The implemented "
"smoothing technique requires the mass matrix inversion "
"which is trivial just for GLL_LAGRANGE_SEM expansions.");
"which is trivial just for GLL_LAGRANGE_SEM and "
"GAUSS_LAGRANGE_SEMexpansions.");
}
......
......@@ -190,7 +190,8 @@ namespace Nektar
LibUtilities::BasisKey TetBc
= expIt->second->m_basisKeyVector[2];
if(TetBa.GetBasisType() == LibUtilities::eGLL_Lagrange)
if(TetBa.GetBasisType() == LibUtilities::eGLL_Lagrange ||
TetBa.GetBasisType() == LibUtilities::eGauss_Lagrange)
{
ASSERTL0(false,"LocalRegions::NodalTetExp is not "
"implemented yet");
......@@ -309,7 +310,8 @@ namespace Nektar
LibUtilities::BasisKey TetBc
= exp->m_basisKeyVector[2];
if(TetBa.GetBasisType() == LibUtilities::eGLL_Lagrange)
if(TetBa.GetBasisType() == LibUtilities::eGLL_Lagrange ||
TetBa.GetBasisType() == LibUtilities::eGauss_Lagrange)
{
ASSERTL0(false,"LocalRegions::NodalTetExp is not "
"implemented yet");
......
......@@ -2352,6 +2352,7 @@ namespace Nektar
(strcmp(LibUtilities::BasisTypeMap[basis[j]], "Modified_B") == 0) ||
(strcmp(LibUtilities::BasisTypeMap[basis[j]], "Modified_C") == 0) ||
(strcmp(LibUtilities::BasisTypeMap[basis[j]], "GLL_Lagrange") == 0) ||
(strcmp(LibUtilities::BasisTypeMap[basis[j]], "Gauss_Lagrange") == 0) ||
(strcmp(LibUtilities::BasisTypeMap[basis[j]], "Fourier") == 0) ||
(strcmp(LibUtilities::BasisTypeMap[basis[j]], "FourierSingleMode") == 0)||
(strcmp(LibUtilities::BasisTypeMap[basis[j]], "FourierHalfModeRe") == 0) ||
......@@ -2984,7 +2985,7 @@ namespace Nektar
break;
default:
{
ASSERTL0(false,"Expansion not defined in switch for this shape");
ASSERTL0(false, "Expansion not defined in switch for this shape");
}
break;
}
......@@ -2997,7 +2998,7 @@ namespace Nektar
{
case eSegment:
{
const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eGaussGaussLegendre);
const LibUtilities::PointsKey pkey(nummodes+1, LibUtilities::eGaussGaussLegendre);
LibUtilities::BasisKey bkey(LibUtilities::eGauss_Lagrange, nummodes, pkey);
returnval.push_back(bkey);
......@@ -3005,7 +3006,7 @@ namespace Nektar
break;
case eQuadrilateral:
{
const LibUtilities::PointsKey pkey(nummodes,LibUtilities::eGaussGaussLegendre);
const LibUtilities::PointsKey pkey(nummodes+1,LibUtilities::eGaussGaussLegendre);
LibUtilities::BasisKey bkey(LibUtilities::eGauss_Lagrange, nummodes, pkey);
returnval.push_back(bkey);
......@@ -3014,7 +3015,7 @@ namespace Nektar
break;
case eHexahedron:
{
const LibUtilities::PointsKey pkey(nummodes,LibUtilities::eGaussGaussLegendre);
const LibUtilities::PointsKey pkey(nummodes+1,LibUtilities::eGaussGaussLegendre);
LibUtilities::BasisKey bkey(LibUtilities::eGauss_Lagrange, nummodes, pkey);
returnval.push_back(bkey);
......@@ -3129,6 +3130,46 @@ namespace Nektar
}
}
break;
case eGauss_Lagrange_SEM:
{
switch (shape)
{
case eSegment:
{
const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eGaussGaussLegendre);
LibUtilities::BasisKey bkey(LibUtilities::eGauss_Lagrange, nummodes, pkey);
returnval.push_back(bkey);
}
break;
case eQuadrilateral:
{
const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eGaussGaussLegendre);
LibUtilities::BasisKey bkey(LibUtilities::eGauss_Lagrange, nummodes, pkey);
returnval.push_back(bkey);
returnval.push_back(bkey);
}
break;
case eHexahedron:
{
const LibUtilities::PointsKey pkey(nummodes, LibUtilities::eGaussGaussLegendre);
LibUtilities::BasisKey bkey(LibUtilities::eGauss_Lagrange, nummodes, pkey);
returnval.push_back(bkey);
returnval.push_back(bkey);
returnval.push_back(bkey);
}
break;
default:
{
ASSERTL0(false, "Expansion not defined in switch for this shape");
}
break;
}
}
break;
case eFourier:
{
......
......@@ -63,6 +63,7 @@ namespace Nektar
eGLL_Lagrange,
eGLL_Lagrange_SEM,
eGauss_Lagrange,
eGauss_Lagrange_SEM,
eFourier,
eFourierSingleMode,
eFourierHalfModeRe,
......@@ -83,7 +84,8 @@ namespace Nektar
"ORTHOGONAL",
"GLL_LAGRANGE",
"GLL_LAGRANGE_SEM",
"Gauss_LAGRANGE",
"GAUSS_LAGRANGE",
"GAUSS_LAGRANGE_SEM",
"FOURIER",
"FOURIERSINGLEMODE",
"FOURIERHALFMODERE",
......
......@@ -715,10 +715,12 @@ namespace Nektar
int StdQuadExp::v_NumDGBndryCoeffs() const
{
ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
GetBasisType(0) == LibUtilities::eGLL_Lagrange,
GetBasisType(0) == LibUtilities::eGLL_Lagrange ||
GetBasisType(0) == LibUtilities::eGauss_Lagrange,
"BasisType is not a boundary interior form");
ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
GetBasisType(1) == LibUtilities::eGLL_Lagrange,
GetBasisType(1) == LibUtilities::eGLL_Lagrange ||
GetBasisType(0) == LibUtilities::eGauss_Lagrange,
"BasisType is not a boundary interior form");
return 2*GetBasisNumModes(0) + 2*GetBasisNumModes(1);
......@@ -1217,7 +1219,8 @@ namespace Nektar
break;
}
}
else if(bType == LibUtilities::eGLL_Lagrange)
else if(bType == LibUtilities::eGLL_Lagrange ||
bType == LibUtilities::eGauss_Lagrange)
{
switch(eid)
{
......
......@@ -370,17 +370,18 @@ namespace Nektar
* \f$ u(\xi_1) \f$
* \param coll_check flag to identify when a Basis->Collocation() call
* should be performed to see if this is a GLL_Lagrange basis with a
* collocation property. (should be set to 0 if taking the inner product
* with respect to the derivative of basis)
* collocation property. (should be set to 0 if taking the inner
* product with respect to the derivative of basis)
* \param outarray the values of the inner product with respect to
* each basis over region will be stored in the array \a outarray as
* output of the function
*/
void StdSegExp::v_IProductWRTBase(const Array<OneD, const NekDouble>& base,
const Array<OneD, const NekDouble>& inarray,
Array<OneD, NekDouble> &outarray,
int coll_check)
void StdSegExp::v_IProductWRTBase(
const Array<OneD, const NekDouble> &base,
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
int coll_check)
{
int nquad = m_base[0]->GetNumPoints();
Array<OneD, NekDouble> tmp(nquad);
......@@ -388,7 +389,8 @@ namespace Nektar
Array<OneD, const NekDouble> w = m_base[0]->GetW();
Vmath::Vmul(nquad, inarray, 1, w, 1, tmp, 1);
/* Comment below was a bug for collocated basis
if(coll_check&&m_base[0]->Collocation())
{
Vmath::Vcopy(nquad, tmp, 1, outarray, 1);
......@@ -397,7 +399,11 @@ namespace Nektar
{
Blas::Dgemv('T',nquad,m_ncoeffs,1.0,base.get(),nquad,
&tmp[0],1,0.0,outarray.get(),1);
}
}*/
// Correct implementation
Blas::Dgemv('T',nquad,m_ncoeffs,1.0,base.get(),nquad,
&tmp[0],1,0.0,outarray.get(),1);
}
/** \brief Inner product of \a inarray over region with respect to the
......@@ -670,7 +676,7 @@ namespace Nektar
switch(Btype)
{
case LibUtilities::eGLL_Lagrange:
case LibUtilities::eGauss_Lagrange:
case LibUtilities::eGauss_Lagrange:
case LibUtilities::eChebyshev:
case LibUtilities::eFourier:
outarray[1]= nummodes-1;
......@@ -697,7 +703,7 @@ namespace Nektar
switch(Btype)
{
case LibUtilities::eGLL_Lagrange:
case LibUtilities::eGauss_Lagrange:
case LibUtilities::eGauss_Lagrange:
case LibUtilities::eChebyshev:
case LibUtilities::eFourier:
for(i = 0 ; i < GetNcoeffs()-2;i++)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment