Commit fba90cb5 authored by Gianmarco Mengaldo's avatar Gianmarco Mengaldo
Browse files

Added a new basis called Gauss_Lagrange.


git-svn-id: https://gforge.sci.utah.edu/svn/nektar/trunk@4007 305cdda6-5ce1-45b3-a98d-dfc68c8b3305
parent f79bf1e3
......@@ -524,14 +524,14 @@ namespace Nektar
}
// define derivative basis
Blas::Dgemm('n','n',numPoints,numModes,numPoints,1.0,D,
numPoints, m_bdata.data(),numPoints,0.0,
m_dbdata.data(),numPoints);
Blas::Dgemm('n', 'n', numPoints, numModes, numPoints, 1.0,
D, numPoints, m_bdata.data(), numPoints, 0.0,
m_dbdata.data(), numPoints);
}//end scope
break;
case eG_Lagrange:
case eGauss_Lagrange:
{
mode = m_bdata.data();
boost::shared_ptr< Points<NekDouble> > m_points = PointsManager()[PointsKey(numModes, eGaussGaussLegendre)];
......@@ -546,9 +546,9 @@ namespace Nektar
}
// define derivative basis
Blas::Dgemm('n','n',numPoints,numModes,numPoints,1.0,D,
numPoints, m_bdata.data(),numPoints,0.0,
m_dbdata.data(),numPoints);
Blas::Dgemm('n', 'n', numPoints, numModes, numPoints, 1.0,
D, numPoints, m_bdata.data(), numPoints, 0.0,
m_dbdata.data(), numPoints);
}//end scope
break;
......@@ -622,43 +622,45 @@ namespace Nektar
case eChebyshev:
mode = m_bdata.data();
for (p=0,scal = 1; p<numModes; ++p,mode += numPoints)
{
Polylib::jacobfd(numPoints, z.data(), mode, NULL, p, -0.5, -0.5);
mode = m_bdata.data();
for(i = 0; i < numPoints; ++i)
for (p=0,scal = 1; p<numModes; ++p,mode += numPoints)
{
mode[i] *= scal;
Polylib::jacobfd(numPoints, z.data(), mode, NULL, p, -0.5, -0.5);
for(i = 0; i < numPoints; ++i)
{
mode[i] *= scal;
}
scal *= 4*(p+1)*(p+1)/(2*p+2)/(2*p+1);
}
scal *= 4*(p+1)*(p+1)/(2*p+2)/(2*p+1);
// Define derivative basis
Blas::Dgemm('n', 'n', numPoints, numModes, numPoints, 1.0,
D, numPoints, m_bdata.data(), numPoints, 0.0,
m_dbdata.data(), numPoints);
}
// define derivative basis
Blas::Dgemm('n','n',numPoints,numModes,numPoints,1.0,D,
numPoints, m_bdata.data(),numPoints,0.0,m_dbdata.data(),
numPoints);
break;
case eMonomial:
{
int P = numModes - 1;
NekDouble *mode = m_bdata.data();
int P = numModes - 1;
NekDouble *mode = m_bdata.data();
for( int p = 0; p <= P; ++p, mode += numPoints )
{
for( int p = 0; p <= P; ++p, mode += numPoints )
{
for( int i = 0; i < numPoints; ++i )
{
mode[i] = pow(z[i], p);
}
}
}
// define derivative basis
Blas::Dgemm('n','n',numPoints,numModes,numPoints,1.0,D,numPoints,
m_bdata.data(),numPoints,0.0,m_dbdata.data(),numPoints);
Blas::Dgemm('n', 'n', numPoints, numModes, numPoints, 1.0,
D, numPoints, m_bdata.data(), numPoints, 0.0,
m_dbdata.data(),numPoints);
}//end scope
break;
......@@ -968,7 +970,8 @@ namespace Nektar
{
return ( m_basistype == eGLL_Lagrange &&
GetPointsType() == eGaussLobattoLegendre &&
GetNumModes() == GetNumPoints() || m_basistype == eG_Lagrange);
GetNumModes() == GetNumPoints() ||
m_basistype == eGauss_Lagrange);
}
// BasisKey compared to BasisKey
......
......@@ -106,7 +106,7 @@ namespace Nektar
case eModified_A:
case eFourier:
case eGLL_Lagrange:
case eG_Lagrange:
case eGauss_Lagrange:
case eLegendre:
case eChebyshev:
case eMonomial:
......
......@@ -59,7 +59,7 @@ namespace Nektar
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$
eG_Lagrange, //!< Lagrange Polynomials using the Gauss points \f$ h_p(z_i) \f$
eGauss_Lagrange, //!< Lagrange Polynomials using the Gauss points \f$ h_p(z_i) \f$
eDG_DG_Left, //!< Derivative of the left correction function for DG FR \f$ dGL_{p}(z_i) \f$
eDG_DG_Right, //!< Derivative of the Right correction function for DG FR \f$ dGR_{p}(z_i) \f$
eDG_SD_Left, //!< Derivative of the left correction function for SD FR \f$ dGL_{p}(z_i) \f$
......@@ -86,7 +86,7 @@ namespace Nektar
"FourierSingleMode",
"FourierHalfModeRe",
"FourierHalfModeIm",
"G_Lagrange",
"Gauss_Lagrange",
"DG_DG_Left",
"DG_DG_Right",
"DG_SD_Left",
......
......@@ -104,21 +104,21 @@ namespace Nektar
Polylib::zwgrjm(m_points[0].data(),m_weights.data(),numpoints,2.0,0.0);
break;
case eGaussKronrodLegendre:
Polylib::zwgk(m_points[0].data(),m_weights.data(),numpoints,0.0,0.0);
break;
case eGaussKronrodLegendre:
Polylib::zwgk(m_points[0].data(),m_weights.data(),numpoints,0.0,0.0);
break;
case eGaussRadauKronrodMLegendre:
Polylib::zwrk(m_points[0].data(),m_weights.data(),numpoints,0.0,0.0);
break;
Polylib::zwrk(m_points[0].data(),m_weights.data(),numpoints,0.0,0.0);
break;
case eGaussRadauKronrodMAlpha1Beta0:
Polylib::zwrk(m_points[0].data(),m_weights.data(),numpoints,1.0,0.0);
break;
Polylib::zwrk(m_points[0].data(),m_weights.data(),numpoints,1.0,0.0);
break;
case eGaussLobattoKronrodLegendre:
Polylib::zwlk(m_points[0].data(),m_weights.data(),numpoints,0.0,0.0);
break;
Polylib::zwlk(m_points[0].data(),m_weights.data(),numpoints,0.0,0.0);
break;
default:
ASSERTL0(false, "Unknown Gauss quadrature point distribution requested");
......@@ -127,12 +127,12 @@ namespace Nektar
void GaussPoints::CalculateWeights()
{
//For Gauss Quadrature, this is done as part of the points computation
// For Gauss Quadrature, this is done as part of the points computation
}
void GaussPoints::CalculateDerivMatrix()
{
// Allocate the derivative matrix.
// Allocate the derivative matrix
PointsBaseType::CalculateDerivMatrix();
int numpoints = m_pointsKey.GetNumPoints();
......@@ -190,25 +190,24 @@ namespace Nektar
Polylib::Dgrjm(dmtemp,m_points[0].data(),numpoints,2.0,0.0);
break;
case eGaussKronrodLegendre:
case eGaussRadauKronrodMLegendre:
case eGaussRadauKronrodMAlpha1Beta0:
case eGaussLobattoKronrodLegendre:
{
for(unsigned int i=0;i<m_pointsKey.GetNumPoints();++i)
{
for(unsigned int j=0;j<m_pointsKey.GetNumPoints();++j)
{
(*m_derivmatrix[0])(i,j) = LagrangePolyDeriv(m_points[0][i],j,m_pointsKey.GetNumPoints(),m_points[0]);
}
}
return;
}
break;
case eGaussKronrodLegendre:
case eGaussRadauKronrodMLegendre:
case eGaussRadauKronrodMAlpha1Beta0:
case eGaussLobattoKronrodLegendre:
{
for(unsigned int i=0;i<m_pointsKey.GetNumPoints();++i)
{
for(unsigned int j=0;j<m_pointsKey.GetNumPoints();++j)
{
(*m_derivmatrix[0])(i,j) = LagrangePolyDeriv(m_points[0][i],j,m_pointsKey.GetNumPoints(),m_points[0]);
}
}
return;
}
break;
default:
ASSERTL0(false, "Unknown Gauss quadrature point distribution requested");
ASSERTL0(false, "Unknown Gauss quadrature point distribution requested");
}
std::copy(dmtemp,dmtemp+totpoints*totpoints,m_derivmatrix[0]->begin());
......@@ -218,68 +217,68 @@ namespace Nektar
{
switch(m_pointsKey.GetPointsType())
{
case eGaussGaussLegendre:
Polylib::Imgj(interp.data(),m_points[0].data(),xpoints.data(),GetNumPoints(),npts,0.0,0.0);
break;
case eGaussGaussLegendre:
Polylib::Imgj(interp.data(),m_points[0].data(),xpoints.data(),GetNumPoints(),npts,0.0,0.0);
break;
case eGaussRadauMLegendre:
Polylib::Imgrjm(interp.data(),m_points[0].data(),xpoints.data(),GetNumPoints(),npts,0.0,0.0);
break;
case eGaussRadauMLegendre:
Polylib::Imgrjm(interp.data(),m_points[0].data(),xpoints.data(),GetNumPoints(),npts,0.0,0.0);
break;
case eGaussRadauPLegendre:
Polylib::Imgrjp(interp.data(),m_points[0].data(),xpoints.data(),GetNumPoints(),npts,0.0,0.0);
break;
case eGaussRadauPLegendre:
Polylib::Imgrjp(interp.data(),m_points[0].data(),xpoints.data(),GetNumPoints(),npts,0.0,0.0);
break;
case eGaussLobattoLegendre:
Polylib::Imglj(interp.data(),m_points[0].data(),xpoints.data(),GetNumPoints(),npts,0.0,0.0);
break;
case eGaussLobattoLegendre:
Polylib::Imglj(interp.data(),m_points[0].data(),xpoints.data(),GetNumPoints(),npts,0.0,0.0);
break;
case eGaussGaussChebyshev:
Polylib::Imgj(interp.data(),m_points[0].data(),xpoints.data(),GetNumPoints(),npts,-0.5,-0.5);
break;
case eGaussGaussChebyshev:
Polylib::Imgj(interp.data(),m_points[0].data(),xpoints.data(),GetNumPoints(),npts,-0.5,-0.5);
break;
case eGaussRadauMChebyshev:
Polylib::Imgrjm(interp.data(),m_points[0].data(),xpoints.data(),GetNumPoints(),npts,-0.5,-0.5);
break;
case eGaussRadauMChebyshev:
Polylib::Imgrjm(interp.data(),m_points[0].data(),xpoints.data(),GetNumPoints(),npts,-0.5,-0.5);
break;
case eGaussRadauPChebyshev:
Polylib::Imgrjp(interp.data(),m_points[0].data(),xpoints.data(),GetNumPoints(),npts,-0.5,-0.5);
break;
case eGaussRadauPChebyshev:
Polylib::Imgrjp(interp.data(),m_points[0].data(),xpoints.data(),GetNumPoints(),npts,-0.5,-0.5);
break;
case eGaussLobattoChebyshev:
Polylib::Imglj(interp.data(),m_points[0].data(),xpoints.data(),GetNumPoints(),npts,-0.5,-0.5);
break;
case eGaussLobattoChebyshev:
Polylib::Imglj(interp.data(),m_points[0].data(),xpoints.data(),GetNumPoints(),npts,-0.5,-0.5);
break;
case eGaussRadauMAlpha0Beta1:
Polylib::Imgrjm(interp.data(),m_points[0].data(),xpoints.data(),GetNumPoints(),npts,0.0,1.0);
break;
case eGaussRadauMAlpha0Beta1:
Polylib::Imgrjm(interp.data(),m_points[0].data(),xpoints.data(),GetNumPoints(),npts,0.0,1.0);
break;
case eGaussRadauMAlpha0Beta2:
Polylib::Imgrjm(interp.data(),m_points[0].data(),xpoints.data(),GetNumPoints(),npts,0.0,2.0);
break;
case eGaussRadauMAlpha0Beta2:
Polylib::Imgrjm(interp.data(),m_points[0].data(),xpoints.data(),GetNumPoints(),npts,0.0,2.0);
break;
case eGaussRadauMAlpha1Beta0:
Polylib::Imgrjm(interp.data(),m_points[0].data(),xpoints.data(),GetNumPoints(),npts,1.0,0.0);
break;
case eGaussRadauMAlpha1Beta0:
Polylib::Imgrjm(interp.data(),m_points[0].data(),xpoints.data(),GetNumPoints(),npts,1.0,0.0);
break;
case eGaussRadauMAlpha2Beta0:
Polylib::Imgrjm(interp.data(),m_points[0].data(),xpoints.data(),GetNumPoints(),npts,2.0,0.0);
break;
case eGaussRadauMAlpha2Beta0:
Polylib::Imgrjm(interp.data(),m_points[0].data(),xpoints.data(),GetNumPoints(),npts,2.0,0.0);
break;
case eGaussKronrodLegendre:
case eGaussRadauKronrodMLegendre:
case eGaussRadauKronrodMAlpha1Beta0:
case eGaussLobattoKronrodLegendre:
{
for(unsigned int i=0;i<npts;++i)
{
for(unsigned int j=0;j<m_pointsKey.GetNumPoints();++j)
{
interp[i + j*npts] = LagrangePoly(xpoints[i],j,m_pointsKey.GetNumPoints(),m_points[0]);
}
}
}
break;
case eGaussKronrodLegendre:
case eGaussRadauKronrodMLegendre:
case eGaussRadauKronrodMAlpha1Beta0:
case eGaussLobattoKronrodLegendre:
{
for(unsigned int i=0;i<npts;++i)
{
for(unsigned int j=0;j<m_pointsKey.GetNumPoints();++j)
{
interp[i + j*npts] = LagrangePoly(xpoints[i],j,m_pointsKey.GetNumPoints(),m_points[0]);
}
}
}
break;
default:
ASSERTL0(false, "Unknown Gauss quadrature point distribution requested");
......@@ -302,7 +301,7 @@ namespace Nektar
PointsManager()[pkey]->GetPoints(xpoints);
/// Delegate to function below.
// Delegate to function below
return GetI(numpoints, xpoints);
}
......@@ -317,7 +316,7 @@ namespace Nektar
{
int numpoints = 1;
/// Delegate to function below.
// Delegate to function below
return GetI(numpoints, x);
}
......
......@@ -1798,7 +1798,7 @@ namespace Nektar
/**
* If one wants to get hold of the underlying data without modifying
* them, rather use the function #GetCoeffs insead.
* them, rather use the function #GetCoeffs instead.
*
* @return (A reference to) the array #m_coeffs.
*/
......
This diff is collapsed.
......@@ -67,6 +67,7 @@ namespace Nektar
eOrthogonal,
eGLL_Lagrange,
eGLL_Lagrange_SEM,
eGauss_Lagrange,
eFourier,
eFourierSingleMode,
eFourierHalfModeRe,
......@@ -87,6 +88,7 @@ namespace Nektar
"ORTHOGONAL",
"GLL_LAGRANGE",
"GLL_LAGRANGE_SEM",
"Gauss_LAGRANGE",
"FOURIER",
"FOURIERSINGLEMODE",
"FOURIERHALFMODERE",
......
......@@ -670,7 +670,7 @@ namespace Nektar
switch(Btype)
{
case LibUtilities::eGLL_Lagrange:
case LibUtilities::eG_Lagrange:
case LibUtilities::eGauss_Lagrange:
case LibUtilities::eChebyshev:
case LibUtilities::eFourier:
outarray[1]= nummodes-1;
......@@ -697,7 +697,7 @@ namespace Nektar
switch(Btype)
{
case LibUtilities::eGLL_Lagrange:
case LibUtilities::eG_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