#include <NodalTriElec.h>

Inheritance diagram for Nektar::LibUtilities::NodalTriElec:

Public Member Functions
virtual	~NodalTriElec ()

	NodalTriElec (const PointsKey &key)

Public Member Functions inherited from Nektar::LibUtilities::Points< NekDouble >
virtual	~Points ()

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 MatrixSharedPtrType &	GetD (Direction dir=xDir) const

const MatrixSharedPtrType	GetI (const PointsKey &key)

const MatrixSharedPtrType	GetI (const Array< OneD, const DataType > &x)

const MatrixSharedPtrType	GetI (unsigned int uint, const Array< OneD, const DataType > &x)

const MatrixSharedPtrType	GetI (const Array< OneD, const DataType > &x, const Array< OneD, const DataType > &y)

const MatrixSharedPtrType	GetI (const Array< OneD, const DataType > &x, const Array< OneD, const DataType > &y, const Array< OneD, const DataType > &z)

const MatrixSharedPtrType	GetGalerkinProjection (const PointsKey &pkey)

Static Public Member Functions
static std::shared_ptr< PointsBaseType >	Create (const PointsKey &key)

Protected Member Functions
virtual const MatrixSharedPtrType	v_GetI (const PointsKey &pkey) override

virtual const MatrixSharedPtrType	v_GetI (const Array< OneD, const NekDouble > &x, const Array< OneD, const NekDouble > &y) override

Protected Member Functions inherited from Nektar::LibUtilities::Points< NekDouble >
virtual void	v_Initialize (void)

virtual void	v_CalculateBaryWeights ()
	This function calculates the barycentric weights used for enhanced interpolation speed. More...

	Points (const PointsKey &key)

virtual const MatrixSharedPtrType	v_GetI (const Array< OneD, const DataType > &x)

virtual const MatrixSharedPtrType	v_GetI (unsigned int, const Array< OneD, const DataType > &x)

virtual const MatrixSharedPtrType	v_GetI (const Array< OneD, const DataType > &x, const Array< OneD, const DataType > &y, const Array< OneD, const DataType > &z)

virtual const MatrixSharedPtrType	v_GetGalerkinProjection (const PointsKey &pkey)

Private Member Functions
	NodalTriElec ()

void	NodalPointReorder2d ()

virtual void	v_CalculatePoints () override

virtual void	v_CalculateWeights () override

virtual void	v_CalculateDerivMatrix () override

void	CalculateInterpMatrix (const Array< OneD, const NekDouble > &xia, const Array< OneD, const NekDouble > &yia, Array< OneD, NekDouble > &interp)

Private Attributes
std::shared_ptr< NodalUtilTriangle >	m_util

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 Attributes inherited from Nektar::LibUtilities::Points< NekDouble >
PointsKey	m_pointsKey
	Points type for this points distributions. More...

Array< OneD, DataType >	m_points [3]
	Storage for the point locations, allowing for up to a 3D points storage. More...

Array< OneD, DataType >	m_weights
	Quadrature weights for the weights. More...

Array< OneD, DataType >	m_bcweights
	Barycentric weights. More...

MatrixSharedPtrType	m_derivmatrix [3]
	Derivative matrices. More...

NekManager< PointsKey, NekMatrix< DataType >, PointsKey::opLess >	m_InterpManager

NekManager< PointsKey, NekMatrix< DataType >, PointsKey::opLess >	m_GalerkinProjectionManager

Detailed Description

Definition at line 52 of file NodalTriElec.h.

Constructor & Destructor Documentation

◆ ~NodalTriElec()

virtual Nektar::LibUtilities::NodalTriElec::~NodalTriElec ( )

inlinevirtual

Definition at line 55 of file NodalTriElec.h.

56 {

57 }

◆ NodalTriElec() [1/2]

Nektar::LibUtilities::NodalTriElec::NodalTriElec ( const PointsKey & key )

inline

Definition at line 62 of file NodalTriElec.h.

                                        : PointsBaseType(key)
     {
     }

◆ NodalTriElec() [2/2]

Nektar::LibUtilities::NodalTriElec::NodalTriElec ( )

inlineprivate

Definition at line 96 of file NodalTriElec.h.

                    : PointsBaseType(NullPointsKey)
     {
     }

Member Function Documentation

◆ CalculateInterpMatrix()

void Nektar::LibUtilities::NodalTriElec::CalculateInterpMatrix	(	const Array< OneD, const NekDouble > &	xia,
		const Array< OneD, const NekDouble > &	yia,
		Array< OneD, NekDouble > &	interp
	)

private

Definition at line 160 of file NodalTriElec.cpp.

 {
     Array<OneD, Array<OneD, NekDouble>> xi(2);
     xi[0] = xia;
     xi[1] = yia;
  
     std::shared_ptr<NekMatrix<NekDouble>> mat =
         m_util->GetInterpolationMatrix(xi);
     Vmath::Vcopy(mat->GetRows() * mat->GetColumns(), mat->GetRawPtr(), 1,
                  &interp[0], 1);
 }

References m_util, and Vmath::Vcopy().

Referenced by v_GetI().

◆ Create()

std::shared_ptr< PointsBaseType > Nektar::LibUtilities::NodalTriElec::Create ( const PointsKey & key )

static

Definition at line 174 of file NodalTriElec.cpp.

 {
     std::shared_ptr<PointsBaseType> returnval(
         MemoryManager<NodalTriElec>::AllocateSharedPtr(key));
     returnval->Initialize();
     return returnval;
 }

◆ NodalPointReorder2d()

void Nektar::LibUtilities::NodalTriElec::NodalPointReorder2d ( )

private

Definition at line 182 of file NodalTriElec.cpp.

 {
     int i, j;
     int cnt;
     int istart, iend;
  
     const int nVerts              = 3;
     const int nEdgeInteriorPoints = GetNumPoints() - 2;
     const int nBoundaryPoints     = 3 * nEdgeInteriorPoints + 3;
  
     if (nEdgeInteriorPoints == 0)
     {
         return;
     }
  
     // group the points of edge 1 together;
     istart = nVerts;
     for (i = cnt = istart; i < nBoundaryPoints; i++)
     {
         if (fabs(m_points[1][i] + 1.0) < NekConstants::kNekZeroTol)
         {
             std::swap(m_points[0][cnt], m_points[0][i]);
             std::swap(m_points[1][cnt], m_points[1][i]);
             cnt++;
         }
     }
  
     // bubble sort edge 1 (counterclockwise numbering)
     iend = istart + nEdgeInteriorPoints;
     for (i = istart; i < iend; i++)
     {
         for (j = istart + 1; j < iend; j++)
         {
             if (m_points[0][j] < m_points[0][j - 1])
             {
                 std::swap(m_points[0][j], m_points[0][j - 1]);
                 std::swap(m_points[1][j], m_points[1][j - 1]);
             }
         }
     }
  
     // group the points of edge 2 together;
     istart = iend;
     for (i = cnt = istart; i < nBoundaryPoints; i++)
     {
         if (fabs(m_points[1][i] + m_points[0][i]) < NekConstants::kNekZeroTol)
         {
             std::swap(m_points[0][cnt], m_points[0][i]);
             std::swap(m_points[1][cnt], m_points[1][i]);
             cnt++;
         }
     }
  
     // bubble sort edge 2 (counterclockwise numbering)
     iend = istart + nEdgeInteriorPoints;
     for (i = istart; i < iend; i++)
     {
         for (j = istart + 1; j < iend; j++)
         {
             if (m_points[1][j] < m_points[1][j - 1])
             {
                 std::swap(m_points[0][j], m_points[0][j - 1]);
                 std::swap(m_points[1][j], m_points[1][j - 1]);
             }
         }
     }
  
     // group the points of edge 3 together;
     istart = iend;
     for (i = cnt = istart; i < nBoundaryPoints; i++)
     {
         if (fabs(m_points[0][i] + 1.0) < NekConstants::kNekZeroTol)
         {
             std::swap(m_points[0][cnt], m_points[0][i]);
             std::swap(m_points[1][cnt], m_points[1][i]);
             cnt++;
         }
     }
     // bubble sort edge 3 (counterclockwise numbering)
     iend = istart + nEdgeInteriorPoints;
     for (i = istart; i < iend; i++)
     {
         for (j = istart + 1; j < iend; j++)
         {
             if (m_points[1][j] > m_points[1][j - 1])
             {
                 std::swap(m_points[0][j], m_points[0][j - 1]);
                 std::swap(m_points[1][j], m_points[1][j - 1]);
             }
         }
     }
  
     if (GetNumPoints() < 5)
     {
         // at numpoints = 4 there is only one interior point so doesnt
         // need sorting
         return;
     }
  
     // someone forgot to finish this piece of code and tell anyone
     // that they didnt
     // face interior nodes needs to be considered
     // make a copy of the unsorted nodes
     // bubble sort by smallest y
     // which will put them into sets of ever decreasing size
     // which can be bubble sorted by x to obtain the distrobution
  
     Array<OneD, NekDouble> xc(m_points[0].size() - iend);
     Array<OneD, NekDouble> yc(m_points[0].size() - iend);
     int ct = 0;
     for (i = iend; i < m_points[0].size(); i++, ct++)
     {
         xc[ct] = m_points[0][i];
         yc[ct] = m_points[1][i];
     }
  
     // sort smallest first
     bool repeat = true;
     while (repeat)
     {
         repeat = false;
         for (i = 0; i < xc.size() - 1; i++)
         {
             if (yc[i] > yc[i + 1])
             {
                 std::swap(xc[i], xc[i + 1]);
                 std::swap(yc[i], yc[i + 1]);
                 repeat = true;
             }
         }
     }
  
     int offset = 0;
     int npl    = GetNumPoints() - 3;
     while (npl > 1)
     {
         repeat = true;
         while (repeat)
         {
             repeat = false;
             for (i = offset; i < offset + npl - 1; i++)
             {
                 if (xc[i] > xc[i + 1])
                 {
                     std::swap(xc[i], xc[i + 1]);
                     std::swap(yc[i], yc[i + 1]);
                     repeat = true;
                 }
             }
         }
         offset += npl;
         npl--;
     }
  
     // copy back in
     ct = 0;
     for (i = iend; i < m_points[0].size(); i++, ct++)
     {
         m_points[0][i] = xc[ct];
         m_points[1][i] = yc[ct];
     }
     return;
 }

References Nektar::LibUtilities::Points< NekDouble >::GetNumPoints(), Nektar::NekConstants::kNekZeroTol, Nektar::LibUtilities::Points< NekDouble >::m_points, and Nektar::Array< OneD, const DataType >::size().

Referenced by v_CalculatePoints().

◆ v_CalculateDerivMatrix()

void Nektar::LibUtilities::NodalTriElec::v_CalculateDerivMatrix ( )

overrideprivatevirtual

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

Definition at line 149 of file NodalTriElec.cpp.

 {
     // Allocate the derivative matrix.
     PointsBaseType::v_CalculateDerivMatrix();
  
     m_derivmatrix[0] = m_util->GetDerivMatrix(0);
     m_derivmatrix[1] = m_util->GetDerivMatrix(1);
 }

References Nektar::LibUtilities::Points< NekDouble >::m_derivmatrix, m_util, and Nektar::LibUtilities::Points< NekDouble >::v_CalculateDerivMatrix().

◆ v_CalculatePoints()

void Nektar::LibUtilities::NodalTriElec::v_CalculatePoints ( )

overrideprivatevirtual

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

Definition at line 49 of file NodalTriElec.cpp.

 {
     // Allocate the storage for points
     PointsBaseType::v_CalculatePoints();
  
     int index = 0, isum = 0;
     const int offset = 3; // offset to match Datafile
     NekDouble b, c;
     unsigned int numPoints = GetNumPoints();
  
     // initialize values
     for (unsigned int i = 0; i < numPoints - 2; ++i)
     {
         index += NodalTriElecNPTS[i];
     }
  
     for (unsigned int i = 0; i < NodalTriElecNPTS[numPoints - 2]; ++i, ++index)
     {
         if (int(NodalTriElecData[index][0]))
         {
             b = NodalTriElecData[index][4];
             c = NodalTriElecData[index][5];
  
             m_points[0][isum] = 2.0 * b - 1.0;
             m_points[1][isum] = 2.0 * c - 1.0;
             isum++;
             continue;
         } // end symmetry1
  
         if (int(NodalTriElecData[index][1]) == 1)
         {
             for (unsigned int j = 0; j < 3; ++j)
             {
                 b = NodalTriElecData[index][offset + perm3A_2d[j][1]];
                 c = NodalTriElecData[index][offset + perm3A_2d[j][2]];
                 m_points[0][isum] = 2.0 * b - 1.0;
                 m_points[1][isum] = 2.0 * c - 1.0;
                 isum++;
             } // end j
             continue;
         } // end symmetry3a
  
         if (int(NodalTriElecData[index][1]) == 2)
         {
             for (unsigned int j = 0; j < 3; ++j)
             {
                 b = NodalTriElecData[index][offset + perm3B_2d[j][1]];
                 c = NodalTriElecData[index][offset + perm3B_2d[j][2]];
                 m_points[0][isum] = 2.0 * b - 1.0;
                 m_points[1][isum] = 2.0 * c - 1.0;
                 isum++;
             } // end j
             continue;
         } // end symmetry3b
  
         if (int(NodalTriElecData[index][2]))
         {
             for (unsigned int j = 0; j < 6; ++j)
             {
                 b = NodalTriElecData[index][offset + perm6_2d[j][1]];
                 c = NodalTriElecData[index][offset + perm6_2d[j][2]];
                 m_points[0][isum] = 2.0 * b - 1.0;
                 m_points[1][isum] = 2.0 * c - 1.0;
                 isum++;
             } // end j
             continue;
         } // end symmetry6
     }     // end npts
  
     //    std::cout << "(x y) = (" << ToVector(m_points[0]) << ", " <<
     //    ToVector(m_points[1]) <<  ")" << std::endl; cout << "numPoints = " <<
     //    numPoints << endl; cout << "NodalTriElecNPTS[numPoints-2] = " <<
     //    NodalTriElecNPTS[numPoints-2] << endl; cout << "isum = " << isum <<
     //    endl; for( int i = 0; i <= numPoints-2; ++i ) {
     //        cout << "NodalTriElecNPTS[" << i << "] = " << NodalTriElecNPTS[i]
     //        << endl;
     //    }
     NodalPointReorder2d();
  
     ASSERTL1((static_cast<unsigned int>(isum) == m_pointsKey.GetTotNumPoints()),
              "sum not equal to npts");
  
     m_util = MemoryManager<NodalUtilTriangle>::AllocateSharedPtr(
         numPoints - 1, m_points[0], m_points[1]);
  
     // exit(0);
 }

◆ v_CalculateWeights()

void Nektar::LibUtilities::NodalTriElec::v_CalculateWeights ( )

overrideprivatevirtual

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

Definition at line 137 of file NodalTriElec.cpp.

 {
     // Allocate the storage for points
     PointsBaseType::v_CalculateWeights();
  
     typedef DataType T;
  
     // Solve the Vandermonde system of integrals for the weight vector
     NekVector<T> w = m_util->GetWeights();
     m_weights      = Array<OneD, T>(w.GetRows(), w.GetPtr());
 }

References Nektar::NekVector< DataType >::GetPtr(), Nektar::NekVector< DataType >::GetRows(), m_util, Nektar::LibUtilities::Points< NekDouble >::m_weights, and Nektar::LibUtilities::Points< NekDouble >::v_CalculateWeights().

◆ v_GetI() [1/2]

virtual const MatrixSharedPtrType Nektar::LibUtilities::NodalTriElec::v_GetI	(	const Array< OneD, const NekDouble > &	x,
		const Array< OneD, const NekDouble > &	y
	)

inlineoverrideprotectedvirtual

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

Definition at line 77 of file NodalTriElec.h.

     {
         size_t numpoints = x.size();
         Array<OneD, NekDouble> interp(GetTotNumPoints() * numpoints);
         CalculateInterpMatrix(x, y, interp);
  
         unsigned int np = GetTotNumPoints();
         NekDouble *d    = interp.data();
         return MemoryManager<NekMatrix<NekDouble>>::AllocateSharedPtr(numpoints,
                                                                       np, d);
     }

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

◆ v_GetI() [2/2]

virtual const MatrixSharedPtrType Nektar::LibUtilities::NodalTriElec::v_GetI ( const PointsKey & pkey )

inlineoverrideprotectedvirtual

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

Definition at line 67 of file NodalTriElec.h.

     {
         ASSERTL0(pkey.GetPointsDim() == 2,
                  "NodalTriElec Points can only interpolate to other 2d "
                  "point distributions");
         Array<OneD, const NekDouble> x, y;
         PointsManager()[pkey]->GetPoints(x, y);
         return GetI(x, y);
     }