Nektar++
Loading...
Searching...
No Matches
Public Member Functions | Protected Member Functions | Private Member Functions | List of all members
Nektar::StdRegions::StdPrismExp Class Reference

Class representing a prismatic element in reference space. More...

#include <StdPrismExp.h>

Inheritance diagram for Nektar::StdRegions::StdPrismExp:
[legend]

Public Member Functions

 StdPrismExp (const LibUtilities::BasisKey &Ba, const LibUtilities::BasisKey &Bb, const LibUtilities::BasisKey &Bc)
 
 StdPrismExp (const LibUtilities::BasisKey &Ba, const LibUtilities::BasisKey &Bb, const LibUtilities::BasisKey &Bc, NekDouble *coeffs, NekDouble *phys)
 
 StdPrismExp ()=default
 
 StdPrismExp (const StdPrismExp &T)=default
 
 ~StdPrismExp () override=default
 
- Public Member Functions inherited from Nektar::StdRegions::StdExpansion3D
 StdExpansion3D (int numcoeffs, const LibUtilities::BasisKey &Ba, const LibUtilities::BasisKey &Bb, const LibUtilities::BasisKey &Bc)
 
 StdExpansion3D ()=default
 
 StdExpansion3D (const StdExpansion3D &T)=default
 
 ~StdExpansion3D () override=default
 
void IProductWRTBaseKernel (const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const Array< OneD, NekDouble > &jac, const bool Deformed, bool CollDir0=false, bool CollDir1=false, bool CollDir2=false)
 
int GetNedges () const
 return the number of edges in 3D expansion
 
int GetEdgeNcoeffs (const int i) const
 This function returns the number of expansion coefficients belonging to the i-th edge.
 
void GetEdgeInteriorToElementMap (const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, Orientation traceOrient=eForwards)
 
- Public Member Functions inherited from Nektar::StdRegions::StdExpansion
 StdExpansion ()
 Default Constructor.
 
 StdExpansion (const int numcoeffs, const int numbases, const LibUtilities::BasisKey &Ba=LibUtilities::NullBasisKey, const LibUtilities::BasisKey &Bb=LibUtilities::NullBasisKey, const LibUtilities::BasisKey &Bc=LibUtilities::NullBasisKey)
 Constructor.
 
 StdExpansion (const StdExpansion &T)
 Copy Constructor.
 
virtual ~StdExpansion ()
 Destructor.
 
int GetNumBases () const
 This function returns the number of 1D bases used in the expansion.
 
const Array< OneD, const LibUtilities::BasisSharedPtr > & GetBase () const
 This function gets the shared point to basis.
 
const LibUtilities::BasisSharedPtrGetBasis (int dir) const
 This function gets the shared point to basis in the dir direction.
 
int GetNcoeffs (void) const
 This function returns the total number of coefficients used in the expansion.
 
int GetTotPoints () const
 This function returns the total number of quadrature points used in the element.
 
LibUtilities::BasisType GetBasisType (const int dir) const
 This function returns the type of basis used in the dir direction.
 
int GetBasisNumModes (const int dir) const
 This function returns the number of expansion modes in the dir direction.
 
int EvalBasisNumModesMax (void) const
 This function returns the maximum number of expansion modes over all local directions.
 
LibUtilities::PointsType GetPointsType (const int dir) const
 This function returns the type of quadrature points used in the dir direction.
 
int GetNumPoints (const int dir) const
 This function returns the number of quadrature points in the dir direction.
 
const Array< OneD, const NekDouble > & GetPoints (const int dir) const
 This function returns a pointer to the array containing the quadrature points in dir direction.
 
int GetNverts () const
 This function returns the number of vertices of the expansion domain.
 
int GetTraceNcoeffs (const int i) const
 This function returns the number of expansion coefficients belonging to the i-th trace.
 
int GetTraceIntNcoeffs (const int i) const
 
int GetTraceNumPoints (const int i) const
 This function returns the number of quadrature points belonging to the i-th trace.
 
const LibUtilities::BasisKey GetTraceBasisKey (const int i, int k=-1, bool UseGLL=false) const
 This function returns the basis key belonging to the i-th trace.
 
LibUtilities::PointsKey GetTracePointsKey (const int i, int k=-1) const
 This function returns the basis key belonging to the i-th trace.
 
int NumBndryCoeffs (void) const
 
int NumDGBndryCoeffs (void) const
 
const LibUtilities::PointsKey GetNodalPointsKey () const
 This function returns the type of expansion Nodal point type if defined.
 
int GetNtraces () const
 Returns the number of trace elements connected to this element.
 
LibUtilities::ShapeType DetShapeType () const
 This function returns the shape of the expansion domain.
 
int GetShapeDimension () const
 
bool IsBoundaryInteriorExpansion () const
 
bool IsNodalNonTensorialExp ()
 
void NodalToModal (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
void BwdTrans (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 This function performs the Backward transformation from coefficient space to physical space.
 
void FwdTrans (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
void FwdTransBndConstrained (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
NekDouble Integral (const Array< OneD, const NekDouble > &inarray)
 This function integrates the specified function over the domain.
 
void FillMode (const int mode, Array< OneD, NekDouble > &outarray)
 This function fills the array outarray with the mode-th mode of the expansion.
 
void IProductWRTBase (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 this function calculates the inner product of a given function f with the different modes of the expansion
 
void IProductWRTDerivBase (const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
void IProductWRTDirectionalDerivBase (const Array< OneD, const NekDouble > &direction, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
int GetElmtId ()
 Get the element id of this expansion when used in a list by returning value of m_elmt_id.
 
void SetElmtId (const int id)
 Set the element id of this expansion when used in a list by returning value of m_elmt_id.
 
void GetCoords (Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2=NullNekDouble1DArray, Array< OneD, NekDouble > &coords_3=NullNekDouble1DArray)
 this function returns the physical coordinates of the quadrature points of the expansion
 
Array< OneD, Array< OneD, NekDouble > > GetCoords ()
 
void GetCoord (const Array< OneD, const NekDouble > &Lcoord, Array< OneD, NekDouble > &coord)
 given the coordinates of a point of the element in the local collapsed coordinate system, this function calculates the physical coordinates of the point
 
DNekMatSharedPtr GetStdMatrix (const StdMatrixKey &mkey)
 
DNekBlkMatSharedPtr GetStdStaticCondMatrix (const StdMatrixKey &mkey)
 
Array< OneD, const NekDoubleGetStdFac (const StdFacKey &mkey)
 
void NormVectorIProductWRTBase (const Array< OneD, const NekDouble > &Fx, Array< OneD, NekDouble > &outarray)
 
void NormVectorIProductWRTBase (const Array< OneD, const NekDouble > &Fx, const Array< OneD, NekDouble > &Fy, Array< OneD, NekDouble > &outarray)
 
void NormVectorIProductWRTBase (const Array< OneD, const NekDouble > &Fx, const Array< OneD, const NekDouble > &Fy, const Array< OneD, const NekDouble > &Fz, Array< OneD, NekDouble > &outarray)
 
void NormVectorIProductWRTBase (const Array< OneD, const Array< OneD, NekDouble > > &Fvec, Array< OneD, NekDouble > &outarray)
 
DNekScalBlkMatSharedPtr GetLocStaticCondMatrix (const LocalRegions::MatrixKey &mkey)
 
void DropLocStaticCondMatrix (const LocalRegions::MatrixKey &mkey)
 
int CalcNumberOfCoefficients (const std::vector< unsigned int > &nummodes, int &modes_offset)
 
NekDouble StdPhysEvaluate (const Array< OneD, const NekDouble > &Lcoord, const Array< OneD, const NekDouble > &physvals)
 
int GetCoordim ()
 
void GetBoundaryMap (Array< OneD, unsigned int > &outarray)
 
void GetInteriorMap (Array< OneD, unsigned int > &outarray)
 
int GetVertexMap (const int localVertexId, bool useCoeffPacking=false)
 
void GetTraceToElementMap (const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, Orientation traceOrient=eForwards, int P=-1, int Q=-1)
 
void GetTraceCoeffMap (const unsigned int traceid, Array< OneD, unsigned int > &maparray)
 
void GetElmtTraceToTraceMap (const unsigned int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, Orientation traceOrient=eForwards, int P=-1, int Q=-1)
 
void GetTraceInteriorToElementMap (const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, const Orientation traceOrient=eForwards)
 
void GetTraceNumModes (const int tid, int &numModes0, int &numModes1, const Orientation traceOrient=eDir1FwdDir1_Dir2FwdDir2)
 
void MultiplyByQuadratureMetric (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
void MultiplyByStdQuadratureMetric (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
DNekMatSharedPtr CreateGeneralMatrix (const StdMatrixKey &mkey)
 this function generates the mass matrix \(\mathbf{M}[i][j] = \int \phi_i(\mathbf{x}) \phi_j(\mathbf{x}) d\mathbf{x}\)
 
void GeneralMatrixOp (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
 
void MassMatrixOp (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
 
void LaplacianMatrixOp (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
 
void ReduceOrderCoeffs (int numMin, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
void SVVLaplacianFilter (Array< OneD, NekDouble > &array, const StdMatrixKey &mkey)
 
void ExponentialFilter (Array< OneD, NekDouble > &array, const NekDouble alpha, const NekDouble exponent, const NekDouble cutoff)
 
void LaplacianMatrixOp (const int k1, const int k2, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
 
void WeakDerivMatrixOp (const int i, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
 
void WeakDirectionalDerivMatrixOp (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
 
void MassLevelCurvatureMatrixOp (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
 
void LinearAdvectionMatrixOp (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
 
void LinearAdvectionDiffusionReactionMatrixOp (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey, bool addDiffusionTerm=true)
 
void HelmholtzMatrixOp (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
 
DNekMatSharedPtr GenMatrix (const StdMatrixKey &mkey)
 
void PhysDeriv (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1=NullNekDouble1DArray, Array< OneD, NekDouble > &out_d2=NullNekDouble1DArray)
 
void PhysDeriv (const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
void PhysDirectionalDeriv (const Array< OneD, const NekDouble > &inarray, const Array< OneD, const NekDouble > &direction, Array< OneD, NekDouble > &outarray)
 
void StdPhysDeriv (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1=NullNekDouble1DArray, Array< OneD, NekDouble > &out_d2=NullNekDouble1DArray)
 
NekDouble PhysEvaluate (const Array< OneD, const NekDouble > &coords, const Array< OneD, const NekDouble > &physvals)
 This function evaluates the expansion at a single (arbitrary) point of the domain.
 
NekDouble PhysEvaluate (const Array< OneD, NekDouble > &coord, const Array< OneD, const NekDouble > &inarray, std::array< NekDouble, 3 > &firstOrderDerivs)
 This function evaluates the first derivative of the expansion at a single (arbitrary) point of the domain.
 
NekDouble PhysEvaluate (const Array< OneD, NekDouble > &coord, const Array< OneD, const NekDouble > &inarray, std::array< NekDouble, 3 > &firstOrderDerivs, std::array< NekDouble, 6 > &secondOrderDerivs)
 
NekDouble PhysEvaluate (const Array< OneD, DNekMatSharedPtr > &I, const Array< OneD, const NekDouble > &physvals)
 This function evaluates the expansion at a single (arbitrary) point of the domain.
 
NekDouble PhysEvaluateBasis (const Array< OneD, const NekDouble > &coords, int mode)
 This function evaluates the basis function mode mode at a point coords of the domain.
 
void ReOrientTracePhysMap (const StdRegions::Orientation orient, Array< OneD, int > &idmap, const int nq0, const int nq1, bool Forwards=true)
 
void LocCoordToLocCollapsed (const Array< OneD, const NekDouble > &xi, Array< OneD, NekDouble > &eta)
 Convert local cartesian coordinate xi into local collapsed coordinates eta.
 
void LocCollapsedToLocCoord (const Array< OneD, const NekDouble > &eta, Array< OneD, NekDouble > &xi)
 Convert local collapsed coordinates eta into local cartesian coordinate xi.
 
void PhysInterp (std::shared_ptr< StdExpansion > fromExp, const Array< OneD, const NekDouble > &fromData, Array< OneD, NekDouble > &toData, bool Transpose=false)
 interpolate from one set of quadrature points available from FromExp to the set of quadrature points in the current expansion. If the points are the same this routine will just copy the data
 
virtual void v_NormVectorIProductWRTBase (const Array< OneD, const NekDouble > &Fx, Array< OneD, NekDouble > &outarray)
 
virtual void v_NormVectorIProductWRTBase (const Array< OneD, const NekDouble > &Fx, const Array< OneD, const NekDouble > &Fy, Array< OneD, NekDouble > &outarray)
 
virtual void v_NormVectorIProductWRTBase (const Array< OneD, const NekDouble > &Fx, const Array< OneD, const NekDouble > &Fy, const Array< OneD, const NekDouble > &Fz, Array< OneD, NekDouble > &outarray)
 
virtual void v_NormVectorIProductWRTBase (const Array< OneD, const Array< OneD, NekDouble > > &Fvec, Array< OneD, NekDouble > &outarray)
 
virtual DNekScalBlkMatSharedPtr v_GetLocStaticCondMatrix (const LocalRegions::MatrixKey &mkey)
 
virtual void v_DropLocStaticCondMatrix (const LocalRegions::MatrixKey &mkey)
 
NekDouble Linf (const Array< OneD, const NekDouble > &phys, const Array< OneD, const NekDouble > &sol=NullNekDouble1DArray)
 Function to evaluate the discrete \( L_\infty\) error \( |\epsilon|_\infty = \max |u - u_{exact}|\) where \( u_{exact}\) is given by the array sol.
 
NekDouble L2 (const Array< OneD, const NekDouble > &phys, const Array< OneD, const NekDouble > &sol=NullNekDouble1DArray)
 Function to evaluate the discrete \( L_2\) error, \( | \epsilon |_{2} = \left [ \int^1_{-1} [u - u_{exact}]^2 dx \right]^{1/2} d\xi_1 \) where \( u_{exact}\) is given by the array sol.
 
NekDouble H1 (const Array< OneD, const NekDouble > &phys, const Array< OneD, const NekDouble > &sol=NullNekDouble1DArray)
 Function to evaluate the discrete \( H^1\) error, \( | \epsilon |^1_{2} = \left [ \int^1_{-1} [u - u_{exact}]^2 + \nabla(u - u_{exact})\cdot\nabla(u - u_{exact})\cdot dx \right]^{1/2} d\xi_1 \) where \( u_{exact}\) is given by the array sol.
 
const LibUtilities::PointsKeyVector GetPointsKeys () const
 
DNekMatSharedPtr BuildInverseTransformationMatrix (const DNekScalMatSharedPtr &m_transformationmatrix)
 
void PhysInterpToSimplexEquiSpaced (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, int npset=-1)
 This function performs an interpolation from the physical space points provided at input into an array of equispaced points which are not the collapsed coordinate. So for a tetrahedron you will only get a tetrahedral number of values.
 
void PhysInterpToGLL (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, int npset=-1)
 
void PhysInterpToPoints (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, int npset, MatrixType distrib)
 
void GetSimplexEquiSpacedConnectivity (Array< OneD, int > &conn, bool standard=true)
 This function provides the connectivity of local simplices (triangles or tets) to connect the equispaced data points provided by PhysInterpToSimplexEquiSpaced.
 
void EquiSpacedToCoeffs (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 This function performs a projection/interpolation from the equispaced points sometimes used in post-processing onto the coefficient space.
 
void EquiSpacedToPhys (const int nequi, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
template<class T >
std::shared_ptr< T > as ()
 
void GenStdMatBwdDeriv (const int dir, DNekMatSharedPtr &mat)
 

Protected Member Functions

void v_StdPhysDeriv (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2) override
 Calculate the derivative of the physical points.
 
void v_BwdTrans (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
 
void v_IProductWRTBaseKernel (const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const Array< OneD, NekDouble > &jac, const bool Deformed, bool CollDir0=false, bool CollDir1=false, bool CollDir2=false) override
 Inner product of inarray over region with respect to the expansion basis (this)->m_base[0] and return in outarray.
 
void v_IProductWRTDerivBase (const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
 Inner product of inarray over region with respect to the object's default expansion basis; output in outarray.
 
void v_LocCoordToLocCollapsed (const Array< OneD, const NekDouble > &xi, Array< OneD, NekDouble > &eta) override
 
void v_LocCollapsedToLocCoord (const Array< OneD, const NekDouble > &eta, Array< OneD, NekDouble > &xi) override
 
void v_GetCoords (Array< OneD, NekDouble > &xi_x, Array< OneD, NekDouble > &xi_y, Array< OneD, NekDouble > &xi_z) override
 
void v_FillMode (const int mode, Array< OneD, NekDouble > &outarray) override
 
NekDouble v_PhysEvaluateBasis (const Array< OneD, const NekDouble > &coords, int mode) final
 
void v_GetTraceNumModes (const int fid, int &numModes0, int &numModes1, Orientation faceOrient=eDir1FwdDir1_Dir2FwdDir2) override
 
NekDouble v_PhysEvalFirstDeriv (const Array< OneD, NekDouble > &coord, const Array< OneD, const NekDouble > &inarray, std::array< NekDouble, 3 > &firstOrderDerivs) override
 
int v_GetNverts () const final
 
int v_GetNedges () const final
 
int v_GetNtraces () const final
 
LibUtilities::ShapeType v_DetShapeType () const override
 Return Shape of region, using ShapeType enum list; i.e. prism.
 
int v_NumBndryCoeffs () const override
 
int v_NumDGBndryCoeffs () const override
 
int v_GetTraceNcoeffs (const int i) const override
 
int v_GetTraceIntNcoeffs (const int i) const override
 
int v_GetTraceNumPoints (const int i) const override
 
int v_GetEdgeNcoeffs (const int i) const override
 
const LibUtilities::BasisKey v_GetTraceBasisKey (const int i, const int k, bool UseGLL=false) const override
 
LibUtilities::PointsKey v_GetTracePointsKey (const int i, const int j) const override
 
int v_CalcNumberOfCoefficients (const std::vector< unsigned int > &nummodes, int &modes_offset) override
 
bool v_IsBoundaryInteriorExpansion () const override
 
int v_GetVertexMap (int localVertexId, bool useCoeffPacking=false) override
 
void v_GetInteriorMap (Array< OneD, unsigned int > &outarray) override
 
void v_GetBoundaryMap (Array< OneD, unsigned int > &outarray) override
 
void v_GetTraceCoeffMap (const unsigned int fid, Array< OneD, unsigned int > &maparray) override
 
void v_GetElmtTraceToTraceMap (const unsigned int fid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, Orientation faceOrient, int P, int Q) override
 
void v_GetEdgeInteriorToElementMap (const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, const Orientation traceOrient=eDir1FwdDir1_Dir2FwdDir2) override
 
void v_GetTraceInteriorToElementMap (const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, const Orientation traceOrient=eDir1FwdDir1_Dir2FwdDir2) override
 
DNekMatSharedPtr v_GenMatrix (const StdMatrixKey &mkey) override
 
DNekMatSharedPtr v_CreateStdMatrix (const StdMatrixKey &mkey) override
 
void v_SVVLaplacianFilter (Array< OneD, NekDouble > &array, const StdMatrixKey &mkey) override
 
void v_ReduceOrderCoeffs (int numMin, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
 
- Protected Member Functions inherited from Nektar::StdRegions::StdExpansion3D
void PhysTensorDeriv (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2)
 Calculate the 3D derivative in the local tensor/collapsed coordinate at the physical points.
 
void v_PhysDeriv (const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
 Calculate the derivative of the physical points in a given direction.
 
NekDouble v_StdPhysEvaluate (const Array< OneD, const NekDouble > &coords, const Array< OneD, const NekDouble > &physvals) override
 This function evaluates the expansion at a single (arbitrary) point of the domain.
 
NekDouble v_PhysEvaluateInterp (const Array< OneD, DNekMatSharedPtr > &I, const Array< OneD, const NekDouble > &physvals) override
 
void v_IProductWRTBase (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
 
void v_MultiplyByStdQuadratureMetric (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
 
void v_LaplacianMatrixOp_MatFree (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
 
void v_HelmholtzMatrixOp_MatFree (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
 
NekDouble BaryTensorDeriv (const Array< OneD, NekDouble > &coord, const Array< OneD, const NekDouble > &inarray, std::array< NekDouble, 3 > &firstOrderDerivs)
 
void v_GetTraceToElementMap (const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, Orientation traceOrient, int P, int Q) override
 
void v_GenStdMatBwdDeriv (const int dir, DNekMatSharedPtr &mat) override
 
void v_PhysInterp (std::shared_ptr< StdExpansion > fromExp, const Array< OneD, const NekDouble > &fromData, Array< OneD, NekDouble > &toData, bool Transpose) override
 
void v_ReOrientTracePhysMap (const StdRegions::Orientation orient, Array< OneD, int > &idmap, const int nq0, const int nq1, bool Forwards) override
 This method produces a mapping.
 
int v_GetShapeDimension () const final
 
bool v_IsCollocatedBasis () const final
 
virtual void v_PhysDeriv (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2, Array< OneD, NekDouble > &out_d3)
 Calculate the derivative of the physical points.
 
virtual void v_PhysDeriv (const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0)
 Calculate the derivative of the physical points in a given direction.
 
- Protected Member Functions inherited from Nektar::StdRegions::StdExpansion
DNekMatSharedPtr CreateStdMatrix (const StdMatrixKey &mkey)
 
std::shared_ptr< Array< OneD, const NekDouble > > CreateStdFac (const StdFacKey &mkey)
 
DNekBlkMatSharedPtr CreateStdStaticCondMatrix (const StdMatrixKey &mkey)
 Create the static condensation of a matrix when using a boundary interior decomposition.
 
void GeneralMatrixOp_MatFree (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
 
void MassMatrixOp_MatFree (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
 
void LaplacianMatrixOp_MatFree (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
 
void LaplacianMatrixOp_MatFree_Kernel (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp)
 
void LaplacianMatrixOp_MatFree_GenericImpl (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
 
void LaplacianMatrixOp_MatFree (const int k1, const int k2, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
 
void WeakDerivMatrixOp_MatFree (const int i, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
 
void WeakDirectionalDerivMatrixOp_MatFree (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
 
void MassLevelCurvatureMatrixOp_MatFree (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
 
void LinearAdvectionMatrixOp_MatFree (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
 
void LinearAdvectionDiffusionReactionMatrixOp_MatFree (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey, bool addDiffusionTerm=true)
 
void HelmholtzMatrixOp_MatFree (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
 
void HelmholtzMatrixOp_MatFree_GenericImpl (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
 
virtual void v_SetCoeffsToOrientation (StdRegions::Orientation dir, Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
template<int DIR, bool DERIV = false, bool DERIV2 = false>
NekDouble BaryEvaluate (const NekDouble &coord, const NekDouble *physvals, NekDouble &deriv, NekDouble &deriv2)
 This function performs the barycentric interpolation of the polynomial stored in coord at a point physvals using barycentric interpolation weights in direction.
 
template<int DIR>
NekDouble BaryEvaluateBasis (const NekDouble &coord, const int &mode)
 
template<int DIR, bool DERIV = false, bool DERIV2 = false>
NekDouble BaryEvaluate (const NekDouble &coord, const NekDouble *physvals)
 Helper function to pass an unused value by reference into BaryEvaluate.
 
template<int DIR, bool DERIV = false, bool DERIV2 = false>
NekDouble BaryEvaluate (const NekDouble &coord, const NekDouble *physvals, NekDouble &deriv)
 
virtual const LibUtilities::PointsKey v_GetNodalPointsKey () const
 
virtual bool v_IsNodalNonTensorialExp ()
 
virtual void v_NodalToModal (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
virtual void v_FwdTrans (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 Transform a given function from physical quadrature space to coefficient space.
 
virtual void v_IProductWRTDirectionalDerivBase (const Array< OneD, const NekDouble > &direction, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
virtual void v_FwdTransBndConstrained (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
virtual void v_PhysDirectionalDeriv (const Array< OneD, const NekDouble > &inarray, const Array< OneD, const NekDouble > &direction, Array< OneD, NekDouble > &outarray)
 Physical derivative along a direction vector.
 
virtual NekDouble v_PhysEvaluate (const Array< OneD, const NekDouble > &coords, const Array< OneD, const NekDouble > &physvals)
 
virtual NekDouble v_PhysEvalFirstSecondDeriv (const Array< OneD, NekDouble > &coord, const Array< OneD, const NekDouble > &inarray, std::array< NekDouble, 3 > &firstOrderDerivs, std::array< NekDouble, 6 > &secondOrderDerivs)
 
virtual void v_GetCoord (const Array< OneD, const NekDouble > &Lcoord, Array< OneD, NekDouble > &coord)
 
virtual int v_GetCoordim () const
 
virtual void v_GetVertexPhysVals (const int vertex, const Array< OneD, const NekDouble > &inarray, NekDouble &outarray)
 
virtual void v_MultiplyByQuadratureMetric (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
virtual void v_MassMatrixOp (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
 
virtual void v_LaplacianMatrixOp (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
 
virtual void v_ExponentialFilter (Array< OneD, NekDouble > &array, const NekDouble alpha, const NekDouble exponent, const NekDouble cutoff)
 
virtual void v_LaplacianMatrixOp (const int k1, const int k2, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
 
virtual void v_WeakDerivMatrixOp (const int i, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
 
virtual void v_WeakDirectionalDerivMatrixOp (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
 
virtual void v_MassLevelCurvatureMatrixOp (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
 
virtual void v_LinearAdvectionMatrixOp (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
 
virtual void v_LinearAdvectionDiffusionReactionMatrixOp (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey, bool addDiffusionTerm=true)
 
virtual void v_HelmholtzMatrixOp (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
 
virtual void v_LaplacianMatrixOp_MatFree_Kernel (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp)
 
virtual DNekMatSharedPtr v_BuildInverseTransformationMatrix (const DNekScalMatSharedPtr &m_transformationmatrix)
 
virtual void v_GetSimplexEquiSpacedConnectivity (Array< OneD, int > &conn, bool standard=true)
 

Private Member Functions

int GetMode (int I, int J, int K)
 Compute the local mode number in the expansion for a particular tensorial combination.
 

Additional Inherited Members

- Protected Attributes inherited from Nektar::StdRegions::StdExpansion
Array< OneD, LibUtilities::BasisSharedPtrm_base
 
int m_elmt_id
 
int m_ncoeffs
 
std::vector< Array< OneD, const NekDouble > > m_weights
 
LibUtilities::NekManager< StdMatrixKey, DNekMat, StdMatrixKey::opLessm_stdMatrixManager
 
LibUtilities::NekManager< StdMatrixKey, DNekBlkMat, StdMatrixKey::opLessm_stdStaticCondMatrixManager
 
LibUtilities::NekManager< StdFacKey, Array< OneD, const NekDouble > > m_stdFacManager
 

Detailed Description

Class representing a prismatic element in reference space.

Definition at line 44 of file StdPrismExp.h.

Constructor & Destructor Documentation

◆ StdPrismExp() [1/4]

Nektar::StdRegions::StdPrismExp::StdPrismExp ( const LibUtilities::BasisKey Ba,
const LibUtilities::BasisKey Bb,
const LibUtilities::BasisKey Bc 
)

Definition at line 52 of file StdPrismExp.cpp.

56 Ba.GetNumModes(), Bb.GetNumModes(), Bc.GetNumModes()),
57 3, Ba, Bb, Bc),
59 Ba.GetNumModes(), Bb.GetNumModes(), Bc.GetNumModes()),
60 Ba, Bb, Bc)
61{
62 ASSERTL0(Ba.GetNumModes() <= Bc.GetNumModes(),
63 "order in 'a' direction is higher than order in 'c' direction");
64
65 // cache integration weights for future use
66 m_weights.push_back(m_base[0]->GetW());
67 m_weights.push_back(m_base[1]->GetW());
68
69 StdFacKey w2key(eWeights2, Bc);
70 // get weights[2] from manager where points are rescaled
71 m_weights.push_back(GetStdFac(w2key));
72}
#define ASSERTL0(condition, msg)
StdExpansion()
Default Constructor.
Array< OneD, const NekDouble > GetStdFac(const StdFacKey &mkey)
Array< OneD, LibUtilities::BasisSharedPtr > m_base
std::vector< Array< OneD, const NekDouble > > m_weights
constexpr int getNumberOfCoefficients(int Na, int Nb, int Nc)

References ASSERTL0, Nektar::StdRegions::eWeights2, Nektar::LibUtilities::BasisKey::GetNumModes(), Nektar::StdRegions::StdExpansion::GetStdFac(), Nektar::StdRegions::StdExpansion::m_base, and Nektar::StdRegions::StdExpansion::m_weights.

◆ StdPrismExp() [2/4]

Nektar::StdRegions::StdPrismExp::StdPrismExp ( const LibUtilities::BasisKey Ba,
const LibUtilities::BasisKey Bb,
const LibUtilities::BasisKey Bc,
NekDouble coeffs,
NekDouble phys 
)

◆ StdPrismExp() [3/4]

Nektar::StdRegions::StdPrismExp::StdPrismExp ( )
default

◆ StdPrismExp() [4/4]

Nektar::StdRegions::StdPrismExp::StdPrismExp ( const StdPrismExp T)
default

◆ ~StdPrismExp()

Nektar::StdRegions::StdPrismExp::~StdPrismExp ( )
overridedefault

Member Function Documentation

◆ GetMode()

int Nektar::StdRegions::StdPrismExp::GetMode ( int  p,
int  q,
int  r 
)
private

Compute the local mode number in the expansion for a particular tensorial combination.

Modes are numbered with the r index travelling fastest, followed by q and then p, and each q-r plane is of size (R+1-p). For example, with P=1, Q=2, R=3, the indexing inside each q-r plane (with r increasing upwards and q to the right) is:

p = 0: p = 1:

3 7 11 2 6 10 14 17 20 1 5 9 13 16 19 0 4 8 12 15 18

Note that in this element, we must have that \( P <= R \).

Definition at line 1947 of file StdPrismExp.cpp.

1948{
1949 int Q = m_base[1]->GetNumModes() - 1;
1950 int R = m_base[2]->GetNumModes() - 1;
1951
1952 return r + // Skip along stacks (r-direction)
1953 q * (R + 1 - p) + // Skip along columns (q-direction)
1954 (Q + 1) * (p * R + 1 -
1955 (p - 2) * (p - 1) / 2); // Skip along rows (p-direction)
1956}
std::vector< double > p(NPUPPER)
std::vector< double > q(NPUPPER *NPUPPER)

References Nektar::StdRegions::StdExpansion::m_base.

Referenced by v_GetBoundaryMap(), v_GetEdgeInteriorToElementMap(), v_GetInteriorMap(), v_GetTraceCoeffMap(), v_GetTraceInteriorToElementMap(), and v_GetVertexMap().

◆ v_BwdTrans()

void Nektar::StdRegions::StdPrismExp::v_BwdTrans ( const Array< OneD, const NekDouble > &  inarray,
Array< OneD, NekDouble > &  outarray 
)
overrideprotectedvirtual
Note
'r' (base[2]) runs fastest in this element.

Perform backwards transformation at the quadrature points:

\( u^{\delta} (\xi_{1i}, \xi_{2j}, \xi_{3k}) = \sum_{m(pqr)} \hat u_{pqr} \phi_{pqr} (\xi_{1i}, \xi_{2j}, \xi_{3k})\)

In the prism this expansion becomes:

\( u (\xi_{1i}, \xi_{2j}, \xi_{3k}) = \sum_{p=0}^{Q_x} \psi_p^a (\xi_{1i}) \lbrace { \sum_{q=0}^{Q_y} \psi_{q}^a (\xi_{2j}) \lbrace { \sum_{r=0}^{Q_z} \hat u_{pqr} \psi_{pr}^b (\xi_{3k}) \rbrace} \rbrace}. \)

And sumfactorizing step of the form is as:\

\( f_{pr} (\xi_{3k}) = \sum_{r=0}^{Q_z} \hat u_{pqr} \psi_{pr}^b (\xi_{3k}),\\ g_{p} (\xi_{2j}, \xi_{3k}) = \sum_{r=0}^{Q_y} \psi_{p}^a (\xi_{2j}) f_{pr} (\xi_{3k}),\ \ u(\xi_{1i}, \xi_{2j}, \xi_{3k}) = \sum_{p=0}^{Q_x} \psi_{p}^a (\xi_{1i}) g_{p} (\xi_{2j}, \xi_{3k}). \)

Implements Nektar::StdRegions::StdExpansion.

Definition at line 186 of file StdPrismExp.cpp.

188{
191 "Basis[1] is not a general tensor type");
192
195 "Basis[2] is not a general tensor type");
196
197 const Array<OneD, const NekDouble> base0 = m_base[0]->GetBdata();
198 const Array<OneD, const NekDouble> base1 = m_base[1]->GetBdata();
199 const Array<OneD, const NekDouble> base2 = m_base[2]->GetBdata();
200
201 int nquad0 = m_base[0]->GetNumPoints();
202 int nquad1 = m_base[1]->GetNumPoints();
203 int nquad2 = m_base[2]->GetNumPoints();
204 int nmodes0 = m_base[0]->GetNumModes();
205 int nmodes1 = m_base[1]->GetNumModes();
206 int nmodes2 = m_base[2]->GetNumModes();
207
208 bool isModified = (m_base[0]->GetBasisType() == LibUtilities::eModified_A);
209
210 std::vector<vec_t, tinysimd::allocator<vec_t>> wsp0(nmodes0 * nmodes1),
211 wsp1(nmodes0);
212
213 // Switch statment using boost_pp and macros. This unfolls intwo a
214 // nested swtich statement where the outer swtich statement runs
215 // from SMIN to SMAX for modal order and the inner switch
216 // statemets run from the outer value of the case to 2*SMAX for
217 // the quadrature order. If you want to see it unwrapped compile
218 // in verbose mode and add --preprocess to the c++ command.
219 // Default case
220#undef BWDTRANS_DEF
221#define BWDTRANS_DEF \
222 BwdTransPrismKernel( \
223 nmodes0, nmodes1, nmodes2, nquad0, nquad1, nquad2, isModified, \
224 (const vec_t *)base0.data(), (const vec_t *)base1.data(), \
225 (const vec_t *)base2.data(), wsp0.data(), wsp1.data(), \
226 (const vec_t *)inarray.data(), (vec_t *)outarray.data())
227
228 // Inner loop case over quarature points
229#undef BWDTRANS_Q
230#define BWDTRANS_Q(r, i) \
231 case NQ(i): \
232 BwdTransPrismKernel( \
233 NM(i), NM(i), NM(i), NQ(i), NQ(i), NQ_M1(i), isModified, \
234 (const vec_t *)base0.data(), (const vec_t *)base1.data(), \
235 (const vec_t *)base2.data(), wsp0.data(), wsp1.data(), \
236 (const vec_t *)inarray.data(), (vec_t *)outarray.data()); \
237 break;
238
239 // outer loop case over modes
240#undef BWDTRANS_M
241#define BWDTRANS_M(r, i) \
242 case NM(i): \
243 { \
244 switch (nquad0) \
245 { \
246 BOOST_PP_FOR_##r((NM(i), NM_P1(i), BOOST_PP_MUL(2, NM(i))), \
247 STDLEV2TEST1, STDLEV2UPDATE1, BWDTRANS_Q) default \
248 : BWDTRANS_DEF; \
249 break; \
250 } \
251 } \
252 break;
253
254 // templated cases on equi-ordered modes and standard quad
255 // usage where quad order goes from mode order to 2(*mode
256 // order)
257 if ((nmodes0 == nmodes1) && (nmodes1 == nmodes2) && (nquad0 == nquad1) &&
258 (nquad1 == nquad2 + 1))
259 {
260 switch (nmodes0)
261 {
262 BOOST_PP_FOR((SMIN, 0, SMAX), STDLEV2TEST, STDLEV2UPDATE,
264 default:
266 break;
267 }
268 }
269 else
270 {
272 }
273}
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
#define BWDTRANS_M(r, i)
#define BWDTRANS_DEF
#define STDLEV2TEST(r, state)
#define STDLEV2UPDATE(r, state)
LibUtilities::BasisType GetBasisType(const int dir) const
This function returns the type of basis used in the dir direction.
@ eModified_B
Principle Modified Functions .
Definition BasisType.h:49
@ eModified_C
Principle Modified Functions .
Definition BasisType.h:50
@ eOrtho_C
Principle Orthogonal Functions .
Definition BasisType.h:46
@ eOrtho_B
Principle Orthogonal Functions .
Definition BasisType.h:44
@ eModified_A
Principle Modified Functions .
Definition BasisType.h:48

References ASSERTL1, BWDTRANS_DEF, BWDTRANS_M, Nektar::LibUtilities::eModified_A, Nektar::LibUtilities::eModified_B, Nektar::LibUtilities::eModified_C, Nektar::LibUtilities::eOrtho_B, Nektar::LibUtilities::eOrtho_C, Nektar::StdRegions::StdExpansion::GetBasisType(), Nektar::StdRegions::StdExpansion::m_base, STDLEV2TEST, and STDLEV2UPDATE.

Referenced by Nektar::StdRegions::StdNodalPrismExp::v_BwdTrans(), and v_FillMode().

◆ v_CalcNumberOfCoefficients()

int Nektar::StdRegions::StdPrismExp::v_CalcNumberOfCoefficients ( const std::vector< unsigned int > &  nummodes,
int &  modes_offset 
)
overrideprotectedvirtual

Reimplemented from Nektar::StdRegions::StdExpansion.

Definition at line 920 of file StdPrismExp.cpp.

922{
924 nummodes[modes_offset], nummodes[modes_offset + 1],
925 nummodes[modes_offset + 2]);
926
927 modes_offset += 3;
928 return nmodes;
929}

References Nektar::LibUtilities::StdPrismData::getNumberOfCoefficients().

◆ v_CreateStdMatrix()

DNekMatSharedPtr Nektar::StdRegions::StdPrismExp::v_CreateStdMatrix ( const StdMatrixKey mkey)
overrideprotectedvirtual

Reimplemented from Nektar::StdRegions::StdExpansion.

Definition at line 1924 of file StdPrismExp.cpp.

1925{
1926 return v_GenMatrix(mkey);
1927}
DNekMatSharedPtr v_GenMatrix(const StdMatrixKey &mkey) override

References v_GenMatrix().

◆ v_DetShapeType()

LibUtilities::ShapeType Nektar::StdRegions::StdPrismExp::v_DetShapeType ( ) const
overrideprotectedvirtual

Return Shape of region, using ShapeType enum list; i.e. prism.

Implements Nektar::StdRegions::StdExpansion.

Definition at line 767 of file StdPrismExp.cpp.

768{
770}

References Nektar::LibUtilities::ePrism.

◆ v_FillMode()

void Nektar::StdRegions::StdPrismExp::v_FillMode ( const int  mode,
Array< OneD, NekDouble > &  outarray 
)
overrideprotectedvirtual

Reimplemented from Nektar::StdRegions::StdExpansion.

Definition at line 681 of file StdPrismExp.cpp.

682{
683 Array<OneD, NekDouble> tmp(m_ncoeffs, 0.0);
684 tmp[mode] = 1.0;
685 StdPrismExp::v_BwdTrans(tmp, outarray);
686}
void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override

References Nektar::StdRegions::StdExpansion::m_ncoeffs, and v_BwdTrans().

Referenced by Nektar::StdRegions::StdNodalPrismExp::GenNBasisTransMatrix().

◆ v_GenMatrix()

DNekMatSharedPtr Nektar::StdRegions::StdPrismExp::v_GenMatrix ( const StdMatrixKey mkey)
overrideprotectedvirtual

Reimplemented from Nektar::StdRegions::StdExpansion.

Definition at line 1769 of file StdPrismExp.cpp.

1770{
1771
1772 MatrixType mtype = mkey.GetMatrixType();
1773
1774 DNekMatSharedPtr Mat;
1775
1776 switch (mtype)
1777 {
1779 {
1780 int nq0 = m_base[0]->GetNumPoints();
1781 int nq1 = m_base[1]->GetNumPoints();
1782 int nq2 = m_base[2]->GetNumPoints();
1783 int nq;
1784
1785 // take definition from key
1786 if (mkey.ConstFactorExists(eFactorConst))
1787 {
1788 nq = (int)mkey.GetConstFactor(eFactorConst);
1789 }
1790 else
1791 {
1792 nq = max(nq0, max(nq1, nq2));
1793 }
1794
1795 int neq =
1797 Array<OneD, Array<OneD, NekDouble>> coords(neq);
1798 Array<OneD, NekDouble> coll(3);
1799 Array<OneD, DNekMatSharedPtr> I(3);
1800 Array<OneD, NekDouble> tmp(nq0);
1801
1802 Mat =
1803 MemoryManager<DNekMat>::AllocateSharedPtr(neq, nq0 * nq1 * nq2);
1804 int cnt = 0;
1805 for (int i = 0; i < nq; ++i)
1806 {
1807 for (int j = 0; j < nq; ++j)
1808 {
1809 for (int k = 0; k < nq - i; ++k, ++cnt)
1810 {
1811 coords[cnt] = Array<OneD, NekDouble>(3);
1812 coords[cnt][0] = -1.0 + 2 * k / (NekDouble)(nq - 1);
1813 coords[cnt][1] = -1.0 + 2 * j / (NekDouble)(nq - 1);
1814 coords[cnt][2] = -1.0 + 2 * i / (NekDouble)(nq - 1);
1815 }
1816 }
1817 }
1818
1819 for (int i = 0; i < neq; ++i)
1820 {
1821 LocCoordToLocCollapsed(coords[i], coll);
1822
1823 I[0] = m_base[0]->GetI(coll);
1824 I[1] = m_base[1]->GetI(coll + 1);
1825 I[2] = m_base[2]->GetI(coll + 2);
1826
1827 // interpolate first coordinate direction
1828 NekDouble fac;
1829 for (int k = 0; k < nq2; ++k)
1830 {
1831 for (int j = 0; j < nq1; ++j)
1832 {
1833
1834 fac = (I[1]->GetPtr())[j] * (I[2]->GetPtr())[k];
1835 Vmath::Smul(nq0, fac, I[0]->GetPtr(), 1, tmp, 1);
1836
1837 Vmath::Vcopy(nq0, &tmp[0], 1,
1838 Mat->GetRawPtr() + k * nq0 * nq1 * neq +
1839 j * nq0 * neq + i,
1840 neq);
1841 }
1842 }
1843 }
1844 }
1845 break;
1846 case ePhysInterpToGLL:
1847 {
1848 int nq0 = m_base[0]->GetNumPoints();
1849 int nq1 = m_base[1]->GetNumPoints();
1850 int nq2 = m_base[2]->GetNumPoints();
1851 int nq;
1852
1853 // take definition from key
1854 if (mkey.ConstFactorExists(eFactorConst))
1855 {
1856 nq = (int)mkey.GetConstFactor(eFactorConst);
1857 }
1858 else
1859 {
1860 nq = max(nq0, max(nq1, nq2));
1861 }
1862
1863 int neq =
1865 Array<OneD, NekDouble> coords(3);
1866 Array<OneD, NekDouble> coll(3);
1867 Array<OneD, DNekMatSharedPtr> I(3);
1868 Array<OneD, NekDouble> tmp(nq0);
1869
1870 Mat =
1871 MemoryManager<DNekMat>::AllocateSharedPtr(neq, nq0 * nq1 * nq2);
1872
1873 const LibUtilities::PointsKey key(nq,
1875
1876 Array<OneD, const NekDouble> x, y, z;
1877 LibUtilities::PointsManager()[key]->GetPoints(x, y, z);
1878
1879 Array<OneD, int> sorted;
1881
1882 for (int i = 0; i < neq; ++i)
1883 {
1884 coords[0] = x[sorted[i]];
1885 coords[1] = y[sorted[i]];
1886 coords[2] = z[sorted[i]];
1887
1888 LocCoordToLocCollapsed(coords, coll);
1889
1890 I[0] = m_base[0]->GetI(coll);
1891 I[1] = m_base[1]->GetI(coll + 1);
1892 I[2] = m_base[2]->GetI(coll + 2);
1893
1894 // interpolate first coordinate direction
1895 NekDouble fac;
1896 for (int k = 0; k < nq2; ++k)
1897 {
1898 for (int j = 0; j < nq1; ++j)
1899 {
1900
1901 fac = (I[1]->GetPtr())[j] * (I[2]->GetPtr())[k];
1902 Vmath::Smul(nq0, fac, I[0]->GetPtr(), 1, tmp, 1);
1903
1904 Vmath::Vcopy(nq0, &tmp[0], 1,
1905 Mat->GetRawPtr() + k * nq0 * nq1 * neq +
1906 j * nq0 * neq + i,
1907 neq);
1908 }
1909 }
1910 }
1911 }
1912 break;
1913
1914 default:
1915 {
1917 }
1918 break;
1919 }
1920
1921 return Mat;
1922}
static void CartesianOrdering(const int nq, Array< OneD, int > &sorted)
Definition NodalUtil.h:363
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
void LocCoordToLocCollapsed(const Array< OneD, const NekDouble > &xi, Array< OneD, NekDouble > &eta)
Convert local cartesian coordinate xi into local collapsed coordinates eta.
DNekMatSharedPtr CreateGeneralMatrix(const StdMatrixKey &mkey)
this function generates the mass matrix
PointsManagerT & PointsManager(void)
@ eNodalPrismElec
3D electrostatically spaced points on a Prism
Definition PointsType.h:87
std::vector< double > z(NPUPPER)
std::shared_ptr< DNekMat > DNekMatSharedPtr
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.hpp:100
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition Vmath.hpp:825
scalarT< T > max(scalarT< T > lhs, scalarT< T > rhs)
Definition scalar.hpp:305

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::LibUtilities::NodalUtilPrism::CartesianOrdering(), Nektar::StdRegions::StdMatrixKey::ConstFactorExists(), Nektar::StdRegions::StdExpansion::CreateGeneralMatrix(), Nektar::StdRegions::eFactorConst, Nektar::LibUtilities::eNodalPrismElec, Nektar::StdRegions::ePhysInterpToEquiSpaced, Nektar::StdRegions::ePhysInterpToGLL, Nektar::StdRegions::StdMatrixKey::GetConstFactor(), Nektar::StdRegions::StdMatrixKey::GetMatrixType(), Nektar::LibUtilities::StdPrismData::getNumberOfCoefficients(), Nektar::StdRegions::StdExpansion::LocCoordToLocCollapsed(), Nektar::StdRegions::StdExpansion::m_base, tinysimd::max(), Nektar::LibUtilities::PointsManager(), Vmath::Smul(), and Vmath::Vcopy().

Referenced by v_CreateStdMatrix().

◆ v_GetBoundaryMap()

void Nektar::StdRegions::StdPrismExp::v_GetBoundaryMap ( Array< OneD, unsigned int > &  outarray)
overrideprotectedvirtual

Reimplemented from Nektar::StdRegions::StdExpansion.

Definition at line 1045 of file StdPrismExp.cpp.

1046{
1049 "BasisType is not a boundary interior form");
1052 "BasisType is not a boundary interior form");
1055 "BasisType is not a boundary interior form");
1056
1057 int P = m_base[0]->GetNumModes() - 1, p;
1058 int Q = m_base[1]->GetNumModes() - 1, q;
1059 int R = m_base[2]->GetNumModes() - 1, r;
1060 int idx = 0;
1061
1062 int nBnd = NumBndryCoeffs();
1063
1064 if (maparray.size() != nBnd)
1065 {
1066 maparray = Array<OneD, unsigned int>(nBnd);
1067 }
1068
1069 // Loop over all boundary modes (in ascending order).
1070 for (p = 0; p <= P; ++p)
1071 {
1072 // First two q-r planes are entirely boundary modes.
1073 if (p <= 1)
1074 {
1075 for (q = 0; q <= Q; ++q)
1076 {
1077 for (r = 0; r <= R - p; ++r)
1078 {
1079 maparray[idx++] = GetMode(p, q, r);
1080 }
1081 }
1082 }
1083 else
1084 {
1085 // Remaining q-r planes contain boundary modes on the two
1086 // left-hand sides and bottom edge.
1087 for (q = 0; q <= Q; ++q)
1088 {
1089 if (q <= 1)
1090 {
1091 for (r = 0; r <= R - p; ++r)
1092 {
1093 maparray[idx++] = GetMode(p, q, r);
1094 }
1095 }
1096 else
1097 {
1098 maparray[idx++] = GetMode(p, q, 0);
1099 }
1100 }
1101 }
1102 }
1103}
int GetMode(int I, int J, int K)
Compute the local mode number in the expansion for a particular tensorial combination.
@ eGLL_Lagrange
Lagrange for SEM basis .
Definition BasisType.h:56

References ASSERTL1, Nektar::LibUtilities::eGLL_Lagrange, Nektar::LibUtilities::eModified_A, Nektar::LibUtilities::eModified_B, Nektar::StdRegions::StdExpansion::GetBasisType(), GetMode(), Nektar::StdRegions::StdExpansion::m_base, Nektar::StdRegions::StdExpansion::NumBndryCoeffs(), and Nektar::LibUtilities::P.

◆ v_GetCoords()

void Nektar::StdRegions::StdPrismExp::v_GetCoords ( Array< OneD, NekDouble > &  xi_x,
Array< OneD, NekDouble > &  xi_y,
Array< OneD, NekDouble > &  xi_z 
)
overrideprotectedvirtual

Reimplemented from Nektar::StdRegions::StdExpansion.

Definition at line 575 of file StdPrismExp.cpp.

578{
579 Array<OneD, const NekDouble> etaBar_x = m_base[0]->GetZ();
580 Array<OneD, const NekDouble> eta_y = m_base[1]->GetZ();
581 Array<OneD, const NekDouble> eta_z = m_base[2]->GetZ();
582 int Qx = GetNumPoints(0);
583 int Qy = GetNumPoints(1);
584 int Qz = GetNumPoints(2);
585
586 // Convert collapsed coordinates into cartesian coordinates: eta --> xi
587 for (int k = 0; k < Qz; ++k)
588 {
589 for (int j = 0; j < Qy; ++j)
590 {
591 for (int i = 0; i < Qx; ++i)
592 {
593 int s = i + Qx * (j + Qy * k);
594 xi_x[s] = (1.0 - eta_z[k]) * (1.0 + etaBar_x[i]) / 2.0 - 1.0;
595 xi_y[s] = eta_y[j];
596 xi_z[s] = eta_z[k];
597 }
598 }
599 }
600}
int GetNumPoints(const int dir) const
This function returns the number of quadrature points in the dir direction.

References Nektar::StdRegions::StdExpansion::GetNumPoints(), and Nektar::StdRegions::StdExpansion::m_base.

◆ v_GetEdgeInteriorToElementMap()

void Nektar::StdRegions::StdPrismExp::v_GetEdgeInteriorToElementMap ( const int  tid,
Array< OneD, unsigned int > &  maparray,
Array< OneD, int > &  signarray,
const Orientation  traceOrient = eDir1FwdDir1_Dir2FwdDir2 
)
overrideprotectedvirtual

Reimplemented from Nektar::StdRegions::StdExpansion3D.

Definition at line 1464 of file StdPrismExp.cpp.

1467{
1468 int i;
1469 bool signChange;
1470 const int P = m_base[0]->GetNumModes() - 1;
1471 const int Q = m_base[1]->GetNumModes() - 1;
1472 const int R = m_base[2]->GetNumModes() - 1;
1473 const int nEdgeIntCoeffs = v_GetEdgeNcoeffs(eid) - 2;
1474
1475 if (maparray.size() != nEdgeIntCoeffs)
1476 {
1477 maparray = Array<OneD, unsigned int>(nEdgeIntCoeffs);
1478 }
1479
1480 if (signarray.size() != nEdgeIntCoeffs)
1481 {
1482 signarray = Array<OneD, int>(nEdgeIntCoeffs, 1);
1483 }
1484 else
1485 {
1486 fill(signarray.data(), signarray.data() + nEdgeIntCoeffs, 1);
1487 }
1488
1489 // If edge is oriented backwards, change sign of modes which have
1490 // degree 2n+1, n >= 1.
1491 signChange = edgeOrient == eBackwards;
1492
1493 switch (eid)
1494 {
1495 case 0:
1496 for (i = 2; i <= P; ++i)
1497 {
1498 maparray[i - 2] = GetMode(i, 0, 0);
1499 }
1500 break;
1501
1502 case 1:
1503 for (i = 2; i <= Q; ++i)
1504 {
1505 maparray[i - 2] = GetMode(1, i, 0);
1506 }
1507 break;
1508
1509 case 2:
1510 // Base quad; reverse direction.
1511 // signChange = !signChange;
1512 for (i = 2; i <= P; ++i)
1513 {
1514 maparray[i - 2] = GetMode(i, 1, 0);
1515 }
1516 break;
1517
1518 case 3:
1519 // Base quad; reverse direction.
1520 // signChange = !signChange;
1521 for (i = 2; i <= Q; ++i)
1522 {
1523 maparray[i - 2] = GetMode(0, i, 0);
1524 }
1525 break;
1526
1527 case 4:
1528 for (i = 2; i <= R; ++i)
1529 {
1530 maparray[i - 2] = GetMode(0, 0, i);
1531 }
1532 break;
1533
1534 case 5:
1535 for (i = 1; i <= R - 1; ++i)
1536 {
1537 maparray[i - 1] = GetMode(1, 0, i);
1538 }
1539 break;
1540
1541 case 6:
1542 for (i = 1; i <= R - 1; ++i)
1543 {
1544 maparray[i - 1] = GetMode(1, 1, i);
1545 }
1546 break;
1547
1548 case 7:
1549 for (i = 2; i <= R; ++i)
1550 {
1551 maparray[i - 2] = GetMode(0, 1, i);
1552 }
1553 break;
1554
1555 case 8:
1556 for (i = 2; i <= Q; ++i)
1557 {
1558 maparray[i - 2] = GetMode(0, i, 1);
1559 }
1560 break;
1561
1562 default:
1563 ASSERTL0(false, "Edge not defined.");
1564 break;
1565 }
1566
1567 if (signChange)
1568 {
1569 for (i = 1; i < nEdgeIntCoeffs; i += 2)
1570 {
1571 signarray[i] = -1;
1572 }
1573 }
1574}
int v_GetEdgeNcoeffs(const int i) const override

References ASSERTL0, Nektar::StdRegions::eBackwards, GetMode(), Nektar::StdRegions::StdExpansion::m_base, Nektar::LibUtilities::P, and v_GetEdgeNcoeffs().

◆ v_GetEdgeNcoeffs()

int Nektar::StdRegions::StdPrismExp::v_GetEdgeNcoeffs ( const int  i) const
overrideprotectedvirtual

Reimplemented from Nektar::StdRegions::StdExpansion3D.

Definition at line 726 of file StdPrismExp.cpp.

727{
728 ASSERTL2(i >= 0 && i <= 8, "edge id is out of range");
729
730 if (i == 0 || i == 2)
731 {
732 return GetBasisNumModes(0);
733 }
734 else if (i == 1 || i == 3 || i == 8)
735 {
736 return GetBasisNumModes(1);
737 }
738 else
739 {
740 return GetBasisNumModes(2);
741 }
742}
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed to...
int GetBasisNumModes(const int dir) const
This function returns the number of expansion modes in the dir direction.

References ASSERTL2, and Nektar::StdRegions::StdExpansion::GetBasisNumModes().

Referenced by v_GetEdgeInteriorToElementMap().

◆ v_GetElmtTraceToTraceMap()

void Nektar::StdRegions::StdPrismExp::v_GetElmtTraceToTraceMap ( const unsigned int  fid,
Array< OneD, unsigned int > &  maparray,
Array< OneD, int > &  signarray,
Orientation  faceOrient,
int  P,
int  Q 
)
overrideprotectedvirtual

Reimplemented from Nektar::StdRegions::StdExpansion.

Definition at line 1206 of file StdPrismExp.cpp.

1210{
1212 "Method only implemented if BasisType is identical"
1213 "in x and y directions");
1216 "Method only implemented for Modified_A BasisType"
1217 "(x and y direction) and Modified_B BasisType (z "
1218 "direction)");
1219
1220 int i, j, k, p, r, nFaceCoeffs, idx = 0;
1221 int nummodesA = 0, nummodesB = 0;
1222
1223 switch (fid)
1224 {
1225 case 0:
1226 nummodesA = m_base[0]->GetNumModes();
1227 nummodesB = m_base[1]->GetNumModes();
1228 break;
1229 case 1:
1230 case 3:
1231 nummodesA = m_base[0]->GetNumModes();
1232 nummodesB = m_base[2]->GetNumModes();
1233 break;
1234 case 2:
1235 case 4:
1236 nummodesA = m_base[1]->GetNumModes();
1237 nummodesB = m_base[2]->GetNumModes();
1238 break;
1239 default:
1240 ASSERTL0(false, "fid must be between 0 and 4");
1241 }
1242
1243 if (P == -1)
1244 {
1245 P = nummodesA;
1246 Q = nummodesB;
1247 nFaceCoeffs = GetTraceNcoeffs(fid);
1248 }
1249 else if (fid == 1 || fid == 3)
1250 {
1251 nFaceCoeffs = P * (2 * Q - P + 1) / 2;
1252 }
1253 else
1254 {
1255 nFaceCoeffs = P * Q;
1256 }
1257
1258 // Allocate the map array and sign array; set sign array to ones (+)
1259 if (maparray.size() != nFaceCoeffs)
1260 {
1261 maparray = Array<OneD, unsigned int>(nFaceCoeffs);
1262 }
1263
1264 if (signarray.size() != nFaceCoeffs)
1265 {
1266 signarray = Array<OneD, int>(nFaceCoeffs, 1);
1267 }
1268 else
1269 {
1270 fill(signarray.data(), signarray.data() + nFaceCoeffs, 1);
1271 }
1272
1273 int minPA = min(nummodesA, P);
1274 int minQB = min(nummodesB, Q);
1275 // triangular faces
1276 if (fid == 1 || fid == 3)
1277 {
1278 // zero signmap and set maparray to zero if elemental
1279 // modes are not as large as face modesl
1280 idx = 0;
1281 int cnt = 0;
1282
1283 for (j = 0; j < minPA; ++j)
1284 {
1285 // set maparray
1286 for (k = 0; k < minQB - j; ++k, ++cnt)
1287 {
1288 maparray[idx++] = cnt;
1289 }
1290
1291 cnt += nummodesB - minQB;
1292
1293 // idx += nummodesB-j;
1294 for (k = nummodesB - j; k < Q - j; ++k)
1295 {
1296 signarray[idx] = 0.0;
1297 maparray[idx++] = maparray[0];
1298 }
1299 }
1300
1301 for (j = nummodesA; j < P; ++j)
1302 {
1303 for (k = 0; k < Q - j; ++k)
1304 {
1305 signarray[idx] = 0.0;
1306 maparray[idx++] = maparray[0];
1307 }
1308 }
1309
1310 // Triangles only have one possible orientation (base
1311 // direction reversed); swap edge modes.
1312 if (faceOrient == eDir1BwdDir1_Dir2FwdDir2)
1313 {
1314 idx = 0;
1315 for (p = 0; p < P; ++p)
1316 {
1317 for (r = 0; r < Q - p; ++r, idx++)
1318 {
1319 if (p > 1)
1320 {
1321 signarray[idx] = p % 2 ? -1 : 1;
1322 }
1323 }
1324 }
1325
1326 swap(maparray[0], maparray[Q]);
1327 for (i = 1; i < Q - 1; ++i)
1328 {
1329 swap(maparray[i + 1], maparray[Q + i]);
1330 }
1331 }
1332 }
1333 else
1334 {
1335 // Set up an array indexing for quads, since the
1336 // ordering may need to be transposed.
1337 Array<OneD, int> arrayindx(nFaceCoeffs, -1);
1338
1339 for (i = 0; i < Q; i++)
1340 {
1341 for (j = 0; j < P; j++)
1342 {
1343 if (faceOrient < eDir1FwdDir2_Dir2FwdDir1)
1344 {
1345 arrayindx[i * P + j] = i * P + j;
1346 }
1347 else
1348 {
1349 arrayindx[i * P + j] = j * Q + i;
1350 }
1351 }
1352 }
1353
1354 // zero signmap and set maparray to zero if elemental
1355 // modes are not as large as face modesl
1356 for (j = 0; j < P; ++j)
1357 {
1358 // set up default maparray
1359 for (k = 0; k < Q; k++)
1360 {
1361 maparray[arrayindx[j + k * P]] = j + k * nummodesA;
1362 }
1363
1364 for (k = nummodesB; k < Q; ++k)
1365 {
1366 signarray[arrayindx[j + k * P]] = 0.0;
1367 maparray[arrayindx[j + k * P]] = maparray[0];
1368 }
1369 }
1370
1371 for (j = nummodesA; j < P; ++j)
1372 {
1373 for (k = 0; k < Q; ++k)
1374 {
1375 signarray[arrayindx[j + k * P]] = 0.0;
1376 maparray[arrayindx[j + k * P]] = maparray[0];
1377 }
1378 }
1379
1380 // The code below is exactly the same as that taken from
1381 // StdHexExp and reverses the 'b' and 'a' directions as
1382 // appropriate (1st and 2nd if statements respectively) in
1383 // quadrilateral faces.
1384 if (faceOrient == eDir1FwdDir1_Dir2BwdDir2 ||
1385 faceOrient == eDir1BwdDir1_Dir2BwdDir2 ||
1386 faceOrient == eDir1BwdDir2_Dir2FwdDir1 ||
1387 faceOrient == eDir1BwdDir2_Dir2BwdDir1)
1388 {
1389 if (faceOrient < eDir1FwdDir2_Dir2FwdDir1)
1390 {
1391 for (i = 3; i < Q; i += 2)
1392 {
1393 for (j = 0; j < P; j++)
1394 {
1395 signarray[arrayindx[i * P + j]] *= -1;
1396 }
1397 }
1398
1399 for (i = 0; i < P; i++)
1400 {
1401 swap(maparray[i], maparray[i + P]);
1402 swap(signarray[i], signarray[i + P]);
1403 }
1404 }
1405 else
1406 {
1407 for (i = 0; i < Q; i++)
1408 {
1409 for (j = 3; j < P; j += 2)
1410 {
1411 signarray[arrayindx[i * P + j]] *= -1;
1412 }
1413 }
1414
1415 for (i = 0; i < Q; i++)
1416 {
1417 swap(maparray[i], maparray[i + Q]);
1418 swap(signarray[i], signarray[i + Q]);
1419 }
1420 }
1421 }
1422
1423 if (faceOrient == eDir1BwdDir1_Dir2FwdDir2 ||
1424 faceOrient == eDir1BwdDir1_Dir2BwdDir2 ||
1425 faceOrient == eDir1FwdDir2_Dir2BwdDir1 ||
1426 faceOrient == eDir1BwdDir2_Dir2BwdDir1)
1427 {
1428 if (faceOrient < eDir1FwdDir2_Dir2FwdDir1)
1429 {
1430 for (i = 0; i < Q; i++)
1431 {
1432 for (j = 3; j < P; j += 2)
1433 {
1434 signarray[arrayindx[i * P + j]] *= -1;
1435 }
1436 }
1437
1438 for (i = 0; i < Q; i++)
1439 {
1440 swap(maparray[i * P], maparray[i * P + 1]);
1441 swap(signarray[i * P], signarray[i * P + 1]);
1442 }
1443 }
1444 else
1445 {
1446 for (i = 3; i < Q; i += 2)
1447 {
1448 for (j = 0; j < P; j++)
1449 {
1450 signarray[arrayindx[i * P + j]] *= -1;
1451 }
1452 }
1453
1454 for (i = 0; i < P; i++)
1455 {
1456 swap(maparray[i * Q], maparray[i * Q + 1]);
1457 swap(signarray[i * Q], signarray[i * Q + 1]);
1458 }
1459 }
1460 }
1461 }
1462}
int GetTraceNcoeffs(const int i) const
This function returns the number of expansion coefficients belonging to the i-th trace.
scalarT< T > min(scalarT< T > lhs, scalarT< T > rhs)
Definition scalar.hpp:300

References ASSERTL0, ASSERTL1, Nektar::StdRegions::eDir1BwdDir1_Dir2BwdDir2, Nektar::StdRegions::eDir1BwdDir1_Dir2FwdDir2, Nektar::StdRegions::eDir1BwdDir2_Dir2BwdDir1, Nektar::StdRegions::eDir1BwdDir2_Dir2FwdDir1, Nektar::StdRegions::eDir1FwdDir1_Dir2BwdDir2, Nektar::StdRegions::eDir1FwdDir2_Dir2BwdDir1, Nektar::StdRegions::eDir1FwdDir2_Dir2FwdDir1, Nektar::LibUtilities::eModified_A, Nektar::LibUtilities::eModified_B, Nektar::StdRegions::StdExpansion::GetBasisType(), Nektar::StdRegions::StdExpansion::GetTraceNcoeffs(), Nektar::StdRegions::StdExpansion::m_base, tinysimd::min(), and Nektar::LibUtilities::P.

◆ v_GetInteriorMap()

void Nektar::StdRegions::StdPrismExp::v_GetInteriorMap ( Array< OneD, unsigned int > &  outarray)
overrideprotectedvirtual

Reimplemented from Nektar::StdRegions::StdExpansion.

Definition at line 1007 of file StdPrismExp.cpp.

1008{
1011 "BasisType is not a boundary interior form");
1014 "BasisType is not a boundary interior form");
1017 "BasisType is not a boundary interior form");
1018
1019 int P = m_base[0]->GetNumModes() - 1, p;
1020 int Q = m_base[1]->GetNumModes() - 1, q;
1021 int R = m_base[2]->GetNumModes() - 1, r;
1022
1023 int nIntCoeffs = m_ncoeffs - NumBndryCoeffs();
1024
1025 if (outarray.size() != nIntCoeffs)
1026 {
1027 outarray = Array<OneD, unsigned int>(nIntCoeffs);
1028 }
1029
1030 int idx = 0;
1031
1032 // Loop over all interior modes.
1033 for (p = 2; p <= P; ++p)
1034 {
1035 for (q = 2; q <= Q; ++q)
1036 {
1037 for (r = 1; r <= R - p; ++r)
1038 {
1039 outarray[idx++] = GetMode(p, q, r);
1040 }
1041 }
1042 }
1043}

References ASSERTL1, Nektar::LibUtilities::eGLL_Lagrange, Nektar::LibUtilities::eModified_A, Nektar::LibUtilities::eModified_B, Nektar::StdRegions::StdExpansion::GetBasisType(), GetMode(), Nektar::StdRegions::StdExpansion::m_base, Nektar::StdRegions::StdExpansion::m_ncoeffs, Nektar::StdRegions::StdExpansion::NumBndryCoeffs(), and Nektar::LibUtilities::P.

◆ v_GetNedges()

int Nektar::StdRegions::StdPrismExp::v_GetNedges ( void  ) const
finalprotectedvirtual

Reimplemented from Nektar::StdRegions::StdExpansion3D.

Definition at line 753 of file StdPrismExp.cpp.

754{
755 return 9;
756}

◆ v_GetNtraces()

int Nektar::StdRegions::StdPrismExp::v_GetNtraces ( ) const
finalprotectedvirtual

Implements Nektar::StdRegions::StdExpansion.

Definition at line 758 of file StdPrismExp.cpp.

759{
760 return 5;
761}

◆ v_GetNverts()

int Nektar::StdRegions::StdPrismExp::v_GetNverts ( ) const
finalprotectedvirtual

Implements Nektar::StdRegions::StdExpansion.

Definition at line 748 of file StdPrismExp.cpp.

749{
750 return 6;
751}

◆ v_GetTraceBasisKey()

const LibUtilities::BasisKey Nektar::StdRegions::StdPrismExp::v_GetTraceBasisKey ( const int  i,
const int  k,
bool  UseGLL = false 
) const
overrideprotectedvirtual

Reimplemented from Nektar::StdRegions::StdExpansion.

Definition at line 890 of file StdPrismExp.cpp.

893{
894 ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
895 ASSERTL2(k >= 0 && k <= 1, "basis key id is out of range");
896
897 switch (i)
898 {
899 case 0:
900 {
901 return EvaluateQuadFaceBasisKey(k, m_base[k]);
902 }
903 case 2:
904 case 4:
905 {
906 return EvaluateQuadFaceBasisKey(k, m_base[k + 1]);
907 }
908 case 1:
909 case 3:
910 {
911 return EvaluateTriFaceBasisKey(k, m_base[2 * k], UseGLL);
912 }
913 break;
914 }
915
916 // Should never get here.
918}
static const BasisKey NullBasisKey(eNoBasisType, 0, NullPointsKey)
Defines a null basis with no type or points.
LibUtilities::BasisKey EvaluateQuadFaceBasisKey(const int facedir, const LibUtilities::BasisSharedPtr &faceDirBasis)
LibUtilities::BasisKey EvaluateTriFaceBasisKey(const int facedir, const LibUtilities::BasisSharedPtr &faceDirBasis, bool UseGLL)

References ASSERTL2, Nektar::StdRegions::EvaluateQuadFaceBasisKey(), Nektar::StdRegions::EvaluateTriFaceBasisKey(), Nektar::StdRegions::StdExpansion::m_base, and Nektar::LibUtilities::NullBasisKey().

◆ v_GetTraceCoeffMap()

void Nektar::StdRegions::StdPrismExp::v_GetTraceCoeffMap ( const unsigned int  fid,
Array< OneD, unsigned int > &  maparray 
)
overrideprotectedvirtual

Reimplemented from Nektar::StdRegions::StdExpansion.

Definition at line 1105 of file StdPrismExp.cpp.

1107{
1109 "Method only implemented if BasisType is identical"
1110 "in x and y directions");
1113 "Method only implemented for Modified_A BasisType"
1114 "(x and y direction) and Modified_B BasisType (z "
1115 "direction)");
1116 int p, q, r, idx = 0;
1117 int P = 0, Q = 0;
1118
1119 switch (fid)
1120 {
1121 case 0:
1122 P = m_base[0]->GetNumModes();
1123 Q = m_base[1]->GetNumModes();
1124 break;
1125 case 1:
1126 case 3:
1127 P = m_base[0]->GetNumModes();
1128 Q = m_base[2]->GetNumModes();
1129 break;
1130 case 2:
1131 case 4:
1132 P = m_base[1]->GetNumModes();
1133 Q = m_base[2]->GetNumModes();
1134 break;
1135 default:
1136 ASSERTL0(false, "fid must be between 0 and 4");
1137 }
1138
1139 if (maparray.size() != P * Q)
1140 {
1141 maparray = Array<OneD, unsigned int>(P * Q);
1142 }
1143
1144 // Set up ordering inside each 2D face. Also for triangular faces,
1145 // populate signarray.
1146 switch (fid)
1147 {
1148 case 0: // Bottom quad
1149 for (q = 0; q < Q; ++q)
1150 {
1151 for (p = 0; p < P; ++p)
1152 {
1153 maparray[q * P + p] = GetMode(p, q, 0);
1154 }
1155 }
1156 break;
1157 case 1: // Left triangle
1158 for (p = 0; p < P; ++p)
1159 {
1160 for (r = 0; r < Q - p; ++r)
1161 {
1162 maparray[idx++] = GetMode(p, 0, r);
1163 }
1164 }
1165 break;
1166 case 2: // Slanted quad
1167 for (q = 0; q < P; ++q)
1168 {
1169 maparray[q] = GetMode(1, q, 0);
1170 }
1171 for (q = 0; q < P; ++q)
1172 {
1173 maparray[P + q] = GetMode(0, q, 1);
1174 }
1175 for (r = 1; r < Q - 1; ++r)
1176 {
1177 for (q = 0; q < P; ++q)
1178 {
1179 maparray[(r + 1) * P + q] = GetMode(1, q, r);
1180 }
1181 }
1182 break;
1183 case 3: // Right triangle
1184 for (p = 0; p < P; ++p)
1185 {
1186 for (r = 0; r < Q - p; ++r)
1187 {
1188 maparray[idx++] = GetMode(p, 1, r);
1189 }
1190 }
1191 break;
1192 case 4: // Rear quad
1193 for (r = 0; r < Q; ++r)
1194 {
1195 for (q = 0; q < P; ++q)
1196 {
1197 maparray[r * P + q] = GetMode(0, q, r);
1198 }
1199 }
1200 break;
1201 default:
1202 ASSERTL0(false, "Face to element map unavailable.");
1203 }
1204}

References ASSERTL0, ASSERTL1, Nektar::LibUtilities::eModified_A, Nektar::LibUtilities::eModified_B, Nektar::StdRegions::StdExpansion::GetBasisType(), GetMode(), Nektar::StdRegions::StdExpansion::m_base, and Nektar::LibUtilities::P.

◆ v_GetTraceInteriorToElementMap()

void Nektar::StdRegions::StdPrismExp::v_GetTraceInteriorToElementMap ( const int  tid,
Array< OneD, unsigned int > &  maparray,
Array< OneD, int > &  signarray,
const Orientation  traceOrient = eDir1FwdDir1_Dir2FwdDir2 
)
overrideprotectedvirtual

Reimplemented from Nektar::StdRegions::StdExpansion.

Definition at line 1576 of file StdPrismExp.cpp.

1579{
1580 const int P = m_base[0]->GetNumModes() - 1;
1581 const int Q = m_base[1]->GetNumModes() - 1;
1582 const int R = m_base[2]->GetNumModes() - 1;
1583 const int nFaceIntCoeffs = v_GetTraceIntNcoeffs(fid);
1584 int p, q, r, idx = 0;
1585 int nummodesA = 0;
1586 int nummodesB = 0;
1587 int i = 0;
1588 int j = 0;
1589
1590 if (maparray.size() != nFaceIntCoeffs)
1591 {
1592 maparray = Array<OneD, unsigned int>(nFaceIntCoeffs);
1593 }
1594
1595 if (signarray.size() != nFaceIntCoeffs)
1596 {
1597 signarray = Array<OneD, int>(nFaceIntCoeffs, 1);
1598 }
1599 else
1600 {
1601 fill(signarray.data(), signarray.data() + nFaceIntCoeffs, 1);
1602 }
1603
1604 // Set up an array indexing for quad faces, since the ordering may
1605 // need to be transposed depending on orientation.
1606 Array<OneD, int> arrayindx(nFaceIntCoeffs);
1607 if (fid != 1 && fid != 3)
1608 {
1609 if (fid == 0) // Base quad
1610 {
1611 nummodesA = P - 1;
1612 nummodesB = Q - 1;
1613 }
1614 else // front and back quad
1615 {
1616 nummodesA = Q - 1;
1617 nummodesB = R - 1;
1618 }
1619
1620 for (i = 0; i < nummodesB; i++)
1621 {
1622 for (j = 0; j < nummodesA; j++)
1623 {
1624 if (faceOrient < eDir1FwdDir2_Dir2FwdDir1)
1625 {
1626 arrayindx[i * nummodesA + j] = i * nummodesA + j;
1627 }
1628 else
1629 {
1630 arrayindx[i * nummodesA + j] = j * nummodesB + i;
1631 }
1632 }
1633 }
1634 }
1635
1636 switch (fid)
1637 {
1638 case 0: // Bottom quad
1639 for (q = 2; q <= Q; ++q)
1640 {
1641 for (p = 2; p <= P; ++p)
1642 {
1643 maparray[arrayindx[(q - 2) * nummodesA + (p - 2)]] =
1644 GetMode(p, q, 0);
1645 }
1646 }
1647 break;
1648
1649 case 1: // Left triangle
1650 for (p = 2; p <= P; ++p)
1651 {
1652 for (r = 1; r <= R - p; ++r)
1653 {
1654 if (faceOrient == eDir1BwdDir1_Dir2FwdDir2)
1655 {
1656 signarray[idx] = p % 2 ? -1 : 1;
1657 }
1658 maparray[idx++] = GetMode(p, 0, r);
1659 }
1660 }
1661 break;
1662
1663 case 2: // Slanted quad
1664 for (r = 1; r <= R - 1; ++r)
1665 {
1666 for (q = 2; q <= Q; ++q)
1667 {
1668 maparray[arrayindx[(r - 1) * nummodesA + (q - 2)]] =
1669 GetMode(1, q, r);
1670 }
1671 }
1672 break;
1673
1674 case 3: // Right triangle
1675 for (p = 2; p <= P; ++p)
1676 {
1677 for (r = 1; r <= R - p; ++r)
1678 {
1679 if (faceOrient == eDir1BwdDir1_Dir2FwdDir2)
1680 {
1681 signarray[idx] = p % 2 ? -1 : 1;
1682 }
1683 maparray[idx++] = GetMode(p, 1, r);
1684 }
1685 }
1686 break;
1687
1688 case 4: // Back quad
1689 for (r = 2; r <= R; ++r)
1690 {
1691 for (q = 2; q <= Q; ++q)
1692 {
1693 maparray[arrayindx[(r - 2) * nummodesA + (q - 2)]] =
1694 GetMode(0, q, r);
1695 }
1696 }
1697 break;
1698
1699 default:
1700 ASSERTL0(false, "Face interior map not available.");
1701 }
1702
1703 // Triangular faces are processed in the above switch loop; for
1704 // remaining quad faces, set up orientation if necessary.
1705 if (fid == 1 || fid == 3)
1706 {
1707 return;
1708 }
1709
1710 if (faceOrient == eDir1FwdDir1_Dir2BwdDir2 ||
1711 faceOrient == eDir1BwdDir1_Dir2BwdDir2 ||
1712 faceOrient == eDir1BwdDir2_Dir2FwdDir1 ||
1713 faceOrient == eDir1BwdDir2_Dir2BwdDir1)
1714 {
1715 if (faceOrient < eDir1FwdDir2_Dir2FwdDir1)
1716 {
1717 for (i = 1; i < nummodesB; i += 2)
1718 {
1719 for (j = 0; j < nummodesA; j++)
1720 {
1721 signarray[arrayindx[i * nummodesA + j]] *= -1;
1722 }
1723 }
1724 }
1725 else
1726 {
1727 for (i = 0; i < nummodesB; i++)
1728 {
1729 for (j = 1; j < nummodesA; j += 2)
1730 {
1731 signarray[arrayindx[i * nummodesA + j]] *= -1;
1732 }
1733 }
1734 }
1735 }
1736
1737 if (faceOrient == eDir1BwdDir1_Dir2FwdDir2 ||
1738 faceOrient == eDir1BwdDir1_Dir2BwdDir2 ||
1739 faceOrient == eDir1FwdDir2_Dir2BwdDir1 ||
1740 faceOrient == eDir1BwdDir2_Dir2BwdDir1)
1741 {
1742 if (faceOrient < eDir1FwdDir2_Dir2FwdDir1)
1743 {
1744 for (i = 0; i < nummodesB; i++)
1745 {
1746 for (j = 1; j < nummodesA; j += 2)
1747 {
1748 signarray[arrayindx[i * nummodesA + j]] *= -1;
1749 }
1750 }
1751 }
1752 else
1753 {
1754 for (i = 1; i < nummodesB; i += 2)
1755 {
1756 for (j = 0; j < nummodesA; j++)
1757 {
1758 signarray[arrayindx[i * nummodesA + j]] *= -1;
1759 }
1760 }
1761 }
1762 }
1763}
int v_GetTraceIntNcoeffs(const int i) const override

References ASSERTL0, Nektar::StdRegions::eDir1BwdDir1_Dir2BwdDir2, Nektar::StdRegions::eDir1BwdDir1_Dir2FwdDir2, Nektar::StdRegions::eDir1BwdDir2_Dir2BwdDir1, Nektar::StdRegions::eDir1BwdDir2_Dir2FwdDir1, Nektar::StdRegions::eDir1FwdDir1_Dir2BwdDir2, Nektar::StdRegions::eDir1FwdDir2_Dir2BwdDir1, Nektar::StdRegions::eDir1FwdDir2_Dir2FwdDir1, GetMode(), Nektar::StdRegions::StdExpansion::m_base, Nektar::LibUtilities::P, and v_GetTraceIntNcoeffs().

◆ v_GetTraceIntNcoeffs()

int Nektar::StdRegions::StdPrismExp::v_GetTraceIntNcoeffs ( const int  i) const
overrideprotectedvirtual

Implements Nektar::StdRegions::StdExpansion.

Definition at line 830 of file StdPrismExp.cpp.

831{
832 ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
833
834 int Pi = GetBasisNumModes(0) - 2;
835 int Qi = GetBasisNumModes(1) - 2;
836 int Ri = GetBasisNumModes(2) - 2;
837
838 if (i == 0)
839 {
840 return Pi * Qi;
841 }
842 else if (i == 1 || i == 3)
843 {
844 return Pi * (2 * Ri - Pi - 1) / 2;
845 }
846 else
847 {
848 return Qi * Ri;
849 }
850}

References ASSERTL2, and Nektar::StdRegions::StdExpansion::GetBasisNumModes().

Referenced by v_GetTraceInteriorToElementMap().

◆ v_GetTraceNcoeffs()

int Nektar::StdRegions::StdPrismExp::v_GetTraceNcoeffs ( const int  i) const
overrideprotectedvirtual

Implements Nektar::StdRegions::StdExpansion.

Definition at line 812 of file StdPrismExp.cpp.

813{
814 ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
815 if (i == 0)
816 {
817 return GetBasisNumModes(0) * GetBasisNumModes(1);
818 }
819 else if (i == 1 || i == 3)
820 {
821 int P = GetBasisNumModes(0) - 1, Q = GetBasisNumModes(2) - 1;
822 return Q + 1 + (P * (1 + 2 * Q - P)) / 2;
823 }
824 else
825 {
826 return GetBasisNumModes(1) * GetBasisNumModes(2);
827 }
828}

References ASSERTL2, Nektar::StdRegions::StdExpansion::GetBasisNumModes(), and Nektar::LibUtilities::P.

◆ v_GetTraceNumModes()

void Nektar::StdRegions::StdPrismExp::v_GetTraceNumModes ( const int  fid,
int &  numModes0,
int &  numModes1,
Orientation  faceOrient = eDir1FwdDir1_Dir2FwdDir2 
)
overrideprotectedvirtual

Reimplemented from Nektar::StdRegions::StdExpansion.

Definition at line 688 of file StdPrismExp.cpp.

690{
691 int nummodes[3] = {m_base[0]->GetNumModes(), m_base[1]->GetNumModes(),
692 m_base[2]->GetNumModes()};
693 switch (fid)
694 {
695 // base quad
696 case 0:
697 {
698 numModes0 = nummodes[0];
699 numModes1 = nummodes[1];
700 }
701 break;
702 // front and back quad
703 case 2:
704 case 4:
705 {
706 numModes0 = nummodes[1];
707 numModes1 = nummodes[2];
708 }
709 break;
710 // triangles
711 case 1:
712 case 3:
713 {
714 numModes0 = nummodes[0];
715 numModes1 = nummodes[2];
716 }
717 break;
718 }
719
720 if (faceOrient >= eDir1FwdDir2_Dir2FwdDir1)
721 {
722 std::swap(numModes0, numModes1);
723 }
724}

References Nektar::StdRegions::eDir1FwdDir2_Dir2FwdDir1, and Nektar::StdRegions::StdExpansion::m_base.

◆ v_GetTraceNumPoints()

int Nektar::StdRegions::StdPrismExp::v_GetTraceNumPoints ( const int  i) const
overrideprotectedvirtual

Implements Nektar::StdRegions::StdExpansion.

Definition at line 852 of file StdPrismExp.cpp.

853{
854 ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
855
856 if (i == 0)
857 {
858 return m_base[0]->GetNumPoints() * m_base[1]->GetNumPoints();
859 }
860 else if (i == 1 || i == 3)
861 {
862 return m_base[0]->GetNumPoints() * m_base[2]->GetNumPoints();
863 }
864 else
865 {
866 return m_base[1]->GetNumPoints() * m_base[2]->GetNumPoints();
867 }
868}

References ASSERTL2, and Nektar::StdRegions::StdExpansion::m_base.

◆ v_GetTracePointsKey()

LibUtilities::PointsKey Nektar::StdRegions::StdPrismExp::v_GetTracePointsKey ( const int  i,
const int  j 
) const
overrideprotectedvirtual

Reimplemented from Nektar::StdRegions::StdExpansion.

Definition at line 870 of file StdPrismExp.cpp.

872{
873 ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
874 ASSERTL2(j == 0 || j == 1, "face direction is out of range");
875
876 if (i == 0)
877 {
878 return m_base[j]->GetPointsKey();
879 }
880 else if (i == 1 || i == 3)
881 {
882 return m_base[2 * j]->GetPointsKey();
883 }
884 else
885 {
886 return m_base[j + 1]->GetPointsKey();
887 }
888}

References ASSERTL2, and Nektar::StdRegions::StdExpansion::m_base.

◆ v_GetVertexMap()

int Nektar::StdRegions::StdPrismExp::v_GetVertexMap ( int  localVertexId,
bool  useCoeffPacking = false 
)
overrideprotectedvirtual

Reimplemented from Nektar::StdRegions::StdExpansion.

Definition at line 942 of file StdPrismExp.cpp.

943{
947 "Mapping not defined for this type of basis");
948
949 int l = 0;
950
951 if (useCoeffPacking == true) // follow packing of coefficients i.e q,r,p
952 {
953 switch (vId)
954 {
955 case 0:
956 l = GetMode(0, 0, 0);
957 break;
958 case 1:
959 l = GetMode(0, 0, 1);
960 break;
961 case 2:
962 l = GetMode(0, 1, 0);
963 break;
964 case 3:
965 l = GetMode(0, 1, 1);
966 break;
967 case 4:
968 l = GetMode(1, 0, 0);
969 break;
970 case 5:
971 l = GetMode(1, 1, 0);
972 break;
973 default:
974 ASSERTL0(false, "local vertex id must be between 0 and 5");
975 }
976 }
977 else
978 {
979 switch (vId)
980 {
981 case 0:
982 l = GetMode(0, 0, 0);
983 break;
984 case 1:
985 l = GetMode(1, 0, 0);
986 break;
987 case 2:
988 l = GetMode(1, 1, 0);
989 break;
990 case 3:
991 l = GetMode(0, 1, 0);
992 break;
993 case 4:
994 l = GetMode(0, 0, 1);
995 break;
996 case 5:
997 l = GetMode(0, 1, 1);
998 break;
999 default:
1000 ASSERTL0(false, "local vertex id must be between 0 and 5");
1001 }
1002 }
1003
1004 return l;
1005}

References ASSERTL0, ASSERTL1, Nektar::LibUtilities::eModified_A, Nektar::LibUtilities::eModified_B, Nektar::StdRegions::StdExpansion::GetBasisType(), and GetMode().

◆ v_IProductWRTBaseKernel()

void Nektar::StdRegions::StdPrismExp::v_IProductWRTBaseKernel ( const Array< OneD, const NekDouble > &  base0,
const Array< OneD, const NekDouble > &  base1,
const Array< OneD, const NekDouble > &  base2,
const Array< OneD, const NekDouble > &  inarray,
Array< OneD, NekDouble > &  outarray,
const Array< OneD, NekDouble > &  jac,
const bool  Deformed,
bool  CollDir0 = false,
bool  CollDir1 = false,
bool  CollDir2 = false 
)
overrideprotectedvirtual

Inner product of inarray over region with respect to the expansion basis (this)->m_base[0] and return in outarray.

Parameters
base0- An array containing the values of the basis in the 0-direction at the quarature poitns
base1- An array containing the values of the basis in the 1-direction at the quarature poitns
base2- An array containing the values of the basis in the 2-direction at the quarature poitns
inarray- Array of values evaluated at the physical quadrature points
outarraythe values of the inner product with respect to each basis over region will be stored in the array outarray as output of the function
jac- An array of size 1 if not deformed or the number of quadrature points if deformed holding the values of the jacobian
Deformed- a bool identifying if the inner product is to be treated as a deformed or regular integration which just relates to how the
jacarray is treated

Implements Nektar::StdRegions::StdExpansion3D.

Definition at line 298 of file StdPrismExp.cpp.

306{
309 "Basis[1] is not a general tensor type");
310
313 "Basis[2] is not a general tensor type");
314
315 int nquad0 = m_base[0]->GetNumPoints();
316 int nquad1 = m_base[1]->GetNumPoints();
317 int nquad2 = m_base[2]->GetNumPoints();
318
319 int order0 = m_base[0]->GetNumModes();
320 int order1 = m_base[1]->GetNumModes();
321 int order2 = m_base[2]->GetNumModes();
322
323 const bool isModified =
324 (m_base[0]->GetBasisType() == LibUtilities::eModified_A);
325
326 std::vector<vec_t, tinysimd::allocator<vec_t>> wsp0(nquad1 * nquad2),
327 wsp1(nquad2), wsp2(order1);
328
329 // Swith statment using boost_pp and macros. This unfolls intwo a
330 // nested swtich statement where the outer swtich statement runs
331 // from SMIN to SMAX for modal order and the inner switch
332 // statemets run from the outer value of the case to 2*SMAX for
333 // the quadrature order. If you want to see it unwrapped compile
334 // in verbose mode and add --preprocess to the c++ command.
335 if (Deformed)
336 {
337 // Default case
338#undef IPRODUCTWRTBASE_DEF
339#define IPRODUCTWRTBASE_DEF \
340 IProductPrismKernel<false, false, true>( \
341 order0, order1, order2, nquad0, nquad1, nquad2, isModified, \
342 (const vec_t *)inarray.data(), (const vec_t *)base0.data(), \
343 (const vec_t *)base1.data(), (const vec_t *)base2.data(), \
344 (const vec_t *)m_weights[0].data(), \
345 (const vec_t *)m_weights[1].data(), \
346 (const vec_t *)m_weights[2].data(), (const vec_t *)jac.data(), \
347 (vec_t *)wsp0.data(), (vec_t *)wsp1.data(), (vec_t *)wsp2.data(), \
348 (vec_t *)outarray.data())
349
350 // Inner loop case over quarature points
351#undef IPRODUCTWRTBASE_Q
352#define IPRODUCTWRTBASE_Q(r, i) \
353 case NQ(i): \
354 IProductPrismKernel<false, false, true>( \
355 NM(i), NM(i), NM(i), NQ(i), NQ(i), NQ_M1(i), isModified, \
356 (const vec_t *)inarray.data(), (const vec_t *)base0.data(), \
357 (const vec_t *)base1.data(), (const vec_t *)base2.data(), \
358 (const vec_t *)m_weights[0].data(), \
359 (const vec_t *)m_weights[1].data(), \
360 (const vec_t *)m_weights[2].data(), (const vec_t *)jac.data(), \
361 (vec_t *)wsp0.data(), (vec_t *)wsp1.data(), (vec_t *)wsp2.data(), \
362 (vec_t *)outarray.data()); \
363 break;
364
365 // outer loop case over modes
366#undef IPRODUCTWRTBASE_M
367#define IPRODUCTWRTBASE_M(r, i) \
368 case NM(i): \
369 { \
370 switch (nquad0) \
371 { \
372 BOOST_PP_FOR_##r((NM(i), NM_P1(i), BOOST_PP_MUL(2, NM(i))), \
373 STDLEV2TEST1, STDLEV2UPDATE1, \
374 IPRODUCTWRTBASE_Q) default : IPRODUCTWRTBASE_DEF; \
375 break; \
376 } \
377 } \
378 break;
379
380 // templated cases on equi-ordered modes and standard quad usage
381 // where quad order goes from mode order to 2(*mode order)
382 if ((order0 == order1) && (order1 == order2) && (nquad0 == nquad1) &&
383 (nquad1 == nquad2 + 1))
384 {
385 switch (order0)
386 {
387 BOOST_PP_FOR((SMIN, 0, SMAX), STDLEV2TEST, STDLEV2UPDATE,
389 default:
391 break;
392 }
393 }
394 else
395 {
397 }
398 }
399 else // non-deformed case
400 {
401 // Default case
402#undef IPRODUCTWRTBASE_DEF
403#define IPRODUCTWRTBASE_DEF \
404 IProductPrismKernel<false, false, false>( \
405 order0, order1, order2, nquad0, nquad1, nquad2, isModified, \
406 (const vec_t *)inarray.data(), (const vec_t *)base0.data(), \
407 (const vec_t *)base1.data(), (const vec_t *)base2.data(), \
408 (const vec_t *)m_weights[0].data(), \
409 (const vec_t *)m_weights[1].data(), \
410 (const vec_t *)m_weights[2].data(), (const vec_t *)jac.data(), \
411 (vec_t *)wsp0.data(), (vec_t *)wsp1.data(), (vec_t *)wsp2.data(), \
412 (vec_t *)outarray.data())
413
414 // Inner loop case over quarature points
415#undef IPRODUCTWRTBASE_Q
416#define IPRODUCTWRTBASE_Q(r, i) \
417 case NQ(i): \
418 IProductPrismKernel<false, false, false>( \
419 NM(i), NM(i), NM(i), NQ(i), NQ(i), NQ_M1(i), isModified, \
420 (const vec_t *)inarray.data(), (const vec_t *)base0.data(), \
421 (const vec_t *)base1.data(), (const vec_t *)base2.data(), \
422 (const vec_t *)m_weights[0].data(), \
423 (const vec_t *)m_weights[1].data(), \
424 (const vec_t *)m_weights[2].data(), (const vec_t *)jac.data(), \
425 (vec_t *)wsp0.data(), (vec_t *)wsp1.data(), (vec_t *)wsp2.data(), \
426 (vec_t *)outarray.data()); \
427 break;
428
429 // outer loop case over modes
430#undef IPRODUCTWRTBASE_M
431#define IPRODUCTWRTBASE_M(r, i) \
432 case NM(i): \
433 { \
434 switch (nquad0) \
435 { \
436 BOOST_PP_FOR_##r((NM(i), NM_P1(i), BOOST_PP_MUL(2, NM(i))), \
437 STDLEV2TEST1, STDLEV2UPDATE1, \
438 IPRODUCTWRTBASE_Q) default : IPRODUCTWRTBASE_DEF; \
439 break; \
440 } \
441 } \
442 break;
443
444 // templated cases on equi-ordered modes and standard quad usage
445 // where quad order goes from mode order to 2(*mode order)
446 if ((order0 == order1) && (order1 == order2) && (nquad0 == nquad1) &&
447 (nquad1 == nquad2 + 1))
448 {
449 switch (order0)
450 {
451 BOOST_PP_FOR((SMIN, 0, SMAX), STDLEV2TEST, STDLEV2UPDATE,
453 default:
455 break;
456 }
457 }
458 else
459 {
461 }
462 }
463}
#define IPRODUCTWRTBASE_DEF
#define IPRODUCTWRTBASE_M(r, i)

References ASSERTL1, Nektar::LibUtilities::eModified_A, Nektar::LibUtilities::eModified_B, Nektar::LibUtilities::eModified_C, Nektar::LibUtilities::eOrtho_B, Nektar::LibUtilities::eOrtho_C, Nektar::StdRegions::StdExpansion::GetBasisType(), IPRODUCTWRTBASE_DEF, IPRODUCTWRTBASE_M, Nektar::StdRegions::StdExpansion::m_base, STDLEV2TEST, and STDLEV2UPDATE.

Referenced by Nektar::LocalRegions::PrismExp::v_IProductWRTDerivBase(), v_IProductWRTDerivBase(), and Nektar::LocalRegions::PrismExp::v_LaplacianMatrixOp_MatFree_Kernel().

◆ v_IProductWRTDerivBase()

void Nektar::StdRegions::StdPrismExp::v_IProductWRTDerivBase ( const int  dir,
const Array< OneD, const NekDouble > &  inarray,
Array< OneD, NekDouble > &  outarray 
)
overrideprotectedvirtual

Inner product of inarray over region with respect to the object's default expansion basis; output in outarray.

Reimplemented from Nektar::StdRegions::StdExpansion.

Definition at line 469 of file StdPrismExp.cpp.

472{
473 ASSERTL0(dir >= 0 && dir <= 2, "input dir is out of range");
474
475 int i;
476 int nquad0 = m_base[0]->GetNumPoints();
477 int nquad1 = m_base[1]->GetNumPoints();
478 int nquad2 = m_base[2]->GetNumPoints();
479
480 Array<OneD, NekDouble> tmp0(nquad0 * nquad1 * nquad2);
481
482 StdFacKey fackey2(eTwoOverOneMinusZ2, m_base[2]->GetBasisKey());
483 Array<OneD, const NekDouble> gfac2 = GetStdFac(fackey2);
484
485 const Array<OneD, const NekDouble> one(1, 1.0);
486
487 // Scale first derivative term by gfac2.
488 if (dir != 1)
489 {
490 for (i = 0; i < nquad2; ++i)
491 {
492 Vmath::Smul(nquad0 * nquad1, gfac2[i],
493 &inarray[0] + i * nquad0 * nquad1, 1,
494 &tmp0[0] + i * nquad0 * nquad1, 1);
495 }
496 }
497
498 switch (dir)
499 {
500 case 0:
501 {
503 m_base[0]->GetDbdata(), m_base[1]->GetBdata(),
504 m_base[2]->GetBdata(), tmp0, outarray, one, false);
505 }
506 break;
507 case 1:
508 {
510 m_base[0]->GetBdata(), m_base[1]->GetDbdata(),
511 m_base[2]->GetBdata(), inarray, outarray, one, false);
512 }
513 break;
514 case 2:
515 {
516 Array<OneD, NekDouble> tmp1(m_ncoeffs);
517
518 StdFacKey fackey0(eHalfMultOnePlusZ0, m_base[0]->GetBasisKey());
519 Array<OneD, const NekDouble> gfac0 = GetStdFac(fackey0);
520
521 // Scale eta_1 derivative with gfac0.
522 for (i = 0; i < nquad1 * nquad2; ++i)
523 {
524 Vmath::Vmul(nquad0, &gfac0[0], 1, &tmp0[0] + i * nquad0, 1,
525 &tmp0[0] + i * nquad0, 1);
526 }
527
529 m_base[0]->GetDbdata(), m_base[1]->GetBdata(),
530 m_base[2]->GetBdata(), tmp0, tmp1, one, false);
531
533 m_base[0]->GetBdata(), m_base[1]->GetBdata(),
534 m_base[2]->GetDbdata(), inarray, outarray, one, false);
535
536 Vmath::Vadd(m_ncoeffs, &tmp1[0], 1, &outarray[0], 1, &outarray[0],
537 1);
538 break;
539 }
540 }
541}
void v_IProductWRTBaseKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const Array< OneD, NekDouble > &jac, const bool Deformed, bool CollDir0=false, bool CollDir1=false, bool CollDir2=false) override
Inner product of inarray over region with respect to the expansion basis (this)->m_base[0] and return...
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.hpp:72
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition Vmath.hpp:180

References ASSERTL0, Nektar::StdRegions::eHalfMultOnePlusZ0, Nektar::StdRegions::eTwoOverOneMinusZ2, Nektar::StdRegions::StdExpansion::GetStdFac(), Nektar::StdRegions::StdExpansion::m_base, Nektar::StdRegions::StdExpansion::m_ncoeffs, Vmath::Smul(), v_IProductWRTBaseKernel(), Vmath::Vadd(), and Vmath::Vmul().

Referenced by Nektar::StdRegions::StdNodalPrismExp::v_IProductWRTDerivBase().

◆ v_IsBoundaryInteriorExpansion()

bool Nektar::StdRegions::StdPrismExp::v_IsBoundaryInteriorExpansion ( ) const
overrideprotectedvirtual

◆ v_LocCollapsedToLocCoord()

void Nektar::StdRegions::StdPrismExp::v_LocCollapsedToLocCoord ( const Array< OneD, const NekDouble > &  eta,
Array< OneD, NekDouble > &  xi 
)
overrideprotectedvirtual

Reimplemented from Nektar::StdRegions::StdExpansion.

Definition at line 567 of file StdPrismExp.cpp.

569{
570 xi[0] = (1.0 + eta[0]) * (1.0 - eta[2]) * 0.5 - 1.0;
571 xi[1] = eta[1];
572 xi[2] = eta[2];
573}

◆ v_LocCoordToLocCollapsed()

void Nektar::StdRegions::StdPrismExp::v_LocCoordToLocCollapsed ( const Array< OneD, const NekDouble > &  xi,
Array< OneD, NekDouble > &  eta 
)
overrideprotectedvirtual

Reimplemented from Nektar::StdRegions::StdExpansion.

Definition at line 547 of file StdPrismExp.cpp.

549{
550 NekDouble d2 = 1.0 - xi[2];
551 if (fabs(d2) < NekConstants::kNekZeroTol)
552 {
553 if (d2 >= 0.)
554 {
556 }
557 else
558 {
560 }
561 }
562 eta[2] = xi[2]; // eta_z = xi_z
563 eta[1] = xi[1]; // eta_y = xi_y
564 eta[0] = 2.0 * (1.0 + xi[0]) / d2 - 1.0;
565}
static const NekDouble kNekZeroTol

References Nektar::NekConstants::kNekZeroTol.

◆ v_NumBndryCoeffs()

int Nektar::StdRegions::StdPrismExp::v_NumBndryCoeffs ( ) const
overrideprotectedvirtual

Implements Nektar::StdRegions::StdExpansion.

Definition at line 772 of file StdPrismExp.cpp.

773{
776 "BasisType is not a boundary interior form");
779 "BasisType is not a boundary interior form");
782 "BasisType is not a boundary interior form");
783
784 int P = m_base[0]->GetNumModes();
785 int Q = m_base[1]->GetNumModes();
786 int R = m_base[2]->GetNumModes();
787
789}
constexpr int getNumberOfBndCoefficients(int Na, int Nb, int Nc)

References ASSERTL1, Nektar::LibUtilities::eGLL_Lagrange, Nektar::LibUtilities::eModified_A, Nektar::LibUtilities::eModified_B, Nektar::StdRegions::StdExpansion::GetBasisType(), Nektar::LibUtilities::StdPrismData::getNumberOfBndCoefficients(), Nektar::StdRegions::StdExpansion::m_base, and Nektar::LibUtilities::P.

◆ v_NumDGBndryCoeffs()

int Nektar::StdRegions::StdPrismExp::v_NumDGBndryCoeffs ( ) const
overrideprotectedvirtual

Implements Nektar::StdRegions::StdExpansion.

Definition at line 791 of file StdPrismExp.cpp.

792{
795 "BasisType is not a boundary interior form");
798 "BasisType is not a boundary interior form");
801 "BasisType is not a boundary interior form");
802
803 int P = m_base[0]->GetNumModes() - 1;
804 int Q = m_base[1]->GetNumModes() - 1;
805 int R = m_base[2]->GetNumModes() - 1;
806
807 return (P + 1) * (Q + 1) // 1 rect. face on base
808 + 2 * (Q + 1) * (R + 1) // other 2 rect. faces
809 + 2 * (R + 1) + P * (1 + 2 * R - P); // 2 tri. faces
810}

References ASSERTL1, Nektar::LibUtilities::eGLL_Lagrange, Nektar::LibUtilities::eModified_A, Nektar::LibUtilities::eModified_B, Nektar::StdRegions::StdExpansion::GetBasisType(), Nektar::StdRegions::StdExpansion::m_base, and Nektar::LibUtilities::P.

◆ v_PhysEvalFirstDeriv()

NekDouble Nektar::StdRegions::StdPrismExp::v_PhysEvalFirstDeriv ( const Array< OneD, NekDouble > &  coord,
const Array< OneD, const NekDouble > &  inarray,
std::array< NekDouble, 3 > &  firstOrderDerivs 
)
overrideprotectedvirtual

Reimplemented from Nektar::StdRegions::StdExpansion.

Definition at line 633 of file StdPrismExp.cpp.

637{
638 // Collapse coordinates
639 Array<OneD, NekDouble> coll(3, 0.0);
640 LocCoordToLocCollapsed(coord, coll);
641
642 // If near singularity do the old interpolation matrix method
643 // @TODO: Dave thinks there might be a way in the Barycentric to
644 // mathematically remove this singularity?
645 if ((1 - coll[2]) < 1e-5)
646 {
647 int totPoints = GetTotPoints();
648 Array<OneD, NekDouble> EphysDeriv0(totPoints), EphysDeriv1(totPoints),
649 EphysDeriv2(totPoints);
650 v_PhysDeriv(inarray, EphysDeriv0, EphysDeriv1, EphysDeriv2);
651
652 Array<OneD, DNekMatSharedPtr> I(3);
653 I[0] = GetBase()[0]->GetI(coll);
654 I[1] = GetBase()[1]->GetI(coll + 1);
655 I[2] = GetBase()[2]->GetI(coll + 2);
656
657 firstOrderDerivs[0] = PhysEvaluate(I, EphysDeriv0);
658 firstOrderDerivs[1] = PhysEvaluate(I, EphysDeriv1);
659 firstOrderDerivs[2] = PhysEvaluate(I, EphysDeriv2);
660 return PhysEvaluate(I, inarray);
661 }
662
663 NekDouble val = BaryTensorDeriv(coll, inarray, firstOrderDerivs);
664
665 NekDouble dEta_bar1 = firstOrderDerivs[0];
666
667 NekDouble fac = 2.0 / (1.0 - coll[2]);
668 firstOrderDerivs[0] = fac * dEta_bar1;
669
670 // divide dEta_Bar1 by (1-eta_z)
671 fac = 1.0 / (1.0 - coll[2]);
672 dEta_bar1 = fac * dEta_bar1;
673
674 // Multiply dEta_Bar1 by (1+eta_x) and add ot out_dxi3
675 fac = 1.0 + coll[0];
676 firstOrderDerivs[2] += fac * dEta_bar1;
677
678 return val;
679}
NekDouble BaryTensorDeriv(const Array< OneD, NekDouble > &coord, const Array< OneD, const NekDouble > &inarray, std::array< NekDouble, 3 > &firstOrderDerivs)
void v_PhysDeriv(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Calculate the derivative of the physical points in a given direction.
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
const Array< OneD, const LibUtilities::BasisSharedPtr > & GetBase() const
This function gets the shared point to basis.
NekDouble PhysEvaluate(const Array< OneD, const NekDouble > &coords, const Array< OneD, const NekDouble > &physvals)
This function evaluates the expansion at a single (arbitrary) point of the domain.

References Nektar::StdRegions::StdExpansion3D::BaryTensorDeriv(), Nektar::StdRegions::StdExpansion::GetBase(), Nektar::StdRegions::StdExpansion::GetTotPoints(), Nektar::StdRegions::StdExpansion::LocCoordToLocCollapsed(), Nektar::StdRegions::StdExpansion::PhysEvaluate(), and Nektar::StdRegions::StdExpansion3D::v_PhysDeriv().

◆ v_PhysEvaluateBasis()

NekDouble Nektar::StdRegions::StdPrismExp::v_PhysEvaluateBasis ( const Array< OneD, const NekDouble > &  coords,
int  mode 
)
finalprotectedvirtual

Reimplemented from Nektar::StdRegions::StdExpansion.

Definition at line 602 of file StdPrismExp.cpp.

604{
605 Array<OneD, NekDouble> coll(3);
606 LocCoordToLocCollapsed(coords, coll);
607
608 const int nm1 = m_base[1]->GetNumModes();
609 const int nm2 = m_base[2]->GetNumModes();
610 const int b = 2 * nm2 + 1;
611
612 const int mode0 = floor(0.5 * (b - sqrt(b * b - 8.0 * mode / nm1)));
613 const int tmp =
614 mode - nm1 * (mode0 * (nm2 - 1) + 1 - (mode0 - 2) * (mode0 - 1) / 2);
615 const int mode1 = tmp / (nm2 - mode0);
616 const int mode2 = tmp % (nm2 - mode0);
617
618 if (mode0 == 0 && mode2 == 1 &&
620 {
621 // handle collapsed top edge to remove mode0 terms
622 return StdExpansion::BaryEvaluateBasis<1>(coll[1], mode1) *
623 StdExpansion::BaryEvaluateBasis<2>(coll[2], mode2);
624 }
625 else
626 {
627 return StdExpansion::BaryEvaluateBasis<0>(coll[0], mode0) *
628 StdExpansion::BaryEvaluateBasis<1>(coll[1], mode1) *
629 StdExpansion::BaryEvaluateBasis<2>(coll[2], mode2);
630 }
631}
scalarT< T > sqrt(scalarT< T > in)
Definition scalar.hpp:290

References Nektar::LibUtilities::eModified_A, Nektar::StdRegions::StdExpansion::GetBasisType(), Nektar::StdRegions::StdExpansion::LocCoordToLocCollapsed(), Nektar::StdRegions::StdExpansion::m_base, and tinysimd::sqrt().

◆ v_ReduceOrderCoeffs()

void Nektar::StdRegions::StdPrismExp::v_ReduceOrderCoeffs ( int  numMin,
const Array< OneD, const NekDouble > &  inarray,
Array< OneD, NekDouble > &  outarray 
)
overrideprotectedvirtual

Reimplemented from Nektar::StdRegions::StdExpansion.

Definition at line 2096 of file StdPrismExp.cpp.

2099{
2100 int nquad0 = m_base[0]->GetNumPoints();
2101 int nquad1 = m_base[1]->GetNumPoints();
2102 int nquad2 = m_base[2]->GetNumPoints();
2103 int nqtot = nquad0 * nquad1 * nquad2;
2104 int nmodes0 = m_base[0]->GetNumModes();
2105 int nmodes1 = m_base[1]->GetNumModes();
2106 int nmodes2 = m_base[2]->GetNumModes();
2107 int numMax = nmodes0;
2108
2109 Array<OneD, NekDouble> coeff(m_ncoeffs);
2110 Array<OneD, NekDouble> coeff_tmp1(m_ncoeffs, 0.0);
2111 Array<OneD, NekDouble> phys_tmp(nqtot, 0.0);
2112 Array<OneD, NekDouble> tmp, tmp2, tmp3, tmp4;
2113
2114 const LibUtilities::PointsKey Pkey0 = m_base[0]->GetPointsKey();
2115 const LibUtilities::PointsKey Pkey1 = m_base[1]->GetPointsKey();
2116 const LibUtilities::PointsKey Pkey2 = m_base[2]->GetPointsKey();
2117
2118 LibUtilities::BasisKey bortho0(LibUtilities::eOrtho_A, nmodes0, Pkey0);
2119 LibUtilities::BasisKey bortho1(LibUtilities::eOrtho_A, nmodes1, Pkey1);
2120 LibUtilities::BasisKey bortho2(LibUtilities::eOrtho_B, nmodes2, Pkey2);
2121
2122 int cnt = 0;
2123 int u = 0;
2124 int i = 0;
2126
2128 bortho0, bortho1, bortho2);
2129
2130 BwdTrans(inarray, phys_tmp);
2131 OrthoPrismExp->FwdTrans(phys_tmp, coeff);
2132
2133 // filtering
2134 for (u = 0; u < numMin; ++u)
2135 {
2136 for (i = 0; i < numMin; ++i)
2137 {
2138 Vmath::Vcopy(numMin - u, tmp = coeff + cnt, 1,
2139 tmp2 = coeff_tmp1 + cnt, 1);
2140 cnt += numMax - u;
2141 }
2142
2143 for (i = numMin; i < numMax; ++i)
2144 {
2145 cnt += numMax - u;
2146 }
2147 }
2148
2149 OrthoPrismExp->BwdTrans(coeff_tmp1, phys_tmp);
2150 StdPrismExp::FwdTrans(phys_tmp, outarray);
2151}
void BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Backward transformation from coefficient space to physical space.
void FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
@ eOrtho_A
Principle Orthogonal Functions .
Definition BasisType.h:42
std::shared_ptr< StdPrismExp > StdPrismExpSharedPtr

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::StdRegions::StdExpansion::BwdTrans(), Nektar::LibUtilities::eOrtho_A, Nektar::LibUtilities::eOrtho_B, Nektar::StdRegions::StdExpansion::FwdTrans(), Nektar::StdRegions::StdExpansion::m_base, Nektar::StdRegions::StdExpansion::m_ncoeffs, and Vmath::Vcopy().

◆ v_StdPhysDeriv()

void Nektar::StdRegions::StdPrismExp::v_StdPhysDeriv ( const Array< OneD, const NekDouble > &  u_physical,
Array< OneD, NekDouble > &  out_dxi1,
Array< OneD, NekDouble > &  out_dxi2,
Array< OneD, NekDouble > &  out_dxi3 
)
overrideprotectedvirtual

Calculate the derivative of the physical points.

The derivative is evaluated at the nodal physical points. Derivatives with respect to the local Cartesian coordinates.

\(\begin{Bmatrix} \frac {\partial} {\partial \xi_1} \\ \frac {\partial} {\partial \xi_2} \\ \frac {\partial} {\partial \xi_3} \end{Bmatrix} = \begin{Bmatrix} \frac 2 {(1-\eta_3)} \frac \partial {\partial \bar \eta_1} \\ \frac {\partial} {\partial \xi_2} \ \ \frac {(1 + \bar \eta_1)} {(1 - \eta_3)} \frac \partial {\partial \bar \eta_1} + \frac {\partial} {\partial \eta_3} \end{Bmatrix}\)

Reimplemented from Nektar::StdRegions::StdExpansion.

Definition at line 91 of file StdPrismExp.cpp.

95{
96 int Qx = m_base[0]->GetNumPoints();
97 int Qy = m_base[1]->GetNumPoints();
98 int Qz = m_base[2]->GetNumPoints();
99 int Qtot = Qx * Qy * Qz;
100
101 Array<OneD, NekDouble> dEta_bar1(Qtot, 0.0);
102
103 Array<OneD, const NekDouble> eta_x, eta_z;
104 eta_x = m_base[0]->GetZ();
105 eta_z = m_base[2]->GetZ();
106
107 int i, k;
108
109 bool Do_1 = (out_dxi1.size() > 0) ? true : false;
110 bool Do_3 = (out_dxi3.size() > 0) ? true : false;
111
112 // out_dXi2 is just a tensor derivative so is just passed through
113 if (Do_3)
114 {
115 PhysTensorDeriv(u_physical, dEta_bar1, out_dxi2, out_dxi3);
116 }
117 else if (Do_1)
118 {
119 PhysTensorDeriv(u_physical, dEta_bar1, out_dxi2, NullNekDouble1DArray);
120 }
121 else // case if just require 2nd direction
122 {
123 PhysTensorDeriv(u_physical, NullNekDouble1DArray, out_dxi2,
125 }
126
127 if (Do_1)
128 {
129 for (k = 0; k < Qz; ++k)
130 {
131 Vmath::Smul(Qx * Qy, 2.0 / (1.0 - eta_z[k]),
132 &dEta_bar1[0] + k * Qx * Qy, 1,
133 &out_dxi1[0] + k * Qx * Qy, 1);
134 }
135 }
136
137 if (Do_3)
138 {
139 // divide dEta_Bar1 by (1-eta_z)
140 for (k = 0; k < Qz; ++k)
141 {
142 Vmath::Smul(Qx * Qy, 1.0 / (1.0 - eta_z[k]),
143 &dEta_bar1[0] + k * Qx * Qy, 1,
144 &dEta_bar1[0] + k * Qx * Qy, 1);
145 }
146
147 // Multiply dEta_Bar1 by (1+eta_x) and add ot out_dxi3
148 for (i = 0; i < Qx; ++i)
149 {
150 Vmath::Svtvp(Qz * Qy, 1.0 + eta_x[i], &dEta_bar1[0] + i, Qx,
151 &out_dxi3[0] + i, Qx, &out_dxi3[0] + i, Qx);
152 }
153 }
154}
void PhysTensorDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2)
Calculate the 3D derivative in the local tensor/collapsed coordinate at the physical points.
static Array< OneD, NekDouble > NullNekDouble1DArray
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Svtvp (scalar times vector plus vector): z = alpha*x + y.
Definition Vmath.hpp:396

References Nektar::StdRegions::StdExpansion::m_base, Nektar::NullNekDouble1DArray, Nektar::StdRegions::StdExpansion3D::PhysTensorDeriv(), Vmath::Smul(), and Vmath::Svtvp().

◆ v_SVVLaplacianFilter()

void Nektar::StdRegions::StdPrismExp::v_SVVLaplacianFilter ( Array< OneD, NekDouble > &  array,
const StdMatrixKey mkey 
)
overrideprotectedvirtual

Reimplemented from Nektar::StdRegions::StdExpansion.

Definition at line 1958 of file StdPrismExp.cpp.

1960{
1961 // Generate an orthonogal expansion
1962 int qa = m_base[0]->GetNumPoints();
1963 int qb = m_base[1]->GetNumPoints();
1964 int qc = m_base[2]->GetNumPoints();
1965 int nmodes_a = m_base[0]->GetNumModes();
1966 int nmodes_b = m_base[1]->GetNumModes();
1967 int nmodes_c = m_base[2]->GetNumModes();
1968 // Declare orthogonal basis.
1969 LibUtilities::PointsKey pa(qa, m_base[0]->GetPointsType());
1970 LibUtilities::PointsKey pb(qb, m_base[1]->GetPointsType());
1971 LibUtilities::PointsKey pc(qc, m_base[2]->GetPointsType());
1972
1973 LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A, nmodes_a, pa);
1974 LibUtilities::BasisKey Bb(LibUtilities::eOrtho_A, nmodes_b, pb);
1975 LibUtilities::BasisKey Bc(LibUtilities::eOrtho_B, nmodes_c, pc);
1976 StdPrismExp OrthoExp(Ba, Bb, Bc);
1977
1978 Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
1979 int i, j, k, cnt = 0;
1980
1981 // project onto modal space.
1982 OrthoExp.FwdTrans(array, orthocoeffs);
1983
1984 if (mkey.ConstFactorExists(eFactorSVVPowerKerDiffCoeff))
1985 {
1986 // Rodrigo's power kernel
1987 NekDouble cutoff = mkey.GetConstFactor(eFactorSVVCutoffRatio);
1988 NekDouble SvvDiffCoeff =
1989 mkey.GetConstFactor(eFactorSVVPowerKerDiffCoeff) *
1990 mkey.GetConstFactor(eFactorSVVDiffCoeff);
1991
1992 for (i = 0; i < nmodes_a; ++i)
1993 {
1994 for (j = 0; j < nmodes_b; ++j)
1995 {
1996 NekDouble fac1 = std::max(
1997 pow((1.0 * i) / (nmodes_a - 1), cutoff * nmodes_a),
1998 pow((1.0 * j) / (nmodes_b - 1), cutoff * nmodes_b));
1999
2000 for (k = 0; k < nmodes_c - i; ++k)
2001 {
2002 NekDouble fac =
2003 std::max(fac1, pow((1.0 * k) / (nmodes_c - 1),
2004 cutoff * nmodes_c));
2005
2006 orthocoeffs[cnt] *= SvvDiffCoeff * fac;
2007 cnt++;
2008 }
2009 }
2010 }
2011 }
2012 else if (mkey.ConstFactorExists(
2013 eFactorSVVDGKerDiffCoeff)) // Rodrigo/Mansoor's DG Kernel
2014 {
2015 NekDouble SvvDiffCoeff = mkey.GetConstFactor(eFactorSVVDGKerDiffCoeff) *
2016 mkey.GetConstFactor(eFactorSVVDiffCoeff);
2017
2018 int max_abc = max(nmodes_a - kSVVDGFiltermodesmin,
2019 nmodes_b - kSVVDGFiltermodesmin);
2020 max_abc = max(max_abc, nmodes_c - kSVVDGFiltermodesmin);
2021 // clamp max_abc
2022 max_abc = max(max_abc, 0);
2023 max_abc = min(max_abc, kSVVDGFiltermodesmax - kSVVDGFiltermodesmin);
2024
2025 for (i = 0; i < nmodes_a; ++i)
2026 {
2027 for (j = 0; j < nmodes_b; ++j)
2028 {
2029 int maxij = max(i, j);
2030
2031 for (k = 0; k < nmodes_c - i; ++k)
2032 {
2033 int maxijk = max(maxij, k);
2034 maxijk = min(maxijk, kSVVDGFiltermodesmax - 1);
2035
2036 orthocoeffs[cnt] *=
2037 SvvDiffCoeff * kSVVDGFilter[max_abc][maxijk];
2038 cnt++;
2039 }
2040 }
2041 }
2042 }
2043 else
2044 {
2045 // SVV filter paramaters (how much added diffusion relative
2046 // to physical one and fraction of modes from which you
2047 // start applying this added diffusion)
2048 //
2049 NekDouble SvvDiffCoeff =
2050 mkey.GetConstFactor(StdRegions::eFactorSVVDiffCoeff);
2051 NekDouble SVVCutOff =
2052 mkey.GetConstFactor(StdRegions::eFactorSVVCutoffRatio);
2053
2054 // Defining the cut of mode
2055 int cutoff_a = (int)(SVVCutOff * nmodes_a);
2056 int cutoff_b = (int)(SVVCutOff * nmodes_b);
2057 int cutoff_c = (int)(SVVCutOff * nmodes_c);
2058 // To avoid the fac[j] from blowing up
2059 NekDouble epsilon = 1;
2060
2061 int nmodes = min(min(nmodes_a, nmodes_b), nmodes_c);
2062 NekDouble cutoff = min(min(cutoff_a, cutoff_b), cutoff_c);
2063
2064 //------"New" Version August 22nd '13--------------------
2065 for (i = 0; i < nmodes_a; ++i) // P
2066 {
2067 for (j = 0; j < nmodes_b; ++j) // Q
2068 {
2069 for (k = 0; k < nmodes_c - i; ++k) // R
2070 {
2071 if (j >= cutoff || i + k >= cutoff)
2072 {
2073 orthocoeffs[cnt] *=
2074 (SvvDiffCoeff *
2075 exp(-(i + k - nmodes) * (i + k - nmodes) /
2076 ((NekDouble)((i + k - cutoff + epsilon) *
2077 (i + k - cutoff + epsilon)))) *
2078 exp(-(j - nmodes) * (j - nmodes) /
2079 ((NekDouble)((j - cutoff + epsilon) *
2080 (j - cutoff + epsilon)))));
2081 }
2082 else
2083 {
2084 orthocoeffs[cnt] *= 0.0;
2085 }
2086 cnt++;
2087 }
2088 }
2089 }
2090 }
2091
2092 // backward transform to physical space
2093 OrthoExp.BwdTrans(orthocoeffs, array);
2094}
LibUtilities::PointsType GetPointsType(const int dir) const
This function returns the type of quadrature points used in the dir direction.
const int kSVVDGFiltermodesmin
const int kSVVDGFiltermodesmax
const NekDouble kSVVDGFilter[9][11]

References Nektar::StdRegions::StdExpansion::BwdTrans(), Nektar::StdRegions::StdMatrixKey::ConstFactorExists(), Nektar::StdRegions::eFactorSVVCutoffRatio, Nektar::StdRegions::eFactorSVVDGKerDiffCoeff, Nektar::StdRegions::eFactorSVVDiffCoeff, Nektar::StdRegions::eFactorSVVPowerKerDiffCoeff, Nektar::LibUtilities::eOrtho_A, Nektar::LibUtilities::eOrtho_B, Nektar::StdRegions::StdExpansion::FwdTrans(), Nektar::StdRegions::StdMatrixKey::GetConstFactor(), Nektar::StdRegions::StdExpansion::GetNcoeffs(), Nektar::StdRegions::StdExpansion::GetPointsType(), Nektar::StdRegions::kSVVDGFilter, Nektar::StdRegions::kSVVDGFiltermodesmax, Nektar::StdRegions::kSVVDGFiltermodesmin, Nektar::StdRegions::StdExpansion::m_base, tinysimd::max(), and tinysimd::min().