A class to assist in the construction of nodal simplex and hybrid elements in two and three dimensions. More...

#include <NodalUtil.h>

Inheritance diagram for Nektar::LibUtilities::NodalUtil:

Collaboration diagram for Nektar::LibUtilities::NodalUtil:

Public Member Functions
NekVector< NekDouble >	GetWeights ()
	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
	NodalUtil (int degree, int dim)
	Set up the NodalUtil object. More...

virtual NekVector< NekDouble >	v_OrthoBasis (const int mode)=0
	Return the values of the orthogonal basis at the nodal points for a given mode. More...

virtual NekVector< NekDouble >	v_OrthoBasisDeriv (const int dir, const int mode)=0
	Return the values of the derivative of the orthogonal basis at the nodal points for a given mode. More...

virtual boost::shared_ptr < NodalUtil >	v_CreateUtil (Array< OneD, Array< OneD, NekDouble > > &xi)=0
	Construct a NodalUtil object of the appropriate element type for a given set of points. More...

virtual NekDouble	v_ModeZeroIntegral ()=0
	Return the value of the integral of the zero-th mode for this element. More...

virtual int	v_NumModes ()=0
	Calculate the number of degrees of freedom for this element. More...

Protected Attributes
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...

Detailed Description

A class to assist in the construction of nodal simplex and hybrid elements in two and three dimensions.

The NodalUtil class and its subclasses are designed to take care of some common issues that arise when considering triangles, tetrahedra and prismatic elements that are equipped with a nodal Lagrangian basis, defined using a set of nodal points $\xi_i$ that we store in the array NodalUtil::m_xi. Since one cannot write this basis analytically, we instead construct the Vandermonde matrix

$\mathbf{V}_{ij} = \psi_j(\xi_i)$

where $\psi_j$ is a basis that spans the polynomial space of the element. The Vandermonde matrix can then be used to construct the integration weights, derivative and interpolation matrices. Although this can be any basis, such as the monomial basis $ x^i y^j z^k $ , in practice this is numerically unstable at high polynomial orders. Elements are therefore expected to use the 'traditional' modal orthogonal basis. See Sherwin & Karniadakis or Hesthaven & Warburton for further details of this basis and the surrounding numerical issues.

This class therefore contains the generic logic needed to construct various matrices, and subclasses override virtual functions that define the orthogonal basis and its derivatives for a particular element type.

Definition at line 84 of file NodalUtil.h.

Constructor & Destructor Documentation

Nektar::LibUtilities::NodalUtil::NodalUtil	(	int	degree,
		int	dim
	)

inlineprotected

Set up the NodalUtil object.

Parameters

dim	Dimension of the element.
degree	Polynomial degree of the element.

Definition at line 101 of file NodalUtil.h.

                                    : m_dim(dim), m_degree(degree), m_xi(dim)
     {
     }

Member Function Documentation

SharedMatrix Nektar::LibUtilities::NodalUtil::GetDerivMatrix ( int dir )

Return the derivative matrix for the nodal distribution.

This routine constructs and returns the derivative matrices $\mathbf{D}_d$ for coordinate directions $0 \leq d \leq 2$ , which can be used to evaluate the derivative of a nodal expansion at the points defined by NodalUtil::m_xi. These are calculated as $\mathbf{D}_d = \mathbf{V}_d \mathbf{V}^{-1}$ , where $\mathbf{V}_d$ is the derivative Vandermonde matrix and $\mathbf{V}$ is the Vandermonde matrix.

Parameters

dir	Coordinate direction in which to evaluate the derivative.

Returns: The derivative matrix for direction dir.

Definition at line 168 of file NodalUtil.cpp.

References GetVandermonde(), and GetVandermondeForDeriv().

 {
     SharedMatrix V  = GetVandermonde();
     SharedMatrix Vd = GetVandermondeForDeriv(dir);
     SharedMatrix D = MemoryManager<NekMatrix<NekDouble> >::AllocateSharedPtr(
         V->GetRows(), V->GetColumns(), 0.0);
 
     V->Invert();
 
     *D = (*Vd) * (*V);
 
     return D;
 }

SharedMatrix Nektar::LibUtilities::NodalUtil::GetInterpolationMatrix ( Array< OneD, Array< OneD, NekDouble > > & xi )

Construct the interpolation matrix used to evaluate the basis at the points xi inside the element.

This routine returns a matrix $\mathbf{I}(\mathbf{a})$ that can be used to evaluate the nodal basis at the points defined by the parameter xi, which is denoted by $\mathbf{a} = (a_1, \dots, a_N)$ and $ N $ is the number of points in xi.

In particular, if the array $\mathbf{u}$ with components $\mathbf{u}_i = u^\delta(\xi_i)$ represents the polynomial approximation of a function $ u $ evaluated at the nodal points NodalUtil::m_xi, then the evaluation of $u^\delta$ evaluated at the input points $\mathbf{a}$ is given by $\mathbf{I}(\mathbf{a})\mathbf{u}$ .

Parameters

xi	An array of first size number of spatial dimensions and secondary size the number of points to interpolate.

Returns: The interpolation matrix for the points xi.

Definition at line 202 of file NodalUtil.cpp.

References GetVandermonde(), and v_CreateUtil().

Referenced by Nektar::Utilities::ProcessVarOpti::BuildDerivUtil().

 {
     boost::shared_ptr<NodalUtil> subUtil = v_CreateUtil(xi);
     SharedMatrix matS = GetVandermonde();
     SharedMatrix matT = subUtil->GetVandermonde();
     SharedMatrix D = MemoryManager<NekMatrix<NekDouble> >::AllocateSharedPtr(
         matT->GetRows(), matS->GetColumns(), 0.0);
 
     matS->Invert();
 
     *D = (*matT) * (*matS);
     return D;
 }

SharedMatrix Nektar::LibUtilities::NodalUtil::GetVandermonde ( )

Return the Vandermonde matrix for the nodal distribution.

This routine constructs and returns the Vandermonde matrix $\mathbf{V}$ , with each entry as $\mathbf{V}_{ij} = (\psi_i(\xi_j))$ where $\psi_i$ is the orthogonal basis obtained through the abstract function NodalUtil::v_OrthoBasis.

Returns: The Vandermonde matrix.

Definition at line 102 of file NodalUtil.cpp.

References m_numPoints, v_NumModes(), and v_OrthoBasis().

Referenced by Nektar::Utilities::ProcessVarOpti::BuildDerivUtil(), GetDerivMatrix(), GetInterpolationMatrix(), GetWeights(), and Nektar::NekMeshUtils::BLMesh::Setup().

 {
     int rows = m_numPoints, cols = v_NumModes();
     SharedMatrix matV = MemoryManager<NekMatrix<NekDouble> >::AllocateSharedPtr(
         rows, cols, 0.0);
 
     for (int j = 0; j < cols; ++j)
     {
         NekVector<NekDouble> col = v_OrthoBasis(j);
         for (int i = 0; i < rows; ++i)
         {
             (*matV)(i, j) = col[i];
         }
     }
 
     return matV;
 }

SharedMatrix Nektar::LibUtilities::NodalUtil::GetVandermondeForDeriv ( int dir )

Return the Vandermonde matrix of the derivative of the basis functions for the nodal distribution.

This routine constructs and returns the Vandermonde matrix for the derivative of the basis functions $\mathbf{V}_d$ for coordinate directions $0 \leq d \leq 2$ , with each entry as $\mathbf{V}_{ij} = (\partial_d \psi_i(\xi_j))$ where $\partial_d\psi_i$ is the derivative of the orthogonal basis obtained through the abstract function NodalUtil::v_OrthoBasisDeriv.

Parameters

dir	Direction of derivative in the standard element.

Returns: Vandermonde matrix corresponding with derivative of the basis functions in direction dir.

Definition at line 135 of file NodalUtil.cpp.

References m_numPoints, v_NumModes(), and v_OrthoBasisDeriv().

Referenced by Nektar::Utilities::ProcessVarOpti::BuildDerivUtil(), GetDerivMatrix(), and Nektar::NekMeshUtils::BLMesh::Setup().

 {
     int rows = m_numPoints, cols = v_NumModes();
     SharedMatrix matV = MemoryManager<NekMatrix<NekDouble> >::AllocateSharedPtr(
         rows, cols, 0.0);
 
     for (int j = 0; j < cols; ++j)
     {
         NekVector<NekDouble> col = v_OrthoBasisDeriv(dir, j);
         for (int i = 0; i < rows; ++i)
         {
             (*matV)(i, j) = col[i];
         }
     }
 
     return matV;
 }

NekVector< NekDouble > Nektar::LibUtilities::NodalUtil::GetWeights ( )

Obtain the integration weights for the given nodal distribution.

This routine constructs the integration weights for the given nodal distribution inside NodalUtil::m_xi. To do this we solve the linear system $\mathbf{V}^\top \mathbf{w} = \mathbf{g}$ where $\mathbf{g}_i = \int_\Omega \psi_i(x)\, dx$ , and we use the fact that under the definition of the orthogonal basis for each element type, $\mathbf{g}_i = 0$ for $ i > 0 $ . We use NodalUtil::v_ModeZeroIntegral to return the analytic value of $\mathbf{g}_0$ .

Returns: Vector of integration weights for the integration points.

Definition at line 66 of file NodalUtil.cpp.

References GetVandermonde(), m_numPoints, Nektar::LinearSystem::SolveTranspose(), v_ModeZeroIntegral(), and v_NumModes().

 {
     // Get number of modes in orthogonal basis
     int numModes = v_NumModes();
 
     // If we have the same number of nodes as modes, then we can solve the
     // linear system V^T w = (1,0,...)
     if (numModes == m_numPoints)
     {
         NekVector<NekDouble> g(m_numPoints, 0.0);
         g(0) = v_ModeZeroIntegral();
 
         SharedMatrix V = GetVandermonde();
 
         // Solve the system V^T w = g to obtain weights.
         LinearSystem matL(V);
         return matL.SolveTranspose(g);
     }
     else
     {
         // System is either over- or under-determined. Need to do least squares
         // here using SVD.
         return NekVector<NekDouble>();
     }
 }

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

protectedpure virtual

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

xi	Distribution of nodal points to create utility with.

Implemented in Nektar::Utilities::NodalUtilTriMonomial, Nektar::LibUtilities::NodalUtilHex, Nektar::LibUtilities::NodalUtilQuad, Nektar::LibUtilities::NodalUtilPrism, Nektar::LibUtilities::NodalUtilTetrahedron, and Nektar::LibUtilities::NodalUtilTriangle.

Referenced by GetInterpolationMatrix().

virtual NekDouble Nektar::LibUtilities::NodalUtil::v_ModeZeroIntegral ( )

protectedpure virtual

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.

Implemented in Nektar::LibUtilities::NodalUtilHex, Nektar::LibUtilities::NodalUtilQuad, Nektar::LibUtilities::NodalUtilPrism, Nektar::LibUtilities::NodalUtilTetrahedron, and Nektar::LibUtilities::NodalUtilTriangle.

Referenced by GetWeights().

virtual int Nektar::LibUtilities::NodalUtil::v_NumModes ( )

protectedpure virtual

Calculate the number of degrees of freedom for this element.

Implemented in Nektar::LibUtilities::NodalUtilHex, Nektar::LibUtilities::NodalUtilQuad, Nektar::LibUtilities::NodalUtilPrism, Nektar::LibUtilities::NodalUtilTetrahedron, and Nektar::LibUtilities::NodalUtilTriangle.

Referenced by GetVandermonde(), GetVandermondeForDeriv(), and GetWeights().

virtual NekVector<NekDouble> Nektar::LibUtilities::NodalUtil::v_OrthoBasis ( const int mode )