Nektar++
Public Member Functions | Static Public Member Functions | Private Member Functions | Static Private Attributes | List of all members
Nektar::LibUtilities::GaussPoints Class Reference

#include <GaussPoints.h>

Inheritance diagram for Nektar::LibUtilities::GaussPoints:
[legend]

Public Member Functions

virtual ~GaussPoints ()
 
std::shared_ptr< NekMatrix< NekDouble > > CreateMatrix (const PointsKey &pkey)
 
const std::shared_ptr< NekMatrix< NekDouble > > GetI (const PointsKey &pkey)
 
const std::shared_ptr< NekMatrix< NekDouble > > GetI (const Array< OneD, const NekDouble > &x)
 
const std::shared_ptr< NekMatrix< NekDouble > > GetI (unsigned int numpoints, const Array< OneD, const NekDouble > &x)
 
std::shared_ptr< NekMatrix< NekDouble > > CreateGPMatrix (const PointsKey &pkey)
 
const std::shared_ptr< NekMatrix< NekDouble > > GetGalerkinProjection (const PointsKey &pkey)
 
 GaussPoints (const PointsKey &pkey)
 
- Public Member Functions inherited from Nektar::LibUtilities::Points< NekDouble >
virtual ~Points ()
 
virtual void Initialize (void)
 
unsigned int GetPointsDim () const
 
unsigned int GetNumPoints () const
 
unsigned int GetTotNumPoints () const
 
PointsType GetPointsType () const
 
const Array< OneD, const DataType > & GetZ () const
 
const Array< OneD, const DataType > & GetW () const
 
void GetZW (Array< OneD, const DataType > &z, Array< OneD, const DataType > &w) const
 
const Array< OneD, const NekDouble > & GetBaryWeights () const
 
void GetPoints (Array< OneD, const DataType > &x) const
 
void GetPoints (Array< OneD, const DataType > &x, Array< OneD, const DataType > &y) const
 
void GetPoints (Array< OneD, const DataType > &x, Array< OneD, const DataType > &y, Array< OneD, const DataType > &z) const
 
const MatrixSharedPtrTypeGetD (Direction dir=xDir) const
 
virtual const MatrixSharedPtrType GetI (const Array< OneD, const DataType > &x, const Array< OneD, const DataType > &y)
 
virtual const MatrixSharedPtrType GetI (const Array< OneD, const DataType > &x, const Array< OneD, const DataType > &y, const Array< OneD, const DataType > &z)
 

Static Public Member Functions

static std::shared_ptr< Points< NekDouble > > Create (const PointsKey &pkey)
 

Private Member Functions

 GaussPoints ()
 These should not be called. All creation is done using the constructor requiring the key, declared above. More...
 
 GaussPoints (const GaussPoints &points)
 
void CalculatePoints ()
 
void CalculateWeights ()
 
void CalculateDerivMatrix ()
 
void CalculateInterpMatrix (unsigned int npts, const Array< OneD, const NekDouble > &xpoints, Array< OneD, NekDouble > &interp)
 
std::shared_ptr< NekMatrix< NekDouble > > CalculateGalerkinProjectionMatrix (const PointsKey &pkey)
 
NekDouble LagrangeInterpolant (NekDouble x, int npts, const Array< OneD, const NekDouble > &xpts, const Array< OneD, const NekDouble > &funcvals)
 functions used by the Kronrod points More...
 
NekDouble LagrangePoly (NekDouble x, int pt, int npts, const Array< OneD, const NekDouble > &xpts)
 
NekDouble LagrangePolyDeriv (NekDouble x, int pt, int npts, const Array< OneD, const NekDouble > &xpts)
 

Static Private Attributes

static bool initPointsManager []
 

Additional Inherited Members

- Public Types inherited from Nektar::LibUtilities::Points< NekDouble >
typedef NekDouble DataType
 
typedef std::shared_ptr< NekMatrix< DataType > > MatrixSharedPtrType
 
- Protected Member Functions inherited from Nektar::LibUtilities::Points< NekDouble >
virtual void CalculateBaryWeights ()
 This function calculates the barycentric weights used for enhanced interpolation speed. More...
 
 Points (const PointsKey &key)
 
- Protected Attributes inherited from Nektar::LibUtilities::Points< NekDouble >
PointsKey m_pointsKey
 Points type for this points distributions. More...
 
Array< OneD, DataTypem_points [3]
 Storage for the point locations, allowing for up to a 3D points storage. More...
 
Array< OneD, DataTypem_weights
 Quadrature weights for the weights. More...
 
Array< OneD, DataTypem_bcweights
 Barycentric weights. More...
 
MatrixSharedPtrType m_derivmatrix [3]
 Derivative matrices. More...
 
NekManager< PointsKey, NekMatrix< DataType >, PointsKey::opLessm_InterpManager
 
NekManager< PointsKey, NekMatrix< DataType >, PointsKey::opLessm_GalerkinProjectionManager
 

Detailed Description

Definition at line 47 of file GaussPoints.h.

Constructor & Destructor Documentation

◆ ~GaussPoints()

virtual Nektar::LibUtilities::GaussPoints::~GaussPoints ( )
inlinevirtual

Definition at line 50 of file GaussPoints.h.

51  {
52  }

◆ GaussPoints() [1/3]

Nektar::LibUtilities::GaussPoints::GaussPoints ( const PointsKey pkey)
inline

Definition at line 73 of file GaussPoints.h.

73  : PointsBaseType(pkey)
74  {
75  namespace pl = std::placeholders;
76  m_InterpManager.RegisterCreator(
77  PointsKey(0, eGaussGaussLegendre),
78  std::bind(&GaussPoints::CreateMatrix, this, pl::_1));
79  m_InterpManager.RegisterCreator(
80  PointsKey(0, eGaussRadauMLegendre),
81  std::bind(&GaussPoints::CreateMatrix, this, pl::_1));
82  m_InterpManager.RegisterCreator(
83  PointsKey(0, eGaussRadauPLegendre),
84  std::bind(&GaussPoints::CreateMatrix, this, pl::_1));
85  m_InterpManager.RegisterCreator(
86  PointsKey(0, eGaussLobattoLegendre),
87  std::bind(&GaussPoints::CreateMatrix, this, pl::_1));
88  m_InterpManager.RegisterCreator(
89  PointsKey(0, eGaussGaussChebyshev),
90  std::bind(&GaussPoints::CreateMatrix, this, pl::_1));
91  m_InterpManager.RegisterCreator(
92  PointsKey(0, eGaussRadauMChebyshev),
93  std::bind(&GaussPoints::CreateMatrix, this, pl::_1));
94  m_InterpManager.RegisterCreator(
95  PointsKey(0, eGaussRadauPChebyshev),
96  std::bind(&GaussPoints::CreateMatrix, this, pl::_1));
97  m_InterpManager.RegisterCreator(
98  PointsKey(0, eGaussLobattoChebyshev),
99  std::bind(&GaussPoints::CreateMatrix, this, pl::_1));
100  m_InterpManager.RegisterCreator(
101  PointsKey(0, eGaussRadauMAlpha0Beta1),
102  std::bind(&GaussPoints::CreateMatrix, this, pl::_1));
103  m_InterpManager.RegisterCreator(
104  PointsKey(0, eGaussRadauMAlpha0Beta2),
105  std::bind(&GaussPoints::CreateMatrix, this, pl::_1));
106  m_InterpManager.RegisterCreator(
107  PointsKey(0, eGaussRadauMAlpha1Beta0),
108  std::bind(&GaussPoints::CreateMatrix, this, pl::_1));
109  m_InterpManager.RegisterCreator(
110  PointsKey(0, eGaussRadauMAlpha2Beta0),
111  std::bind(&GaussPoints::CreateMatrix, this, pl::_1));
112  m_InterpManager.RegisterCreator(
113  PointsKey(0, eGaussKronrodLegendre),
114  std::bind(&GaussPoints::CreateMatrix, this, pl::_1));
115  m_InterpManager.RegisterCreator(
116  PointsKey(0, eGaussRadauKronrodMLegendre),
117  std::bind(&GaussPoints::CreateMatrix, this, pl::_1));
118  m_InterpManager.RegisterCreator(
119  PointsKey(0, eGaussRadauKronrodMAlpha1Beta0),
120  std::bind(&GaussPoints::CreateMatrix, this, pl::_1));
121  m_InterpManager.RegisterCreator(
122  PointsKey(0, eGaussLobattoKronrodLegendre),
123  std::bind(&GaussPoints::CreateMatrix, this, pl::_1));
124  m_InterpManager.RegisterCreator(
125  PointsKey(0, eFourierEvenlySpaced),
126  std::bind(&GaussPoints::CreateMatrix, this, pl::_1));
127  m_InterpManager.RegisterCreator(
128  PointsKey(0, ePolyEvenlySpaced),
129  std::bind(&GaussPoints::CreateMatrix, this, pl::_1));
130  m_InterpManager.RegisterCreator(
131  PointsKey(0, eBoundaryLayerPoints),
132  std::bind(&GaussPoints::CreateMatrix, this, pl::_1));
133  m_InterpManager.RegisterCreator(
134  PointsKey(0, eBoundaryLayerPointsRev),
135  std::bind(&GaussPoints::CreateMatrix, this, pl::_1));
136  m_GalerkinProjectionManager.RegisterCreator(
137  PointsKey(0, eGaussGaussLegendre),
138  std::bind(&GaussPoints::CreateGPMatrix, this, pl::_1));
139  m_GalerkinProjectionManager.RegisterCreator(
140  PointsKey(0, eGaussRadauMLegendre),
141  std::bind(&GaussPoints::CreateGPMatrix, this, pl::_1));
142  m_GalerkinProjectionManager.RegisterCreator(
143  PointsKey(0, eGaussRadauPLegendre),
144  std::bind(&GaussPoints::CreateGPMatrix, this, pl::_1));
145  m_GalerkinProjectionManager.RegisterCreator(
146  PointsKey(0, eGaussLobattoLegendre),
147  std::bind(&GaussPoints::CreateGPMatrix, this, pl::_1));
148  m_GalerkinProjectionManager.RegisterCreator(
149  PointsKey(0, eGaussGaussChebyshev),
150  std::bind(&GaussPoints::CreateGPMatrix, this, pl::_1));
151  m_GalerkinProjectionManager.RegisterCreator(
152  PointsKey(0, eGaussRadauMChebyshev),
153  std::bind(&GaussPoints::CreateGPMatrix, this, pl::_1));
154  m_GalerkinProjectionManager.RegisterCreator(
155  PointsKey(0, eGaussRadauPChebyshev),
156  std::bind(&GaussPoints::CreateGPMatrix, this, pl::_1));
157  m_GalerkinProjectionManager.RegisterCreator(
158  PointsKey(0, eGaussLobattoChebyshev),
159  std::bind(&GaussPoints::CreateGPMatrix, this, pl::_1));
160  m_GalerkinProjectionManager.RegisterCreator(
161  PointsKey(0, eGaussRadauMAlpha0Beta1),
162  std::bind(&GaussPoints::CreateGPMatrix, this, pl::_1));
163  m_GalerkinProjectionManager.RegisterCreator(
164  PointsKey(0, eGaussRadauMAlpha0Beta2),
165  std::bind(&GaussPoints::CreateGPMatrix, this, pl::_1));
166  m_GalerkinProjectionManager.RegisterCreator(
167  PointsKey(0, eGaussRadauMAlpha1Beta0),
168  std::bind(&GaussPoints::CreateGPMatrix, this, pl::_1));
169  m_GalerkinProjectionManager.RegisterCreator(
170  PointsKey(0, eGaussRadauMAlpha2Beta0),
171  std::bind(&GaussPoints::CreateGPMatrix, this, pl::_1));
172  m_GalerkinProjectionManager.RegisterCreator(
173  PointsKey(0, eGaussKronrodLegendre),
174  std::bind(&GaussPoints::CreateGPMatrix, this, pl::_1));
175  m_GalerkinProjectionManager.RegisterCreator(
176  PointsKey(0, eGaussRadauKronrodMLegendre),
177  std::bind(&GaussPoints::CreateGPMatrix, this, pl::_1));
178  m_GalerkinProjectionManager.RegisterCreator(
179  PointsKey(0, eGaussRadauKronrodMAlpha1Beta0),
180  std::bind(&GaussPoints::CreateGPMatrix, this, pl::_1));
181  m_GalerkinProjectionManager.RegisterCreator(
182  PointsKey(0, eGaussLobattoKronrodLegendre),
183  std::bind(&GaussPoints::CreateGPMatrix, this, pl::_1));
184  m_GalerkinProjectionManager.RegisterCreator(
185  PointsKey(0, eFourierEvenlySpaced),
186  std::bind(&GaussPoints::CreateGPMatrix, this, pl::_1));
187  m_GalerkinProjectionManager.RegisterCreator(
188  PointsKey(0, ePolyEvenlySpaced),
189  std::bind(&GaussPoints::CreateGPMatrix, this, pl::_1));
190  m_GalerkinProjectionManager.RegisterCreator(
191  PointsKey(0, eBoundaryLayerPoints),
192  std::bind(&GaussPoints::CreateGPMatrix, this, pl::_1));
193  m_GalerkinProjectionManager.RegisterCreator(
194  PointsKey(0, eBoundaryLayerPointsRev),
195  std::bind(&GaussPoints::CreateGPMatrix, this, pl::_1));
196  }
std::shared_ptr< NekMatrix< NekDouble > > CreateGPMatrix(const PointsKey &pkey)
std::shared_ptr< NekMatrix< NekDouble > > CreateMatrix(const PointsKey &pkey)
NekManager< PointsKey, NekMatrix< DataType >, PointsKey::opLess > m_InterpManager
Definition: Points.h:411
NekManager< PointsKey, NekMatrix< DataType >, PointsKey::opLess > m_GalerkinProjectionManager
Definition: Points.h:413
Points< NekDouble > PointsBaseType
@ eGaussRadauMLegendre
1D Gauss-Radau-Legendre quadrature points, pinned at x=-1
Definition: PointsType.h:49
@ eBoundaryLayerPoints
1D power law distribution for boundary layer points
Definition: PointsType.h:79
@ eGaussRadauMAlpha0Beta1
Gauss Radau pinned at x=-1,.
Definition: PointsType.h:60
@ eGaussRadauKronrodMLegendre
1D Gauss-Radau-Kronrod-Legendre quadrature points, pinned at x=-1
Definition: PointsType.h:69
@ eGaussLobattoChebyshev
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:59
@ eFourierEvenlySpaced
1D Evenly-spaced points using Fourier Fit
Definition: PointsType.h:76
@ eGaussRadauPChebyshev
1D Gauss-Radau-Chebyshev quadrature points, pinned at x=1
Definition: PointsType.h:57
@ eGaussRadauMChebyshev
1D Gauss-Radau-Chebyshev quadrature points, pinned at x=-1
Definition: PointsType.h:55
@ eGaussRadauKronrodMAlpha1Beta0
1D Gauss-Radau-Kronrod-Legendre pinned at x=-1,
Definition: PointsType.h:71
@ eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:53
@ eGaussGaussChebyshev
1D Gauss-Gauss-Chebyshev quadrature points
Definition: PointsType.h:54
@ eBoundaryLayerPointsRev
1D power law distribution for boundary layer points
Definition: PointsType.h:81
@ ePolyEvenlySpaced
1D Evenly-spaced points using Lagrange polynomial
Definition: PointsType.h:75
@ eGaussLobattoKronrodLegendre
1D Lobatto Kronrod quadrature points
Definition: PointsType.h:74
@ eGaussGaussLegendre
1D Gauss-Gauss-Legendre quadrature points
Definition: PointsType.h:48
@ eGaussRadauPLegendre
1D Gauss-Radau-Legendre quadrature points, pinned at x=1
Definition: PointsType.h:51

References CreateGPMatrix(), CreateMatrix(), Nektar::LibUtilities::eBoundaryLayerPoints, Nektar::LibUtilities::eBoundaryLayerPointsRev, Nektar::LibUtilities::eFourierEvenlySpaced, Nektar::LibUtilities::eGaussGaussChebyshev, Nektar::LibUtilities::eGaussGaussLegendre, Nektar::LibUtilities::eGaussLobattoChebyshev, Nektar::LibUtilities::eGaussLobattoKronrodLegendre, Nektar::LibUtilities::eGaussLobattoLegendre, Nektar::LibUtilities::eGaussRadauKronrodMAlpha1Beta0, Nektar::LibUtilities::eGaussRadauKronrodMLegendre, Nektar::LibUtilities::eGaussRadauMAlpha0Beta1, Nektar::LibUtilities::eGaussRadauMChebyshev, Nektar::LibUtilities::eGaussRadauMLegendre, Nektar::LibUtilities::eGaussRadauPChebyshev, Nektar::LibUtilities::eGaussRadauPLegendre, Nektar::LibUtilities::ePolyEvenlySpaced, Nektar::LibUtilities::Points< NekDouble >::m_GalerkinProjectionManager, and Nektar::LibUtilities::Points< NekDouble >::m_InterpManager.

◆ GaussPoints() [2/3]

Nektar::LibUtilities::GaussPoints::GaussPoints ( )
private

These should not be called. All creation is done using the constructor requiring the key, declared above.

◆ GaussPoints() [3/3]

Nektar::LibUtilities::GaussPoints::GaussPoints ( const GaussPoints points)
private

Member Function Documentation

◆ CalculateDerivMatrix()

void Nektar::LibUtilities::GaussPoints::CalculateDerivMatrix ( )
privatevirtual

Reimplemented from Nektar::LibUtilities::Points< NekDouble >.

Definition at line 184 of file GaussPoints.cpp.

185 {
186  // Allocate the derivative matrix
188 
189  int numpoints = m_pointsKey.GetNumPoints();
190  int totpoints = m_pointsKey.GetTotNumPoints();
191  Array<OneD, NekDouble> dmtempSharedArray =
192  Array<OneD, NekDouble>(totpoints * totpoints);
193  NekDouble *dmtemp = dmtempSharedArray.data();
194 
195  switch (m_pointsKey.GetPointsType())
196  {
197  case eGaussGaussLegendre:
198  Polylib::Dgj(dmtemp, m_points[0].data(), numpoints, 0.0, 0.0);
199  break;
200 
202  Polylib::Dgrjm(dmtemp, m_points[0].data(), numpoints, 0.0, 0.0);
203  break;
204 
206  Polylib::Dgrjp(dmtemp, m_points[0].data(), numpoints, 0.0, 0.0);
207  break;
208 
210  Polylib::Dglj(dmtemp, m_points[0].data(), numpoints, 0.0, 0.0);
211  break;
212 
214  Polylib::Dgj(dmtemp, m_points[0].data(), numpoints, -0.5, -0.5);
215  break;
216 
218  Polylib::Dgrjm(dmtemp, m_points[0].data(), numpoints, -0.5, -0.5);
219  break;
220 
222  Polylib::Dgrjp(dmtemp, m_points[0].data(), numpoints, -0.5, -0.5);
223  break;
224 
226  Polylib::Dglj(dmtemp, m_points[0].data(), numpoints, -0.5, -0.5);
227  break;
228 
230  Polylib::Dgrjm(dmtemp, m_points[0].data(), numpoints, 0.0, 1.0);
231  break;
232 
233  case eGaussRadauMAlpha0Beta2:
234  Polylib::Dgrjm(dmtemp, m_points[0].data(), numpoints, 0.0, 2.0);
235  break;
236 
237  case eGaussRadauMAlpha1Beta0:
238  Polylib::Dgrjm(dmtemp, m_points[0].data(), numpoints, 1.0, 0.0);
239  break;
240 
241  case eGaussRadauMAlpha2Beta0:
242  Polylib::Dgrjm(dmtemp, m_points[0].data(), numpoints, 2.0, 0.0);
243  break;
244 
245  case eGaussKronrodLegendre:
249  {
250  for (unsigned int i = 0; i < m_pointsKey.GetNumPoints(); ++i)
251  {
252  for (unsigned int j = 0; j < m_pointsKey.GetNumPoints(); ++j)
253  {
254  (*m_derivmatrix[0])(i, j) = LagrangePolyDeriv(
255  m_points[0][i], j, m_pointsKey.GetNumPoints(),
256  m_points[0]);
257  }
258  }
259  return;
260  }
261  break;
262 
263  default:
265  "Unknown Gauss quadrature point distribution requested");
266  }
267 
268  std::copy(dmtemp, dmtemp + totpoints * totpoints,
269  m_derivmatrix[0]->begin());
270 }
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition: ErrorUtil.hpp:209
NekDouble LagrangePolyDeriv(NekDouble x, int pt, int npts, const Array< OneD, const NekDouble > &xpts)
Array< OneD, DataType > m_points[3]
Storage for the point locations, allowing for up to a 3D points storage.
Definition: Points.h:403
MatrixSharedPtrType m_derivmatrix[3]
Derivative matrices.
Definition: Points.h:409
PointsKey m_pointsKey
Points type for this points distributions.
Definition: Points.h:400
unsigned int GetTotNumPoints() const
Definition: Points.h:177
PointsType GetPointsType() const
Definition: Points.h:109
unsigned int GetNumPoints() const
Definition: Points.h:104
def copy(self)
Definition: pycml.py:2663
double NekDouble
void Dgj(double *D, const double *z, const int np, const double alpha, const double beta)
Compute the Derivative Matrix and its transpose associated with the Gauss-Jacobi zeros.
Definition: Polylib.cpp:558
void Dglj(double *D, const double *z, const int np, const double alpha, const double beta)
Compute the Derivative Matrix associated with the Gauss-Lobatto-Jacobi zeros.
Definition: Polylib.cpp:718
void Dgrjm(double *D, const double *z, const int np, const double alpha, const double beta)
Compute the Derivative Matrix and its transpose associated with the Gauss-Radau-Jacobi zeros with a z...
Definition: Polylib.cpp:605
void Dgrjp(double *D, const double *z, const int np, const double alpha, const double beta)
Compute the Derivative Matrix associated with the Gauss-Radau-Jacobi zeros with a zero at z=1.
Definition: Polylib.cpp:661

References Nektar::LibUtilities::Points< NekDouble >::CalculateDerivMatrix(), CellMLToNektar.pycml::copy(), Polylib::Dgj(), Polylib::Dglj(), Polylib::Dgrjm(), Polylib::Dgrjp(), Nektar::ErrorUtil::efatal, Nektar::LibUtilities::eGaussGaussChebyshev, Nektar::LibUtilities::eGaussGaussLegendre, Nektar::LibUtilities::eGaussLobattoChebyshev, Nektar::LibUtilities::eGaussLobattoKronrodLegendre, Nektar::LibUtilities::eGaussLobattoLegendre, Nektar::LibUtilities::eGaussRadauKronrodMAlpha1Beta0, Nektar::LibUtilities::eGaussRadauKronrodMLegendre, Nektar::LibUtilities::eGaussRadauMAlpha0Beta1, Nektar::LibUtilities::eGaussRadauMChebyshev, Nektar::LibUtilities::eGaussRadauMLegendre, Nektar::LibUtilities::eGaussRadauPChebyshev, Nektar::LibUtilities::eGaussRadauPLegendre, Nektar::LibUtilities::PointsKey::GetNumPoints(), Nektar::LibUtilities::PointsKey::GetPointsType(), Nektar::LibUtilities::PointsKey::GetTotNumPoints(), LagrangePolyDeriv(), Nektar::LibUtilities::Points< NekDouble >::m_derivmatrix, Nektar::LibUtilities::Points< NekDouble >::m_points, Nektar::LibUtilities::Points< NekDouble >::m_pointsKey, and NEKERROR.

◆ CalculateGalerkinProjectionMatrix()

std::shared_ptr< NekMatrix< NekDouble > > Nektar::LibUtilities::GaussPoints::CalculateGalerkinProjectionMatrix ( const PointsKey pkey)
private

Definition at line 491 of file GaussPoints.cpp.

493 {
494  int numpointsfrom = pkey.GetNumPoints();
495  int numpointsto = GetNumPoints();
496 
497  Array<OneD, const NekDouble> weightsfrom;
498 
499  weightsfrom = PointsManager()[pkey]->GetW();
500 
501  std::shared_ptr<NekMatrix<NekDouble>> Interp = GetI(pkey);
502 
503  Array<OneD, NekDouble> GalProj(numpointsfrom * numpointsto);
504 
505  // set up inner product matrix and multiply by inverse of
506  // diagaonal mass matrix
507  for (int i = 0; i < numpointsto; ++i)
508  {
509  Vmath::Vmul(numpointsfrom, Interp->GetPtr().get() + i * numpointsfrom,
510  1, &weightsfrom[0], 1, &GalProj[0] + i, numpointsto);
511  Vmath::Smul(numpointsfrom, 1.0 / m_weights[i], &GalProj[0] + i,
512  numpointsto, &GalProj[0] + i, numpointsto);
513  }
514 
515  NekDouble *t = GalProj.data();
516  std::shared_ptr<NekMatrix<NekDouble>> returnval(
517  MemoryManager<NekMatrix<NekDouble>>::AllocateSharedPtr(
518  numpointsto, numpointsfrom, t));
519 
520  return returnval;
521 }
const std::shared_ptr< NekMatrix< NekDouble > > GetI(const PointsKey &pkey)
unsigned int GetNumPoints() const
Definition: Points.h:276
Array< OneD, DataType > m_weights
Quadrature weights for the weights.
Definition: Points.h:405
PointsManagerT & PointsManager(void)
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:209
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:248

References GetI(), Nektar::LibUtilities::PointsKey::GetNumPoints(), Nektar::LibUtilities::Points< NekDouble >::GetNumPoints(), Nektar::LibUtilities::Points< NekDouble >::m_weights, Nektar::LibUtilities::PointsManager(), Vmath::Smul(), and Vmath::Vmul().

Referenced by CreateGPMatrix().

◆ CalculateInterpMatrix()

void Nektar::LibUtilities::GaussPoints::CalculateInterpMatrix ( unsigned int  npts,
const Array< OneD, const NekDouble > &  xpoints,
Array< OneD, NekDouble > &  interp 
)
private

Definition at line 272 of file GaussPoints.cpp.

275 {
276  switch (m_pointsKey.GetPointsType())
277  {
278  case eGaussGaussLegendre:
279  Polylib::Imgj(interp.data(), m_points[0].data(), xpoints.data(),
280  GetNumPoints(), npts, 0.0, 0.0);
281  break;
282 
284  Polylib::Imgrjm(interp.data(), m_points[0].data(), xpoints.data(),
285  GetNumPoints(), npts, 0.0, 0.0);
286  break;
287 
289  Polylib::Imgrjp(interp.data(), m_points[0].data(), xpoints.data(),
290  GetNumPoints(), npts, 0.0, 0.0);
291  break;
292 
294  Polylib::Imglj(interp.data(), m_points[0].data(), xpoints.data(),
295  GetNumPoints(), npts, 0.0, 0.0);
296  break;
297 
299  Polylib::Imgj(interp.data(), m_points[0].data(), xpoints.data(),
300  GetNumPoints(), npts, -0.5, -0.5);
301  break;
302 
304  Polylib::Imgrjm(interp.data(), m_points[0].data(), xpoints.data(),
305  GetNumPoints(), npts, -0.5, -0.5);
306  break;
307 
309  Polylib::Imgrjp(interp.data(), m_points[0].data(), xpoints.data(),
310  GetNumPoints(), npts, -0.5, -0.5);
311  break;
312 
314  Polylib::Imglj(interp.data(), m_points[0].data(), xpoints.data(),
315  GetNumPoints(), npts, -0.5, -0.5);
316  break;
317 
319  Polylib::Imgrjm(interp.data(), m_points[0].data(), xpoints.data(),
320  GetNumPoints(), npts, 0.0, 1.0);
321  break;
322 
323  case eGaussRadauMAlpha0Beta2:
324  Polylib::Imgrjm(interp.data(), m_points[0].data(), xpoints.data(),
325  GetNumPoints(), npts, 0.0, 2.0);
326  break;
327 
328  case eGaussRadauMAlpha1Beta0:
329  Polylib::Imgrjm(interp.data(), m_points[0].data(), xpoints.data(),
330  GetNumPoints(), npts, 1.0, 0.0);
331  break;
332 
333  case eGaussRadauMAlpha2Beta0:
334  Polylib::Imgrjm(interp.data(), m_points[0].data(), xpoints.data(),
335  GetNumPoints(), npts, 2.0, 0.0);
336  break;
337 
338  case eGaussKronrodLegendre:
342  {
343  for (unsigned int i = 0; i < npts; ++i)
344  {
345  for (unsigned int j = 0; j < m_pointsKey.GetNumPoints(); ++j)
346  {
347  interp[i + j * npts] = LagrangePoly(
348  xpoints[i], j, m_pointsKey.GetNumPoints(), m_points[0]);
349  }
350  }
351  }
352  break;
353 
354  default:
356  "Unknown Gauss quadrature point distribution requested");
357  }
358 }
NekDouble LagrangePoly(NekDouble x, int pt, int npts, const Array< OneD, const NekDouble > &xpts)
void Imgrjp(double *im, const double *zgrj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
Interpolation Operator from Gauss-Radau-Jacobi points (including z=1) to an arbitrary distrubtion at ...
Definition: Polylib.cpp:985
void Imgj(double *im, const double *zgj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
Interpolation Operator from Gauss-Jacobi points to an arbitrary distribution at points zm.
Definition: Polylib.cpp:923
void Imgrjm(double *im, const double *zgrj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
Interpolation Operator from Gauss-Radau-Jacobi points (including z=-1) to an arbitrary distrubtion at...
Definition: Polylib.cpp:954
void Imglj(double *im, const double *zglj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
Interpolation Operator from Gauss-Lobatto-Jacobi points to an arbitrary distrubtion at points zm.
Definition: Polylib.cpp:1016

References Nektar::Array< OneD, DataType >::data(), Nektar::ErrorUtil::efatal, Nektar::LibUtilities::eGaussGaussChebyshev, Nektar::LibUtilities::eGaussGaussLegendre, Nektar::LibUtilities::eGaussLobattoChebyshev, Nektar::LibUtilities::eGaussLobattoKronrodLegendre, Nektar::LibUtilities::eGaussLobattoLegendre, Nektar::LibUtilities::eGaussRadauKronrodMAlpha1Beta0, Nektar::LibUtilities::eGaussRadauKronrodMLegendre, Nektar::LibUtilities::eGaussRadauMAlpha0Beta1, Nektar::LibUtilities::eGaussRadauMChebyshev, Nektar::LibUtilities::eGaussRadauMLegendre, Nektar::LibUtilities::eGaussRadauPChebyshev, Nektar::LibUtilities::eGaussRadauPLegendre, Nektar::LibUtilities::PointsKey::GetNumPoints(), Nektar::LibUtilities::Points< NekDouble >::GetNumPoints(), Nektar::LibUtilities::PointsKey::GetPointsType(), Polylib::Imgj(), Polylib::Imglj(), Polylib::Imgrjm(), Polylib::Imgrjp(), LagrangePoly(), Nektar::LibUtilities::Points< NekDouble >::m_points, Nektar::LibUtilities::Points< NekDouble >::m_pointsKey, and NEKERROR.

Referenced by GetI().

◆ CalculatePoints()

void Nektar::LibUtilities::GaussPoints::CalculatePoints ( )
privatevirtual

Reimplemented from Nektar::LibUtilities::Points< NekDouble >.

Definition at line 83 of file GaussPoints.cpp.

84 {
85  // Allocate the storage for points and weights
88 
89  int numpoints = m_pointsKey.GetNumPoints();
90 
91  switch (m_pointsKey.GetPointsType())
92  {
94  Polylib::zwgj(m_points[0].data(), m_weights.data(), numpoints, 0.0,
95  0.0);
96  break;
97 
99  Polylib::zwgrjm(m_points[0].data(), m_weights.data(), numpoints,
100  0.0, 0.0);
101  break;
102 
104  Polylib::zwgrjp(m_points[0].data(), m_weights.data(), numpoints,
105  0.0, 0.0);
106  break;
107 
109  Polylib::zwglj(m_points[0].data(), m_weights.data(), numpoints, 0.0,
110  0.0);
111  break;
112 
114  Polylib::zwgj(m_points[0].data(), m_weights.data(), numpoints, -0.5,
115  -0.5);
116  break;
117 
119  Polylib::zwgrjm(m_points[0].data(), m_weights.data(), numpoints,
120  -0.5, -0.5);
121  break;
122 
124  Polylib::zwgrjp(m_points[0].data(), m_weights.data(), numpoints,
125  -0.5, -0.5);
126  break;
127 
129  Polylib::zwglj(m_points[0].data(), m_weights.data(), numpoints,
130  -0.5, -0.5);
131  break;
132 
134  Polylib::zwgrjm(m_points[0].data(), m_weights.data(), numpoints,
135  0.0, 1.0);
136  break;
137 
138  case eGaussRadauMAlpha0Beta2:
139  Polylib::zwgrjm(m_points[0].data(), m_weights.data(), numpoints,
140  0.0, 2.0);
141  break;
142 
143  case eGaussRadauMAlpha1Beta0:
144  Polylib::zwgrjm(m_points[0].data(), m_weights.data(), numpoints,
145  1.0, 0.0);
146  break;
147 
148  case eGaussRadauMAlpha2Beta0:
149  Polylib::zwgrjm(m_points[0].data(), m_weights.data(), numpoints,
150  2.0, 0.0);
151  break;
152 
153  case eGaussKronrodLegendre:
154  Polylib::zwgk(m_points[0].data(), m_weights.data(), numpoints, 0.0,
155  0.0);
156  break;
157 
159  Polylib::zwrk(m_points[0].data(), m_weights.data(), numpoints, 0.0,
160  0.0);
161  break;
162 
164  Polylib::zwrk(m_points[0].data(), m_weights.data(), numpoints, 1.0,
165  0.0);
166  break;
167 
169  Polylib::zwlk(m_points[0].data(), m_weights.data(), numpoints, 0.0,
170  0.0);
171  break;
172 
173  default:
175  "Unknown Gauss quadrature point distribution requested");
176  }
177 }
void zwgrjm(double *z, double *w, const int np, const double alpha, const double beta)
Gauss-Radau-Jacobi zeros and weights with end point at z=-1.
Definition: Polylib.cpp:127
void zwgrjp(double *z, double *w, const int np, const double alpha, const double beta)
Gauss-Radau-Jacobi zeros and weights with end point at z=1.
Definition: Polylib.cpp:168
void zwglj(double *z, double *w, const int np, const double alpha, const double beta)
Gauss-Lobatto-Jacobi zeros and weights with end point at z=-1,1.
Definition: Polylib.cpp:207
void zwgj(double *z, double *w, const int np, const double alpha, const double beta)
Gauss-Jacobi zeros and weights.
Definition: Polylib.cpp:99
void zwgk(double *z, double *w, const int npt, const double alpha, const double beta)
Gauss-Kronrod-Jacobi zeros and weights.
Definition: Polylib.cpp:256
void zwrk(double *z, double *w, const int npt, const double alpha, const double beta)
Gauss-Radau-Kronrod-Jacobi zeros and weights.
Definition: Polylib.cpp:328
void zwlk(double *z, double *w, const int npt, const double alpha, const double beta)
Gauss-Lobatto-Kronrod-Jacobi zeros and weights.
Definition: Polylib.cpp:434

References Nektar::LibUtilities::Points< NekDouble >::CalculatePoints(), Nektar::LibUtilities::Points< NekDouble >::CalculateWeights(), Nektar::Array< OneD, DataType >::data(), Nektar::ErrorUtil::efatal, Nektar::LibUtilities::eGaussGaussChebyshev, Nektar::LibUtilities::eGaussGaussLegendre, Nektar::LibUtilities::eGaussLobattoChebyshev, Nektar::LibUtilities::eGaussLobattoKronrodLegendre, Nektar::LibUtilities::eGaussLobattoLegendre, Nektar::LibUtilities::eGaussRadauKronrodMAlpha1Beta0, Nektar::LibUtilities::eGaussRadauKronrodMLegendre, Nektar::LibUtilities::eGaussRadauMAlpha0Beta1, Nektar::LibUtilities::eGaussRadauMChebyshev, Nektar::LibUtilities::eGaussRadauMLegendre, Nektar::LibUtilities::eGaussRadauPChebyshev, Nektar::LibUtilities::eGaussRadauPLegendre, Nektar::LibUtilities::PointsKey::GetNumPoints(), Nektar::LibUtilities::PointsKey::GetPointsType(), Nektar::LibUtilities::Points< NekDouble >::m_points, Nektar::LibUtilities::Points< NekDouble >::m_pointsKey, Nektar::LibUtilities::Points< NekDouble >::m_weights, NEKERROR, Polylib::zwgj(), Polylib::zwgk(), Polylib::zwglj(), Polylib::zwgrjm(), Polylib::zwgrjp(), Polylib::zwlk(), and Polylib::zwrk().

◆ CalculateWeights()

void Nektar::LibUtilities::GaussPoints::CalculateWeights ( )
privatevirtual

Reimplemented from Nektar::LibUtilities::Points< NekDouble >.

Definition at line 179 of file GaussPoints.cpp.

180 {
181  // For Gauss Quadrature, this is done as part of the points computation
182 }

◆ Create()

std::shared_ptr< Points< NekDouble > > Nektar::LibUtilities::GaussPoints::Create ( const PointsKey pkey)
static

Definition at line 360 of file GaussPoints.cpp.

361 {
362  std::shared_ptr<Points<NekDouble>> returnval(
364 
365  returnval->Initialize();
366 
367  return returnval;
368 }
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.

◆ CreateGPMatrix()

std::shared_ptr< NekMatrix< NekDouble > > Nektar::LibUtilities::GaussPoints::CreateGPMatrix ( const PointsKey pkey)

Definition at line 481 of file GaussPoints.cpp.

483 {
484  std::shared_ptr<NekMatrix<NekDouble>> returnval =
486 
487  // Delegate to function below
488  return returnval;
489 }
std::shared_ptr< NekMatrix< NekDouble > > CalculateGalerkinProjectionMatrix(const PointsKey &pkey)

References CalculateGalerkinProjectionMatrix().

Referenced by GaussPoints().

◆ CreateMatrix()

std::shared_ptr< NekMatrix< NekDouble > > Nektar::LibUtilities::GaussPoints::CreateMatrix ( const PointsKey pkey)

Definition at line 370 of file GaussPoints.cpp.

372 {
373  int numpoints = pkey.GetNumPoints();
374  Array<OneD, const NekDouble> xpoints;
375 
376  PointsManager()[pkey]->GetPoints(xpoints);
377 
378  // Delegate to function below
379  return GetI(numpoints, xpoints);
380 }

References GetI(), Nektar::LibUtilities::PointsKey::GetNumPoints(), and Nektar::LibUtilities::PointsManager().

Referenced by GaussPoints().

◆ GetGalerkinProjection()

const std::shared_ptr< NekMatrix< NekDouble > > Nektar::LibUtilities::GaussPoints::GetGalerkinProjection ( const PointsKey pkey)
virtual

◆ GetI() [1/3]

const std::shared_ptr< NekMatrix< NekDouble > > Nektar::LibUtilities::GaussPoints::GetI ( const Array< OneD, const NekDouble > &  x)
virtual

Reimplemented from Nektar::LibUtilities::Points< NekDouble >.

Definition at line 391 of file GaussPoints.cpp.

393 {
394  int numpoints = 1;
395 
396  // Delegate to function below
397  return GetI(numpoints, x);
398 }

References GetI().

◆ GetI() [2/3]

const std::shared_ptr< NekMatrix< NekDouble > > Nektar::LibUtilities::GaussPoints::GetI ( const PointsKey pkey)
virtual

Reimplemented from Nektar::LibUtilities::Points< NekDouble >.

Definition at line 382 of file GaussPoints.cpp.

384 {
385  ASSERTL0(pkey.GetPointsDim() == 1,
386  "Gauss Points can only interp to other 1d point distributions");
387 
388  return m_InterpManager[pkey];
389 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215

References ASSERTL0, Nektar::LibUtilities::PointsKey::GetPointsDim(), and Nektar::LibUtilities::Points< NekDouble >::m_InterpManager.

Referenced by CalculateGalerkinProjectionMatrix(), CreateMatrix(), and GetI().

◆ GetI() [3/3]

const std::shared_ptr< NekMatrix< NekDouble > > Nektar::LibUtilities::GaussPoints::GetI ( unsigned int  numpoints,
const Array< OneD, const NekDouble > &  x 
)
virtual

Reimplemented from Nektar::LibUtilities::Points< NekDouble >.

Definition at line 400 of file GaussPoints.cpp.

402 {
403  Array<OneD, NekDouble> interp(GetNumPoints() * numpoints);
404 
405  CalculateInterpMatrix(numpoints, x, interp);
406 
407  NekDouble *t = interp.data();
408  unsigned int np = GetNumPoints();
409  std::shared_ptr<NekMatrix<NekDouble>> returnval(
410  MemoryManager<NekMatrix<NekDouble>>::AllocateSharedPtr(numpoints, np,
411  t));
412 
413  return returnval;
414 }
void CalculateInterpMatrix(unsigned int npts, const Array< OneD, const NekDouble > &xpoints, Array< OneD, NekDouble > &interp)

References CalculateInterpMatrix(), and Nektar::LibUtilities::Points< NekDouble >::GetNumPoints().

◆ LagrangeInterpolant()

NekDouble Nektar::LibUtilities::GaussPoints::LagrangeInterpolant ( NekDouble  x,
int  npts,
const Array< OneD, const NekDouble > &  xpts,
const Array< OneD, const NekDouble > &  funcvals 
)
private

functions used by the Kronrod points

Definition at line 416 of file GaussPoints.cpp.

419 {
420  NekDouble sum = 0.0;
421 
422  for (int i = 0; i < npts; ++i)
423  {
424  sum += funcvals[i] * LagrangePoly(x, i, npts, xpts);
425  }
426  return sum;
427 }

References LagrangePoly().

◆ LagrangePoly()

NekDouble Nektar::LibUtilities::GaussPoints::LagrangePoly ( NekDouble  x,
int  pt,
int  npts,
const Array< OneD, const NekDouble > &  xpts 
)
private

Definition at line 429 of file GaussPoints.cpp.

431 {
432  NekDouble h = 1.0;
433 
434  for (int i = 0; i < pt; ++i)
435  {
436  h = h * (x - xpts[i]) / (xpts[pt] - xpts[i]);
437  }
438 
439  for (int i = pt + 1; i < npts; ++i)
440  {
441  h = h * (x - xpts[i]) / (xpts[pt] - xpts[i]);
442  }
443 
444  return h;
445 }

Referenced by CalculateInterpMatrix(), and LagrangeInterpolant().

◆ LagrangePolyDeriv()

NekDouble Nektar::LibUtilities::GaussPoints::LagrangePolyDeriv ( NekDouble  x,
int  pt,
int  npts,
const Array< OneD, const NekDouble > &  xpts 
)
private

Definition at line 447 of file GaussPoints.cpp.

449 {
450  NekDouble h;
451  NekDouble y = 0.0;
452 
453  for (int j = 0; j < npts; ++j)
454  {
455  if (j != pt)
456  {
457  h = 1.0;
458  for (int i = 0; i < npts; ++i)
459  {
460  if (i != pt)
461  {
462  if (i != j)
463  {
464  h = h * (x - xpts[i]);
465  }
466  h = h / (xpts[pt] - xpts[i]);
467  }
468  }
469  y = y + h;
470  }
471  }
472  return y;
473 }

Referenced by CalculateDerivMatrix().

Member Data Documentation

◆ initPointsManager

bool Nektar::LibUtilities::GaussPoints::initPointsManager
staticprivate

Definition at line 199 of file GaussPoints.h.