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

Specialisation of the NodalUtil class to support nodal tetrahedral elements. More...

#include <NodalUtil.h>

Inheritance diagram for Nektar::LibUtilities::NodalUtilTetrahedron:
Inheritance graph
[legend]
Collaboration diagram for Nektar::LibUtilities::NodalUtilTetrahedron:
Collaboration graph
[legend]

Public Member Functions

 NodalUtilTetrahedron (int degree, Array< OneD, NekDouble > r, Array< OneD, NekDouble > s, Array< OneD, NekDouble > t)
 Construct the nodal utility class for a tetrahedron. More...
 
virtual ~NodalUtilTetrahedron ()
 
- Public Member Functions inherited from Nektar::LibUtilities::NodalUtil
NekVector< NekDoubleGetWeights ()
 Obtain the integration weights for the given nodal distribution. More...
 
SharedMatrix GetVandermonde ()
 Return the Vandermonde matrix for the nodal distribution. More...
 
SharedMatrix GetVandermondeForDeriv (int dir)
 Return the Vandermonde matrix of the derivative of the basis functions for the nodal distribution. More...
 
SharedMatrix GetDerivMatrix (int dir)
 Return the derivative matrix for the nodal distribution. More...
 
SharedMatrix GetInterpolationMatrix (Array< OneD, Array< OneD, NekDouble > > &xi)
 Construct the interpolation matrix used to evaluate the basis at the points xi inside the element. More...
 

Protected Member Functions

virtual NekVector< NekDoublev_OrthoBasis (const int mode)
 Return the value of the modal functions for the tetrahedral element at the nodal points m_xi for a given mode. More...
 
virtual NekVector< NekDoublev_OrthoBasisDeriv (const int dir, const int mode)
 Return the value of the derivative of the modal functions for the tetrahedral element at the nodal points m_xi for a given mode. More...
 
virtual boost::shared_ptr
< NodalUtil
v_CreateUtil (Array< OneD, Array< OneD, NekDouble > > &xi)
 Construct a NodalUtil object of the appropriate element type for a given set of points. More...
 
virtual NekDouble v_ModeZeroIntegral ()
 Return the value of the integral of the zero-th mode for this element. More...
 
virtual int v_NumModes ()
 Calculate the number of degrees of freedom for this element. More...
 
- Protected Member Functions inherited from Nektar::LibUtilities::NodalUtil
 NodalUtil (int degree, int dim)
 Set up the NodalUtil object. More...
 

Protected Attributes

std::vector< Modem_ordering
 Mapping from the $ (i,j,k) $ indexing of the basis to a continuous ordering. More...
 
Array< OneD, Array< OneD,
NekDouble > > 
m_eta
 Collapsed coordinates $ (\eta_1, \eta_2, \eta_3) $ of the nodal points. More...
 
- Protected Attributes inherited from Nektar::LibUtilities::NodalUtil
int m_dim
 Dimension of the nodal element. More...
 
int m_degree
 Degree of the nodal element. More...
 
int m_numPoints
 Total number of nodal points. More...
 
Array< OneD, Array< OneD,
NekDouble > > 
m_xi
 Coordinates of the nodal points defining the basis. More...
 

Private Types

typedef boost::tuple< int, int,
int > 
Mode
 

Detailed Description

Specialisation of the NodalUtil class to support nodal tetrahedral elements.

Definition at line 215 of file NodalUtil.h.

Member Typedef Documentation

typedef boost::tuple<int, int, int> Nektar::LibUtilities::NodalUtilTetrahedron::Mode
private

Definition at line 217 of file NodalUtil.h.

Constructor & Destructor Documentation

Nektar::LibUtilities::NodalUtilTetrahedron::NodalUtilTetrahedron ( int  degree,
Array< OneD, NekDouble r,
Array< OneD, NekDouble s,
Array< OneD, NekDouble t 
)

Construct the nodal utility class for a tetrahedron.

The constructor of this class sets up two member variables used in the evaluation of the orthogonal basis:

  • NodalUtilTetrahedron::m_eta is used to construct the collapsed coordinate locations of the nodal points $ (\eta_1, \eta_2, \eta_3) $ inside the cube $[-1,1]^3$ on which the orthogonal basis functions are defined.
  • NodalUtilTetrahedron::m_ordering constructs a mapping from the index set $ I = \{ (i,j,k)\ |\ 0\leq i,j,k \leq P, i+j \leq P, i+j+k \leq P \}$ to an ordering $ 0 \leq m(ijk) \leq (P+1)(P+2)(P+3)/6 $ that defines the monomials $ \xi_1^i \xi_2^j \xi_3^k $ that span the tetrahedral space. This is then used to calculate which $ (i,j,k) $ triple (represented as a boost tuple) corresponding to a column of the Vandermonde matrix when calculating the orthogonal polynomials.
Parameters
degreePolynomial order of this nodal tetrahedron
r$ \xi_1 $-coordinates of nodal points in the standard element.
s$ \xi_2 $-coordinates of nodal points in the standard element.
t$ \xi_3 $-coordinates of nodal points in the standard element.

Definition at line 420 of file NodalUtil.cpp.

References Nektar::NekConstants::kNekZeroTol, Nektar::LibUtilities::NodalUtil::m_degree, m_eta, Nektar::LibUtilities::NodalUtil::m_numPoints, m_ordering, and Nektar::LibUtilities::NodalUtil::m_xi.

424  : NodalUtil(degree, 3), m_eta(3)
425 {
426  m_numPoints = r.num_elements();
427  m_xi[0] = r;
428  m_xi[1] = s;
429  m_xi[2] = t;
430 
431  for (int i = 0; i <= m_degree; ++i)
432  {
433  for (int j = 0; j <= m_degree - i; ++j)
434  {
435  for (int k = 0; k <= m_degree - i - j; ++k)
436  {
437  m_ordering.push_back(Mode(i, j, k));
438  }
439  }
440  }
441 
442  // Calculate collapsed coordinates from r/s values
443  m_eta[0] = Array<OneD, NekDouble>(m_numPoints);
444  m_eta[1] = Array<OneD, NekDouble>(m_numPoints);
445  m_eta[2] = Array<OneD, NekDouble>(m_numPoints);
446 
447  for (int i = 0; i < m_numPoints; ++i)
448  {
449  if (fabs(m_xi[2][i] - 1.0) < NekConstants::kNekZeroTol)
450  {
451  // Very top point of the tetrahedron
452  m_eta[0][i] = -1.0;
453  m_eta[1][i] = -1.0;
454  m_eta[2][i] = m_xi[2][i];
455  }
456  else
457  {
458  if (fabs(m_xi[1][i] - 1.0) < NekConstants::kNekZeroTol)
459  {
460  // Distant diagonal edge shared by all eta_x coordinate planes:
461  // the xi_y == -xi_z line
462  m_eta[0][i] = -1.0;
463  }
464  else if (fabs(m_xi[1][i] + m_xi[2][i]) < NekConstants::kNekZeroTol)
465  {
466  m_eta[0][i] = -1.0;
467  }
468  else
469  {
470  m_eta[0][i] = 2.0 * (1.0 + m_xi[0][i]) /
471  (-m_xi[1][i] - m_xi[2][i]) - 1.0;
472  }
473  m_eta[1][i] = 2.0 * (1.0 + m_xi[1][i]) / (1.0 - m_xi[2][i]) - 1.0;
474  m_eta[2][i] = m_xi[2][i];
475  }
476  }
477 }
std::vector< Mode > m_ordering
Mapping from the indexing of the basis to a continuous ordering.
Definition: NodalUtil.h:232
Array< OneD, Array< OneD, NekDouble > > m_eta
Collapsed coordinates of the nodal points.
Definition: NodalUtil.h:236
boost::tuple< int, int, int > Mode
Definition: NodalUtil.h:217
int m_degree
Degree of the nodal element.
Definition: NodalUtil.h:108
static const NekDouble kNekZeroTol
NodalUtil(int degree, int dim)
Set up the NodalUtil object.
Definition: NodalUtil.h:101
Array< OneD, Array< OneD, NekDouble > > m_xi
Coordinates of the nodal points defining the basis.
Definition: NodalUtil.h:112
int m_numPoints
Total number of nodal points.
Definition: NodalUtil.h:110
virtual Nektar::LibUtilities::NodalUtilTetrahedron::~NodalUtilTetrahedron ( )
inlinevirtual

Definition at line 225 of file NodalUtil.h.

226  {
227  }

Member Function Documentation

virtual boost::shared_ptr<NodalUtil> Nektar::LibUtilities::NodalUtilTetrahedron::v_CreateUtil ( Array< OneD, Array< OneD, NekDouble > > &  xi)
inlineprotectedvirtual

Construct a NodalUtil object of the appropriate element type for a given set of points.

This function is used inside NodalUtil::GetInterpolationMatrix so that the (potentially non-square) Vandermonde matrix can be constructed to create the interpolation matrix at an arbitrary set of points in the domain.

Parameters
xiDistribution of nodal points to create utility with.

Implements Nektar::LibUtilities::NodalUtil.

Definition at line 242 of file NodalUtil.h.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), and Nektar::LibUtilities::NodalUtil::m_degree.

244  {
246  m_degree, xi[0], xi[1], xi[2]);
247  }
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
int m_degree
Degree of the nodal element.
Definition: NodalUtil.h:108
virtual NekDouble Nektar::LibUtilities::NodalUtilTetrahedron::v_ModeZeroIntegral ( )
inlineprotectedvirtual

Return the value of the integral of the zero-th mode for this element.

Note that for the orthogonal basis under consideration, all modes integrate to zero asides from the zero-th mode. This function is used in NodalUtil::GetWeights to determine integration weights.

Implements Nektar::LibUtilities::NodalUtil.

Definition at line 249 of file NodalUtil.h.

250  {
251  return 8.0 * sqrt(2.0) / 3.0;
252  }
virtual int Nektar::LibUtilities::NodalUtilTetrahedron::v_NumModes ( )
inlineprotectedvirtual

Calculate the number of degrees of freedom for this element.

Implements Nektar::LibUtilities::NodalUtil.

Definition at line 254 of file NodalUtil.h.

References Nektar::LibUtilities::NodalUtil::m_degree.

255  {
256  return (m_degree + 1) * (m_degree + 2) * (m_degree + 3) / 6;
257  }
int m_degree
Degree of the nodal element.
Definition: NodalUtil.h:108
NekVector< NekDouble > Nektar::LibUtilities::NodalUtilTetrahedron::v_OrthoBasis ( const int  mode)
protectedvirtual

Return the value of the modal functions for the tetrahedral element at the nodal points m_xi for a given mode.

In a tetrahedron, we use the orthogonal basis

\[ \psi_{m(ijk)} = \sqrt{8} P^{(0,0)}_i(\xi_1) P_j^{(2i+1,0)}(\xi_2) P_k^{(2i+2j+2,0)}(\xi_3) (1-\xi_2)^i (1-\xi_3)^{i+j} \]

where $ m(ijk) $ is the mapping defined in m_ordering and $ J_n^{(\alpha,\beta)}(z) $ denotes the standard Jacobi polynomial.

Parameters
modeThe mode of the orthogonal basis to evaluate.
Returns
Vector containing orthogonal basis evaluated at the points m_xi.

Implements Nektar::LibUtilities::NodalUtil.

Definition at line 495 of file NodalUtil.cpp.

References Polylib::jacobfd(), m_eta, Nektar::LibUtilities::NodalUtil::m_numPoints, m_ordering, and CG_Iterations::modes.

496 {
497  std::vector<NekDouble> jacA(m_numPoints), jacB(m_numPoints);
498  std::vector<NekDouble> jacC(m_numPoints);
499  Mode modes = m_ordering[mode];
500 
501  const int I = modes.get<0>(), J = modes.get<1>(), K = modes.get<2>();
502 
503  // Calculate Jacobi polynomials
505  m_numPoints, &m_eta[0][0], &jacA[0], NULL, I, 0.0, 0.0);
507  m_numPoints, &m_eta[1][0], &jacB[0], NULL, J, 2.0 * I + 1.0, 0.0);
509  m_numPoints, &m_eta[2][0], &jacC[0], NULL, K, 2.0 * (I+J) + 2.0, 0.0);
510 
511  NekVector<NekDouble> ret(m_numPoints);
512  NekDouble sqrt8 = sqrt(8.0);
513 
514  for (int i = 0; i < m_numPoints; ++i)
515  {
516  ret[i] = sqrt8 * jacA[i] * jacB[i] * jacC[i] *
517  pow(1.0 - m_eta[1][i], I) * pow(1.0 - m_eta[2][i], I + J);
518  }
519 
520  return ret;
521 }
std::vector< Mode > m_ordering
Mapping from the indexing of the basis to a continuous ordering.
Definition: NodalUtil.h:232
Array< OneD, Array< OneD, NekDouble > > m_eta
Collapsed coordinates of the nodal points.
Definition: NodalUtil.h:236
boost::tuple< int, int, int > Mode
Definition: NodalUtil.h:217
double NekDouble
int m_numPoints
Total number of nodal points.
Definition: NodalUtil.h:110
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:1920
NekVector< NekDouble > Nektar::LibUtilities::NodalUtilTetrahedron::v_OrthoBasisDeriv ( const int  dir,
const int  mode 
)
protectedvirtual

Return the value of the derivative of the modal functions for the tetrahedral element at the nodal points m_xi for a given mode.

Note that this routine must use the chain rule combined with the collapsed coordinate derivatives as described in Sherwin & Karniadakis (2nd edition), pg 152.

Parameters
dirCoordinate direction in which to evaluate the derivative.
modeThe mode of the orthogonal basis to evaluate.
Returns
Vector containing the derivative of the orthogonal basis evaluated at the points m_xi.

Implements Nektar::LibUtilities::NodalUtil.

Definition at line 537 of file NodalUtil.cpp.

References Polylib::jacobd(), Polylib::jacobfd(), m_eta, Nektar::LibUtilities::NodalUtil::m_numPoints, m_ordering, and CG_Iterations::modes.

539 {
540  std::vector<NekDouble> jacA(m_numPoints), jacB(m_numPoints);
541  std::vector<NekDouble> jacC(m_numPoints);
542  std::vector<NekDouble> jacDerivA(m_numPoints), jacDerivB(m_numPoints);
543  std::vector<NekDouble> jacDerivC(m_numPoints);
544  Mode modes = m_ordering[mode];
545 
546  const int I = modes.get<0>(), J = modes.get<1>(), K = modes.get<2>();
547 
548  // Calculate Jacobi polynomials
550  m_numPoints, &m_eta[0][0], &jacA[0], NULL, I, 0.0, 0.0);
552  m_numPoints, &m_eta[1][0], &jacB[0], NULL, J, 2.0 * I + 1.0, 0.0);
554  m_numPoints, &m_eta[2][0], &jacC[0], NULL, K, 2.0 * (I+J) + 2.0, 0.0);
556  m_numPoints, &m_eta[0][0], &jacDerivA[0], I, 0.0, 0.0);
558  m_numPoints, &m_eta[1][0], &jacDerivB[0], J, 2.0 * I + 1.0, 0.0);
560  m_numPoints, &m_eta[2][0], &jacDerivC[0], K, 2.0 * (I+J) + 2.0, 0.0);
561 
562  NekVector<NekDouble> ret(m_numPoints);
563  NekDouble sqrt8 = sqrt(8.0);
564 
565  // Always compute x-derivative since this term appears in the latter two
566  // terms.
567  for (int i = 0; i < m_numPoints; ++i)
568  {
569  ret[i] = 4.0 * sqrt8 * jacDerivA[i] * jacB[i] * jacC[i];
570 
571  if (I > 0)
572  {
573  ret[i] *= pow(1 - m_eta[1][i], I - 1);
574  }
575 
576  if (I + J > 0)
577  {
578  ret[i] *= pow(1 - m_eta[2][i], I + J - 1);
579  }
580  }
581 
582  if (dir >= 1)
583  {
584  // Multiply by (1+a)/2
585  NekVector<NekDouble> tmp(m_numPoints);
586 
587  for (int i = 0; i < m_numPoints; ++i)
588  {
589  ret[i] *= 0.5 * (m_eta[0][i] + 1.0);
590 
591  tmp[i] = 2.0 * sqrt8 * jacA[i] * jacC[i];
592  if (I + J > 0)
593  {
594  tmp[i] *= pow(1.0 - m_eta[2][i], I + J - 1);
595  }
596 
597  NekDouble tmp2 = jacDerivB[i] * pow(1.0 - m_eta[1][i], I);
598  if (I > 0)
599  {
600  tmp2 -= I * jacB[i] * pow(1.0 - m_eta[1][i], I - 1);
601  }
602 
603  tmp[i] *= tmp2;
604  }
605 
606  if (dir == 1)
607  {
608  ret += tmp;
609  return ret;
610  }
611 
612  for (int i = 0; i < m_numPoints; ++i)
613  {
614  ret[i] += 0.5 * (1.0 + m_eta[1][i]) * tmp[i];
615 
616  NekDouble tmp2 = jacDerivC[i] * pow(1.0 - m_eta[2][i], I + J);
617  if (I + J > 0)
618  {
619  tmp2 -= jacC[i] * (I + J) * pow(1.0 - m_eta[2][i], I + J - 1);
620  }
621 
622  ret[i] += sqrt8 * jacA[i] * jacB[i] * pow(1.0 - m_eta[1][i], I) *
623  tmp2;
624  }
625  }
626 
627  return ret;
628 }
std::vector< Mode > m_ordering
Mapping from the indexing of the basis to a continuous ordering.
Definition: NodalUtil.h:232
Array< OneD, Array< OneD, NekDouble > > m_eta
Collapsed coordinates of the nodal points.
Definition: NodalUtil.h:236
boost::tuple< int, int, int > Mode
Definition: NodalUtil.h:217
void jacobd(const int np, const double *z, double *polyd, const int n, const double alpha, const double beta)
Calculate the derivative of Jacobi polynomials.
Definition: Polylib.cpp:2120
double NekDouble
int m_numPoints
Total number of nodal points.
Definition: NodalUtil.h:110
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:1920

Member Data Documentation

Array<OneD, Array<OneD, NekDouble> > Nektar::LibUtilities::NodalUtilTetrahedron::m_eta
protected

Collapsed coordinates $ (\eta_1, \eta_2, \eta_3) $ of the nodal points.

Definition at line 236 of file NodalUtil.h.

Referenced by NodalUtilTetrahedron(), v_OrthoBasis(), and v_OrthoBasisDeriv().

std::vector<Mode> Nektar::LibUtilities::NodalUtilTetrahedron::m_ordering
protected

Mapping from the $ (i,j,k) $ indexing of the basis to a continuous ordering.

Definition at line 232 of file NodalUtil.h.

Referenced by NodalUtilTetrahedron(), v_OrthoBasis(), and v_OrthoBasisDeriv().