Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | Static Public Member Functions | Protected Attributes | Private Member Functions | List of all members
Nektar::LibUtilities::Basis Class Reference

Represents a basis of a given type. More...

#include <Basis.h>

Collaboration diagram for Nektar::LibUtilities::Basis:
Collaboration graph
[legend]

Public Member Functions

virtual ~Basis ()
 Destructor. More...
 
int GetNumModes () const
 Return order of basis from the basis specification. More...
 
int GetTotNumModes () const
 Return total number of modes from the basis specification. More...
 
int GetNumPoints () const
 Return the number of points from the basis specification. More...
 
int GetTotNumPoints () const
 Return total number of points from the basis specification. More...
 
BasisType GetBasisType () const
 Return the type of expansion basis. More...
 
PointsKey GetPointsKey () const
 Return the points specification for the basis. More...
 
PointsType GetPointsType () const
 Return the type of quadrature. More...
 
const Array< OneD, const
NekDouble > & 
GetZ () const
 
const Array< OneD, const
NekDouble > & 
GetW () const
 
void GetZW (Array< OneD, const NekDouble > &z, Array< OneD, const NekDouble > &w) const
 
const boost::shared_ptr
< NekMatrix< NekDouble > > & 
GetD (Direction dir=xDir) const
 
const boost::shared_ptr
< NekMatrix< NekDouble > > 
GetI (const Array< OneD, const NekDouble > &x)
 
const boost::shared_ptr
< NekMatrix< NekDouble > > 
GetI (const BasisKey &bkey)
 
bool ExactIprodInt () const
 Determine if basis has exact integration for inner product. More...
 
bool Collocation () const
 Determine if basis has collocation properties. More...
 
const Array< OneD, const
NekDouble > & 
GetBdata () const
 Return basis definition array m_bdata. More...
 
const Array< OneD, const
NekDouble > & 
GetDbdata () const
 Return basis definition array m_dbdata. More...
 
const BasisKey GetBasisKey () const
 Returns the specification for the Basis. More...
 
virtual void Initialize ()
 

Static Public Member Functions

static boost::shared_ptr< BasisCreate (const BasisKey &bkey)
 Returns a new instance of a Basis with given BasisKey. More...
 

Protected Attributes

BasisKey m_basisKey
 Basis specification. More...
 
PointsSharedPtr m_points
 Set of points. More...
 
Array< OneD, NekDoublem_bdata
 Basis definition. More...
 
Array< OneD, NekDoublem_dbdata
 Derivative Basis definition. More...
 
NekManager< BasisKey,
NekMatrix< NekDouble >
, BasisKey::opLess
m_InterpManager
 

Private Member Functions

 Basis (const BasisKey &bkey)
 Private constructor with BasisKey. More...
 
 Basis ()
 Private default constructor. More...
 
boost::shared_ptr< NekMatrix
< NekDouble > > 
CalculateInterpMatrix (const BasisKey &tbasis0)
 Calculate the interpolation Matrix for coefficient from one base (m_basisKey) to another (tbasis0) More...
 
void GenBasis ()
 Generate appropriate basis and their derivatives. More...
 

Detailed Description

Represents a basis of a given type.

Definition at line 206 of file Basis.h.

Constructor & Destructor Documentation

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

Destructor.

Definition at line 213 of file Basis.h.

214  {
215  };
Nektar::LibUtilities::Basis::Basis ( const BasisKey bkey)
private

Private constructor with BasisKey.

Definition at line 93 of file Basis.cpp.

References CalculateInterpMatrix(), and m_InterpManager.

93  :
94  m_basisKey(bkey),
95  m_points(PointsManager()[bkey.GetPointsKey()]),
96  m_bdata(bkey.GetTotNumModes()*bkey.GetTotNumPoints()),
97  m_dbdata(bkey.GetTotNumModes()*bkey.GetTotNumPoints())
98  {
99  m_InterpManager.RegisterGlobalCreator(boost::bind(&Basis::CalculateInterpMatrix,this,_1));
100  }
BasisKey m_basisKey
Basis specification.
Definition: Basis.h:328
Array< OneD, NekDouble > m_bdata
Basis definition.
Definition: Basis.h:330
Array< OneD, NekDouble > m_dbdata
Derivative Basis definition.
Definition: Basis.h:331
PointsSharedPtr m_points
Set of points.
Definition: Basis.h:329
PointsManagerT & PointsManager(void)
boost::shared_ptr< NekMatrix< NekDouble > > CalculateInterpMatrix(const BasisKey &tbasis0)
Calculate the interpolation Matrix for coefficient from one base (m_basisKey) to another (tbasis0) ...
Definition: Basis.cpp:122
NekManager< BasisKey, NekMatrix< NekDouble >, BasisKey::opLess > m_InterpManager
Definition: Basis.h:333
Nektar::LibUtilities::Basis::Basis ( )
inlineprivate

Private default constructor.

Definition at line 340 of file Basis.h.

References ErrorUtil::efatal, and NEKERROR.

Referenced by Create().

341  {
343  "Default Constructor for Basis should not be called");
344  }
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
Definition: ErrorUtil.hpp:158
BasisKey m_basisKey
Basis specification.
Definition: Basis.h:328
static const BasisKey NullBasisKey(eNoBasisType, 0, NullPointsKey)
Defines a null basis with no type or points.

Member Function Documentation

boost::shared_ptr< NekMatrix< NekDouble > > Nektar::LibUtilities::Basis::CalculateInterpMatrix ( const BasisKey tbasis0)
private

Calculate the interpolation Matrix for coefficient from one base (m_basisKey) to another (tbasis0)

Definition at line 122 of file Basis.cpp.

References Nektar::LibUtilities::BasisManager(), Nektar::LibUtilities::eGaussLobattoLegendre, Nektar::LibUtilities::BasisKey::GetBasisType(), Nektar::LibUtilities::BasisKey::GetNumModes(), and m_basisKey.

Referenced by Basis().

123  {
124  int dim = m_basisKey.GetNumModes();
125  const PointsKey pkey(dim,LibUtilities::eGaussLobattoLegendre);
126  BasisKey fbkey(m_basisKey.GetBasisType(),dim,pkey);
127  BasisKey tbkey(tbasis0.GetBasisType(),dim,pkey);
128 
129  // "Constructur" of the basis
130  BasisSharedPtr fbasis = BasisManager()[fbkey];
131  BasisSharedPtr tbasis = BasisManager()[tbkey];
132 
133  // Get B Matrices
134  Array<OneD, NekDouble> fB_data = fbasis->GetBdata();
135  Array<OneD, NekDouble> tB_data = tbasis->GetBdata();
136 
137  // Convert to a NekMatrix
138  NekMatrix<NekDouble> fB(dim,dim,fB_data);
139  NekMatrix<NekDouble> tB(dim,dim,tB_data);
140 
141  // Invert the "to" matrix: tu = tB^(-1)*fB fu = ftB fu
142  tB.Invert();
143 
144  // Compute transformation matrix
145  Array<OneD, NekDouble> zero1D(dim*dim,0.0);
146  boost::shared_ptr< NekMatrix<NekDouble> > ftB(MemoryManager<NekMatrix<NekDouble> >::AllocateSharedPtr(dim,dim,zero1D));
147  (*ftB) = tB*fB;
148 
149  return ftB;
150  }
BasisKey m_basisKey
Basis specification.
Definition: Basis.h:328
BasisType GetBasisType() const
Return type of expansion basis.
Definition: Basis.h:139
BasisManagerT & BasisManager(void)
int GetNumModes() const
Returns the order of the basis.
Definition: Basis.h:84
boost::shared_ptr< Basis > BasisSharedPtr
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:50
bool Nektar::LibUtilities::Basis::Collocation ( ) const
inline

Determine if basis has collocation properties.

Definition at line 302 of file Basis.h.

References Nektar::LibUtilities::BasisKey::Collocation(), and m_basisKey.

303  {
304  return m_basisKey.Collocation();
305  }
BasisKey m_basisKey
Basis specification.
Definition: Basis.h:328
bool Collocation() const
Determine if basis has collocation properties.
Definition: Basis.cpp:692
boost::shared_ptr< Basis > Nektar::LibUtilities::Basis::Create ( const BasisKey bkey)
static

Returns a new instance of a Basis with given BasisKey.

Definition at line 102 of file Basis.cpp.

References Basis().

103  {
104  boost::shared_ptr<Basis> returnval(new Basis(bkey));
105  returnval->Initialize();
106 
107  return returnval;
108  }
Basis()
Private default constructor.
Definition: Basis.h:340
bool Nektar::LibUtilities::Basis::ExactIprodInt ( void  ) const
inline

Determine if basis has exact integration for inner product.

Definition at line 296 of file Basis.h.

References Nektar::LibUtilities::BasisKey::ExactIprodInt(), and m_basisKey.

297  {
298  return m_basisKey.ExactIprodInt();
299  }
BasisKey m_basisKey
Basis specification.
Definition: Basis.h:328
bool ExactIprodInt() const
Determine if basis has exact integration for inner product.
Definition: Basis.cpp:665
void Nektar::LibUtilities::Basis::GenBasis ( )
private

Generate appropriate basis and their derivatives.

The following expansions are generated depending on the enum type defined in m_basisKey.m_basistype:

NOTE: This definition does not follow the order in the Karniadakis & Sherwin book since this leads to a more compact hierarchical pattern for implementation purposes. The order of these modes dictates the ordering of the expansion coefficients.

In the following m_numModes = P

eModified_A:

m_bdata[i + j*m_numpoints] = $ \phi^a_i(z_j) = \left \{ \begin{array}{ll} \left ( \frac{1-z_j}{2}\right ) & i = 0 \\ \\ \left ( \frac{1+z_j}{2}\right ) & i = 1 \\ \\ \left ( \frac{1-z_j}{2}\right )\left ( \frac{1+z_j}{2}\right ) P^{1,1}_{i-2}(z_j) & 2\leq i < P\\ \end{array} \right . $

eModified_B:

m_bdata[n(i,j) + k*m_numpoints] = $ \phi^b_{ij}(z_k) = \left \{ \begin{array}{lll} \phi^a_j(z_k) & i = 0, & 0\leq j < P \\ \\ \left ( \frac{1-z_k}{2}\right )^{i} & 1 \leq i < P,& j = 0 \\ \\ \left ( \frac{1-z_k}{2}\right )^{i} \left ( \frac{1+z_k}{2}\right ) P^{2i-1,1}_{j-1}(z_k) & 1 \leq i < P,\ & 1\leq j < P-i\ \\ \end{array} \right . , $

where $ n(i,j) $ is a consecutive ordering of the triangular indices $ 0 \leq i, i+j < P $ where j runs fastest.

eModified_C:

m_bdata[n(i,j,k) + l*m_numpoints] = $ \phi^c_{ij,k}(z_l) = \phi^b_{i+j,k}(z_l) = \left \{ \begin{array}{llll} \phi^b_{j,k}(z_l) & i = 0, & 0\leq j < P & 0\leq k < P-j\\ \\ \left ( \frac{1-z_l}{2}\right )^{i+j} & 1\leq i < P,\ & 0\leq j < P-i,\ & k = 0 \\ \\ \left ( \frac{1-z_l}{2}\right )^{i+j} \left ( \frac{1+z_l}{2}\right ) P^{2i+2j-1,1}_{k-1}(z_k) & 1\leq i < P,& 0\leq j < P-i& 1\leq k < P-i-j \\ \\ \end{array} \right . , $

where $ n(i,j,k) $ is a consecutive ordering of the triangular indices $ 0 \leq i, i+j, i+j+k < P $ where k runs fastest, then j and finally i.

Orthogonal basis A

$\tilde \psi_p^a (\eta_1) = L_p(\eta_1) = P_p^{0,0}(\eta_1)$

Orthogonal basis B

$\tilde \psi_{pq}^b(\eta_2) = \left ( {1 - \eta_2} \over 2 \right)^p P_q^{2p+1,0}(\eta_2)$ \

Orthogonal basis C

$\tilde \psi_{pqr}^c = \left ( {1 - \eta_3} \over 2 \right)^{p+q} P_r^{2p+2q+2, 0}(\eta_3)$ \

Definition at line 215 of file Basis.cpp.

References ASSERTL0, Nektar::LibUtilities::BasisManager(), Nektar::LibUtilities::eChebyshev, Nektar::LibUtilities::eFourier, Nektar::LibUtilities::eFourierHalfModeIm, Nektar::LibUtilities::eFourierHalfModeRe, Nektar::LibUtilities::eFourierSingleMode, Nektar::LibUtilities::eGauss_Lagrange, Nektar::LibUtilities::eGaussGaussLegendre, Nektar::LibUtilities::eGaussLobattoLegendre, Nektar::LibUtilities::eGLL_Lagrange, Nektar::LibUtilities::eLegendre, Nektar::LibUtilities::eModified_A, Nektar::LibUtilities::eModified_B, Nektar::LibUtilities::eModified_C, Nektar::LibUtilities::eMonomial, Nektar::LibUtilities::eOrtho_A, Nektar::LibUtilities::eOrtho_B, Nektar::LibUtilities::eOrtho_C, GetBasisType(), Nektar::LibUtilities::BasisKey::GetNumModes(), GetNumModes(), GetNumPoints(), Nektar::LibUtilities::BasisKey::GetPointsKey(), Polylib::hgj(), Polylib::hglj(), Polylib::jacobfd(), m_basisKey, m_bdata, m_dbdata, m_points, Nektar::LibUtilities::PointsManager(), and Vmath::Vcopy().

Referenced by Initialize().

216  {
217  int i,p,q;
218  NekDouble scal;
219  Array<OneD, NekDouble> modeSharedArray;
220  NekDouble *mode;
221  Array<OneD, const NekDouble> z;
222  Array<OneD, const NekDouble> w;
223  const NekDouble *D;
224 
225  m_points->GetZW(z,w);
226 
227  D = &(m_points->GetD()->GetPtr())[0];
228  int numModes = GetNumModes();
229  int numPoints = GetNumPoints();
230 
231  switch(GetBasisType())
232  {
233 
234  /** \brief Orthogonal basis A
235 
236  \f$\tilde \psi_p^a (\eta_1) = L_p(\eta_1) = P_p^{0,0}(\eta_1)\f$
237 
238  */
239  case eOrtho_A:
240  case eLegendre:
241  mode = m_bdata.data();
242 
243  for (p=0; p<numModes; ++p, mode += numPoints)
244  {
245  Polylib::jacobfd(numPoints, z.data(), mode, NULL, p, 0.0, 0.0);
246  // normalise
247  scal = sqrt(0.5*(2.0*p+1.0));
248  for(i = 0; i < numPoints; ++i)
249  {
250  mode[i] *= scal;
251  }
252  }
253  // define derivative basis
254  Blas::Dgemm('n','n',numPoints,numModes,numPoints,1.0,D,numPoints,
255  m_bdata.data(),numPoints,0.0,m_dbdata.data(),numPoints);
256  break;
257 
258  /** \brief Orthogonal basis B
259 
260  \f$\tilde \psi_{pq}^b(\eta_2) = \left ( {1 - \eta_2} \over 2 \right)^p P_q^{2p+1,0}(\eta_2)\f$ \\
261 
262  */
263 
264  // This is tilde psi_pq in Spen's book, page 105
265  // The 3-dimensional array is laid out in memory such that
266  // 1) Eta_y values are the changing the fastest, then q and p.
267  // 2) q index increases by the stride of numPoints.
268  case eOrtho_B:
269  {
270  NekDouble *mode = m_bdata.data();
271 
272  for( int p = 0; p < numModes; ++p )
273  {
274  for( int q = 0; q < numModes - p; ++q, mode += numPoints )
275  {
276  Polylib::jacobfd(numPoints, z.data(), mode, NULL, q, 2*p + 1.0, 0.0);
277  for( int j = 0; j < numPoints; ++j )
278  {
279  mode[j] *= sqrt(p+q+1.0)*pow(0.5*(1.0 - z[j]), p);
280  }
281  }
282  }
283 
284  // define derivative basis
285  Blas::Dgemm('n','n',numPoints,numModes*(numModes+1)/2,numPoints,1.0,D,numPoints,
286  m_bdata.data(),numPoints,0.0,m_dbdata.data(),numPoints);
287  }
288  break;
289 
290  /** \brief Orthogonal basis C
291 
292  \f$\tilde \psi_{pqr}^c = \left ( {1 - \eta_3} \over 2 \right)^{p+q} P_r^{2p+2q+2, 0}(\eta_3)\f$ \\
293 
294  */
295 
296  // This is tilde psi_pqr in Spen's book, page 105
297  // The 4-dimensional array is laid out in memory such that
298  // 1) Eta_z values are the changing the fastest, then r, q, and finally p.
299  // 2) r index increases by the stride of numPoints.
300  case eOrtho_C:
301  {
302  int P = numModes - 1, Q = numModes - 1, R = numModes - 1;
303  NekDouble *mode = m_bdata.data();
304 
305  for( int p = 0; p <= P; ++p )
306  {
307  for( int q = 0; q <= Q - p; ++q )
308  {
309  for( int r = 0; r <= R - p - q; ++r, mode += numPoints )
310  {
311  Polylib::jacobfd(numPoints, z.data(), mode, NULL, r, 2*p + 2*q + 2.0, 0.0);
312  for( int k = 0; k < numPoints; ++k )
313  {
314  // Note factor of 0.5 is part of normalisation
315  mode[k] *= pow(0.5*(1.0 - z[k]), p+q);
316 
317  // finish normalisation
318  mode[k] *= sqrt(r+p+q+1.5);
319  }
320  }
321  }
322  }
323 
324  // Define derivative basis
325  Blas::Dgemm('n','n',numPoints,numModes*(numModes+1)*
326  (numModes+2)/6,numPoints,1.0, D, numPoints,
327  m_bdata.data(),numPoints,0.0,m_dbdata.data(),numPoints);
328  }
329  break;
330 
331  case eModified_A:
332 
333  // Note the following packing deviates from the
334  // definition in the Book by Karniadakis in that we
335  // put the vertex degrees of freedom at the lower
336  // index range to follow a more hierarchic structure.
337 
338  for(i = 0; i < numPoints; ++i)
339  {
340  m_bdata[i] = 0.5*(1-z[i]);
341  m_bdata[numPoints + i] = 0.5*(1+z[i]);
342  }
343 
344  mode = m_bdata.data() + 2*numPoints;
345 
346  for(p = 2; p < numModes; ++p, mode += numPoints)
347  {
348  Polylib::jacobfd(numPoints, z.data(), mode, NULL, p-2,1.0,1.0);
349 
350  for(i = 0; i < numPoints; ++i)
351  {
352  mode[i] *= m_bdata[i]*m_bdata[numPoints+i];
353  }
354  }
355 
356  // define derivative basis
357  Blas::Dgemm('n','n',numPoints,numModes,numPoints,1.0,D,
358  numPoints,m_bdata.data(),numPoints,0.0,m_dbdata.data(),
359  numPoints);
360  break;
361 
362  case eModified_B:
363  {
364 
365  // Note the following packing deviates from the
366  // definition in the Book by Karniadakis in two
367  // ways. 1) We put the vertex degrees of freedom
368  // at the lower index range to follow a more
369  // hierarchic structure. 2) We do not duplicate
370  // the singular vertex definition so that only a
371  // triangular number (i.e. (modes)*(modes+1)/2) of
372  // modes are required consistent with the
373  // orthogonal basis.
374 
375  // In the current structure the q index runs
376  // faster than the p index so that the matrix has
377  // a more compact structure
378 
379  const NekDouble *one_m_z_pow, *one_p_z;
380 
381  // bdata should be of size order*(order+1)/2*zorder
382 
383  // first fow
384  for(i = 0; i < numPoints; ++i)
385  {
386  m_bdata[0*numPoints + i] = 0.5*(1-z[i]);
387  m_bdata[1*numPoints + i] = 0.5*(1+z[i]);
388  }
389 
390  mode = m_bdata.data() + 2*numPoints;
391 
392  for(q = 2; q < numModes; ++q, mode+=numPoints)
393  {
394  Polylib::jacobfd(numPoints, z.data(), mode, NULL, q-2,1.0,1.0);
395 
396  for(i = 0; i < numPoints; ++i)
397  {
398  mode[i] *= m_bdata[i]*m_bdata[numPoints+i];
399  }
400  }
401 
402  // second row
403  for(i = 0; i < numPoints; ++i)
404  {
405  mode[i] = 0.5*(1-z[i]);
406  }
407 
408  mode += numPoints;
409 
410  for(q = 2; q < numModes; ++q, mode+=numPoints)
411  {
412  Polylib::jacobfd(numPoints, z.data(), mode, NULL, q-2,1.0,1.0);
413 
414  for(i = 0; i < numPoints; ++i)
415  {
416  mode[i] *= m_bdata[i]*m_bdata[numPoints+i];
417  }
418  }
419 
420  // third and higher rows
421  one_m_z_pow = m_bdata.data();
422  one_p_z = m_bdata.data()+numPoints;
423 
424  for(p = 2; p < numModes; ++p)
425  {
426  for(i = 0; i < numPoints; ++i)
427  {
428  mode[i] = m_bdata[i]*one_m_z_pow[i];
429  }
430 
431  one_m_z_pow = mode;
432  mode += numPoints;
433 
434  for(q = 1; q < numModes-p; ++q, mode+=numPoints)
435  {
436  Polylib::jacobfd(numPoints,z.data(),mode,NULL,q-1,2*p-1,1.0);
437 
438  for(i = 0; i < numPoints; ++i)
439  {
440  mode[i] *= one_m_z_pow[i]*one_p_z[i];
441  }
442  }
443  }
444 
445  Blas::Dgemm('n','n',numPoints,numModes*(numModes+1)/2,
446  numPoints,1.0,D,numPoints,
447  m_bdata.data(),numPoints,0.0,m_dbdata.data(),numPoints);
448  }
449  break;
450 
451 
452  case eModified_C:
453  {
454  // Note the following packing deviates from the
455  // definition in the Book by Karniadakis in two
456  // ways. 1) We put the vertex degrees of freedom
457  // at the lower index range to follow a more
458  // hierarchic structure. 2) We do not duplicate
459  // the singular vertex definition (or the
460  // duplicated face information in the book ) so
461  // that only a tetrahedral number
462  // (i.e. (modes)*(modes+1)*(modes+2)/6) of modes
463  // are required consistent with the orthogonal
464  // basis.
465 
466  // In the current structure the r index runs
467  // fastest rollowed by q and than the p index so
468  // that the matrix has a more compact structure
469 
470  // Note that eModified_C is a re-organisation/
471  // duplication of eModified_B so will get a
472  // temporary Modified_B expansion and copy the
473  // correct components.
474 
475  // Generate Modified_B basis;
476  BasisKey ModBKey(eModified_B,m_basisKey.GetNumModes(),
478  BasisSharedPtr ModB = BasisManager()[ModBKey];
479 
480  Array<OneD, const NekDouble> ModB_data = ModB->GetBdata();
481 
482  // Copy Modified_B basis into first
483  // (numModes*(numModes+1)/2)*numPoints entires of
484  // bdata. This fills in the complete (r,p) face.
485 
486  // Set up \phi^c_{p,q,r} = \phi^b_{p+q,r}
487 
488  int N;
489  int B_offset = 0;
490  int offset = 0;
491  for(p = 0; p < numModes; ++p)
492  {
493  N = numPoints*(numModes-p)*(numModes-p+1)/2;
494  Vmath::Vcopy(N, &ModB_data[0]+B_offset,1,&m_bdata[0] + offset,1);
495  B_offset += numPoints*(numModes-p);
496  offset += N;
497  }
498 
499  // set up derivative of basis.
500  Blas::Dgemm('n','n',numPoints,
501  numModes*(numModes+1)*(numModes+2)/6,
502  numPoints,1.0,D,numPoints,
503  m_bdata.data(),numPoints,0.0,
504  m_dbdata.data(),numPoints);
505  }
506  break;
507 
508  case eGLL_Lagrange:
509  {
510  mode = m_bdata.data();
511  boost::shared_ptr< Points<NekDouble> > m_points = PointsManager()[PointsKey(numModes, eGaussLobattoLegendre)];
512  const Array<OneD, const NekDouble>& zp(m_points->GetZ());
513 
514  for (p=0; p<numModes; ++p, mode += numPoints)
515  {
516  for(q = 0; q < numPoints; ++q)
517  {
518  mode[q] = Polylib::hglj(p, z[q], zp.data(), numModes, 0.0, 0.0);
519  }
520  }
521 
522  // define derivative basis
523  Blas::Dgemm('n', 'n', numPoints, numModes, numPoints, 1.0,
524  D, numPoints, m_bdata.data(), numPoints, 0.0,
525  m_dbdata.data(), numPoints);
526 
527  }//end scope
528  break;
529  case eGauss_Lagrange:
530  {
531  mode = m_bdata.data();
532  boost::shared_ptr< Points<NekDouble> > m_points = PointsManager()[PointsKey(numModes, eGaussGaussLegendre)];
533  const Array<OneD, const NekDouble>& zp(m_points->GetZ());
534 
535  for (p=0; p<numModes; ++p,mode += numPoints)
536  {
537  for(q = 0; q < numPoints; ++q)
538  {
539  mode[q] = Polylib::hgj(p, z[q], zp.data(), numModes, 0.0, 0.0);
540  }
541  }
542 
543  // define derivative basis
544  Blas::Dgemm('n', 'n', numPoints, numModes, numPoints, 1.0,
545  D, numPoints, m_bdata.data(), numPoints, 0.0,
546  m_dbdata.data(), numPoints);
547 
548  }//end scope
549  break;
550  case eFourier:
551 
552  ASSERTL0(numModes%2==0, "Fourier modes should be a factor of 2");
553 
554  for(i = 0; i < numPoints; ++i)
555  {
556  m_bdata[i] = 1.0;
557  m_bdata[numPoints+i] = 0.0;
558 
559  m_dbdata[i] = m_dbdata[numPoints+i] = 0.0;
560  }
561 
562  for (p=1; p < numModes/2; ++p)
563  {
564  for(i = 0; i < numPoints; ++i)
565  {
566  m_bdata[ 2*p *numPoints+i] = cos(p*M_PI* (z[i]+1) );
567  m_bdata[(2*p+1)*numPoints+i] = -sin(p*M_PI* (z[i]+1) );
568 
569  m_dbdata[ 2*p *numPoints+i] = -p*M_PI*sin(p*M_PI* (z[i]+1) );
570  m_dbdata[(2*p+1)*numPoints+i] = -p*M_PI*cos(p*M_PI* (z[i]+1) );
571  }
572  }
573 
574  break;
575 
576 
577  // Fourier Single Mode (1st mode)
578  case eFourierSingleMode:
579 
580  for(i = 0; i < numPoints; ++i)
581  {
582  m_bdata[i] = cos(M_PI* (z[i]+1) );
583  m_bdata[numPoints+i] = -sin(M_PI* (z[i]+1) );
584 
585  m_dbdata[i] = -M_PI*sin(M_PI* (z[i]+1) );
586  m_dbdata[numPoints+i] = -M_PI*cos(M_PI* (z[i]+1) );
587  }
588 
589  for (p=1; p < numModes/2; ++p)
590  {
591  for(i = 0; i < numPoints; ++i)
592  {
593  m_bdata[ 2*p *numPoints+i] = 0.;
594  m_bdata[(2*p+1)*numPoints+i] = 0.;
595 
596  m_dbdata[ 2*p *numPoints+i] = 0.;
597  m_dbdata[(2*p+1)*numPoints+i] = 0.;
598  }
599  }
600  break;
601 
602  //Fourier Real Half Mode
603  case eFourierHalfModeRe:
604  m_bdata[0] = cos(M_PI*z[0]);
605  m_dbdata[0] = -M_PI*sin(M_PI*z[0]);
606  break;
607 
608  //Fourier Imaginary Half Mode
609  case eFourierHalfModeIm:
610  m_bdata[0] = -sin(M_PI*z[0]);
611  m_dbdata[0] = -M_PI*cos(M_PI*z[0]);
612  break;
613 
614  case eChebyshev:
615  {
616  mode = m_bdata.data();
617 
618  for (p=0,scal = 1; p<numModes; ++p,mode += numPoints)
619  {
620  Polylib::jacobfd(numPoints, z.data(), mode, NULL, p, -0.5, -0.5);
621 
622  for(i = 0; i < numPoints; ++i)
623  {
624  mode[i] *= scal;
625  }
626 
627  scal *= 4*(p+1)*(p+1)/(2*p+2)/(2*p+1);
628  }
629 
630  // Define derivative basis
631  Blas::Dgemm('n', 'n', numPoints, numModes, numPoints, 1.0,
632  D, numPoints, m_bdata.data(), numPoints, 0.0,
633  m_dbdata.data(), numPoints);
634  }
635  break;
636 
637  case eMonomial:
638  {
639  int P = numModes - 1;
640  NekDouble *mode = m_bdata.data();
641 
642  for( int p = 0; p <= P; ++p, mode += numPoints )
643  {
644  for( int i = 0; i < numPoints; ++i )
645  {
646  mode[i] = pow(z[i], p);
647  }
648  }
649 
650  // define derivative basis
651  Blas::Dgemm('n', 'n', numPoints, numModes, numPoints, 1.0,
652  D, numPoints, m_bdata.data(), numPoints, 0.0,
653  m_dbdata.data(),numPoints);
654  }//end scope
655  break;
656  default:
657  ASSERTL0(false, "Basis Type not known or "
658  "not implemented at this time.");
659  }
660  }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
BasisKey m_basisKey
Basis specification.
Definition: Basis.h:328
Principle Modified Functions .
Definition: BasisType.h:51
double hgj(const int i, const double z, const double *zgj, const int np, const double alpha, const double beta)
Compute the value of the i th Lagrangian interpolant through.
Definition: Polylib.cpp:1340
Array< OneD, NekDouble > m_bdata
Basis definition.
Definition: Basis.h:330
Principle Modified Functions .
Definition: BasisType.h:49
Lagrange Polynomials using the Gauss points .
Definition: BasisType.h:54
Fourier Expansion .
Definition: BasisType.h:52
Chebyshev Polynomials .
Definition: BasisType.h:56
1D Gauss-Gauss-Legendre quadrature points
Definition: PointsType.h:47
Principle Orthogonal Functions .
Definition: BasisType.h:47
BasisManagerT & BasisManager(void)
Array< OneD, NekDouble > m_dbdata
Derivative Basis definition.
Definition: Basis.h:331
int GetNumModes() const
Return order of basis from the basis specification.
Definition: Basis.h:218
Fourier Modified expansions with just the real part of the first mode .
Definition: BasisType.h:59
Principle Modified Functions .
Definition: BasisType.h:50
PointsSharedPtr m_points
Set of points.
Definition: Basis.h:329
Principle Orthogonal Functions .
Definition: BasisType.h:48
PointsManagerT & PointsManager(void)
Principle Orthogonal Functions .
Definition: BasisType.h:46
double NekDouble
Fourier Modified expansions with just the imaginary part of the first mode .
Definition: BasisType.h:60
PointsKey GetPointsKey() const
Return distribution of points.
Definition: Basis.h:145
Fourier ModifiedExpansion with just the first mode .
Definition: BasisType.h:58
int GetNumPoints() const
Return the number of points from the basis specification.
Definition: Basis.h:230
Legendre Polynomials . Same as Ortho_A.
Definition: BasisType.h:55
double hglj(const int i, const double z, const double *zglj, const int np, const double alpha, const double beta)
Compute the value of the i th Lagrangian interpolant through the.
Definition: Polylib.cpp:1584
int GetNumModes() const
Returns the order of the basis.
Definition: Basis.h:84
boost::shared_ptr< Basis > BasisSharedPtr
Lagrange for SEM basis .
Definition: BasisType.h:53
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
BasisType GetBasisType() const
Return the type of expansion basis.
Definition: Basis.h:242
void jacobfd(const int np, const double *z, double *poly_in, double *polyd, const int n, const double alpha, const double beta)
Routine to calculate Jacobi polynomials, , and their first derivative, .
Definition: Polylib.cpp:1946
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:50
Monomial polynomials .
Definition: BasisType.h:57
const BasisKey Nektar::LibUtilities::Basis::GetBasisKey ( ) const
inline

Returns the specification for the Basis.

Definition at line 320 of file Basis.h.

References m_basisKey.

321  {
322  return m_basisKey;
323  }
BasisKey m_basisKey
Basis specification.
Definition: Basis.h:328
BasisType Nektar::LibUtilities::Basis::GetBasisType ( ) const
inline

Return the type of expansion basis.

Definition at line 242 of file Basis.h.

References Nektar::LibUtilities::BasisKey::GetBasisType(), and m_basisKey.

Referenced by GenBasis().

243  {
244  return m_basisKey.GetBasisType();
245  }
BasisKey m_basisKey
Basis specification.
Definition: Basis.h:328
BasisType GetBasisType() const
Return type of expansion basis.
Definition: Basis.h:139
const Array<OneD, const NekDouble>& Nektar::LibUtilities::Basis::GetBdata ( ) const
inline

Return basis definition array m_bdata.

Definition at line 308 of file Basis.h.

References m_bdata.

309  {
310  return m_bdata;
311  }
Array< OneD, NekDouble > m_bdata
Basis definition.
Definition: Basis.h:330
const boost::shared_ptr<NekMatrix<NekDouble> >& Nektar::LibUtilities::Basis::GetD ( Direction  dir = xDir) const
inline

Definition at line 275 of file Basis.h.

References m_points.

277  {
278  return m_points->GetD(dir);
279  }
PointsSharedPtr m_points
Set of points.
Definition: Basis.h:329
const Array<OneD, const NekDouble>& Nektar::LibUtilities::Basis::GetDbdata ( ) const
inline

Return basis definition array m_dbdata.

Definition at line 314 of file Basis.h.

References m_dbdata.

315  {
316  return m_dbdata;
317  }
Array< OneD, NekDouble > m_dbdata
Derivative Basis definition.
Definition: Basis.h:331
const boost::shared_ptr<NekMatrix<NekDouble> > Nektar::LibUtilities::Basis::GetI ( const Array< OneD, const NekDouble > &  x)
inline

Definition at line 281 of file Basis.h.

References m_points.

283  {
284  return m_points->GetI(x);
285  }
PointsSharedPtr m_points
Set of points.
Definition: Basis.h:329
const boost::shared_ptr<NekMatrix<NekDouble> > Nektar::LibUtilities::Basis::GetI ( const BasisKey bkey)
inline

Definition at line 287 of file Basis.h.

References ASSERTL0, Nektar::LibUtilities::PointsKey::GetPointsDim(), Nektar::LibUtilities::BasisKey::GetPointsKey(), and m_InterpManager.

289  {
290  ASSERTL0(bkey.GetPointsKey().GetPointsDim()==1,
291  "Interpolation only to other 1d basis");
292  return m_InterpManager[bkey];
293  }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
NekManager< BasisKey, NekMatrix< NekDouble >, BasisKey::opLess > m_InterpManager
Definition: Basis.h:333
int Nektar::LibUtilities::Basis::GetNumModes ( ) const
inline

Return order of basis from the basis specification.

Definition at line 218 of file Basis.h.

References Nektar::LibUtilities::BasisKey::GetNumModes(), and m_basisKey.

Referenced by GenBasis(), and Initialize().

219  {
220  return m_basisKey.GetNumModes();
221  }
BasisKey m_basisKey
Basis specification.
Definition: Basis.h:328
int GetNumModes() const
Returns the order of the basis.
Definition: Basis.h:84
int Nektar::LibUtilities::Basis::GetNumPoints ( ) const
inline

Return the number of points from the basis specification.

Definition at line 230 of file Basis.h.

References Nektar::LibUtilities::BasisKey::GetNumPoints(), and m_basisKey.

Referenced by GenBasis().

231  {
232  return m_basisKey.GetNumPoints();
233  }
BasisKey m_basisKey
Basis specification.
Definition: Basis.h:328
int GetNumPoints() const
Return points order at which basis is defined.
Definition: Basis.h:128
PointsKey Nektar::LibUtilities::Basis::GetPointsKey ( ) const
inline

Return the points specification for the basis.

Definition at line 248 of file Basis.h.

References Nektar::LibUtilities::BasisKey::GetPointsKey(), and m_basisKey.

249  {
250  return m_basisKey.GetPointsKey();
251  }
BasisKey m_basisKey
Basis specification.
Definition: Basis.h:328
PointsKey GetPointsKey() const
Return distribution of points.
Definition: Basis.h:145
PointsType Nektar::LibUtilities::Basis::GetPointsType ( ) const
inline

Return the type of quadrature.

Definition at line 254 of file Basis.h.

References Nektar::LibUtilities::BasisKey::GetPointsType(), and m_basisKey.

255  {
256  return m_basisKey.GetPointsType();
257  }
BasisKey m_basisKey
Basis specification.
Definition: Basis.h:328
PointsType GetPointsType() const
Return type of quadrature.
Definition: Basis.h:151
int Nektar::LibUtilities::Basis::GetTotNumModes ( ) const
inline

Return total number of modes from the basis specification.

Definition at line 224 of file Basis.h.

References Nektar::LibUtilities::BasisKey::GetTotNumModes(), and m_basisKey.

225  {
226  return m_basisKey.GetTotNumModes();
227  }
BasisKey m_basisKey
Basis specification.
Definition: Basis.h:328
int GetTotNumModes() const
Definition: Basis.h:90
int Nektar::LibUtilities::Basis::GetTotNumPoints ( ) const
inline

Return total number of points from the basis specification.

Definition at line 236 of file Basis.h.

References Nektar::LibUtilities::BasisKey::GetTotNumPoints(), and m_basisKey.

Referenced by Initialize().

237  {
238  return m_basisKey.GetTotNumPoints();
239  }
BasisKey m_basisKey
Basis specification.
Definition: Basis.h:328
int GetTotNumPoints() const
Definition: Basis.h:133
const Array<OneD, const NekDouble>& Nektar::LibUtilities::Basis::GetW ( ) const
inline

Definition at line 264 of file Basis.h.

References m_points.

265  {
266  return m_points->GetW();
267  }
PointsSharedPtr m_points
Set of points.
Definition: Basis.h:329
const Array<OneD, const NekDouble>& Nektar::LibUtilities::Basis::GetZ ( ) const
inline

Definition at line 259 of file Basis.h.

References m_points.

260  {
261  return m_points->GetZ();
262  }
PointsSharedPtr m_points
Set of points.
Definition: Basis.h:329
void Nektar::LibUtilities::Basis::GetZW ( Array< OneD, const NekDouble > &  z,
Array< OneD, const NekDouble > &  w 
) const
inline

Definition at line 269 of file Basis.h.

References m_points.

271  {
272  m_points->GetZW(z,w);
273  }
PointsSharedPtr m_points
Set of points.
Definition: Basis.h:329
void Nektar::LibUtilities::Basis::Initialize ( )
virtual

Definition at line 110 of file Basis.cpp.

References ASSERTL0, GenBasis(), GetNumModes(), and GetTotNumPoints().

112  {
113  ASSERTL0(GetNumModes()>0, "Cannot call Basis initialisation with zero or negative order");
114  ASSERTL0(GetTotNumPoints()>0, "Cannot call Basis initialisation with zero or negative numbers of points");
115 
116  GenBasis();
117  };
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
int GetNumModes() const
Return order of basis from the basis specification.
Definition: Basis.h:218
int GetTotNumPoints() const
Return total number of points from the basis specification.
Definition: Basis.h:236
void GenBasis()
Generate appropriate basis and their derivatives.
Definition: Basis.cpp:215

Member Data Documentation

BasisKey Nektar::LibUtilities::Basis::m_basisKey
protected
Array<OneD, NekDouble> Nektar::LibUtilities::Basis::m_bdata
protected

Basis definition.

Definition at line 330 of file Basis.h.

Referenced by GenBasis(), and GetBdata().

Array<OneD, NekDouble> Nektar::LibUtilities::Basis::m_dbdata
protected

Derivative Basis definition.

Definition at line 331 of file Basis.h.

Referenced by GenBasis(), and GetDbdata().

NekManager<BasisKey, NekMatrix<NekDouble>, BasisKey::opLess> Nektar::LibUtilities::Basis::m_InterpManager
protected

Definition at line 333 of file Basis.h.

Referenced by Basis(), and GetI().

PointsSharedPtr Nektar::LibUtilities::Basis::m_points
protected

Set of points.

Definition at line 329 of file Basis.h.

Referenced by GenBasis(), GetD(), GetI(), GetW(), GetZ(), and GetZW().