Nektar++
Public Member Functions | Protected Member Functions | Protected Attributes | Private Member Functions | Private Attributes | List of all members
Nektar::LocalRegions::Expansion3D Class Reference

#include <Expansion3D.h>

Inheritance diagram for Nektar::LocalRegions::Expansion3D:
[legend]

Public Member Functions

 Expansion3D (SpatialDomains::Geometry3DSharedPtr pGeom)
 
 ~Expansion3D () override=default
 
void SetTraceToGeomOrientation (Array< OneD, NekDouble > &inout)
 Align trace orientation with the geometry orientation. More...
 
void SetFaceToGeomOrientation (const int face, Array< OneD, NekDouble > &inout)
 Align face orientation with the geometry orientation. More...
 
void AddHDGHelmholtzFaceTerms (const NekDouble tau, const int edge, Array< OneD, NekDouble > &facePhys, const StdRegions::VarCoeffMap &dirForcing, Array< OneD, NekDouble > &outarray)
 
void AddNormTraceInt (const int dir, Array< OneD, ExpansionSharedPtr > &FaceExp, Array< OneD, Array< OneD, NekDouble > > &faceCoeffs, Array< OneD, NekDouble > &outarray)
 
void AddNormTraceInt (const int dir, Array< OneD, const NekDouble > &inarray, Array< OneD, ExpansionSharedPtr > &FaceExp, Array< OneD, NekDouble > &outarray, const StdRegions::VarCoeffMap &varcoeffs)
 
void AddFaceBoundaryInt (const int face, ExpansionSharedPtr &FaceExp, Array< OneD, NekDouble > &facePhys, Array< OneD, NekDouble > &outarray, const StdRegions::VarCoeffMap &varcoeffs=StdRegions::NullVarCoeffMap)
 
SpatialDomains::Geometry3DSharedPtr GetGeom3D () const
 
void v_ReOrientTracePhysMap (const StdRegions::Orientation orient, Array< OneD, int > &idmap, const int nq0, const int nq1) override
 
void v_NormVectorIProductWRTBase (const Array< OneD, const Array< OneD, NekDouble > > &Fvec, Array< OneD, NekDouble > &outarray) override
 
Array< OneD, unsigned int > GetEdgeInverseBoundaryMap (int eid)
 
Array< OneD, unsigned int > GetTraceInverseBoundaryMap (int fid, StdRegions::Orientation faceOrient=StdRegions::eNoOrientation, int P1=-1, int P2=-1)
 
void GetInverseBoundaryMaps (Array< OneD, unsigned int > &vmap, Array< OneD, Array< OneD, unsigned int > > &emap, Array< OneD, Array< OneD, unsigned int > > &fmap)
 
DNekScalMatSharedPtr CreateMatrix (const MatrixKey &mkey)
 
- Public Member Functions inherited from Nektar::LocalRegions::Expansion
 Expansion (SpatialDomains::GeometrySharedPtr pGeom)
 
 Expansion (const Expansion &pSrc)
 
 ~Expansion () override
 
void SetTraceExp (const int traceid, ExpansionSharedPtr &f)
 
ExpansionSharedPtr GetTraceExp (const int traceid)
 
DNekScalMatSharedPtr GetLocMatrix (const LocalRegions::MatrixKey &mkey)
 
void DropLocMatrix (const LocalRegions::MatrixKey &mkey)
 
DNekScalMatSharedPtr GetLocMatrix (const StdRegions::MatrixType mtype, const StdRegions::ConstFactorMap &factors=StdRegions::NullConstFactorMap, const StdRegions::VarCoeffMap &varcoeffs=StdRegions::NullVarCoeffMap)
 
SpatialDomains::GeometrySharedPtr GetGeom () const
 
void Reset ()
 
IndexMapValuesSharedPtr CreateIndexMap (const IndexMapKey &ikey)
 
DNekScalBlkMatSharedPtr CreateStaticCondMatrix (const MatrixKey &mkey)
 
const SpatialDomains::GeomFactorsSharedPtrGetMetricInfo () const
 
DNekMatSharedPtr BuildTransformationMatrix (const DNekScalMatSharedPtr &r_bnd, const StdRegions::MatrixType matrixType)
 
DNekMatSharedPtr BuildVertexMatrix (const DNekScalMatSharedPtr &r_bnd)
 
void ExtractDataToCoeffs (const NekDouble *data, const std::vector< unsigned int > &nummodes, const int nmodes_offset, NekDouble *coeffs, std::vector< LibUtilities::BasisType > &fromType)
 
void AddEdgeNormBoundaryInt (const int edge, const std::shared_ptr< Expansion > &EdgeExp, const Array< OneD, const NekDouble > &Fx, const Array< OneD, const NekDouble > &Fy, Array< OneD, NekDouble > &outarray)
 
void AddEdgeNormBoundaryInt (const int edge, const std::shared_ptr< Expansion > &EdgeExp, const Array< OneD, const NekDouble > &Fn, Array< OneD, NekDouble > &outarray)
 
void AddFaceNormBoundaryInt (const int face, const std::shared_ptr< Expansion > &FaceExp, const Array< OneD, const NekDouble > &Fn, Array< OneD, NekDouble > &outarray)
 
void DGDeriv (const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, ExpansionSharedPtr > &EdgeExp, Array< OneD, Array< OneD, NekDouble > > &coeffs, Array< OneD, NekDouble > &outarray)
 
NekDouble VectorFlux (const Array< OneD, Array< OneD, NekDouble > > &vec)
 
void NormalTraceDerivFactors (Array< OneD, Array< OneD, NekDouble > > &factors, Array< OneD, Array< OneD, NekDouble > > &d0factors, Array< OneD, Array< OneD, NekDouble > > &d1factors)
 
IndexMapValuesSharedPtr GetIndexMap (const IndexMapKey &ikey)
 
void AlignVectorToCollapsedDir (const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
 
ExpansionSharedPtr GetLeftAdjacentElementExp () const
 
ExpansionSharedPtr GetRightAdjacentElementExp () const
 
int GetLeftAdjacentElementTrace () const
 
int GetRightAdjacentElementTrace () const
 
void SetAdjacentElementExp (int traceid, ExpansionSharedPtr &e)
 
StdRegions::Orientation GetTraceOrient (int trace)
 
void SetCoeffsToOrientation (StdRegions::Orientation dir, Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
void DivideByQuadratureMetric (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 Divided by the metric jacobi and quadrature weights. More...
 
void GetTraceQFactors (const int trace, Array< OneD, NekDouble > &outarray)
 Extract the metric factors to compute the contravariant fluxes along edge edge and stores them into outarray following the local edge orientation (i.e. anticlockwise convention). More...
 
void GetTracePhysVals (const int trace, const StdRegions::StdExpansionSharedPtr &TraceExp, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, StdRegions::Orientation orient=StdRegions::eNoOrientation)
 
void GetTracePhysMap (const int edge, Array< OneD, int > &outarray)
 
void ReOrientTracePhysMap (const StdRegions::Orientation orient, Array< OneD, int > &idmap, const int nq0, const int nq1)
 
const NormalVectorGetTraceNormal (const int id)
 
void ComputeTraceNormal (const int id)
 
const Array< OneD, const NekDouble > & GetPhysNormals (void)
 
void SetPhysNormals (Array< OneD, const NekDouble > &normal)
 
void SetUpPhysNormals (const int trace)
 
void AddRobinMassMatrix (const int traceid, const Array< OneD, const NekDouble > &primCoeffs, DNekMatSharedPtr &inoutmat)
 
void TraceNormLen (const int traceid, NekDouble &h, NekDouble &p)
 
void AddRobinTraceContribution (const int traceid, const Array< OneD, const NekDouble > &primCoeffs, const Array< OneD, NekDouble > &incoeffs, Array< OneD, NekDouble > &coeffs)
 
const Array< OneD, const NekDouble > & GetElmtBndNormDirElmtLen (const int nbnd) const
 
void StdDerivBaseOnTraceMat (Array< OneD, DNekMatSharedPtr > &DerivMat)
 
- Public Member Functions inherited from Nektar::StdRegions::StdExpansion
 StdExpansion ()
 Default Constructor. More...
 
 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. More...
 
 StdExpansion (const StdExpansion &T)
 Copy Constructor. More...
 
virtual ~StdExpansion ()
 Destructor. More...
 
int GetNumBases () const
 This function returns the number of 1D bases used in the expansion. More...
 
const Array< OneD, const LibUtilities::BasisSharedPtr > & GetBase () const
 This function gets the shared point to basis. More...
 
const LibUtilities::BasisSharedPtrGetBasis (int dir) const
 This function gets the shared point to basis in the dir direction. More...
 
int GetNcoeffs (void) const
 This function returns the total number of coefficients used in the expansion. More...
 
int GetTotPoints () const
 This function returns the total number of quadrature points used in the element. More...
 
LibUtilities::BasisType GetBasisType (const int dir) const
 This function returns the type of basis used in the dir direction. More...
 
int GetBasisNumModes (const int dir) const
 This function returns the number of expansion modes in the dir direction. More...
 
int EvalBasisNumModesMax (void) const
 This function returns the maximum number of expansion modes over all local directions. More...
 
LibUtilities::PointsType GetPointsType (const int dir) const
 This function returns the type of quadrature points used in the dir direction. More...
 
int GetNumPoints (const int dir) const
 This function returns the number of quadrature points in the dir direction. More...
 
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. More...
 
int GetNverts () const
 This function returns the number of vertices of the expansion domain. More...
 
int GetTraceNcoeffs (const int i) const
 This function returns the number of expansion coefficients belonging to the i-th trace. More...
 
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. More...
 
const LibUtilities::BasisKey GetTraceBasisKey (const int i, int k=-1) const
 This function returns the basis key belonging to the i-th trace. More...
 
LibUtilities::PointsKey GetTracePointsKey (const int i, int k=-1) const
 This function returns the basis key belonging to the i-th trace. More...
 
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. More...
 
int GetNtraces () const
 Returns the number of trace elements connected to this element. More...
 
LibUtilities::ShapeType DetShapeType () const
 This function returns the shape of the expansion domain. More...
 
std::shared_ptr< StdExpansionGetStdExp () const
 
std::shared_ptr< StdExpansionGetLinStdExp (void) const
 
int GetShapeDimension () const
 
bool IsBoundaryInteriorExpansion () const
 
bool IsNodalNonTensorialExp ()
 
void BwdTrans (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 This function performs the Backward transformation from coefficient space to physical space. More...
 
void FwdTrans (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 This function performs the Forward transformation from physical space to coefficient space. More...
 
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. More...
 
void FillMode (const int mode, Array< OneD, NekDouble > &outarray)
 This function fills the array outarray with the mode-th mode of the expansion. More...
 
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 More...
 
void IProductWRTBase (const Array< OneD, const NekDouble > &base, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, int coll_check)
 
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. More...
 
void SetElmtId (const int id)
 Set the element id of this expansion when used in a list by returning value of m_elmt_id. More...
 
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 More...
 
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 More...
 
DNekMatSharedPtr GetStdMatrix (const StdMatrixKey &mkey)
 
DNekBlkMatSharedPtr GetStdStaticCondMatrix (const StdMatrixKey &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}\) More...
 
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 PhysDeriv_s (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_ds)
 
void PhysDeriv_n (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_dn)
 
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)
 
void StdPhysDeriv (const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
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. More...
 
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. More...
 
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. More...
 
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. More...
 
void LocCoordToLocCollapsed (const Array< OneD, const NekDouble > &xi, Array< OneD, NekDouble > &eta)
 Convert local cartesian coordinate xi into local collapsed coordinates eta. More...
 
void LocCollapsedToLocCoord (const Array< OneD, const NekDouble > &eta, Array< OneD, NekDouble > &xi)
 Convert local collapsed coordinates eta into local cartesian coordinate xi. More...
 
virtual int v_CalcNumberOfCoefficients (const std::vector< unsigned int > &nummodes, int &modes_offset)
 
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. More...
 
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. More...
 
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. More...
 
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. More...
 
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. More...
 
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. More...
 
template<class T >
std::shared_ptr< T > as ()
 
void IProductWRTBase_SumFac (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
 
void GenStdMatBwdDeriv (const int dir, DNekMatSharedPtr &mat)
 
- 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 PhysTensorDeriv (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray_d1, Array< OneD, NekDouble > &outarray_d2, Array< OneD, NekDouble > &outarray_d3)
 Calculate the 3D derivative in the local tensor/collapsed coordinate at the physical points. More...
 
void BwdTrans_SumFacKernel (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, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0, bool doCheckCollDir1, bool doCheckCollDir2)
 
void IProductWRTBase_SumFacKernel (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, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0, bool doCheckCollDir1, bool doCheckCollDir2)
 
int GetNedges () const
 return the number of edges in 3D expansion More...
 
int GetEdgeNcoeffs (const int i) const
 This function returns the number of expansion coefficients belonging to the i-th edge. More...
 
void GetEdgeInteriorToElementMap (const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, Orientation traceOrient=eForwards)
 

Protected Member Functions

void v_DGDeriv (const int dir, const Array< OneD, const NekDouble > &incoeffs, Array< OneD, ExpansionSharedPtr > &FaceExp, Array< OneD, Array< OneD, NekDouble > > &faceCoeffs, Array< OneD, NekDouble > &out_d) override
 Evaluate coefficients of weak deriviative in the direction dir given the input coefficicents incoeffs and the imposed boundary values in EdgeExp (which will have its phys space updated). More...
 
DNekMatSharedPtr v_GenMatrix (const StdRegions::StdMatrixKey &mkey) override
 
void v_AddFaceNormBoundaryInt (const int face, const ExpansionSharedPtr &FaceExp, const Array< OneD, const NekDouble > &Fn, Array< OneD, NekDouble > &outarray) override
 
void v_AddRobinMassMatrix (const int face, const Array< OneD, const NekDouble > &primCoeffs, DNekMatSharedPtr &inoutmat) override
 
StdRegions::Orientation v_GetTraceOrient (int face) override
 
void v_GetTracePhysVals (const int face, const StdRegions::StdExpansionSharedPtr &FaceExp, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, StdRegions::Orientation orient) override
 Extract the physical values along face face from inarray into outarray following the local face orientation and point distribution defined by defined in FaceExp. More...
 
void v_GenTraceExp (const int traceid, ExpansionSharedPtr &exp) override
 
void GetPhysFaceVarCoeffsFromElement (const int face, ExpansionSharedPtr &FaceExp, const Array< OneD, const NekDouble > &varcoeff, Array< OneD, NekDouble > &outarray)
 
DNekMatSharedPtr v_BuildTransformationMatrix (const DNekScalMatSharedPtr &r_bnd, const StdRegions::MatrixType matrixType) override
 
DNekMatSharedPtr v_BuildInverseTransformationMatrix (const DNekScalMatSharedPtr &transformationmatrix) override
 Build inverse and inverse transposed transformation matrix: \(\mathbf{R^{-1}}\) and \(\mathbf{R^{-T}}\). More...
 
DNekMatSharedPtr v_BuildVertexMatrix (const DNekScalMatSharedPtr &r_bnd) override
 
void v_TraceNormLen (const int traceid, NekDouble &h, NekDouble &p) override
 
- Protected Member Functions inherited from Nektar::LocalRegions::Expansion
void ComputeLaplacianMetric ()
 
void ComputeQuadratureMetric ()
 
void ComputeGmatcdotMF (const Array< TwoD, const NekDouble > &df, const Array< OneD, const NekDouble > &direction, Array< OneD, Array< OneD, NekDouble > > &dfdir)
 
Array< OneD, NekDoubleGetMF (const int dir, const int shapedim, const StdRegions::VarCoeffMap &varcoeffs)
 
Array< OneD, NekDoubleGetMFDiv (const int dir, const StdRegions::VarCoeffMap &varcoeffs)
 
Array< OneD, NekDoubleGetMFMag (const int dir, const StdRegions::VarCoeffMap &varcoeffs)
 
void v_MultiplyByQuadratureMetric (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
 
virtual void v_DivideByQuadratureMetric (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
virtual void v_ComputeLaplacianMetric ()
 
int v_GetCoordim () const override
 
void v_GetCoords (Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3) override
 
virtual DNekScalMatSharedPtr v_GetLocMatrix (const LocalRegions::MatrixKey &mkey)
 
virtual void v_DropLocMatrix (const LocalRegions::MatrixKey &mkey)
 
virtual DNekMatSharedPtr v_BuildTransformationMatrix (const DNekScalMatSharedPtr &r_bnd, const StdRegions::MatrixType matrixType)
 
virtual DNekMatSharedPtr v_BuildVertexMatrix (const DNekScalMatSharedPtr &r_bnd)
 
virtual void v_ExtractDataToCoeffs (const NekDouble *data, const std::vector< unsigned int > &nummodes, const int nmodes_offset, NekDouble *coeffs, std::vector< LibUtilities::BasisType > &fromType)
 
virtual void v_AddEdgeNormBoundaryInt (const int edge, const std::shared_ptr< Expansion > &EdgeExp, const Array< OneD, const NekDouble > &Fx, const Array< OneD, const NekDouble > &Fy, Array< OneD, NekDouble > &outarray)
 
virtual void v_AddEdgeNormBoundaryInt (const int edge, const std::shared_ptr< Expansion > &EdgeExp, const Array< OneD, const NekDouble > &Fn, Array< OneD, NekDouble > &outarray)
 
virtual void v_AddFaceNormBoundaryInt (const int face, const std::shared_ptr< Expansion > &FaceExp, const Array< OneD, const NekDouble > &Fn, Array< OneD, NekDouble > &outarray)
 
virtual void v_DGDeriv (const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, ExpansionSharedPtr > &EdgeExp, Array< OneD, Array< OneD, NekDouble > > &coeffs, Array< OneD, NekDouble > &outarray)
 
virtual NekDouble v_VectorFlux (const Array< OneD, Array< OneD, NekDouble > > &vec)
 
virtual void v_NormalTraceDerivFactors (Array< OneD, Array< OneD, NekDouble > > &factors, Array< OneD, Array< OneD, NekDouble > > &d0factors, Array< OneD, Array< OneD, NekDouble > > &d1factors)
 
virtual void v_AlignVectorToCollapsedDir (const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
 
virtual StdRegions::Orientation v_GetTraceOrient (int trace)
 
void v_SetCoeffsToOrientation (StdRegions::Orientation dir, Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
 
virtual void v_GetTraceQFactors (const int trace, Array< OneD, NekDouble > &outarray)
 
virtual void v_GetTracePhysVals (const int trace, const StdRegions::StdExpansionSharedPtr &TraceExp, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, StdRegions::Orientation orient)
 
virtual void v_GetTracePhysMap (const int edge, Array< OneD, int > &outarray)
 
virtual void v_ReOrientTracePhysMap (const StdRegions::Orientation orient, Array< OneD, int > &idmap, const int nq0, const int nq1=-1)
 
virtual void v_ComputeTraceNormal (const int id)
 
virtual const Array< OneD, const NekDouble > & v_GetPhysNormals ()
 
virtual void v_SetPhysNormals (Array< OneD, const NekDouble > &normal)
 
virtual void v_SetUpPhysNormals (const int id)
 
virtual void v_AddRobinMassMatrix (const int face, const Array< OneD, const NekDouble > &primCoeffs, DNekMatSharedPtr &inoutmat)
 
virtual void v_AddRobinTraceContribution (const int traceid, const Array< OneD, const NekDouble > &primCoeffs, const Array< OneD, NekDouble > &incoeffs, Array< OneD, NekDouble > &coeffs)
 
virtual void v_TraceNormLen (const int traceid, NekDouble &h, NekDouble &p)
 
virtual void v_GenTraceExp (const int traceid, ExpansionSharedPtr &exp)
 
- Protected Member Functions inherited from Nektar::StdRegions::StdExpansion
DNekMatSharedPtr CreateStdMatrix (const StdMatrixKey &mkey)
 
DNekBlkMatSharedPtr CreateStdStaticCondMatrix (const StdMatrixKey &mkey)
 Create the static condensation of a matrix when using a boundary interior decomposition. More...
 
void BwdTrans_SumFac (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
void IProductWRTDerivBase_SumFac (const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
void IProductWRTDirectionalDerivBase_SumFac (const Array< OneD, const NekDouble > &direction, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
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)
 
virtual void v_SetCoeffsToOrientation (Array< OneD, NekDouble > &coeffs, StdRegions::Orientation dir)
 
virtual NekDouble v_StdPhysEvaluate (const Array< OneD, const NekDouble > &Lcoord, const Array< OneD, const NekDouble > &physvals)
 
virtual void v_GenStdMatBwdDeriv (const int dir, DNekMatSharedPtr &mat)
 
virtual void v_MultiplyByStdQuadratureMetric (const 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. More...
 
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. More...
 
template<int DIR, bool DERIV = false, bool DERIV2 = false>
NekDouble BaryEvaluate (const NekDouble &coord, const NekDouble *physvals, NekDouble &deriv)
 
- Protected Member Functions inherited from Nektar::StdRegions::StdExpansion3D
NekDouble v_PhysEvaluate (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. More...
 
NekDouble v_PhysEvaluate (const Array< OneD, DNekMatSharedPtr > &I, const Array< OneD, const NekDouble > &physvals) override
 
NekDouble v_PhysEvaluate (const Array< OneD, NekDouble > &coord, const Array< OneD, const NekDouble > &inarray, std::array< NekDouble, 3 > &firstOrderDerivs) override
 
virtual void v_BwdTrans_SumFacKernel (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, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0, bool doCheckCollDir1, bool doCheckCollDir2)=0
 
virtual void v_IProductWRTBase_SumFacKernel (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, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0, bool doCheckCollDir1, bool doCheckCollDir2)=0
 
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 v_Integral (const Array< OneD, const NekDouble > &inarray) override
 Integrates the specified function over the domain. More...
 
virtual int v_GetNedges (void) const
 
virtual int v_GetEdgeNcoeffs (const int i) const
 
NekDouble BaryTensorDeriv (const Array< OneD, NekDouble > &coord, const Array< OneD, const NekDouble > &inarray, std::array< NekDouble, 3 > &firstOrderDerivs)
 
virtual void v_GetEdgeInteriorToElementMap (const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, Orientation traceOrient=eForwards)
 
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
 

Protected Attributes

std::map< int, NormalVectorm_faceNormals
 
- Protected Attributes inherited from Nektar::LocalRegions::Expansion
LibUtilities::NekManager< IndexMapKey, IndexMapValues, IndexMapKey::opLessm_indexMapManager
 
std::map< int, ExpansionWeakPtrm_traceExp
 
SpatialDomains::GeometrySharedPtr m_geom
 
SpatialDomains::GeomFactorsSharedPtr m_metricinfo
 
MetricMap m_metrics
 
std::map< int, NormalVectorm_traceNormals
 
ExpansionWeakPtr m_elementLeft
 
ExpansionWeakPtr m_elementRight
 
int m_elementTraceLeft = -1
 
int m_elementTraceRight = -1
 
std::map< int, Array< OneD, NekDouble > > m_elmtBndNormDirElmtLen
 the element length in each element boundary(Vertex, edge or face) normal direction calculated based on the local m_metricinfo times the standard element length (which is 2.0) More...
 
- Protected Attributes inherited from Nektar::StdRegions::StdExpansion
Array< OneD, LibUtilities::BasisSharedPtrm_base
 
int m_elmt_id
 
int m_ncoeffs
 
LibUtilities::NekManager< StdMatrixKey, DNekMat, StdMatrixKey::opLessm_stdMatrixManager
 
LibUtilities::NekManager< StdMatrixKey, DNekBlkMat, StdMatrixKey::opLessm_stdStaticCondMatrixManager
 

Private Member Functions

Array< OneD, NekDoubleGetnFacecdotMF (const int dir, const int face, ExpansionSharedPtr &FaceExp_f, const Array< OneD, const Array< OneD, NekDouble > > &normals, const StdRegions::VarCoeffMap &varcoeffs)
 

Private Attributes

std::vector< bool > m_requireNeg
 

Detailed Description

Definition at line 55 of file Expansion3D.h.

Constructor & Destructor Documentation

◆ Expansion3D()

Nektar::LocalRegions::Expansion3D::Expansion3D ( SpatialDomains::Geometry3DSharedPtr  pGeom)
inline

Definition at line 59 of file Expansion3D.h.

61 {
62 }
std::vector< bool > m_requireNeg
Definition: Expansion3D.h:168
Expansion(SpatialDomains::GeometrySharedPtr pGeom)
Definition: Expansion.cpp:43

◆ ~Expansion3D()

Nektar::LocalRegions::Expansion3D::~Expansion3D ( )
overridedefault

Member Function Documentation

◆ AddFaceBoundaryInt()

void Nektar::LocalRegions::Expansion3D::AddFaceBoundaryInt ( const int  face,
ExpansionSharedPtr FaceExp,
Array< OneD, NekDouble > &  facePhys,
Array< OneD, NekDouble > &  outarray,
const StdRegions::VarCoeffMap varcoeffs = StdRegions::NullVarCoeffMap 
)
inline

For a given face add the \tilde{F}_1j contributions

Definition at line 322 of file Expansion3D.cpp.

326{
327 int i;
328 int order_f = FaceExp->GetNcoeffs();
329 Array<OneD, NekDouble> coeff(order_f);
330
331 IndexMapKey ikey(eFaceToElement, DetShapeType(), GetBasisNumModes(0),
333 GetTraceOrient(face));
335
336 FaceExp->IProductWRTBase(facePhys, coeff);
337
338 // add data to out array
339 for (i = 0; i < order_f; ++i)
340 {
341 outarray[(*map)[i].index] += (*map)[i].sign * coeff[i];
342 }
343}
StdRegions::Orientation GetTraceOrient(int trace)
Definition: Expansion.h:168
IndexMapValuesSharedPtr GetIndexMap(const IndexMapKey &ikey)
Definition: Expansion.h:146
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:367
int GetBasisNumModes(const int dir) const
This function returns the number of expansion modes in the dir direction.
Definition: StdExpansion.h:169
std::shared_ptr< IndexMapValues > IndexMapValuesSharedPtr
Definition: IndexMapKey.h:126

References Nektar::StdRegions::StdExpansion::DetShapeType(), Nektar::LocalRegions::eFaceToElement, Nektar::StdRegions::StdExpansion::GetBasisNumModes(), Nektar::LocalRegions::Expansion::GetIndexMap(), and Nektar::LocalRegions::Expansion::GetTraceOrient().

Referenced by AddNormTraceInt().

◆ AddHDGHelmholtzFaceTerms()

void Nektar::LocalRegions::Expansion3D::AddHDGHelmholtzFaceTerms ( const NekDouble  tau,
const int  edge,
Array< OneD, NekDouble > &  facePhys,
const StdRegions::VarCoeffMap dirForcing,
Array< OneD, NekDouble > &  outarray 
)
inline

Definition at line 50 of file Expansion3D.cpp.

53{
54 ExpansionSharedPtr FaceExp = GetTraceExp(face);
55 int i, j, n;
56 int nquad_f = FaceExp->GetNumPoints(0) * FaceExp->GetNumPoints(1);
57 int order_f = FaceExp->GetNcoeffs();
58 int coordim = GetCoordim();
59 int ncoeffs = GetNcoeffs();
60 bool mmf = (varcoeffs.find(StdRegions::eVarCoeffMF1x) != varcoeffs.end());
61
62 Array<OneD, NekDouble> inval(nquad_f);
63 Array<OneD, NekDouble> outcoeff(order_f);
64 Array<OneD, NekDouble> tmpcoeff(ncoeffs);
65
66 const Array<OneD, const Array<OneD, NekDouble>> &normals =
67 GetTraceNormal(face);
68
70
71 DNekVec Coeffs(ncoeffs, outarray, eWrapper);
72 DNekVec Tmpcoeff(ncoeffs, tmpcoeff, eWrapper);
73
74 IndexMapKey ikey(eFaceToElement, DetShapeType(), GetBasisNumModes(0),
76 GetTraceOrient(face));
78
82
83 // @TODO Variable coefficients
84 /*
85 StdRegions::VarCoeffType VarCoeff[3] = {StdRegions::eVarCoeffD00,
86 StdRegions::eVarCoeffD11,
87 StdRegions::eVarCoeffD22};
88 Array<OneD, NekDouble> varcoeff_work(nquad_f);
89 StdRegions::VarCoeffMap::const_iterator x;
90 ///// @TODO: What direction to use here??
91 if ((x = varcoeffs.find(VarCoeff[0])) != varcoeffs.end())
92 {
93 GetPhysFaceVarCoeffsFromElement(face,FaceExp,x->second,varcoeff_work);
94 Vmath::Vmul(nquad_f,varcoeff_work,1,FaceExp->GetPhys(),1,FaceExp->UpdatePhys(),1);
95 }
96 */
97
98 //================================================================
99 // Add F = \tau <phi_i,in_phys>
100 // Fill face and take inner product
101 FaceExp->IProductWRTBase(facePhys, outcoeff);
102
103 for (i = 0; i < order_f; ++i)
104 {
105 outarray[(*map)[i].index] += (*map)[i].sign * tau * outcoeff[i];
106 }
107 //================================================================
108
109 //===============================================================
110 // Add -\sum_i D_i^T M^{-1} G_i + E_i M^{-1} G_i =
111 // \sum_i D_i M^{-1} G_i term
112
113 // Three independent direction
114 for (n = 0; n < coordim; ++n)
115 {
116 if (mmf)
117 {
119 Weight[StdRegions::eVarCoeffMass] = GetMFMag(n, varcoeffs);
120
121 MatrixKey invMasskey(StdRegions::eInvMass, DetShapeType(), *this,
123
124 invMass = *GetLocMatrix(invMasskey);
125
126 Array<OneD, NekDouble> ncdotMF_f =
127 GetnFacecdotMF(n, face, FaceExp, normals, varcoeffs);
128
129 Vmath::Vmul(nquad_f, ncdotMF_f, 1, facePhys, 1, inval, 1);
130 }
131 else
132 {
133 Vmath::Vmul(nquad_f, normals[n], 1, facePhys, 1, inval, 1);
134 }
135
136 NekDouble scale = invMass.Scale();
137 const NekDouble *data = invMass.GetRawPtr();
138
139 // @TODO Multiply by variable coefficients
140 // @TODO: Document this (probably not needed)
141 /*
142 StdRegions::VarCoeffMap::const_iterator x;
143 if ((x = varcoeffs.find(VarCoeff[n])) != varcoeffs.end())
144 {
145 GetPhysEdgeVarCoeffsFromElement(edge,FaceExp,x->second,varcoeff_work);
146 Vmath::Vmul(nquad_f,varcoeff_work,1,FaceExp->GetPhys(),1,FaceExp->UpdatePhys(),1);
147 }
148 */
149
150 FaceExp->IProductWRTBase(inval, outcoeff);
151
152 // M^{-1} G
153 for (i = 0; i < ncoeffs; ++i)
154 {
155 tmpcoeff[i] = 0;
156 for (j = 0; j < order_f; ++j)
157 {
158 tmpcoeff[i] += scale * data[i + (*map)[j].index * ncoeffs] *
159 (*map)[j].sign * outcoeff[j];
160 }
161 }
162
163 if (mmf)
164 {
165 StdRegions::VarCoeffMap VarCoeffDirDeriv;
166 VarCoeffDirDeriv[StdRegions::eVarCoeffMF] =
167 GetMF(n, coordim, varcoeffs);
168 VarCoeffDirDeriv[StdRegions::eVarCoeffMFDiv] =
169 GetMFDiv(n, varcoeffs);
170
173 VarCoeffDirDeriv);
174
175 DNekScalMat &Dmat = *GetLocMatrix(Dmatkey);
176
177 Coeffs = Coeffs + Dmat * Tmpcoeff;
178 }
179 else
180 {
181 DNekScalMat &Dmat = *GetLocMatrix(DerivType[n]);
182 Coeffs = Coeffs + Dmat * Tmpcoeff;
183 }
184
185 /*
186 if(varcoeffs.find(VarCoeff[n]) != varcoeffs.end())
187 {
188 MatrixKey mkey(DerivType[n], DetExpansionType(), *this,
189 StdRegions::NullConstFactorMap, varcoeffs); DNekScalMat &Dmat =
190 *GetLocMatrix(mkey); Coeffs = Coeffs + Dmat*Tmpcoeff;
191 }
192
193 else
194 {
195 DNekScalMat &Dmat = *GetLocMatrix(DerivType[n]);
196 Coeffs = Coeffs + Dmat*Tmpcoeff;
197 }
198 */
199 }
200}
Array< OneD, NekDouble > GetnFacecdotMF(const int dir, const int face, ExpansionSharedPtr &FaceExp_f, const Array< OneD, const Array< OneD, NekDouble > > &normals, const StdRegions::VarCoeffMap &varcoeffs)
Array< OneD, NekDouble > GetMFMag(const int dir, const StdRegions::VarCoeffMap &varcoeffs)
Definition: Expansion.cpp:709
ExpansionSharedPtr GetTraceExp(const int traceid)
Definition: Expansion.h:414
DNekScalMatSharedPtr GetLocMatrix(const LocalRegions::MatrixKey &mkey)
Definition: Expansion.cpp:84
Array< OneD, NekDouble > GetMFDiv(const int dir, const StdRegions::VarCoeffMap &varcoeffs)
Definition: Expansion.cpp:686
const NormalVector & GetTraceNormal(const int id)
Definition: Expansion.cpp:251
Array< OneD, NekDouble > GetMF(const int dir, const int shapedim, const StdRegions::VarCoeffMap &varcoeffs)
Definition: Expansion.cpp:633
int GetNcoeffs(void) const
This function returns the total number of coefficients used in the expansion.
Definition: StdExpansion.h:124
std::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:66
static ConstFactorMap NullConstFactorMap
Definition: StdRegions.hpp:431
std::map< StdRegions::VarCoeffType, VarCoeffEntry > VarCoeffMap
Definition: StdRegions.hpp:375
NekMatrix< NekMatrix< NekDouble, StandardMatrixTag >, ScaledMatrixTag > DNekScalMat
NekVector< NekDouble > DNekVec
Definition: NekTypeDefs.hpp:48
double NekDouble
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

References Nektar::StdRegions::StdExpansion::DetShapeType(), Nektar::LocalRegions::eFaceToElement, Nektar::StdRegions::eInvMass, Nektar::StdRegions::eVarCoeffMass, Nektar::StdRegions::eVarCoeffMF, Nektar::StdRegions::eVarCoeffMF1x, Nektar::StdRegions::eVarCoeffMFDiv, Nektar::StdRegions::eWeakDeriv0, Nektar::StdRegions::eWeakDeriv1, Nektar::StdRegions::eWeakDeriv2, Nektar::StdRegions::eWeakDirectionalDeriv, Nektar::eWrapper, Nektar::StdRegions::StdExpansion::GetBasisNumModes(), Nektar::StdRegions::StdExpansion::GetCoordim(), Nektar::LocalRegions::Expansion::GetIndexMap(), Nektar::LocalRegions::Expansion::GetLocMatrix(), Nektar::LocalRegions::Expansion::GetMF(), Nektar::LocalRegions::Expansion::GetMFDiv(), Nektar::LocalRegions::Expansion::GetMFMag(), Nektar::StdRegions::StdExpansion::GetNcoeffs(), GetnFacecdotMF(), Nektar::LocalRegions::Expansion::GetTraceExp(), Nektar::LocalRegions::Expansion::GetTraceNormal(), Nektar::LocalRegions::Expansion::GetTraceOrient(), Nektar::StdRegions::NullConstFactorMap, and Vmath::Vmul().

Referenced by v_GenMatrix().

◆ AddNormTraceInt() [1/2]

void Nektar::LocalRegions::Expansion3D::AddNormTraceInt ( const int  dir,
Array< OneD, const NekDouble > &  inarray,
Array< OneD, ExpansionSharedPtr > &  FaceExp,
Array< OneD, NekDouble > &  outarray,
const StdRegions::VarCoeffMap varcoeffs 
)
inline

Computes the C matrix entries due to the presence of the identity matrix in Eqn. 32.

@TODO: Document this

Definition at line 232 of file Expansion3D.cpp.

237{
238 int i, f, cnt;
239 int order_f, nquad_f;
240 int nfaces = GetNtraces();
241
242 cnt = 0;
243 for (f = 0; f < nfaces; ++f)
244 {
245 order_f = FaceExp[f]->GetNcoeffs();
246 nquad_f = FaceExp[f]->GetNumPoints(0) * FaceExp[f]->GetNumPoints(1);
247
248 const Array<OneD, const Array<OneD, NekDouble>> &normals =
250 Array<OneD, NekDouble> faceCoeffs(order_f);
251 Array<OneD, NekDouble> facePhys(nquad_f);
252
253 for (i = 0; i < order_f; ++i)
254 {
255 faceCoeffs[i] = inarray[i + cnt];
256 }
257 cnt += order_f;
258
259 FaceExp[f]->BwdTrans(faceCoeffs, facePhys);
260
261 // Multiply by variable coefficient
262 /// @TODO: Document this
263 // StdRegions::VarCoeffType VarCoeff[3] =
264 // {StdRegions::eVarCoeffD00,
265 // StdRegions::eVarCoeffD11,
266 // StdRegions::eVarCoeffD22};
267 // StdRegions::VarCoeffMap::const_iterator x;
268 // Array<OneD, NekDouble> varcoeff_work(nquad_e);
269
270 // if ((x = varcoeffs.find(VarCoeff[dir])) !=
271 // varcoeffs.end())
272 // {
273 // GetPhysEdgeVarCoeffsFromElement(e,EdgeExp[e],x->second,varcoeff_work);
274 // Vmath::Vmul(nquad_e,varcoeff_work,1,EdgeExp[e]->GetPhys(),1,EdgeExp[e]->UpdatePhys(),1);
275 // }
276 StdRegions::VarCoeffMap::const_iterator x;
277 if ((x = varcoeffs.find(StdRegions::eVarCoeffMF1x)) != varcoeffs.end())
278 {
279 Array<OneD, NekDouble> ncdotMF_f =
280 GetnFacecdotMF(dir, f, FaceExp[f], normals, varcoeffs);
281
282 Vmath::Vmul(nquad_f, ncdotMF_f, 1, facePhys, 1, facePhys, 1);
283 }
284 else
285 {
286 Vmath::Vmul(nquad_f, normals[dir], 1, facePhys, 1, facePhys, 1);
287 }
288
289 AddFaceBoundaryInt(f, FaceExp[f], facePhys, outarray, varcoeffs);
290 }
291}
void AddFaceBoundaryInt(const int face, ExpansionSharedPtr &FaceExp, Array< OneD, NekDouble > &facePhys, Array< OneD, NekDouble > &outarray, const StdRegions::VarCoeffMap &varcoeffs=StdRegions::NullVarCoeffMap)
int GetNtraces() const
Returns the number of trace elements connected to this element.
Definition: StdExpansion.h:351

References AddFaceBoundaryInt(), Nektar::StdRegions::eVarCoeffMF1x, GetnFacecdotMF(), Nektar::StdRegions::StdExpansion::GetNtraces(), Nektar::LocalRegions::Expansion::GetTraceNormal(), and Vmath::Vmul().

◆ AddNormTraceInt() [2/2]

void Nektar::LocalRegions::Expansion3D::AddNormTraceInt ( const int  dir,
Array< OneD, ExpansionSharedPtr > &  FaceExp,
Array< OneD, Array< OneD, NekDouble > > &  faceCoeffs,
Array< OneD, NekDouble > &  outarray 
)
inline

Definition at line 294 of file Expansion3D.cpp.

298{
299 int f;
300 int nquad_f;
301 int nfaces = GetNtraces();
302
303 for (f = 0; f < nfaces; ++f)
304 {
305 nquad_f = FaceExp[f]->GetNumPoints(0) * FaceExp[f]->GetNumPoints(1);
306
307 const Array<OneD, const Array<OneD, NekDouble>> &normals =
309 Array<OneD, NekDouble> facePhys(nquad_f);
310
311 FaceExp[f]->BwdTrans(faceCoeffs[f], facePhys);
312
313 Vmath::Vmul(nquad_f, normals[dir], 1, facePhys, 1, facePhys, 1);
314
315 AddFaceBoundaryInt(f, FaceExp[f], facePhys, outarray);
316 }
317}

References AddFaceBoundaryInt(), Nektar::StdRegions::StdExpansion::GetNtraces(), Nektar::LocalRegions::Expansion::GetTraceNormal(), and Vmath::Vmul().

Referenced by v_DGDeriv(), and v_GenMatrix().

◆ CreateMatrix()

DNekScalMatSharedPtr Nektar::LocalRegions::Expansion3D::CreateMatrix ( const MatrixKey mkey)

Definition at line 407 of file Expansion3D.cpp.

408{
409 DNekScalMatSharedPtr returnval;
411
413 "Geometric information is not set up");
414
415 switch (mkey.GetMatrixType())
416 {
418 {
419 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
420 mkey.GetNVarCoeff())
421 {
422 NekDouble one = 1.0;
423 DNekMatSharedPtr mat = GenMatrix(mkey);
424 returnval =
426 }
427 else
428 {
429 NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
430 DNekMatSharedPtr mat = GetStdMatrix(mkey);
431 returnval =
433 }
434 }
435 break;
437 {
438 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
439 {
440 NekDouble one = 1.0;
441 StdRegions::StdMatrixKey masskey(StdRegions::eMass,
442 DetShapeType(), *this);
443 DNekMatSharedPtr mat = GenMatrix(masskey);
444 mat->Invert();
445 returnval =
447 }
448 else
449 {
450 NekDouble fac = 1.0 / (m_metricinfo->GetJac(ptsKeys))[0];
451 DNekMatSharedPtr mat = GetStdMatrix(mkey);
452 returnval =
454 }
455 }
456 break;
460 {
461 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
462 mkey.GetNVarCoeff())
463 {
464 NekDouble one = 1.0;
465 DNekMatSharedPtr mat = GenMatrix(mkey);
466
467 returnval =
469 }
470 else
471 {
472 NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
473 Array<TwoD, const NekDouble> df =
474 m_metricinfo->GetDerivFactors(ptsKeys);
475 int dir = 0;
476 if (mkey.GetMatrixType() == StdRegions::eWeakDeriv0)
477 {
478 dir = 0;
479 }
480 if (mkey.GetMatrixType() == StdRegions::eWeakDeriv1)
481 {
482 dir = 1;
483 }
484 if (mkey.GetMatrixType() == StdRegions::eWeakDeriv2)
485 {
486 dir = 2;
487 }
488
489 MatrixKey deriv0key(StdRegions::eWeakDeriv0,
490 mkey.GetShapeType(), *this);
491 MatrixKey deriv1key(StdRegions::eWeakDeriv1,
492 mkey.GetShapeType(), *this);
493 MatrixKey deriv2key(StdRegions::eWeakDeriv2,
494 mkey.GetShapeType(), *this);
495
496 DNekMat &deriv0 = *GetStdMatrix(deriv0key);
497 DNekMat &deriv1 = *GetStdMatrix(deriv1key);
498 DNekMat &deriv2 = *GetStdMatrix(deriv2key);
499
500 int rows = deriv0.GetRows();
501 int cols = deriv1.GetColumns();
502
503 DNekMatSharedPtr WeakDeriv =
505 (*WeakDeriv) = df[3 * dir][0] * deriv0 +
506 df[3 * dir + 1][0] * deriv1 +
507 df[3 * dir + 2][0] * deriv2;
508
510 jac, WeakDeriv);
511 }
512 }
513 break;
515 {
516 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
517 (mkey.HasVarCoeff(StdRegions::eVarCoeffLaplacian)) ||
518 (mkey.HasVarCoeff(StdRegions::eVarCoeffD00)) ||
519 (mkey.HasVarCoeff(StdRegions::eVarCoeffD01)) ||
520 (mkey.HasVarCoeff(StdRegions::eVarCoeffD10)) ||
521 (mkey.HasVarCoeff(StdRegions::eVarCoeffD02)) ||
522 (mkey.HasVarCoeff(StdRegions::eVarCoeffD20)) ||
523 (mkey.HasVarCoeff(StdRegions::eVarCoeffD11)) ||
524 (mkey.HasVarCoeff(StdRegions::eVarCoeffD12)) ||
525 (mkey.HasVarCoeff(StdRegions::eVarCoeffD21)) ||
526 (mkey.HasVarCoeff(StdRegions::eVarCoeffD22)) ||
527 (mkey.ConstFactorExists(StdRegions::eFactorSVVCutoffRatio)))
528 {
529 NekDouble one = 1.0;
530 DNekMatSharedPtr mat = GenMatrix(mkey);
531
532 returnval =
534 }
535 else
536 {
537 MatrixKey lap00key(StdRegions::eLaplacian00,
538 mkey.GetShapeType(), *this);
539 MatrixKey lap01key(StdRegions::eLaplacian01,
540 mkey.GetShapeType(), *this);
541 MatrixKey lap02key(StdRegions::eLaplacian02,
542 mkey.GetShapeType(), *this);
543 MatrixKey lap11key(StdRegions::eLaplacian11,
544 mkey.GetShapeType(), *this);
545 MatrixKey lap12key(StdRegions::eLaplacian12,
546 mkey.GetShapeType(), *this);
547 MatrixKey lap22key(StdRegions::eLaplacian22,
548 mkey.GetShapeType(), *this);
549
550 DNekMat &lap00 = *GetStdMatrix(lap00key);
551 DNekMat &lap01 = *GetStdMatrix(lap01key);
552 DNekMat &lap02 = *GetStdMatrix(lap02key);
553 DNekMat &lap11 = *GetStdMatrix(lap11key);
554 DNekMat &lap12 = *GetStdMatrix(lap12key);
555 DNekMat &lap22 = *GetStdMatrix(lap22key);
556
557 NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
558 Array<TwoD, const NekDouble> gmat =
559 m_metricinfo->GetGmat(ptsKeys);
560
561 int rows = lap00.GetRows();
562 int cols = lap00.GetColumns();
563
564 DNekMatSharedPtr lap =
566
567 (*lap) = gmat[0][0] * lap00 + gmat[4][0] * lap11 +
568 gmat[8][0] * lap22 +
569 gmat[3][0] * (lap01 + Transpose(lap01)) +
570 gmat[6][0] * (lap02 + Transpose(lap02)) +
571 gmat[7][0] * (lap12 + Transpose(lap12));
572
573 returnval =
575 }
576 }
577 break;
579 {
580 NekDouble factor = mkey.GetConstFactor(StdRegions::eFactorLambda);
581 MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this);
582 DNekScalMat &MassMat = *GetLocMatrix(masskey);
583 MatrixKey lapkey(StdRegions::eLaplacian, mkey.GetShapeType(), *this,
584 mkey.GetConstFactors(), mkey.GetVarCoeffs());
585 DNekScalMat &LapMat = *GetLocMatrix(lapkey);
586
587 int rows = LapMat.GetRows();
588 int cols = LapMat.GetColumns();
589
590 DNekMatSharedPtr helm =
592
593 NekDouble one = 1.0;
594 (*helm) = LapMat + factor * MassMat;
595
596 returnval =
598 }
599 break;
601 {
602 MatrixKey helmkey(mkey, StdRegions::eHelmholtz);
603 DNekScalMat &HelmMat = *GetLocMatrix(helmkey);
604
605 // Generate a local copy of traceMat
606 MatrixKey key(mkey, StdRegions::eNormDerivOnTrace);
608
609 ASSERTL1(mkey.ConstFactorExists(StdRegions::eFactorGJP),
610 "Need to specify eFactorGJP to construct "
611 "a HelmholtzGJP matrix");
612
613 NekDouble factor = mkey.GetConstFactor(StdRegions::eFactorGJP);
614
615 factor /= HelmMat.Scale();
616
617 int ntot = HelmMat.GetRows() * HelmMat.GetColumns();
618
619 Vmath::Svtvp(ntot, factor, &NDTraceMat->GetPtr()[0], 1,
620 HelmMat.GetRawPtr(), 1, &NDTraceMat->GetPtr()[0], 1);
621
623 HelmMat.Scale(), NDTraceMat);
624 }
625 break;
627 {
628 NekDouble lambda = mkey.GetConstFactor(StdRegions::eFactorLambda);
629
630 // Construct mass matrix
631 // Check for mass-specific varcoeffs to avoid unncessary
632 // re-computation of the elemental matrix every time step
634 if (mkey.HasVarCoeff(StdRegions::eVarCoeffMass))
635 {
636 massVarcoeffs[StdRegions::eVarCoeffMass] =
637 mkey.GetVarCoeff(StdRegions::eVarCoeffMass);
638 }
639 MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this,
640 mkey.GetConstFactors(), massVarcoeffs);
641 DNekScalMat &MassMat = *GetLocMatrix(masskey);
642
643 // Construct laplacian matrix (Check for varcoeffs)
644 // Take all varcoeffs if one or more are detected
645 // TODO We might want to have a map
646 // from MatrixType to Vector of Varcoeffs and vice-versa
648 if ((mkey.HasVarCoeff(StdRegions::eVarCoeffLaplacian)) ||
649 (mkey.HasVarCoeff(StdRegions::eVarCoeffD00)) ||
650 (mkey.HasVarCoeff(StdRegions::eVarCoeffD01)) ||
651 (mkey.HasVarCoeff(StdRegions::eVarCoeffD10)) ||
652 (mkey.HasVarCoeff(StdRegions::eVarCoeffD02)) ||
653 (mkey.HasVarCoeff(StdRegions::eVarCoeffD20)) ||
654 (mkey.HasVarCoeff(StdRegions::eVarCoeffD11)) ||
655 (mkey.HasVarCoeff(StdRegions::eVarCoeffD12)) ||
656 (mkey.HasVarCoeff(StdRegions::eVarCoeffD21)) ||
657 (mkey.HasVarCoeff(StdRegions::eVarCoeffD22)))
658 {
659 lapVarcoeffs = mkey.GetVarCoeffs();
660 }
661 MatrixKey lapkey(StdRegions::eLaplacian, mkey.GetShapeType(), *this,
662 mkey.GetConstFactors(), lapVarcoeffs);
663 DNekScalMat &LapMat = *GetLocMatrix(lapkey);
664
665 // Construct advection matrix
666 // Check for varcoeffs not required;
667 // assume advection velocity is always time-dependent
668 MatrixKey advkey(mkey, StdRegions::eLinearAdvection);
669 DNekScalMat &AdvMat = *GetLocMatrix(advkey);
670
671 int rows = LapMat.GetRows();
672 int cols = LapMat.GetColumns();
673
674 DNekMatSharedPtr adr =
676
677 NekDouble one = 1.0;
678 (*adr) = LapMat - lambda * MassMat + AdvMat;
679
681
682 // Clear memory for time-dependent matrices
683 DropLocMatrix(advkey);
684 if (!massVarcoeffs.empty())
685 {
686 DropLocMatrix(masskey);
687 }
688 if (!lapVarcoeffs.empty())
689 {
690 DropLocMatrix(lapkey);
691 }
692 }
693 break;
695 {
696 // Copied mostly from ADR solve to have fine-grain control
697 // over updating only advection matrix, relevant for performance!
698 NekDouble lambda = mkey.GetConstFactor(StdRegions::eFactorLambda);
699
700 // Construct mass matrix (Check for varcoeffs)
702 if (mkey.HasVarCoeff(StdRegions::eVarCoeffMass))
703 {
704 massVarcoeffs[StdRegions::eVarCoeffMass] =
705 mkey.GetVarCoeff(StdRegions::eVarCoeffMass);
706 }
707 MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this,
708 mkey.GetConstFactors(), massVarcoeffs);
709 DNekScalMat &MassMat = *GetLocMatrix(masskey);
710
711 // Construct laplacian matrix (Check for varcoeffs)
713 if ((mkey.HasVarCoeff(StdRegions::eVarCoeffLaplacian)) ||
714 (mkey.HasVarCoeff(StdRegions::eVarCoeffD00)) ||
715 (mkey.HasVarCoeff(StdRegions::eVarCoeffD01)) ||
716 (mkey.HasVarCoeff(StdRegions::eVarCoeffD10)) ||
717 (mkey.HasVarCoeff(StdRegions::eVarCoeffD02)) ||
718 (mkey.HasVarCoeff(StdRegions::eVarCoeffD20)) ||
719 (mkey.HasVarCoeff(StdRegions::eVarCoeffD11)) ||
720 (mkey.HasVarCoeff(StdRegions::eVarCoeffD12)) ||
721 (mkey.HasVarCoeff(StdRegions::eVarCoeffD21)) ||
722 (mkey.HasVarCoeff(StdRegions::eVarCoeffD22)))
723 {
724 lapVarcoeffs = mkey.GetVarCoeffs();
725 }
726 MatrixKey lapkey(StdRegions::eLaplacian, mkey.GetShapeType(), *this,
727 mkey.GetConstFactors(), lapVarcoeffs);
728 DNekScalMat &LapMat = *GetLocMatrix(lapkey);
729
730 // Construct advection matrix
731 // (assume advection velocity defined and non-zero)
732 // Could check L2(AdvectionVelocity) or HasVarCoeff
733 MatrixKey advkey(mkey, StdRegions::eLinearAdvection);
734 DNekScalMat &AdvMat = *GetLocMatrix(advkey);
735
736 // Generate a local copy of traceMat
737 MatrixKey gjpkey(StdRegions::eNormDerivOnTrace, mkey.GetShapeType(),
738 *this, mkey.GetConstFactors());
739 DNekScalMat &NDTraceMat = *GetLocMatrix(gjpkey);
740
741 NekDouble gjpfactor = mkey.GetConstFactor(StdRegions::eFactorGJP);
742 ASSERTL1(mkey.ConstFactorExists(StdRegions::eFactorGJP),
743 "Need to specify eFactorGJP to construct "
744 "a LinearAdvectionDiffusionReactionGJP matrix");
745
746 int rows = LapMat.GetRows();
747 int cols = LapMat.GetColumns();
748
749 DNekMatSharedPtr adr =
751
752 NekDouble one = 1.0;
753 (*adr) =
754 LapMat - lambda * MassMat + AdvMat + gjpfactor * NDTraceMat;
755
757
758 // Clear memory
759 DropLocMatrix(advkey);
760 DropLocMatrix(masskey);
761 DropLocMatrix(lapkey);
762 DropLocMatrix(gjpkey);
763 }
764 break;
766 {
767 NekDouble one = 1.0;
769
771 }
772 break;
774 {
775 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
776 {
777 NekDouble one = 1.0;
778 DNekMatSharedPtr mat = GenMatrix(mkey);
779 returnval =
781 }
782 else
783 {
784 NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
785 DNekMatSharedPtr mat = GetStdMatrix(mkey);
786 returnval =
788 }
789 }
790 break;
792 {
793 NekDouble one = 1.0;
794
796 *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
797 DNekMatSharedPtr mat = GenMatrix(hkey);
798
799 mat->Invert();
801 }
802 break;
804 {
805 NekDouble one = 1.0;
806 MatrixKey helmkey(StdRegions::eHelmholtz, mkey.GetShapeType(),
807 *this, mkey.GetConstFactors(),
808 mkey.GetVarCoeffs());
809 DNekScalBlkMatSharedPtr helmStatCond =
810 GetLocStaticCondMatrix(helmkey);
811 DNekScalMatSharedPtr A = helmStatCond->GetBlock(0, 0);
813
815 }
816 break;
818 {
819 NekDouble one = 1.0;
820 MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this);
821 DNekScalBlkMatSharedPtr massStatCond =
822 GetLocStaticCondMatrix(masskey);
823 DNekScalMatSharedPtr A = massStatCond->GetBlock(0, 0);
825
827 }
828 break;
830 {
831 NekDouble one = 1.0;
832 MatrixKey helmkey(StdRegions::eHelmholtz, mkey.GetShapeType(),
833 *this, mkey.GetConstFactors(),
834 mkey.GetVarCoeffs());
835 DNekScalBlkMatSharedPtr helmStatCond =
836 GetLocStaticCondMatrix(helmkey);
837 DNekScalMatSharedPtr A = helmStatCond->GetBlock(0, 0);
838
841 BuildTransformationMatrix(A, mkey.GetMatrixType());
842
844 }
845 break;
847 {
848 NekDouble one = 1.0;
849 MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this);
851 DNekScalMatSharedPtr A = StatCond->GetBlock(0, 0);
852
855 BuildTransformationMatrix(A, mkey.GetMatrixType());
856
858 }
859 break;
860 default:
861 {
862 NekDouble one = 1.0;
863 DNekMatSharedPtr mat = GenMatrix(mkey);
864
866 }
867 break;
868 }
869
870 return returnval;
871}
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:242
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed to...
Definition: ErrorUtil.hpp:265
DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey) override
void DropLocMatrix(const LocalRegions::MatrixKey &mkey)
Definition: Expansion.cpp:90
DNekMatSharedPtr BuildVertexMatrix(const DNekScalMatSharedPtr &r_bnd)
Definition: Expansion.cpp:101
SpatialDomains::GeomFactorsSharedPtr m_metricinfo
Definition: Expansion.h:274
DNekMatSharedPtr BuildTransformationMatrix(const DNekScalMatSharedPtr &r_bnd, const StdRegions::MatrixType matrixType)
Definition: Expansion.cpp:95
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:603
const LibUtilities::PointsKeyVector GetPointsKeys() const
DNekScalBlkMatSharedPtr GetLocStaticCondMatrix(const LocalRegions::MatrixKey &mkey)
Definition: StdExpansion.h:641
DNekMatSharedPtr GenMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:844
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:231
@ eNoGeomType
No type defined.
@ eDeformed
Geometry is curved or has non-constant factors.
@ eLinearAdvectionDiffusionReactionGJP
Definition: StdRegions.hpp:105
static VarCoeffMap NullVarCoeffMap
Definition: StdRegions.hpp:376
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
NekMatrix< NekDouble, StandardMatrixTag > DNekMat
Definition: NekTypeDefs.hpp:50
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:79
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:75
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::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL1, ASSERTL2, Nektar::LocalRegions::Expansion::BuildTransformationMatrix(), Nektar::LocalRegions::Expansion::BuildVertexMatrix(), Nektar::StdRegions::StdMatrixKey::ConstFactorExists(), Nektar::StdRegions::StdExpansion::DetShapeType(), Nektar::LocalRegions::Expansion::DropLocMatrix(), Nektar::SpatialDomains::eDeformed, Nektar::StdRegions::eFactorGJP, Nektar::StdRegions::eFactorLambda, Nektar::StdRegions::eFactorSVVCutoffRatio, Nektar::StdRegions::eHelmholtz, Nektar::StdRegions::eHelmholtzGJP, Nektar::StdRegions::eHybridDGHelmholtz, Nektar::StdRegions::eInvHybridDGHelmholtz, Nektar::StdRegions::eInvMass, Nektar::StdRegions::eIProductWRTBase, Nektar::StdRegions::eLaplacian, Nektar::StdRegions::eLaplacian00, Nektar::StdRegions::eLaplacian01, Nektar::StdRegions::eLaplacian02, Nektar::StdRegions::eLaplacian11, Nektar::StdRegions::eLaplacian12, Nektar::StdRegions::eLaplacian22, Nektar::StdRegions::eLinearAdvection, Nektar::StdRegions::eLinearAdvectionDiffusionReaction, Nektar::StdRegions::eLinearAdvectionDiffusionReactionGJP, Nektar::StdRegions::eMass, Nektar::SpatialDomains::eNoGeomType, Nektar::StdRegions::eNormDerivOnTrace, Nektar::StdRegions::ePreconLinearSpace, Nektar::StdRegions::ePreconLinearSpaceMass, Nektar::StdRegions::ePreconR, Nektar::StdRegions::ePreconRMass, Nektar::StdRegions::eVarCoeffD00, Nektar::StdRegions::eVarCoeffD01, Nektar::StdRegions::eVarCoeffD02, Nektar::StdRegions::eVarCoeffD10, Nektar::StdRegions::eVarCoeffD11, Nektar::StdRegions::eVarCoeffD12, Nektar::StdRegions::eVarCoeffD20, Nektar::StdRegions::eVarCoeffD21, Nektar::StdRegions::eVarCoeffD22, Nektar::StdRegions::eVarCoeffLaplacian, Nektar::StdRegions::eVarCoeffMass, Nektar::StdRegions::eWeakDeriv0, Nektar::StdRegions::eWeakDeriv1, Nektar::StdRegions::eWeakDeriv2, Nektar::StdRegions::StdExpansion::GenMatrix(), Nektar::StdRegions::StdMatrixKey::GetConstFactor(), Nektar::StdRegions::StdMatrixKey::GetConstFactors(), Nektar::LocalRegions::Expansion::GetLocMatrix(), Nektar::StdRegions::StdExpansion::GetLocStaticCondMatrix(), Nektar::StdRegions::StdMatrixKey::GetMatrixType(), Nektar::StdRegions::StdMatrixKey::GetNVarCoeff(), Nektar::StdRegions::StdExpansion::GetPointsKeys(), Nektar::StdRegions::StdMatrixKey::GetShapeType(), Nektar::StdRegions::StdExpansion::GetStdMatrix(), Nektar::StdRegions::StdMatrixKey::GetVarCoeff(), Nektar::StdRegions::StdMatrixKey::GetVarCoeffs(), Nektar::StdRegions::StdMatrixKey::HasVarCoeff(), Nektar::LocalRegions::Expansion::m_metricinfo, Nektar::StdRegions::NullVarCoeffMap, Vmath::Svtvp(), Nektar::Transpose(), and v_GenMatrix().

◆ GetEdgeInverseBoundaryMap()

Array< OneD, unsigned int > Nektar::LocalRegions::Expansion3D::GetEdgeInverseBoundaryMap ( int  eid)

Definition at line 2489 of file Expansion3D.cpp.

2490{
2491 int n, j;
2492 int nEdgeCoeffs;
2493 int nBndCoeffs = NumBndryCoeffs();
2494
2495 Array<OneD, unsigned int> bmap(nBndCoeffs);
2496 GetBoundaryMap(bmap);
2497
2498 // Map from full system to statically condensed system (i.e reverse
2499 // GetBoundaryMap)
2500 map<int, int> invmap;
2501 for (j = 0; j < nBndCoeffs; ++j)
2502 {
2503 invmap[bmap[j]] = j;
2504 }
2505
2506 // Number of interior edge coefficients
2507 nEdgeCoeffs = GetEdgeNcoeffs(eid) - 2;
2508
2510
2511 Array<OneD, unsigned int> edgemaparray(nEdgeCoeffs);
2512 StdRegions::Orientation eOrient = geom->GetEorient(eid);
2513 Array<OneD, unsigned int> maparray = Array<OneD, unsigned int>(nEdgeCoeffs);
2514 Array<OneD, int> signarray = Array<OneD, int>(nEdgeCoeffs, 1);
2515
2516 // maparray is the location of the edge within the matrix
2517 GetEdgeInteriorToElementMap(eid, maparray, signarray, eOrient);
2518
2519 for (n = 0; n < nEdgeCoeffs; ++n)
2520 {
2521 edgemaparray[n] = invmap[maparray[n]];
2522 }
2523
2524 return edgemaparray;
2525}
SpatialDomains::Geometry3DSharedPtr GetGeom3D() const
Definition: Expansion3D.h:176
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)
void GetBoundaryMap(Array< OneD, unsigned int > &outarray)
Definition: StdExpansion.h:669
std::shared_ptr< Geometry3D > Geometry3DSharedPtr
Definition: Geometry3D.h:50

References Nektar::StdRegions::StdExpansion::GetBoundaryMap(), Nektar::StdRegions::StdExpansion3D::GetEdgeInteriorToElementMap(), Nektar::StdRegions::StdExpansion3D::GetEdgeNcoeffs(), GetGeom3D(), and Nektar::StdRegions::StdExpansion::NumBndryCoeffs().

Referenced by v_BuildInverseTransformationMatrix(), and v_BuildTransformationMatrix().

◆ GetGeom3D()

SpatialDomains::Geometry3DSharedPtr Nektar::LocalRegions::Expansion3D::GetGeom3D ( ) const
inline

Definition at line 176 of file Expansion3D.h.

177{
178 return std::dynamic_pointer_cast<SpatialDomains::Geometry3D>(m_geom);
179}
SpatialDomains::GeometrySharedPtr m_geom
Definition: Expansion.h:273

References Nektar::LocalRegions::Expansion::m_geom.

Referenced by GetEdgeInverseBoundaryMap(), GetTraceInverseBoundaryMap(), v_BuildTransformationMatrix(), and v_GetTracePhysVals().

◆ GetInverseBoundaryMaps()

void Nektar::LocalRegions::Expansion3D::GetInverseBoundaryMaps ( Array< OneD, unsigned int > &  vmap,
Array< OneD, Array< OneD, unsigned int > > &  emap,
Array< OneD, Array< OneD, unsigned int > > &  fmap 
)

Definition at line 2641 of file Expansion3D.cpp.

2645{
2646 int n, j;
2647 int nEdgeCoeffs;
2648 int nFaceCoeffs;
2649
2650 int nBndCoeffs = NumBndryCoeffs();
2651
2652 Array<OneD, unsigned int> bmap(nBndCoeffs);
2653 GetBoundaryMap(bmap);
2654
2655 // Map from full system to statically condensed system (i.e reverse
2656 // GetBoundaryMap)
2657 map<int, int> reversemap;
2658 for (j = 0; j < bmap.size(); ++j)
2659 {
2660 reversemap[bmap[j]] = j;
2661 }
2662
2663 int nverts = GetNverts();
2664 vmap = Array<OneD, unsigned int>(nverts);
2665 for (n = 0; n < nverts; ++n)
2666 {
2667 int id = GetVertexMap(n);
2668 vmap[n] = reversemap[id]; // not sure what should be true here.
2669 }
2670
2671 int nedges = GetNedges();
2672 emap = Array<OneD, Array<OneD, unsigned int>>(nedges);
2673
2674 for (int eid = 0; eid < nedges; ++eid)
2675 {
2676 // Number of interior edge coefficients
2677 nEdgeCoeffs = GetEdgeNcoeffs(eid) - 2;
2678
2679 Array<OneD, unsigned int> edgemaparray(nEdgeCoeffs);
2680 Array<OneD, unsigned int> maparray =
2681 Array<OneD, unsigned int>(nEdgeCoeffs);
2682 Array<OneD, int> signarray = Array<OneD, int>(nEdgeCoeffs, 1);
2683
2684 // maparray is the location of the edge within the matrix
2685 GetEdgeInteriorToElementMap(eid, maparray, signarray,
2687
2688 for (n = 0; n < nEdgeCoeffs; ++n)
2689 {
2690 edgemaparray[n] = reversemap[maparray[n]];
2691 }
2692 emap[eid] = edgemaparray;
2693 }
2694
2695 int nfaces = GetNtraces();
2696 fmap = Array<OneD, Array<OneD, unsigned int>>(nfaces);
2697
2698 for (int fid = 0; fid < nfaces; ++fid)
2699 {
2700 // Number of interior face coefficients
2701 nFaceCoeffs = GetTraceIntNcoeffs(fid);
2702
2703 Array<OneD, unsigned int> facemaparray(nFaceCoeffs);
2704 Array<OneD, unsigned int> maparray =
2705 Array<OneD, unsigned int>(nFaceCoeffs);
2706 Array<OneD, int> signarray = Array<OneD, int>(nFaceCoeffs, 1);
2707
2708 // maparray is the location of the face within the matrix
2709 GetTraceInteriorToElementMap(fid, maparray, signarray,
2711
2712 for (n = 0; n < nFaceCoeffs; ++n)
2713 {
2714 facemaparray[n] = reversemap[maparray[n]];
2715 }
2716
2717 fmap[fid] = facemaparray;
2718 }
2719}
int GetNedges() const
return the number of edges in 3D expansion
int GetTraceIntNcoeffs(const int i) const
Definition: StdExpansion.h:266
int GetVertexMap(const int localVertexId, bool useCoeffPacking=false)
Definition: StdExpansion.h:679
int GetNverts() const
This function returns the number of vertices of the expansion domain.
Definition: StdExpansion.h:246
void GetTraceInteriorToElementMap(const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, const Orientation traceOrient=eForwards)
Definition: StdExpansion.h:708

References Nektar::StdRegions::eDir1FwdDir1_Dir2FwdDir2, Nektar::StdRegions::eForwards, Nektar::StdRegions::StdExpansion::GetBoundaryMap(), Nektar::StdRegions::StdExpansion3D::GetEdgeInteriorToElementMap(), Nektar::StdRegions::StdExpansion3D::GetEdgeNcoeffs(), Nektar::StdRegions::StdExpansion3D::GetNedges(), Nektar::StdRegions::StdExpansion::GetNtraces(), Nektar::StdRegions::StdExpansion::GetNverts(), Nektar::StdRegions::StdExpansion::GetTraceInteriorToElementMap(), Nektar::StdRegions::StdExpansion::GetTraceIntNcoeffs(), Nektar::StdRegions::StdExpansion::GetVertexMap(), and Nektar::StdRegions::StdExpansion::NumBndryCoeffs().

◆ GetnFacecdotMF()

Array< OneD, NekDouble > Nektar::LocalRegions::Expansion3D::GetnFacecdotMF ( const int  dir,
const int  face,
ExpansionSharedPtr FaceExp_f,
const Array< OneD, const Array< OneD, NekDouble > > &  normals,
const StdRegions::VarCoeffMap varcoeffs 
)
private

Definition at line 2946 of file Expansion3D.cpp.

2950{
2951 int nquad_f = FaceExp_f->GetNumPoints(0) * FaceExp_f->GetNumPoints(1);
2952 int coordim = GetCoordim();
2953
2954 int nquad0 = m_base[0]->GetNumPoints();
2955 int nquad1 = m_base[1]->GetNumPoints();
2956 int nquad2 = m_base[2]->GetNumPoints();
2957 int nqtot = nquad0 * nquad1 * nquad2;
2958
2959 StdRegions::VarCoeffType MMFCoeffs[15] = {
2968
2969 StdRegions::VarCoeffMap::const_iterator MFdir;
2970
2971 Array<OneD, NekDouble> nFacecdotMF(nqtot, 0.0);
2972 Array<OneD, NekDouble> tmp(nqtot);
2973 Array<OneD, NekDouble> tmp_f(nquad_f);
2974 for (int k = 0; k < coordim; k++)
2975 {
2976 MFdir = varcoeffs.find(MMFCoeffs[dir * 5 + k]);
2977 tmp = MFdir->second.GetValue();
2978
2979 GetPhysFaceVarCoeffsFromElement(face, FaceExp_f, tmp, tmp_f);
2980
2981 Vmath::Vvtvp(nquad_f, &tmp_f[0], 1, &normals[k][0], 1, &nFacecdotMF[0],
2982 1, &nFacecdotMF[0], 1);
2983 }
2984
2985 return nFacecdotMF;
2986}
void GetPhysFaceVarCoeffsFromElement(const int face, ExpansionSharedPtr &FaceExp, const Array< OneD, const NekDouble > &varcoeff, Array< OneD, NekDouble > &outarray)
Array< OneD, LibUtilities::BasisSharedPtr > m_base
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.hpp:366

References Nektar::StdRegions::eVarCoeffMF1Div, Nektar::StdRegions::eVarCoeffMF1Mag, Nektar::StdRegions::eVarCoeffMF1x, Nektar::StdRegions::eVarCoeffMF1y, Nektar::StdRegions::eVarCoeffMF1z, Nektar::StdRegions::eVarCoeffMF2Div, Nektar::StdRegions::eVarCoeffMF2Mag, Nektar::StdRegions::eVarCoeffMF2x, Nektar::StdRegions::eVarCoeffMF2y, Nektar::StdRegions::eVarCoeffMF2z, Nektar::StdRegions::eVarCoeffMF3Div, Nektar::StdRegions::eVarCoeffMF3Mag, Nektar::StdRegions::eVarCoeffMF3x, Nektar::StdRegions::eVarCoeffMF3y, Nektar::StdRegions::eVarCoeffMF3z, Nektar::StdRegions::StdExpansion::GetCoordim(), GetPhysFaceVarCoeffsFromElement(), Nektar::StdRegions::StdExpansion::m_base, and Vmath::Vvtvp().

Referenced by AddHDGHelmholtzFaceTerms(), AddNormTraceInt(), and v_GenMatrix().

◆ GetPhysFaceVarCoeffsFromElement()

void Nektar::LocalRegions::Expansion3D::GetPhysFaceVarCoeffsFromElement ( const int  face,
ExpansionSharedPtr FaceExp,
const Array< OneD, const NekDouble > &  varcoeff,
Array< OneD, NekDouble > &  outarray 
)
protected

Definition at line 202 of file Expansion3D.cpp.

206{
207 Array<OneD, NekDouble> tmp(GetNcoeffs());
208 Array<OneD, NekDouble> facetmp(FaceExp->GetNcoeffs());
209
210 // FwdTrans varcoeffs
211 FwdTrans(varcoeff, tmp);
212
213 // Map to edge
214 Array<OneD, unsigned int> emap;
215 Array<OneD, int> sign;
216
217 GetTraceToElementMap(face, emap, sign, GetTraceOrient(face));
218
219 for (unsigned int i = 0; i < FaceExp->GetNcoeffs(); ++i)
220 {
221 facetmp[i] = tmp[emap[i]];
222 }
223
224 // BwdTrans
225 FaceExp->BwdTrans(facetmp, outarray);
226}
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:47
void GetTraceToElementMap(const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, Orientation traceOrient=eForwards, int P=-1, int Q=-1)
Definition: StdExpansion.h:684
void FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Forward transformation from physical space to coefficient space.

References Nektar::StdRegions::StdExpansion::FwdTrans(), Nektar::StdRegions::StdExpansion::GetNcoeffs(), Nektar::LocalRegions::Expansion::GetTraceOrient(), Nektar::StdRegions::StdExpansion::GetTraceToElementMap(), and sign.

Referenced by GetnFacecdotMF().

◆ GetTraceInverseBoundaryMap()

Array< OneD, unsigned int > Nektar::LocalRegions::Expansion3D::GetTraceInverseBoundaryMap ( int  fid,
StdRegions::Orientation  faceOrient = StdRegions::eNoOrientation,
int  P1 = -1,
int  P2 = -1 
)

Definition at line 2527 of file Expansion3D.cpp.

2529{
2530 int n, j;
2531 int nFaceCoeffs;
2532
2533 int nBndCoeffs = NumBndryCoeffs();
2534
2535 Array<OneD, unsigned int> bmap(nBndCoeffs);
2536 GetBoundaryMap(bmap);
2537
2538 // Map from full system to statically condensed system (i.e reverse
2539 // GetBoundaryMap)
2540 map<int, int> reversemap;
2541 for (j = 0; j < bmap.size(); ++j)
2542 {
2543 reversemap[bmap[j]] = j;
2544 }
2545
2546 // Number of interior face coefficients
2547 nFaceCoeffs = GetTraceIntNcoeffs(fid);
2548
2550 Array<OneD, unsigned int> maparray = Array<OneD, unsigned int>(nFaceCoeffs);
2551 Array<OneD, int> signarray = Array<OneD, int>(nFaceCoeffs, 1);
2552
2553 if (faceOrient == StdRegions::eNoOrientation)
2554 {
2555 fOrient = GetTraceOrient(fid);
2556 }
2557 else
2558 {
2559 fOrient = faceOrient;
2560 }
2561
2562 // maparray is the location of the face within the matrix
2563 GetTraceInteriorToElementMap(fid, maparray, signarray, fOrient);
2564
2565 Array<OneD, unsigned int> facemaparray;
2566 int locP1, locP2;
2567 GetTraceNumModes(fid, locP1, locP2, fOrient);
2568
2569 if (P1 == -1)
2570 {
2571 P1 = locP1;
2572 }
2573 else
2574 {
2575 ASSERTL1(P1 <= locP1, "Expect value of passed P1 to "
2576 "be lower or equal to face num modes");
2577 }
2578
2579 if (P2 == -1)
2580 {
2581 P2 = locP2;
2582 }
2583 else
2584 {
2585 ASSERTL1(P2 <= locP2, "Expect value of passed P2 to "
2586 "be lower or equal to face num modes");
2587 }
2588
2589 switch (GetGeom3D()->GetFace(fid)->GetShapeType())
2590 {
2592 {
2593 if (((P1 - 3) > 0) && ((P2 - 3) > 0))
2594 {
2595 facemaparray = Array<OneD, unsigned int>(
2597 P2 - 3));
2598 int cnt = 0;
2599 int cnt1 = 0;
2600 for (n = 0; n < P1 - 3; ++n)
2601 {
2602 for (int m = 0; m < P2 - 3 - n; ++m, ++cnt)
2603 {
2604 facemaparray[cnt] = reversemap[maparray[cnt1 + m]];
2605 }
2606 cnt1 += locP2 - 3 - n;
2607 }
2608 }
2609 }
2610 break;
2612 {
2613 if (((P1 - 2) > 0) && ((P2 - 2) > 0))
2614 {
2615 facemaparray = Array<OneD, unsigned int>(
2617 P2 - 2));
2618 int cnt = 0;
2619 int cnt1 = 0;
2620 for (n = 0; n < P2 - 2; ++n)
2621 {
2622 for (int m = 0; m < P1 - 2; ++m, ++cnt)
2623 {
2624 facemaparray[cnt] = reversemap[maparray[cnt1 + m]];
2625 }
2626 cnt1 += locP1 - 2;
2627 }
2628 }
2629 }
2630 break;
2631 default:
2632 {
2633 ASSERTL0(false, "Invalid shape type.");
2634 }
2635 break;
2636 }
2637
2638 return facemaparray;
2639}
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
void GetTraceNumModes(const int tid, int &numModes0, int &numModes1, const Orientation traceOrient=eDir1FwdDir1_Dir2FwdDir2)
Definition: StdExpansion.h:716
int getNumberOfCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:133
int getNumberOfCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:109

References ASSERTL0, ASSERTL1, Nektar::StdRegions::eNoOrientation, Nektar::LibUtilities::eQuadrilateral, Nektar::LibUtilities::eTriangle, Nektar::StdRegions::StdExpansion::GetBoundaryMap(), GetGeom3D(), Nektar::LibUtilities::StdTriData::getNumberOfCoefficients(), Nektar::LibUtilities::StdQuadData::getNumberOfCoefficients(), Nektar::StdRegions::StdExpansion::GetTraceInteriorToElementMap(), Nektar::StdRegions::StdExpansion::GetTraceIntNcoeffs(), Nektar::StdRegions::StdExpansion::GetTraceNumModes(), Nektar::LocalRegions::Expansion::GetTraceOrient(), and Nektar::StdRegions::StdExpansion::NumBndryCoeffs().

Referenced by v_BuildInverseTransformationMatrix(), and v_BuildTransformationMatrix().

◆ SetFaceToGeomOrientation()

void Nektar::LocalRegions::Expansion3D::SetFaceToGeomOrientation ( const int  face,
Array< OneD, NekDouble > &  inout 
)

Align face orientation with the geometry orientation.

Definition at line 348 of file Expansion3D.cpp.

350{
351 int j, k;
352 int nface = GetTraceNcoeffs(face);
353 Array<OneD, NekDouble> f_in(nface);
354 Vmath::Vcopy(nface, &inout[0], 1, &f_in[0], 1);
355
356 // retreiving face to element map for standard face orientation and
357 // for actual face orientation
358 IndexMapKey ikey1(eFaceToElement, DetShapeType(), GetBasisNumModes(0),
362 IndexMapKey ikey2(eFaceToElement, DetShapeType(), GetBasisNumModes(0),
364 GetTraceOrient(face));
366
367 ASSERTL1((*map1).size() == (*map2).size(),
368 "There is an error with the GetTraceToElementMap");
369
370 for (j = 0; j < (*map1).size(); ++j)
371 {
372 // j = index in the standard orientation
373 for (k = 0; k < (*map2).size(); ++k)
374 {
375 // k = index in the actual orientation
376 if ((*map1)[j].index == (*map2)[k].index && k != j)
377 {
378 inout[k] = f_in[j];
379 // checking if sign is changing
380 if ((*map1)[j].sign != (*map2)[k].sign)
381 {
382 inout[k] *= -1.0;
383 }
384 break;
385 }
386 }
387 }
388}
int GetTraceNcoeffs(const int i) const
This function returns the number of expansion coefficients belonging to the i-th trace.
Definition: StdExpansion.h:261
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825

References ASSERTL1, Nektar::StdRegions::StdExpansion::DetShapeType(), Nektar::StdRegions::eDir1FwdDir1_Dir2FwdDir2, Nektar::LocalRegions::eFaceToElement, Nektar::StdRegions::StdExpansion::GetBasisNumModes(), Nektar::LocalRegions::Expansion::GetIndexMap(), Nektar::StdRegions::StdExpansion::GetTraceNcoeffs(), Nektar::LocalRegions::Expansion::GetTraceOrient(), sign, and Vmath::Vcopy().

Referenced by SetTraceToGeomOrientation(), and v_GenMatrix().

◆ SetTraceToGeomOrientation()

void Nektar::LocalRegions::Expansion3D::SetTraceToGeomOrientation ( Array< OneD, NekDouble > &  inout)

Align trace orientation with the geometry orientation.

Definition at line 393 of file Expansion3D.cpp.

394{
395 int i, cnt = 0;
396 int nfaces = GetNtraces();
397
398 Array<OneD, NekDouble> f_tmp;
399
400 for (i = 0; i < nfaces; ++i)
401 {
402 SetFaceToGeomOrientation(i, f_tmp = inout + cnt);
403 cnt += GetTraceNcoeffs(i);
404 }
405}
void SetFaceToGeomOrientation(const int face, Array< OneD, NekDouble > &inout)
Align face orientation with the geometry orientation.

References Nektar::StdRegions::StdExpansion::GetNtraces(), Nektar::StdRegions::StdExpansion::GetTraceNcoeffs(), and SetFaceToGeomOrientation().

Referenced by v_GenMatrix().

◆ v_AddFaceNormBoundaryInt()

void Nektar::LocalRegions::Expansion3D::v_AddFaceNormBoundaryInt ( const int  face,
const ExpansionSharedPtr FaceExp,
const Array< OneD, const NekDouble > &  Fn,
Array< OneD, NekDouble > &  outarray 
)
overrideprotected

Definition at line 1649 of file Expansion3D.cpp.

1652{
1653 int i, j;
1654
1655 /*
1656 * Coming into this routine, the velocity V will have been
1657 * multiplied by the trace normals to give the input vector Vn. By
1658 * convention, these normals are inwards facing for elements which
1659 * have FaceExp as their right-adjacent face. This conditional
1660 * statement therefore determines whether the normals must be
1661 * negated, since the integral being performed here requires an
1662 * outwards facing normal.
1663 */
1664 if (m_requireNeg.size() == 0)
1665 {
1666 m_requireNeg.resize(GetNtraces());
1667
1668 for (i = 0; i < GetNtraces(); ++i)
1669 {
1670 m_requireNeg[i] = false;
1671
1672 ExpansionSharedPtr faceExp = m_traceExp[i].lock();
1673
1674 if (faceExp->GetRightAdjacentElementExp())
1675 {
1676 if (faceExp->GetRightAdjacentElementExp()
1677 ->GetGeom()
1678 ->GetGlobalID() == GetGeom()->GetGlobalID())
1679 {
1680 m_requireNeg[i] = true;
1681 }
1682 }
1683 }
1684 }
1685
1686 IndexMapKey ikey(eFaceToElement, DetShapeType(), GetBasisNumModes(0),
1688 GetTraceOrient(face));
1690
1691 int order_e = (*map).size(); // Order of the element
1692 int n_coeffs = FaceExp->GetNcoeffs();
1693
1694 Array<OneD, NekDouble> faceCoeffs(n_coeffs);
1695
1696 if (n_coeffs != order_e) // Going to orthogonal space
1697 {
1698 FaceExp->FwdTrans(Fn, faceCoeffs);
1699
1700 int NumModesElementMax = FaceExp->GetBasis(0)->GetNumModes();
1701 int NumModesElementMin = m_base[0]->GetNumModes();
1702
1703 FaceExp->ReduceOrderCoeffs(NumModesElementMin, faceCoeffs, faceCoeffs);
1704
1705 StdRegions::StdMatrixKey masskey(StdRegions::eMass,
1706 FaceExp->DetShapeType(), *FaceExp);
1707 FaceExp->MassMatrixOp(faceCoeffs, faceCoeffs, masskey);
1708
1709 // Reorder coefficients for the lower degree face.
1710 int offset1 = 0, offset2 = 0;
1711
1712 if (FaceExp->DetShapeType() == LibUtilities::eQuadrilateral)
1713 {
1714 for (i = 0; i < NumModesElementMin; ++i)
1715 {
1716 for (j = 0; j < NumModesElementMin; ++j)
1717 {
1718 faceCoeffs[offset1 + j] = faceCoeffs[offset2 + j];
1719 }
1720 offset1 += NumModesElementMin;
1721 offset2 += NumModesElementMax;
1722 }
1723
1724 // Extract lower degree modes. TODO: Check this is correct.
1725 for (i = NumModesElementMin; i < NumModesElementMax; ++i)
1726 {
1727 for (j = NumModesElementMin; j < NumModesElementMax; ++j)
1728 {
1729 faceCoeffs[i * NumModesElementMax + j] = 0.0;
1730 }
1731 }
1732 }
1733
1734 if (FaceExp->DetShapeType() == LibUtilities::eTriangle)
1735 {
1736
1737 // Reorder coefficients for the lower degree face.
1738 int offset1 = 0, offset2 = 0;
1739
1740 for (i = 0; i < NumModesElementMin; ++i)
1741 {
1742 for (j = 0; j < NumModesElementMin - i; ++j)
1743 {
1744 faceCoeffs[offset1 + j] = faceCoeffs[offset2 + j];
1745 }
1746 offset1 += NumModesElementMin - i;
1747 offset2 += NumModesElementMax - i;
1748 }
1749 }
1750 }
1751 else
1752 {
1753 FaceExp->IProductWRTBase(Fn, faceCoeffs);
1754 }
1755
1756 if (m_requireNeg[face])
1757 {
1758 for (i = 0; i < order_e; ++i)
1759 {
1760 outarray[(*map)[i].index] -= (*map)[i].sign * faceCoeffs[i];
1761 }
1762 }
1763 else
1764 {
1765 for (i = 0; i < order_e; ++i)
1766 {
1767 outarray[(*map)[i].index] += (*map)[i].sign * faceCoeffs[i];
1768 }
1769 }
1770}
SpatialDomains::GeometrySharedPtr GetGeom() const
Definition: Expansion.cpp:167
std::map< int, ExpansionWeakPtr > m_traceExp
Definition: Expansion.h:272

References Nektar::StdRegions::StdExpansion::DetShapeType(), Nektar::LocalRegions::eFaceToElement, Nektar::StdRegions::eMass, Nektar::LibUtilities::eQuadrilateral, Nektar::LibUtilities::eTriangle, Nektar::StdRegions::StdExpansion::GetBasisNumModes(), Nektar::LocalRegions::Expansion::GetGeom(), Nektar::LocalRegions::Expansion::GetIndexMap(), Nektar::StdRegions::StdExpansion::GetNtraces(), Nektar::LocalRegions::Expansion::GetTraceOrient(), Nektar::StdRegions::StdExpansion::m_base, m_requireNeg, and Nektar::LocalRegions::Expansion::m_traceExp.

◆ v_AddRobinMassMatrix()

void Nektar::LocalRegions::Expansion3D::v_AddRobinMassMatrix ( const int  face,
const Array< OneD, const NekDouble > &  primCoeffs,
DNekMatSharedPtr inoutmat 
)
overrideprotectedvirtual

Reimplemented from Nektar::LocalRegions::Expansion.

Definition at line 1806 of file Expansion3D.cpp.

1809{
1811 "Not set up for non boundary-interior expansions");
1812 ASSERTL1(inoutmat->GetRows() == inoutmat->GetColumns(),
1813 "Assuming that input matrix was square");
1814
1815 int i, j;
1816 int id1, id2;
1817 ExpansionSharedPtr faceExp = m_traceExp[face].lock();
1818 int order_f = faceExp->GetNcoeffs();
1819
1820 Array<OneD, unsigned int> map;
1821 Array<OneD, int> sign;
1822
1823 StdRegions::VarCoeffMap varcoeffs;
1824 varcoeffs[StdRegions::eVarCoeffMass] = primCoeffs;
1825
1826 LibUtilities::ShapeType shapeType = faceExp->DetShapeType();
1827
1828 LocalRegions::MatrixKey mkey(StdRegions::eMass, shapeType, *faceExp,
1830
1831 DNekScalMat &facemat = *faceExp->GetLocMatrix(mkey);
1832
1833 // Now need to identify a map which takes the local face
1834 // mass matrix to the matrix stored in inoutmat;
1835 // This can currently be deduced from the size of the matrix
1836
1837 // - if inoutmat.m_rows() == v_NCoeffs() it is a full
1838 // matrix system
1839
1840 // - if inoutmat.m_rows() == v_GetNverts() it is a vertex space
1841 // preconditioner.
1842
1843 // - if inoutmat.m_rows() == v_NumBndCoeffs() it is a
1844 // boundary CG system
1845
1846 // - if inoutmat.m_rows() == v_NumDGBndCoeffs() it is a
1847 // trace DG system; still needs implementing.
1848 int rows = inoutmat->GetRows();
1849
1850 if (rows == GetNcoeffs())
1851 {
1852 GetTraceToElementMap(face, map, sign, GetTraceOrient(face));
1853 }
1854 else if (rows == GetNverts())
1855 {
1856 int nfvert = faceExp->GetNverts();
1857
1858 // Need to find where linear vertices are in facemat
1859 Array<OneD, unsigned int> linmap;
1860 Array<OneD, int> linsign;
1861
1862 // Use a linear expansion to get correct mapping
1863 GetLinStdExp()->GetTraceToElementMap(face, linmap, linsign,
1864 GetTraceOrient(face));
1865
1866 // zero out sign map to remove all other modes
1867 sign = Array<OneD, int>(order_f, 0);
1868 map = Array<OneD, unsigned int>(order_f, (unsigned int)0);
1869
1870 int fmap;
1871 // Reset sign map to only have contribution from vertices
1872 for (i = 0; i < nfvert; ++i)
1873 {
1874 fmap = faceExp->GetVertexMap(i, true);
1875 sign[fmap] = 1;
1876
1877 // need to reset map
1878 map[fmap] = linmap[i];
1879 }
1880 }
1881 else if (rows == NumBndryCoeffs())
1882 {
1883 int nbndry = NumBndryCoeffs();
1884 Array<OneD, unsigned int> bmap(nbndry);
1885
1886 GetTraceToElementMap(face, map, sign, GetTraceOrient(face));
1887 GetBoundaryMap(bmap);
1888
1889 for (i = 0; i < order_f; ++i)
1890 {
1891 for (j = 0; j < nbndry; ++j)
1892 {
1893 if (map[i] == bmap[j])
1894 {
1895 map[i] = j;
1896 break;
1897 }
1898 }
1899 ASSERTL1(j != nbndry, "Did not find number in map");
1900 }
1901 }
1902 else if (rows == NumDGBndryCoeffs())
1903 {
1904 // possibly this should be a separate method
1905 int cnt = 0;
1906 map = Array<OneD, unsigned int>(order_f);
1907 sign = Array<OneD, int>(order_f, 1);
1908
1909 IndexMapKey ikey1(eFaceToElement, DetShapeType(), GetBasisNumModes(0),
1911 GetTraceOrient(face));
1913 IndexMapKey ikey2(eFaceToElement, DetShapeType(), GetBasisNumModes(0),
1917
1918 ASSERTL1((*map1).size() == (*map2).size(),
1919 "There is an error with the GetTraceToElementMap");
1920
1921 for (i = 0; i < face; ++i)
1922 {
1923 cnt += GetTraceNcoeffs(i);
1924 }
1925
1926 for (i = 0; i < (*map1).size(); ++i)
1927 {
1928 int idx = -1;
1929
1930 for (j = 0; j < (*map2).size(); ++j)
1931 {
1932 if ((*map1)[i].index == (*map2)[j].index)
1933 {
1934 idx = j;
1935 break;
1936 }
1937 }
1938
1939 ASSERTL2(idx >= 0, "Index not found");
1940 map[i] = idx + cnt;
1941 sign[i] = (*map2)[idx].sign;
1942 }
1943 }
1944 else
1945 {
1946 ASSERTL0(false, "Could not identify matrix type from dimension");
1947 }
1948
1949 for (i = 0; i < order_f; ++i)
1950 {
1951 id1 = map[i];
1952 for (j = 0; j < order_f; ++j)
1953 {
1954 id2 = map[j];
1955 (*inoutmat)(id1, id2) += facemat(i, j) * sign[i] * sign[j];
1956 }
1957 }
1958}
std::shared_ptr< StdExpansion > GetLinStdExp(void) const
Definition: StdExpansion.h:377

References ASSERTL0, ASSERTL1, ASSERTL2, Nektar::StdRegions::StdExpansion::DetShapeType(), Nektar::StdRegions::eDir1FwdDir1_Dir2FwdDir2, Nektar::LocalRegions::eFaceToElement, Nektar::StdRegions::eMass, Nektar::StdRegions::eVarCoeffMass, Nektar::StdRegions::StdExpansion::GetBasisNumModes(), Nektar::StdRegions::StdExpansion::GetBoundaryMap(), Nektar::LocalRegions::Expansion::GetIndexMap(), Nektar::StdRegions::StdExpansion::GetLinStdExp(), Nektar::StdRegions::StdExpansion::GetNcoeffs(), Nektar::StdRegions::StdExpansion::GetNverts(), Nektar::StdRegions::StdExpansion::GetTraceNcoeffs(), Nektar::LocalRegions::Expansion::GetTraceOrient(), Nektar::StdRegions::StdExpansion::GetTraceToElementMap(), Nektar::StdRegions::StdExpansion::IsBoundaryInteriorExpansion(), Nektar::LocalRegions::Expansion::m_traceExp, Nektar::StdRegions::NullConstFactorMap, Nektar::StdRegions::StdExpansion::NumBndryCoeffs(), Nektar::StdRegions::StdExpansion::NumDGBndryCoeffs(), and sign.

◆ v_BuildInverseTransformationMatrix()

DNekMatSharedPtr Nektar::LocalRegions::Expansion3D::v_BuildInverseTransformationMatrix ( const DNekScalMatSharedPtr transformationmatrix)
overrideprotectedvirtual

Build inverse and inverse transposed transformation matrix: \(\mathbf{R^{-1}}\) and \(\mathbf{R^{-T}}\).

\f\mathbf{R^{-T}}=[\left[\begin{array}{ccc} \mathbf{I} & -\mathbf{R}_{ef} & -\mathbf{R}_{ve}+\mathbf{R}_{ve}\mathbf{R}_{vf} \ 0 & \mathbf{I} & \mathbf{R}_{ef} \ 0 & 0 & \mathbf{I}} \end{array}\right]\f]

Reimplemented from Nektar::StdRegions::StdExpansion.

Definition at line 2372 of file Expansion3D.cpp.

2374{
2375 int i, j, n, eid = 0, fid = 0;
2376 int nCoeffs = NumBndryCoeffs();
2377 NekDouble MatrixValue;
2378 DNekScalMat &R = (*transformationmatrix);
2379
2380 // Define storage for vertex transpose matrix and zero all entries
2381 MatrixStorage storage = eFULL;
2382
2383 // Inverse transformation matrix
2384 DNekMatSharedPtr inversetransformationmatrix =
2385 MemoryManager<DNekMat>::AllocateSharedPtr(nCoeffs, nCoeffs, 0.0,
2386 storage);
2387 DNekMat &InvR = (*inversetransformationmatrix);
2388
2389 int nVerts = GetNverts();
2390 int nEdges = GetNedges();
2391 int nFaces = GetNtraces();
2392
2393 int nedgemodes = 0;
2394 int nfacemodes = 0;
2395 int nedgemodestotal = 0;
2396 int nfacemodestotal = 0;
2397
2398 for (eid = 0; eid < nEdges; ++eid)
2399 {
2400 nedgemodes = GetEdgeNcoeffs(eid) - 2;
2401 nedgemodestotal += nedgemodes;
2402 }
2403
2404 for (fid = 0; fid < nFaces; ++fid)
2405 {
2406 nfacemodes = GetTraceIntNcoeffs(fid);
2407 nfacemodestotal += nfacemodes;
2408 }
2409
2410 Array<OneD, unsigned int> edgemodearray(nedgemodestotal);
2411 Array<OneD, unsigned int> facemodearray(nfacemodestotal);
2412
2413 int offset = 0;
2414
2415 // Create array of edge modes
2416 for (eid = 0; eid < nEdges; ++eid)
2417 {
2418 Array<OneD, unsigned int> edgearray = GetEdgeInverseBoundaryMap(eid);
2419 nedgemodes = GetEdgeNcoeffs(eid) - 2;
2420
2421 // Only copy if there are edge modes
2422 if (nedgemodes)
2423 {
2424 Vmath::Vcopy(nedgemodes, &edgearray[0], 1, &edgemodearray[offset],
2425 1);
2426 }
2427
2428 offset += nedgemodes;
2429 }
2430
2431 offset = 0;
2432
2433 // Create array of face modes
2434 for (fid = 0; fid < nFaces; ++fid)
2435 {
2436 Array<OneD, unsigned int> facearray = GetTraceInverseBoundaryMap(fid);
2437 nfacemodes = GetTraceIntNcoeffs(fid);
2438
2439 // Only copy if there are face modes
2440 if (nfacemodes)
2441 {
2442 Vmath::Vcopy(nfacemodes, &facearray[0], 1, &facemodearray[offset],
2443 1);
2444 }
2445
2446 offset += nfacemodes;
2447 }
2448
2449 // Vertex-edge/face
2450 for (i = 0; i < nVerts; ++i)
2451 {
2452 for (j = 0; j < nedgemodestotal; ++j)
2453 {
2454 InvR.SetValue(GetVertexMap(i), edgemodearray[j],
2455 -R(GetVertexMap(i), edgemodearray[j]));
2456 }
2457 for (j = 0; j < nfacemodestotal; ++j)
2458 {
2459 InvR.SetValue(GetVertexMap(i), facemodearray[j],
2460 -R(GetVertexMap(i), facemodearray[j]));
2461 for (n = 0; n < nedgemodestotal; ++n)
2462 {
2463 MatrixValue = InvR.GetValue(GetVertexMap(i), facemodearray[j]) +
2464 R(GetVertexMap(i), edgemodearray[n]) *
2465 R(edgemodearray[n], facemodearray[j]);
2466 InvR.SetValue(GetVertexMap(i), facemodearray[j], MatrixValue);
2467 }
2468 }
2469 }
2470
2471 // Edge-face contributions
2472 for (i = 0; i < nedgemodestotal; ++i)
2473 {
2474 for (j = 0; j < nfacemodestotal; ++j)
2475 {
2476 InvR.SetValue(edgemodearray[i], facemodearray[j],
2477 -R(edgemodearray[i], facemodearray[j]));
2478 }
2479 }
2480
2481 for (i = 0; i < nCoeffs; ++i)
2482 {
2483 InvR.SetValue(i, i, 1.0);
2484 }
2485
2486 return inversetransformationmatrix;
2487}
Array< OneD, unsigned int > GetTraceInverseBoundaryMap(int fid, StdRegions::Orientation faceOrient=StdRegions::eNoOrientation, int P1=-1, int P2=-1)
Array< OneD, unsigned int > GetEdgeInverseBoundaryMap(int eid)

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::eFULL, GetEdgeInverseBoundaryMap(), Nektar::StdRegions::StdExpansion3D::GetEdgeNcoeffs(), Nektar::StdRegions::StdExpansion3D::GetNedges(), Nektar::StdRegions::StdExpansion::GetNtraces(), Nektar::StdRegions::StdExpansion::GetNverts(), Nektar::StdRegions::StdExpansion::GetTraceIntNcoeffs(), GetTraceInverseBoundaryMap(), Nektar::StdRegions::StdExpansion::GetVertexMap(), Nektar::StdRegions::StdExpansion::NumBndryCoeffs(), and Vmath::Vcopy().

◆ v_BuildTransformationMatrix()

DNekMatSharedPtr Nektar::LocalRegions::Expansion3D::v_BuildTransformationMatrix ( const DNekScalMatSharedPtr r_bnd,
const StdRegions::MatrixType  matrixType 
)
overrideprotectedvirtual

The matrix component of \(\mathbf{R}\) is given by

\[ \mathbf{R^{T}_{v}}= -\mathbf{S}^{-1}_{ef,ef}\mathbf{S}^{T}_{v,ef}\]

For every vertex mode we extract the submatrices from statically condensed matrix \(\mathbf{S}\) corresponding to the coupling between the attached edges and faces of a vertex ( \(\mathbf{S_{ef,ef}}\)). This matrix is then inverted and multiplied by the submatrix representing the coupling between a vertex and the attached edges and faces ( \(\mathbf{S_{v,ef}}\)).

Reimplemented from Nektar::LocalRegions::Expansion.

Definition at line 1990 of file Expansion3D.cpp.

1992{
1993 int nVerts, nEdges;
1994 int eid, fid, vid, n, i;
1995
1996 int nBndCoeffs = NumBndryCoeffs();
1997
1999
2000 // Get geometric information about this element
2001 nVerts = GetNverts();
2002 nEdges = GetNedges();
2003
2004 /*************************************/
2005 /* Vetex-edge & vertex-face matrices */
2006 /*************************************/
2007
2008 /**
2009 * The matrix component of \f$\mathbf{R}\f$ is given by \f[
2010 * \mathbf{R^{T}_{v}}=
2011 * -\mathbf{S}^{-1}_{ef,ef}\mathbf{S}^{T}_{v,ef}\f]
2012 *
2013 * For every vertex mode we extract the submatrices from statically
2014 * condensed matrix \f$\mathbf{S}\f$ corresponding to the coupling
2015 * between the attached edges and faces of a vertex
2016 * (\f$\mathbf{S_{ef,ef}}\f$). This matrix is then inverted and
2017 * multiplied by the submatrix representing the coupling between a
2018 * vertex and the attached edges and faces
2019 * (\f$\mathbf{S_{v,ef}}\f$).
2020 */
2021
2022 int nmodes;
2023 int m;
2024 NekDouble VertexEdgeFaceValue;
2025
2026 // The number of connected edges/faces is 3 (for all elements)
2027 int nConnectedEdges = 3;
2028 int nConnectedFaces = 3;
2029
2030 // Location in the matrix
2031 Array<OneD, Array<OneD, unsigned int>> MatEdgeLocation(nConnectedEdges);
2032 Array<OneD, Array<OneD, unsigned int>> MatFaceLocation(nConnectedFaces);
2033
2034 // Define storage for vertex transpose matrix and zero all entries
2035 MatrixStorage storage = eFULL;
2036 DNekMatSharedPtr transformationmatrix;
2037
2038 transformationmatrix = MemoryManager<DNekMat>::AllocateSharedPtr(
2039 nBndCoeffs, nBndCoeffs, 0.0, storage);
2040
2041 DNekMat &R = (*transformationmatrix);
2042
2043 // Build the vertex-edge/face transform matrix: This matrix is
2044 // constructed from the submatrices corresponding to the couping
2045 // between each vertex and the attached edges/faces
2046 for (vid = 0; vid < nVerts; ++vid)
2047 {
2048 // Row and column size of the vertex-edge/face matrix
2049 int efRow = GetEdgeNcoeffs(geom->GetVertexEdgeMap(vid, 0)) +
2050 GetEdgeNcoeffs(geom->GetVertexEdgeMap(vid, 1)) +
2051 GetEdgeNcoeffs(geom->GetVertexEdgeMap(vid, 2)) +
2052 GetTraceIntNcoeffs(geom->GetVertexFaceMap(vid, 0)) +
2053 GetTraceIntNcoeffs(geom->GetVertexFaceMap(vid, 1)) +
2054 GetTraceIntNcoeffs(geom->GetVertexFaceMap(vid, 2)) - 6;
2055
2056 int nedgemodesconnected =
2057 GetEdgeNcoeffs(geom->GetVertexEdgeMap(vid, 0)) +
2058 GetEdgeNcoeffs(geom->GetVertexEdgeMap(vid, 1)) +
2059 GetEdgeNcoeffs(geom->GetVertexEdgeMap(vid, 2)) - 6;
2060 Array<OneD, unsigned int> edgemodearray(nedgemodesconnected);
2061
2062 int nfacemodesconnected =
2063 GetTraceIntNcoeffs(geom->GetVertexFaceMap(vid, 0)) +
2064 GetTraceIntNcoeffs(geom->GetVertexFaceMap(vid, 1)) +
2065 GetTraceIntNcoeffs(geom->GetVertexFaceMap(vid, 2));
2066 Array<OneD, unsigned int> facemodearray(nfacemodesconnected);
2067
2068 int offset = 0;
2069
2070 // Create array of edge modes
2071 for (eid = 0; eid < nConnectedEdges; ++eid)
2072 {
2073 MatEdgeLocation[eid] =
2074 GetEdgeInverseBoundaryMap(geom->GetVertexEdgeMap(vid, eid));
2075 nmodes = MatEdgeLocation[eid].size();
2076
2077 if (nmodes)
2078 {
2079 Vmath::Vcopy(nmodes, &MatEdgeLocation[eid][0], 1,
2080 &edgemodearray[offset], 1);
2081 }
2082
2083 offset += nmodes;
2084 }
2085
2086 offset = 0;
2087
2088 // Create array of face modes
2089 for (fid = 0; fid < nConnectedFaces; ++fid)
2090 {
2091 MatFaceLocation[fid] =
2092 GetTraceInverseBoundaryMap(geom->GetVertexFaceMap(vid, fid));
2093 nmodes = MatFaceLocation[fid].size();
2094
2095 if (nmodes)
2096 {
2097 Vmath::Vcopy(nmodes, &MatFaceLocation[fid][0], 1,
2098 &facemodearray[offset], 1);
2099 }
2100 offset += nmodes;
2101 }
2102
2103 DNekMatSharedPtr vertexedgefacetransformmatrix =
2104 MemoryManager<DNekMat>::AllocateSharedPtr(1, efRow, 0.0, storage);
2105 DNekMat &Sveft = (*vertexedgefacetransformmatrix);
2106
2107 DNekMatSharedPtr vertexedgefacecoupling =
2108 MemoryManager<DNekMat>::AllocateSharedPtr(1, efRow, 0.0, storage);
2109 DNekMat &Svef = (*vertexedgefacecoupling);
2110
2111 // Vertex-edge coupling
2112 for (n = 0; n < nedgemodesconnected; ++n)
2113 {
2114 // Matrix value for each coefficient location
2115 VertexEdgeFaceValue = (*r_bnd)(GetVertexMap(vid), edgemodearray[n]);
2116
2117 // Set the value in the vertex edge/face matrix
2118 Svef.SetValue(0, n, VertexEdgeFaceValue);
2119 }
2120
2121 // Vertex-face coupling
2122 for (n = 0; n < nfacemodesconnected; ++n)
2123 {
2124 // Matrix value for each coefficient location
2125 VertexEdgeFaceValue = (*r_bnd)(GetVertexMap(vid), facemodearray[n]);
2126
2127 // Set the value in the vertex edge/face matrix
2128 Svef.SetValue(0, n + nedgemodesconnected, VertexEdgeFaceValue);
2129 }
2130
2131 /*
2132 * Build the edge-face transform matrix: This matrix is
2133 * constructed from the submatrices corresponding to the couping
2134 * between the edges and faces on the attached faces/edges of a
2135 * vertex.
2136 */
2137
2138 // Allocation of matrix to store edge/face-edge/face coupling
2139 DNekMatSharedPtr edgefacecoupling =
2141 storage);
2142 DNekMat &Sefef = (*edgefacecoupling);
2143
2144 NekDouble EdgeEdgeValue, FaceFaceValue;
2145
2146 // Edge-edge coupling (S_{ee})
2147 for (m = 0; m < nedgemodesconnected; ++m)
2148 {
2149 for (n = 0; n < nedgemodesconnected; ++n)
2150 {
2151 // Matrix value for each coefficient location
2152 EdgeEdgeValue = (*r_bnd)(edgemodearray[n], edgemodearray[m]);
2153
2154 // Set the value in the vertex edge/face matrix
2155 Sefef.SetValue(n, m, EdgeEdgeValue);
2156 }
2157 }
2158
2159 // Face-face coupling (S_{ff})
2160 for (n = 0; n < nfacemodesconnected; ++n)
2161 {
2162 for (m = 0; m < nfacemodesconnected; ++m)
2163 {
2164 // Matrix value for each coefficient location
2165 FaceFaceValue = (*r_bnd)(facemodearray[n], facemodearray[m]);
2166 // Set the value in the vertex edge/face matrix
2167 Sefef.SetValue(nedgemodesconnected + n, nedgemodesconnected + m,
2168 FaceFaceValue);
2169 }
2170 }
2171
2172 // Edge-face coupling (S_{ef} and trans(S_{ef}))
2173 for (n = 0; n < nedgemodesconnected; ++n)
2174 {
2175 for (m = 0; m < nfacemodesconnected; ++m)
2176 {
2177 // Matrix value for each coefficient location
2178 FaceFaceValue = (*r_bnd)(edgemodearray[n], facemodearray[m]);
2179
2180 // Set the value in the vertex edge/face matrix
2181 Sefef.SetValue(n, nedgemodesconnected + m, FaceFaceValue);
2182
2183 FaceFaceValue = (*r_bnd)(facemodearray[m], edgemodearray[n]);
2184
2185 // and transpose
2186 Sefef.SetValue(nedgemodesconnected + m, n, FaceFaceValue);
2187 }
2188 }
2189
2190 // Invert edge-face coupling matrix
2191 if (efRow)
2192 {
2193 Sefef.Invert();
2194
2195 // R_{v}=-S_{v,ef}inv(S_{ef,ef})
2196 Sveft = -Svef * Sefef;
2197 }
2198
2199 // Populate R with R_{ve} components
2200 for (n = 0; n < edgemodearray.size(); ++n)
2201 {
2202 R.SetValue(GetVertexMap(vid), edgemodearray[n], Sveft(0, n));
2203 }
2204
2205 // Populate R with R_{vf} components
2206 for (n = 0; n < facemodearray.size(); ++n)
2207 {
2208 R.SetValue(GetVertexMap(vid), facemodearray[n],
2209 Sveft(0, n + nedgemodesconnected));
2210 }
2211 }
2212
2213 /********************/
2214 /* edge-face matrix */
2215 /********************/
2216
2217 /*
2218 * The matrix component of \f$\mathbf{R}\f$ is given by \f[
2219 * \mathbf{R^{T}_{ef}}=-\mathbf{S}^{-1}_{ff}\mathbf{S}^{T}_{ef}\f]
2220 *
2221 * For each edge extract the submatrices from statically condensed
2222 * matrix \f$\mathbf{S}\f$ corresponding to inner products of modes
2223 * on the two attached faces within themselves as well as the
2224 * coupling matrix between the two faces
2225 * (\f$\mathbf{S}_{ff}\f$). This matrix of face coupling is then
2226 * inverted and multiplied by the submatrices of corresponding to
2227 * the coupling between the edge and attached faces
2228 * (\f$\mathbf{S}_{ef}\f$).
2229 */
2230
2231 NekDouble EdgeFaceValue, FaceFaceValue;
2232 int efCol, efRow, nedgemodes;
2233
2234 // Number of attached faces is always 2
2235 nConnectedFaces = 2;
2236
2237 // Location in the matrix
2238 MatEdgeLocation = Array<OneD, Array<OneD, unsigned int>>(nEdges);
2239 MatFaceLocation = Array<OneD, Array<OneD, unsigned int>>(nConnectedFaces);
2240
2241 // Build the edge/face transform matrix: This matrix is constructed
2242 // from the submatrices corresponding to the couping between a
2243 // specific edge and the two attached faces.
2244 for (eid = 0; eid < nEdges; ++eid)
2245 {
2246 // Row and column size of the vertex-edge/face matrix
2247 efCol = GetTraceIntNcoeffs(geom->GetEdgeFaceMap(eid, 0)) +
2248 GetTraceIntNcoeffs(geom->GetEdgeFaceMap(eid, 1));
2249 efRow = GetEdgeNcoeffs(eid) - 2;
2250
2251 // Edge-face coupling matrix
2252 DNekMatSharedPtr efedgefacecoupling =
2254 storage);
2255 DNekMat &Mef = (*efedgefacecoupling);
2256
2257 // Face-face coupling matrix
2258 DNekMatSharedPtr effacefacecoupling =
2260 storage);
2261 DNekMat &Meff = (*effacefacecoupling);
2262
2263 // Edge-face transformation matrix
2264 DNekMatSharedPtr edgefacetransformmatrix =
2266 storage);
2267 DNekMat &Meft = (*edgefacetransformmatrix);
2268
2269 int nfacemodesconnected =
2270 GetTraceIntNcoeffs(geom->GetEdgeFaceMap(eid, 0)) +
2271 GetTraceIntNcoeffs(geom->GetEdgeFaceMap(eid, 1));
2272 Array<OneD, unsigned int> facemodearray(nfacemodesconnected);
2273
2274 // Create array of edge modes
2275 Array<OneD, unsigned int> inedgearray = GetEdgeInverseBoundaryMap(eid);
2276 nedgemodes = GetEdgeNcoeffs(eid) - 2;
2277 Array<OneD, unsigned int> edgemodearray(nedgemodes);
2278
2279 if (nedgemodes)
2280 {
2281 Vmath::Vcopy(nedgemodes, &inedgearray[0], 1, &edgemodearray[0], 1);
2282 }
2283
2284 int offset = 0;
2285
2286 // Create array of face modes
2287 for (fid = 0; fid < nConnectedFaces; ++fid)
2288 {
2289 MatFaceLocation[fid] =
2290 GetTraceInverseBoundaryMap(geom->GetEdgeFaceMap(eid, fid));
2291 nmodes = MatFaceLocation[fid].size();
2292
2293 if (nmodes)
2294 {
2295 Vmath::Vcopy(nmodes, &MatFaceLocation[fid][0], 1,
2296 &facemodearray[offset], 1);
2297 }
2298 offset += nmodes;
2299 }
2300
2301 // Edge-face coupling
2302 for (n = 0; n < nedgemodes; ++n)
2303 {
2304 for (m = 0; m < nfacemodesconnected; ++m)
2305 {
2306 // Matrix value for each coefficient location
2307 EdgeFaceValue = (*r_bnd)(edgemodearray[n], facemodearray[m]);
2308
2309 // Set the value in the edge/face matrix
2310 Mef.SetValue(n, m, EdgeFaceValue);
2311 }
2312 }
2313
2314 // Face-face coupling
2315 for (n = 0; n < nfacemodesconnected; ++n)
2316 {
2317 for (m = 0; m < nfacemodesconnected; ++m)
2318 {
2319 // Matrix value for each coefficient location
2320 FaceFaceValue = (*r_bnd)(facemodearray[n], facemodearray[m]);
2321
2322 // Set the value in the vertex edge/face matrix
2323 Meff.SetValue(n, m, FaceFaceValue);
2324 }
2325 }
2326
2327 if (efCol)
2328 {
2329 // Invert edge-face coupling matrix
2330 Meff.Invert();
2331
2332 // trans(R_{ef})=-S_{ef}*(inv(S_{ff})
2333 Meft = -Mef * Meff;
2334 }
2335
2336 // Populate transformation matrix with Meft
2337 for (n = 0; n < Meft.GetRows(); ++n)
2338 {
2339 for (m = 0; m < Meft.GetColumns(); ++m)
2340 {
2341 R.SetValue(edgemodearray[n], facemodearray[m], Meft(n, m));
2342 }
2343 }
2344 }
2345
2346 for (i = 0; i < R.GetRows(); ++i)
2347 {
2348 R.SetValue(i, i, 1.0);
2349 }
2350
2351 if ((matrixType == StdRegions::ePreconR) ||
2352 (matrixType == StdRegions::ePreconRMass))
2353 {
2354 return transformationmatrix;
2355 }
2356 else
2357 {
2358 NEKERROR(ErrorUtil::efatal, "unkown matrix type");
2359 return NullDNekMatSharedPtr;
2360 }
2361}
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition: ErrorUtil.hpp:202
static DNekMatSharedPtr NullDNekMatSharedPtr
Definition: NekTypeDefs.hpp:83

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::ErrorUtil::efatal, Nektar::eFULL, Nektar::StdRegions::ePreconR, Nektar::StdRegions::ePreconRMass, GetEdgeInverseBoundaryMap(), Nektar::StdRegions::StdExpansion3D::GetEdgeNcoeffs(), GetGeom3D(), Nektar::StdRegions::StdExpansion3D::GetNedges(), Nektar::StdRegions::StdExpansion::GetNverts(), Nektar::StdRegions::StdExpansion::GetTraceIntNcoeffs(), GetTraceInverseBoundaryMap(), Nektar::StdRegions::StdExpansion::GetVertexMap(), NEKERROR, Nektar::NullDNekMatSharedPtr, Nektar::StdRegions::StdExpansion::NumBndryCoeffs(), and Vmath::Vcopy().

◆ v_BuildVertexMatrix()

DNekMatSharedPtr Nektar::LocalRegions::Expansion3D::v_BuildVertexMatrix ( const DNekScalMatSharedPtr r_bnd)
overrideprotectedvirtual

Reimplemented from Nektar::LocalRegions::Expansion.

Definition at line 1960 of file Expansion3D.cpp.

1962{
1963 MatrixStorage storage = eFULL;
1964 DNekMatSharedPtr vertexmatrix;
1965
1966 int nVerts, vid1, vid2, vMap1, vMap2;
1967 NekDouble VertexValue;
1968
1969 nVerts = GetNverts();
1970
1971 vertexmatrix =
1972 MemoryManager<DNekMat>::AllocateSharedPtr(nVerts, nVerts, 0.0, storage);
1973 DNekMat &VertexMat = (*vertexmatrix);
1974
1975 for (vid1 = 0; vid1 < nVerts; ++vid1)
1976 {
1977 vMap1 = GetVertexMap(vid1, true);
1978
1979 for (vid2 = 0; vid2 < nVerts; ++vid2)
1980 {
1981 vMap2 = GetVertexMap(vid2, true);
1982 VertexValue = (*r_bnd)(vMap1, vMap2);
1983 VertexMat.SetValue(vid1, vid2, VertexValue);
1984 }
1985 }
1986
1987 return vertexmatrix;
1988}

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::eFULL, Nektar::StdRegions::StdExpansion::GetNverts(), and Nektar::StdRegions::StdExpansion::GetVertexMap().

◆ v_DGDeriv()

void Nektar::LocalRegions::Expansion3D::v_DGDeriv ( const int  dir,
const Array< OneD, const NekDouble > &  incoeffs,
Array< OneD, ExpansionSharedPtr > &  FaceExp,
Array< OneD, Array< OneD, NekDouble > > &  faceCoeffs,
Array< OneD, NekDouble > &  out_d 
)
overrideprotectedvirtual

Evaluate coefficients of weak deriviative in the direction dir given the input coefficicents incoeffs and the imposed boundary values in EdgeExp (which will have its phys space updated).

Reimplemented from Nektar::LocalRegions::Expansion.

Definition at line 1777 of file Expansion3D.cpp.

1782{
1783 int ncoeffs = GetNcoeffs();
1787
1789 DNekScalMat &Dmat = *GetLocMatrix(DerivType[dir]);
1790
1791 Array<OneD, NekDouble> coeffs = incoeffs;
1792 DNekVec Coeffs(ncoeffs, coeffs, eWrapper);
1793
1794 Coeffs = Transpose(Dmat) * Coeffs;
1795 Vmath::Neg(ncoeffs, coeffs, 1);
1796
1797 // Add the boundary integral including the relevant part of
1798 // the normal
1799 AddNormTraceInt(dir, FaceExp, faceCoeffs, coeffs);
1800
1801 DNekVec Out_d(ncoeffs, out_d, eWrapper);
1802
1803 Out_d = InvMass * Coeffs;
1804}
void AddNormTraceInt(const int dir, Array< OneD, ExpansionSharedPtr > &FaceExp, Array< OneD, Array< OneD, NekDouble > > &faceCoeffs, Array< OneD, NekDouble > &outarray)
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.hpp:292

References AddNormTraceInt(), Nektar::StdRegions::eInvMass, Nektar::StdRegions::eWeakDeriv0, Nektar::StdRegions::eWeakDeriv1, Nektar::StdRegions::eWeakDeriv2, Nektar::eWrapper, Nektar::LocalRegions::Expansion::GetLocMatrix(), Nektar::StdRegions::StdExpansion::GetNcoeffs(), Vmath::Neg(), and Nektar::Transpose().

◆ v_GenMatrix()

DNekMatSharedPtr Nektar::LocalRegions::Expansion3D::v_GenMatrix ( const StdRegions::StdMatrixKey mkey)
overrideprotectedvirtual

Computes matrices needed for the HDG formulation. References to equations relate to the following paper (with a suitable changes in formulation to adapt to 3D): R. M. Kirby, S. J. Sherwin, B. Cockburn, To CG or to HDG: A Comparative Study, J. Sci. Comp P1-30 DOI 10.1007/s10915-011-9501-7 NOTE: VARIABLE COEFFICIENTS CASE IS NOT IMPLEMENTED

Reimplemented from Nektar::StdRegions::StdExpansion.

Reimplemented in Nektar::LocalRegions::HexExp, Nektar::LocalRegions::PrismExp, Nektar::LocalRegions::PyrExp, and Nektar::LocalRegions::TetExp.

Definition at line 882 of file Expansion3D.cpp.

883{
884 // Variable coefficients are not implemented/////////
885 ASSERTL1(!mkey.HasVarCoeff(StdRegions::eVarCoeffD00),
886 "Matrix construction is not implemented for variable "
887 "coefficients at the moment");
888 ////////////////////////////////////////////////////
889
890 DNekMatSharedPtr returnval;
891
892 switch (mkey.GetMatrixType())
893 {
894 // (Z^e)^{-1} (Eqn. 33, P22)
896 {
898 "HybridDGHelmholtz matrix not set up "
899 "for non boundary-interior expansions");
900
901 int i, j, k;
902 NekDouble lambdaval =
903 mkey.GetConstFactor(StdRegions::eFactorLambda);
904 NekDouble tau = mkey.GetConstFactor(StdRegions::eFactorTau);
905 int ncoeffs = GetNcoeffs();
906 int nfaces = GetNtraces();
907
908 Array<OneD, unsigned int> fmap;
909 Array<OneD, int> sign;
910 ExpansionSharedPtr FaceExp;
911 ExpansionSharedPtr FaceExp2;
912
913 int order_f, coordim = GetCoordim();
918
919 returnval =
921 DNekMat &Mat = *returnval;
922 Vmath::Zero(ncoeffs * ncoeffs, Mat.GetPtr(), 1);
923
924 // StdRegions::VarCoeffType Coeffs[3] = {StdRegions::eVarCoeffD00,
925 // StdRegions::eVarCoeffD11,
926 // StdRegions::eVarCoeffD22};
927 StdRegions::VarCoeffMap::const_iterator x;
928 const StdRegions::VarCoeffMap &varcoeffs = mkey.GetVarCoeffs();
929
930 for (i = 0; i < coordim; ++i)
931 {
932 if ((x = varcoeffs.find(StdRegions::eVarCoeffMF1x)) !=
933 varcoeffs.end())
934 {
935 StdRegions::VarCoeffMap VarCoeffDirDeriv;
936 VarCoeffDirDeriv[StdRegions::eVarCoeffMF] =
937 GetMF(i, coordim, varcoeffs);
938 VarCoeffDirDeriv[StdRegions::eVarCoeffMFDiv] =
939 GetMFDiv(i, varcoeffs);
940
941 MatrixKey Dmatkey(StdRegions::eWeakDirectionalDeriv,
942 DetShapeType(), *this,
944 VarCoeffDirDeriv);
945
946 DNekScalMat &Dmat = *GetLocMatrix(Dmatkey);
947
950 GetMFMag(i, mkey.GetVarCoeffs());
951
952 MatrixKey invMasskey(StdRegions::eInvMass, DetShapeType(),
954 Weight);
955
956 DNekScalMat &invMass = *GetLocMatrix(invMasskey);
957
958 Mat = Mat + Dmat * invMass * Transpose(Dmat);
959 }
960 else
961 {
962 DNekScalMat &Dmat = *GetLocMatrix(DerivType[i]);
963 Mat = Mat + Dmat * invMass * Transpose(Dmat);
964 }
965
966 /*
967 if(mkey.HasVarCoeff(Coeffs[i]))
968 {
969 MatrixKey DmatkeyL(DerivType[i], DetExpansionType(), *this,
970 StdRegions::NullConstFactorMap,
971 mkey.GetVarCoeffAsMap(Coeffs[i]));
972 MatrixKey DmatkeyR(DerivType[i], DetExpansionType(), *this);
973
974 DNekScalMat &DmatL = *GetLocMatrix(DmatkeyL);
975 DNekScalMat &DmatR = *GetLocMatrix(DmatkeyR);
976 Mat = Mat + DmatL*invMass*Transpose(DmatR);
977 }
978 else
979 {
980 DNekScalMat &Dmat = *GetLocMatrix(DerivType[i]);
981 Mat = Mat + Dmat*invMass*Transpose(Dmat);
982 }
983 */
984 }
985
986 // Add Mass Matrix Contribution for Helmholtz problem
988 Mat = Mat + lambdaval * Mass;
989
990 // Add tau*E_l using elemental mass matrices on each edge
991 for (i = 0; i < nfaces; ++i)
992 {
993 FaceExp = GetTraceExp(i);
994 order_f = FaceExp->GetNcoeffs();
995
996 IndexMapKey ikey(eFaceToElement, DetShapeType(),
1000
1001 // @TODO: Document
1002 /*
1003 StdRegions::VarCoeffMap edgeVarCoeffs;
1004 if (mkey.HasVarCoeff(StdRegions::eVarCoeffD00))
1005 {
1006 Array<OneD, NekDouble> mu(nq);
1007 GetPhysEdgeVarCoeffsFromElement(
1008 i, EdgeExp2,
1009 mkey.GetVarCoeff(StdRegions::eVarCoeffD00), mu);
1010 edgeVarCoeffs[StdRegions::eVarCoeffMass] = mu;
1011 }
1012 DNekScalMat &eMass = *EdgeExp->GetLocMatrix(
1013 StdRegions::eMass,
1014 StdRegions::NullConstFactorMap, edgeVarCoeffs);
1015 */
1016
1017 DNekScalMat &eMass = *FaceExp->GetLocMatrix(StdRegions::eMass);
1018
1019 for (j = 0; j < order_f; ++j)
1020 {
1021 for (k = 0; k < order_f; ++k)
1022 {
1023 Mat((*map)[j].index, (*map)[k].index) +=
1024 tau * (*map)[j].sign * (*map)[k].sign * eMass(j, k);
1025 }
1026 }
1027 }
1028 break;
1029 }
1030
1031 // U^e (P22)
1033 {
1034 int i, j, k;
1035 int nbndry = NumDGBndryCoeffs();
1036 int ncoeffs = GetNcoeffs();
1037 int nfaces = GetNtraces();
1038 NekDouble tau = mkey.GetConstFactor(StdRegions::eFactorTau);
1039
1040 Array<OneD, NekDouble> lambda(nbndry);
1041 DNekVec Lambda(nbndry, lambda, eWrapper);
1042 Array<OneD, NekDouble> ulam(ncoeffs);
1043 DNekVec Ulam(ncoeffs, ulam, eWrapper);
1044 Array<OneD, NekDouble> f(ncoeffs);
1045 DNekVec F(ncoeffs, f, eWrapper);
1046
1047 ExpansionSharedPtr FaceExp;
1048 // declare matrix space
1049 returnval =
1051 DNekMat &Umat = *returnval;
1052
1053 // Z^e matrix
1055 *this, mkey.GetConstFactors(),
1056 mkey.GetVarCoeffs());
1057 DNekScalMat &invHmat = *GetLocMatrix(newkey);
1058
1059 Array<OneD, unsigned int> fmap;
1060 Array<OneD, int> sign;
1061
1062 // alternative way to add boundary terms contribution
1063 int bndry_cnt = 0;
1064 for (i = 0; i < nfaces; ++i)
1065 {
1066 FaceExp = GetTraceExp(
1067 i); // temporary, need to rewrite AddHDGHelmholtzFaceTerms
1068 int nface = GetTraceNcoeffs(i);
1069 Array<OneD, NekDouble> face_lambda(nface);
1070
1071 const Array<OneD, const Array<OneD, NekDouble>> normals =
1072 GetTraceNormal(i);
1073
1074 for (j = 0; j < nface; ++j)
1075 {
1076 Vmath::Zero(nface, &face_lambda[0], 1);
1077 Vmath::Zero(ncoeffs, &f[0], 1);
1078 face_lambda[j] = 1.0;
1079
1080 SetFaceToGeomOrientation(i, face_lambda);
1081
1082 Array<OneD, NekDouble> tmp(FaceExp->GetTotPoints());
1083 FaceExp->BwdTrans(face_lambda, tmp);
1084 AddHDGHelmholtzFaceTerms(tau, i, tmp, mkey.GetVarCoeffs(),
1085 f);
1086
1087 Ulam = invHmat * F; // generate Ulam from lambda
1088
1089 // fill column of matrix
1090 for (k = 0; k < ncoeffs; ++k)
1091 {
1092 Umat(k, bndry_cnt) = Ulam[k];
1093 }
1094
1095 ++bndry_cnt;
1096 }
1097 }
1098
1099 //// Set up face expansions from local geom info
1100 // for(i = 0; i < nfaces; ++i)
1101 //{
1102 // FaceExp[i] = GetTraceExp(i);
1103 //}
1104 //
1105 //// for each degree of freedom of the lambda space
1106 //// calculate Umat entry
1107 //// Generate Lambda to U_lambda matrix
1108 // for(j = 0; j < nbndry; ++j)
1109 //{
1110 // // standard basis vectors e_j
1111 // Vmath::Zero(nbndry,&lambda[0],1);
1112 // Vmath::Zero(ncoeffs,&f[0],1);
1113 // lambda[j] = 1.0;
1114 //
1115 // //cout << Lambda;
1116 // SetTraceToGeomOrientation(lambda);
1117 // //cout << Lambda << endl;
1118 //
1119 // // Compute F = [I D_1 M^{-1} D_2 M^{-1}] C e_j
1120 // AddHDGHelmholtzTraceTerms(tau, lambda, FaceExp,
1121 // mkey.GetVarCoeffs(), f);
1122 //
1123 // // Compute U^e_j
1124 // Ulam = invHmat*F; // generate Ulam from lambda
1125 //
1126 // // fill column of matrix
1127 // for(k = 0; k < ncoeffs; ++k)
1128 // {
1129 // Umat(k,j) = Ulam[k];
1130 // }
1131 //}
1132 }
1133 break;
1134 // Q_0, Q_1, Q_2 matrices (P23)
1135 // Each are a product of a row of Eqn 32 with the C matrix.
1136 // Rather than explicitly computing all of Eqn 32, we note each
1137 // row is almost a multiple of U^e, so use that as our starting
1138 // point.
1142 {
1143 int i = 0;
1144 int j = 0;
1145 int k = 0;
1146 int dir = 0;
1147 int nbndry = NumDGBndryCoeffs();
1148 int coordim = GetCoordim();
1149 int ncoeffs = GetNcoeffs();
1150 int nfaces = GetNtraces();
1151
1152 Array<OneD, NekDouble> lambda(nbndry);
1153 DNekVec Lambda(nbndry, lambda, eWrapper);
1154 Array<OneD, ExpansionSharedPtr> FaceExp(nfaces);
1155
1156 Array<OneD, NekDouble> ulam(ncoeffs);
1157 DNekVec Ulam(ncoeffs, ulam, eWrapper);
1158 Array<OneD, NekDouble> f(ncoeffs);
1159 DNekVec F(ncoeffs, f, eWrapper);
1160
1161 // declare matrix space
1162 returnval =
1164 DNekMat &Qmat = *returnval;
1165
1166 // Lambda to U matrix
1167 MatrixKey lamToUkey(StdRegions::eHybridDGLamToU, DetShapeType(),
1168 *this, mkey.GetConstFactors(),
1169 mkey.GetVarCoeffs());
1170 DNekScalMat &lamToU = *GetLocMatrix(lamToUkey);
1171
1172 // Inverse mass matrix
1174
1175 for (i = 0; i < nfaces; ++i)
1176 {
1177 FaceExp[i] = GetTraceExp(i);
1178 }
1179
1180 // Weak Derivative matrix
1182 switch (mkey.GetMatrixType())
1183 {
1185 dir = 0;
1187 break;
1189 dir = 1;
1191 break;
1193 dir = 2;
1195 break;
1196 default:
1197 ASSERTL0(false, "Direction not known");
1198 break;
1199 }
1200
1201 // DNekScalMatSharedPtr Dmat;
1202 // DNekScalMatSharedPtr &invMass;
1203 StdRegions::VarCoeffMap::const_iterator x;
1204 const StdRegions::VarCoeffMap &varcoeffs = mkey.GetVarCoeffs();
1205 if ((x = varcoeffs.find(StdRegions::eVarCoeffMF1x)) !=
1206 varcoeffs.end())
1207 {
1208 StdRegions::VarCoeffMap VarCoeffDirDeriv;
1209 VarCoeffDirDeriv[StdRegions::eVarCoeffMF] =
1210 GetMF(dir, coordim, varcoeffs);
1211 VarCoeffDirDeriv[StdRegions::eVarCoeffMFDiv] =
1212 GetMFDiv(dir, varcoeffs);
1213
1214 MatrixKey Dmatkey(
1216 StdRegions::NullConstFactorMap, VarCoeffDirDeriv);
1217
1218 Dmat = GetLocMatrix(Dmatkey);
1219
1222 GetMFMag(dir, mkey.GetVarCoeffs());
1223
1224 MatrixKey invMasskey(StdRegions::eInvMass, DetShapeType(),
1226 Weight);
1227
1228 invMass = *GetLocMatrix(invMasskey);
1229 }
1230 else
1231 {
1232 // Inverse mass matrix
1234 }
1235
1236 // for each degree of freedom of the lambda space
1237 // calculate Qmat entry
1238 // Generate Lambda to Q_lambda matrix
1239 for (j = 0; j < nbndry; ++j)
1240 {
1241 Vmath::Zero(nbndry, &lambda[0], 1);
1242 lambda[j] = 1.0;
1243
1244 // for lambda[j] = 1 this is the solution to ulam
1245 for (k = 0; k < ncoeffs; ++k)
1246 {
1247 Ulam[k] = lamToU(k, j);
1248 }
1249
1250 // -D^T ulam
1251 Vmath::Neg(ncoeffs, &ulam[0], 1);
1252 F = Transpose(*Dmat) * Ulam;
1253
1255
1256 // Add the C terms resulting from the I's on the
1257 // diagonals of Eqn 32
1258 AddNormTraceInt(dir, lambda, FaceExp, f, mkey.GetVarCoeffs());
1259
1260 // finally multiply by inverse mass matrix
1261 Ulam = invMass * F;
1262
1263 // fill column of matrix (Qmat is in column major format)
1264 Vmath::Vcopy(ncoeffs, &ulam[0], 1,
1265 &(Qmat.GetPtr())[0] + j * ncoeffs, 1);
1266 }
1267 }
1268 break;
1269 // Matrix K (P23)
1271 {
1272 int i, j, f, cnt;
1273 int order_f, nquad_f;
1274 int nbndry = NumDGBndryCoeffs();
1275 int nfaces = GetNtraces();
1276 NekDouble tau = mkey.GetConstFactor(StdRegions::eFactorTau);
1277
1278 Array<OneD, NekDouble> work, varcoeff_work;
1279 Array<OneD, const Array<OneD, NekDouble>> normals;
1280 Array<OneD, ExpansionSharedPtr> FaceExp(nfaces);
1281 Array<OneD, NekDouble> lam(nbndry);
1282
1283 Array<OneD, unsigned int> fmap;
1284 Array<OneD, int> sign;
1285
1286 // declare matrix space
1287 returnval =
1289 DNekMat &BndMat = *returnval;
1290
1291 DNekScalMatSharedPtr LamToQ[3];
1292
1293 // Matrix to map Lambda to U
1294 MatrixKey LamToUkey(StdRegions::eHybridDGLamToU, DetShapeType(),
1295 *this, mkey.GetConstFactors(),
1296 mkey.GetVarCoeffs());
1297 DNekScalMat &LamToU = *GetLocMatrix(LamToUkey);
1298
1299 // Matrix to map Lambda to Q0
1300 MatrixKey LamToQ0key(StdRegions::eHybridDGLamToQ0, DetShapeType(),
1301 *this, mkey.GetConstFactors(),
1302 mkey.GetVarCoeffs());
1303 LamToQ[0] = GetLocMatrix(LamToQ0key);
1304
1305 // Matrix to map Lambda to Q1
1306 MatrixKey LamToQ1key(StdRegions::eHybridDGLamToQ1, DetShapeType(),
1307 *this, mkey.GetConstFactors(),
1308 mkey.GetVarCoeffs());
1309 LamToQ[1] = GetLocMatrix(LamToQ1key);
1310
1311 // Matrix to map Lambda to Q2
1312 MatrixKey LamToQ2key(StdRegions::eHybridDGLamToQ2, DetShapeType(),
1313 *this, mkey.GetConstFactors(),
1314 mkey.GetVarCoeffs());
1315 LamToQ[2] = GetLocMatrix(LamToQ2key);
1316
1317 // Set up edge segment expansions from local geom info
1318 const StdRegions::VarCoeffMap &varcoeffs = mkey.GetVarCoeffs();
1319 for (i = 0; i < nfaces; ++i)
1320 {
1321 FaceExp[i] = GetTraceExp(i);
1322 }
1323
1324 // Set up matrix derived from <mu, Q_lam.n - \tau (U_lam - Lam) >
1325 for (i = 0; i < nbndry; ++i)
1326 {
1327 cnt = 0;
1328
1329 Vmath::Zero(nbndry, lam, 1);
1330 lam[i] = 1.0;
1332
1333 for (f = 0; f < nfaces; ++f)
1334 {
1335 order_f = FaceExp[f]->GetNcoeffs();
1336 nquad_f = FaceExp[f]->GetNumPoints(0) *
1337 FaceExp[f]->GetNumPoints(1);
1338 normals = GetTraceNormal(f);
1339
1340 work = Array<OneD, NekDouble>(nquad_f);
1341 varcoeff_work = Array<OneD, NekDouble>(nquad_f);
1342
1343 IndexMapKey ikey(eFaceToElement, DetShapeType(),
1347
1348 // @TODO Variable coefficients
1349 /*
1350 StdRegions::VarCoeffType VarCoeff[3] =
1351 {StdRegions::eVarCoeffD00, StdRegions::eVarCoeffD11,
1352 StdRegions::eVarCoeffD22};
1353 const StdRegions::VarCoeffMap &varcoeffs =
1354 mkey.GetVarCoeffs(); StdRegions::VarCoeffMap::const_iterator
1355 x;
1356 */
1357
1358 // Q0 * n0 (BQ_0 terms)
1359 Array<OneD, NekDouble> faceCoeffs(order_f);
1360 Array<OneD, NekDouble> facePhys(nquad_f);
1361 for (j = 0; j < order_f; ++j)
1362 {
1363 faceCoeffs[j] =
1364 (*map)[j].sign * (*LamToQ[0])((*map)[j].index, i);
1365 }
1366
1367 FaceExp[f]->BwdTrans(faceCoeffs, facePhys);
1368
1369 // @TODO Variable coefficients
1370 // Multiply by variable coefficient
1371 /*
1372 if ((x = varcoeffs.find(VarCoeff[0])) != varcoeffs.end())
1373 {
1374 GetPhysEdgeVarCoeffsFromElement(e,EdgeExp[e],x->second,varcoeff_work);
1375 Vmath::Vmul(nquad_e,varcoeff_work,1,EdgeExp[e]->GetPhys(),1,EdgeExp[e]->UpdatePhys(),1);
1376 }
1377 */
1378
1379 if (varcoeffs.find(StdRegions::eVarCoeffMF1x) !=
1380 varcoeffs.end())
1381 {
1382 Array<OneD, NekDouble> ncdotMF = GetnFacecdotMF(
1383 0, f, FaceExp[f], normals, varcoeffs);
1384
1385 Vmath::Vmul(nquad_f, ncdotMF, 1, facePhys, 1, work, 1);
1386 }
1387 else
1388 {
1389 Vmath::Vmul(nquad_f, normals[0], 1, facePhys, 1, work,
1390 1);
1391 }
1392
1393 // Q1 * n1 (BQ_1 terms)
1394 for (j = 0; j < order_f; ++j)
1395 {
1396 faceCoeffs[j] =
1397 (*map)[j].sign * (*LamToQ[1])((*map)[j].index, i);
1398 }
1399
1400 FaceExp[f]->BwdTrans(faceCoeffs, facePhys);
1401
1402 // @TODO Variable coefficients
1403 // Multiply by variable coefficients
1404 /*
1405 if ((x = varcoeffs.find(VarCoeff[1])) != varcoeffs.end())
1406 {
1407 GetPhysEdgeVarCoeffsFromElement(e,EdgeExp[e],x->second,varcoeff_work);
1408 Vmath::Vmul(nquad_e,varcoeff_work,1,EdgeExp[e]->GetPhys(),1,EdgeExp[e]->UpdatePhys(),1);
1409 }
1410 */
1411
1412 if ((varcoeffs.find(StdRegions::eVarCoeffMF1x)) !=
1413 varcoeffs.end())
1414 {
1415 Array<OneD, NekDouble> ncdotMF = GetnFacecdotMF(
1416 1, f, FaceExp[f], normals, varcoeffs);
1417
1418 Vmath::Vvtvp(nquad_f, ncdotMF, 1, facePhys, 1, work, 1,
1419 work, 1);
1420 }
1421 else
1422 {
1423 Vmath::Vvtvp(nquad_f, normals[1], 1, facePhys, 1, work,
1424 1, work, 1);
1425 }
1426
1427 // Q2 * n2 (BQ_2 terms)
1428 for (j = 0; j < order_f; ++j)
1429 {
1430 faceCoeffs[j] =
1431 (*map)[j].sign * (*LamToQ[2])((*map)[j].index, i);
1432 }
1433
1434 FaceExp[f]->BwdTrans(faceCoeffs, facePhys);
1435
1436 // @TODO Variable coefficients
1437 // Multiply by variable coefficients
1438 /*
1439 if ((x = varcoeffs.find(VarCoeff[2])) != varcoeffs.end())
1440 {
1441 GetPhysEdgeVarCoeffsFromElement(e,EdgeExp[e],x->second,varcoeff_work);
1442 Vmath::Vmul(nquad_e,varcoeff_work,1,EdgeExp[e]->GetPhys(),1,EdgeExp[e]->UpdatePhys(),1);
1443 }
1444 */
1445
1446 if (varcoeffs.find(StdRegions::eVarCoeffMF1x) !=
1447 varcoeffs.end())
1448 {
1449 Array<OneD, NekDouble> ncdotMF = GetnFacecdotMF(
1450 2, f, FaceExp[f], normals, varcoeffs);
1451
1452 Vmath::Vvtvp(nquad_f, ncdotMF, 1, facePhys, 1, work, 1,
1453 work, 1);
1454 }
1455 else
1456 {
1457 Vmath::Vvtvp(nquad_f, normals[2], 1, facePhys, 1, work,
1458 1, work, 1);
1459 }
1460
1461 // - tau (ulam - lam)
1462 // Corresponds to the G and BU terms.
1463 for (j = 0; j < order_f; ++j)
1464 {
1465 faceCoeffs[j] =
1466 (*map)[j].sign * LamToU((*map)[j].index, i) -
1467 lam[cnt + j];
1468 }
1469
1470 FaceExp[f]->BwdTrans(faceCoeffs, facePhys);
1471
1472 // @TODO Variable coefficients
1473 // Multiply by variable coefficients
1474 /*
1475 if ((x = varcoeffs.find(VarCoeff[0])) != varcoeffs.end())
1476 {
1477 GetPhysEdgeVarCoeffsFromElement(e,FaceExp[f],x->second,varcoeff_work);
1478 Vmath::Vmul(nquad_f,varcoeff_work,1,FaceExp[f]->GetPhys(),1,FaceExp[f]->UpdatePhys(),1);
1479 }
1480 */
1481
1482 Vmath::Svtvp(nquad_f, -tau, facePhys, 1, work, 1, work, 1);
1483
1484 // @TODO Add variable coefficients
1485 FaceExp[f]->IProductWRTBase(work, faceCoeffs);
1486
1487 SetFaceToGeomOrientation(f, faceCoeffs);
1488
1489 for (j = 0; j < order_f; ++j)
1490 {
1491 BndMat(cnt + j, i) = faceCoeffs[j];
1492 }
1493
1494 cnt += order_f;
1495 }
1496 }
1497 }
1498 break;
1499 // HDG postprocessing
1501 {
1502 MatrixKey lapkey(StdRegions::eLaplacian, DetShapeType(), *this,
1503 mkey.GetConstFactors(), mkey.GetVarCoeffs());
1504 DNekScalMat &LapMat = *GetLocMatrix(lapkey);
1505
1507 LapMat.GetRows(), LapMat.GetColumns());
1508 DNekMatSharedPtr lmat = returnval;
1509
1510 (*lmat) = LapMat;
1511
1512 // replace first column with inner product wrt 1
1513 int nq = GetTotPoints();
1514 Array<OneD, NekDouble> tmp(nq);
1515 Array<OneD, NekDouble> outarray(m_ncoeffs);
1516 Vmath::Fill(nq, 1.0, tmp, 1);
1517 IProductWRTBase(tmp, outarray);
1518
1519 Vmath::Vcopy(m_ncoeffs, &outarray[0], 1, &(lmat->GetPtr())[0], 1);
1520
1521 lmat->Invert();
1522 }
1523 break;
1525 {
1526 int ntraces = GetNtraces();
1527 int ncoords = GetCoordim();
1528 int nphys = GetTotPoints();
1529 Array<OneD, const Array<OneD, NekDouble>> normals;
1530 Array<OneD, NekDouble> phys(nphys);
1531 returnval =
1533 DNekMat &Mat = *returnval;
1534 Vmath::Zero(m_ncoeffs * m_ncoeffs, Mat.GetPtr(), 1);
1535
1536 Array<OneD, Array<OneD, NekDouble>> Deriv(3, NullNekDouble1DArray);
1537
1538 for (int d = 0; d < ncoords; ++d)
1539 {
1540 Deriv[d] = Array<OneD, NekDouble>(nphys);
1541 }
1542
1543 Array<OneD, int> tracepts(ntraces);
1544 Array<OneD, ExpansionSharedPtr> traceExp(ntraces);
1545 int maxtpts = 0;
1546 for (int t = 0; t < ntraces; ++t)
1547 {
1548 traceExp[t] = GetTraceExp(t);
1549 tracepts[t] = traceExp[t]->GetTotPoints();
1550 maxtpts = (maxtpts > tracepts[t]) ? maxtpts : tracepts[t];
1551 }
1552
1553 Array<OneD, NekDouble> val(maxtpts), tmp, tmp1;
1554
1555 Array<OneD, Array<OneD, NekDouble>> dphidn(ntraces);
1556 for (int t = 0; t < ntraces; ++t)
1557 {
1558 dphidn[t] =
1559 Array<OneD, NekDouble>(m_ncoeffs * tracepts[t], 0.0);
1560 }
1561
1562 for (int i = 0; i < m_ncoeffs; ++i)
1563 {
1564 FillMode(i, phys);
1565 PhysDeriv(phys, Deriv[0], Deriv[1], Deriv[2]);
1566
1567 for (int t = 0; t < ntraces; ++t)
1568 {
1569 const NormalVector norm = GetTraceNormal(t);
1570
1571 LibUtilities::BasisKey fromkey0 = GetTraceBasisKey(t, 0);
1572 LibUtilities::BasisKey fromkey1 = GetTraceBasisKey(t, 1);
1573 LibUtilities::BasisKey tokey0 =
1574 traceExp[t]->GetBasis(0)->GetBasisKey();
1575 LibUtilities::BasisKey tokey1 =
1576 traceExp[t]->GetBasis(1)->GetBasisKey();
1577 bool DoInterp =
1578 (fromkey0 != tokey0) || (fromkey1 != tokey1);
1579
1580 // for variable p need add check and
1581 // interpolation here.
1582
1583 Array<OneD, NekDouble> n(tracepts[t]);
1584 ;
1585 for (int d = 0; d < ncoords; ++d)
1586 {
1587 // if variable p may need to interpolate
1588 if (DoInterp)
1589 {
1590 LibUtilities::Interp2D(fromkey0, fromkey1, norm[d],
1591 tokey0, tokey1, n);
1592 }
1593 else
1594 {
1595 n = norm[d];
1596 }
1597
1598 GetTracePhysVals(t, traceExp[t], Deriv[d], val,
1599 v_GetTraceOrient(t));
1600
1601 Vmath::Vvtvp(tracepts[t], n, 1, val, 1,
1602 tmp = dphidn[t] + i * tracepts[t], 1,
1603 tmp1 = dphidn[t] + i * tracepts[t], 1);
1604 }
1605 }
1606 }
1607
1608 for (int t = 0; t < ntraces; ++t)
1609 {
1610 int nt = tracepts[t];
1611 NekDouble h, p;
1612 TraceNormLen(t, h, p);
1613
1614 // scaling from GJP paper
1615 NekDouble scale =
1616 (p == 1) ? 0.02 * h * h : 0.8 * pow(p + 1, -4.0) * h * h;
1617
1618 for (int i = 0; i < m_ncoeffs; ++i)
1619 {
1620 for (int j = i; j < m_ncoeffs; ++j)
1621 {
1622 Vmath::Vmul(nt, dphidn[t] + i * nt, 1,
1623 dphidn[t] + j * nt, 1, val, 1);
1624 Mat(i, j) =
1625 Mat(i, j) + scale * traceExp[t]->Integral(val);
1626 }
1627 }
1628 }
1629
1630 // fill in symmetric components.
1631 for (int i = 0; i < m_ncoeffs; ++i)
1632 {
1633 for (int j = 0; j < i; ++j)
1634 {
1635 Mat(i, j) = Mat(j, i);
1636 }
1637 }
1638 }
1639 break;
1640 default:
1641 ASSERTL0(false,
1642 "This matrix type cannot be generated from this class");
1643 break;
1644 }
1645
1646 return returnval;
1647}
void SetTraceToGeomOrientation(Array< OneD, NekDouble > &inout)
Align trace orientation with the geometry orientation.
StdRegions::Orientation v_GetTraceOrient(int face) override
void AddHDGHelmholtzFaceTerms(const NekDouble tau, const int edge, Array< OneD, NekDouble > &facePhys, const StdRegions::VarCoeffMap &dirForcing, Array< OneD, NekDouble > &outarray)
Definition: Expansion3D.cpp:50
void GetTracePhysVals(const int trace, const StdRegions::StdExpansionSharedPtr &TraceExp, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, StdRegions::Orientation orient=StdRegions::eNoOrientation)
Definition: Expansion.h:200
void TraceNormLen(const int traceid, NekDouble &h, NekDouble &p)
Definition: Expansion.h:250
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:134
void FillMode(const int mode, Array< OneD, NekDouble > &outarray)
This function fills the array outarray with the mode-th mode of the expansion.
Definition: StdExpansion.h:491
const LibUtilities::BasisKey GetTraceBasisKey(const int i, int k=-1) const
This function returns the basis key belonging to the i-th trace.
Definition: StdExpansion.h:299
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 expa...
Definition: StdExpansion.h:528
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)
Definition: StdExpansion.h:849
void Interp2D(const BasisKey &fbasis0, const BasisKey &fbasis1, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, Array< OneD, NekDouble > &to)
this function interpolates a 2D function evaluated at the quadrature points of the 2D basis,...
Definition: Interp.cpp:101
Array< OneD, Array< OneD, NekDouble > > NormalVector
Definition: Expansion.h:53
std::vector< double > d(NPUPPER *NPUPPER)
static Array< OneD, NekDouble > NullNekDouble1DArray
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.hpp:273
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.hpp:54

References AddHDGHelmholtzFaceTerms(), AddNormTraceInt(), Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, ASSERTL1, Nektar::UnitTests::d(), Nektar::StdRegions::StdExpansion::DetShapeType(), Nektar::LocalRegions::eFaceToElement, Nektar::StdRegions::eFactorLambda, Nektar::StdRegions::eFactorTau, Nektar::StdRegions::eHybridDGHelmBndLam, Nektar::StdRegions::eHybridDGHelmholtz, Nektar::StdRegions::eHybridDGLamToQ0, Nektar::StdRegions::eHybridDGLamToQ1, Nektar::StdRegions::eHybridDGLamToQ2, Nektar::StdRegions::eHybridDGLamToU, Nektar::StdRegions::eInvHybridDGHelmholtz, Nektar::StdRegions::eInvLaplacianWithUnityMean, Nektar::StdRegions::eInvMass, Nektar::StdRegions::eLaplacian, Nektar::StdRegions::eMass, Nektar::StdRegions::eNormDerivOnTrace, Nektar::StdRegions::eVarCoeffD00, Nektar::StdRegions::eVarCoeffMass, Nektar::StdRegions::eVarCoeffMF, Nektar::StdRegions::eVarCoeffMF1x, Nektar::StdRegions::eVarCoeffMFDiv, Nektar::StdRegions::eWeakDeriv0, Nektar::StdRegions::eWeakDeriv1, Nektar::StdRegions::eWeakDeriv2, Nektar::StdRegions::eWeakDirectionalDeriv, Nektar::eWrapper, Vmath::Fill(), Nektar::StdRegions::StdExpansion::FillMode(), Nektar::StdRegions::StdExpansion::GetBasisNumModes(), Nektar::StdRegions::StdMatrixKey::GetConstFactor(), Nektar::StdRegions::StdMatrixKey::GetConstFactors(), Nektar::StdRegions::StdExpansion::GetCoordim(), Nektar::LocalRegions::Expansion::GetIndexMap(), Nektar::LocalRegions::Expansion::GetLocMatrix(), Nektar::StdRegions::StdMatrixKey::GetMatrixType(), Nektar::LocalRegions::Expansion::GetMF(), Nektar::LocalRegions::Expansion::GetMFDiv(), Nektar::LocalRegions::Expansion::GetMFMag(), Nektar::StdRegions::StdExpansion::GetNcoeffs(), GetnFacecdotMF(), Nektar::StdRegions::StdExpansion::GetNtraces(), Nektar::StdRegions::StdExpansion::GetTotPoints(), Nektar::StdRegions::StdExpansion::GetTraceBasisKey(), Nektar::LocalRegions::Expansion::GetTraceExp(), Nektar::StdRegions::StdExpansion::GetTraceNcoeffs(), Nektar::LocalRegions::Expansion::GetTraceNormal(), Nektar::LocalRegions::Expansion::GetTraceOrient(), Nektar::LocalRegions::Expansion::GetTracePhysVals(), Nektar::StdRegions::StdMatrixKey::GetVarCoeffs(), Nektar::StdRegions::StdMatrixKey::HasVarCoeff(), Nektar::LibUtilities::Interp2D(), Nektar::StdRegions::StdExpansion::IProductWRTBase(), Nektar::StdRegions::StdExpansion::IsBoundaryInteriorExpansion(), Nektar::StdRegions::StdExpansion::m_ncoeffs, Vmath::Neg(), Nektar::StdRegions::NullConstFactorMap, Nektar::NullNekDouble1DArray, Nektar::StdRegions::StdExpansion::NumDGBndryCoeffs(), CellMLToNektar.cellml_metadata::p, Nektar::StdRegions::StdExpansion::PhysDeriv(), SetFaceToGeomOrientation(), SetTraceToGeomOrientation(), sign, Vmath::Svtvp(), Nektar::LocalRegions::Expansion::TraceNormLen(), Nektar::Transpose(), v_GetTraceOrient(), Vmath::Vcopy(), Vmath::Vmul(), Vmath::Vvtvp(), and Vmath::Zero().

Referenced by CreateMatrix(), Nektar::LocalRegions::HexExp::v_GenMatrix(), Nektar::LocalRegions::PrismExp::v_GenMatrix(), Nektar::LocalRegions::PyrExp::v_GenMatrix(), and Nektar::LocalRegions::TetExp::v_GenMatrix().

◆ v_GenTraceExp()

void Nektar::LocalRegions::Expansion3D::v_GenTraceExp ( const int  traceid,
ExpansionSharedPtr exp 
)
overrideprotectedvirtual

Reimplemented from Nektar::LocalRegions::Expansion.

Definition at line 2783 of file Expansion3D.cpp.

2784{
2785 SpatialDomains::GeometrySharedPtr faceGeom = m_geom->GetFace(traceid);
2786 if (faceGeom->GetNumVerts() == 3)
2787 {
2789 GetTraceBasisKey(traceid, 0), GetTraceBasisKey(traceid, 1),
2790 m_geom->GetFace(traceid));
2791 }
2792 else
2793 {
2795 GetTraceBasisKey(traceid, 0), GetTraceBasisKey(traceid, 1),
2796 m_geom->GetFace(traceid));
2797 }
2798}
std::shared_ptr< Geometry > GeometrySharedPtr
Definition: Geometry.h:51

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::StdRegions::StdExpansion::GetTraceBasisKey(), and Nektar::LocalRegions::Expansion::m_geom.

◆ v_GetTraceOrient()

StdRegions::Orientation Nektar::LocalRegions::Expansion3D::v_GetTraceOrient ( int  face)
overrideprotectedvirtual

Reimplemented from Nektar::LocalRegions::Expansion.

Definition at line 2721 of file Expansion3D.cpp.

2722{
2723 return m_geom->GetForient(face);
2724}

References Nektar::LocalRegions::Expansion::m_geom.

Referenced by v_GenMatrix().

◆ v_GetTracePhysVals()

void Nektar::LocalRegions::Expansion3D::v_GetTracePhysVals ( const int  face,
const StdRegions::StdExpansionSharedPtr FaceExp,
const Array< OneD, const NekDouble > &  inarray,
Array< OneD, NekDouble > &  outarray,
StdRegions::Orientation  orient 
)
overrideprotectedvirtual

Extract the physical values along face face from inarray into outarray following the local face orientation and point distribution defined by defined in FaceExp.

Reimplemented from Nektar::LocalRegions::Expansion.

Definition at line 2731 of file Expansion3D.cpp.

2735{
2736 if (orient == StdRegions::eNoOrientation)
2737 {
2738 orient = GetTraceOrient(face);
2739 }
2740
2741 int nq0 = FaceExp->GetNumPoints(0);
2742 int nq1 = FaceExp->GetNumPoints(1);
2743
2744 int nfacepts = GetTraceNumPoints(face);
2745 int dir0 = GetGeom3D()->GetDir(face, 0);
2746 int dir1 = GetGeom3D()->GetDir(face, 1);
2747
2748 Array<OneD, NekDouble> o_tmp(nfacepts);
2749 Array<OneD, NekDouble> o_tmp2(FaceExp->GetTotPoints());
2750 Array<OneD, int> faceids;
2751
2752 // Get local face pts and put into o_tmp
2753 GetTracePhysMap(face, faceids);
2754 // The static cast is necessary because faceids should be
2755 // Array<OneD, size_t> faceids ... or at leaste the same type as
2756 // faceids.size() ...
2757 Vmath::Gathr(static_cast<int>(faceids.size()), inarray, faceids, o_tmp);
2758
2759 int to_id0, to_id1;
2760
2762 {
2763 to_id0 = 0;
2764 to_id1 = 1;
2765 }
2766 else // transpose points key evaluation
2767 {
2768 to_id0 = 1;
2769 to_id1 = 0;
2770 }
2771
2772 // interpolate to points distrbution given in FaceExp
2774 m_base[dir0]->GetPointsKey(), m_base[dir1]->GetPointsKey(), o_tmp.get(),
2775 FaceExp->GetBasis(to_id0)->GetPointsKey(),
2776 FaceExp->GetBasis(to_id1)->GetPointsKey(), o_tmp2.get());
2777
2778 // Reshuffule points as required and put into outarray.
2779 v_ReOrientTracePhysMap(orient, faceids, nq0, nq1);
2780 Vmath::Scatr(nq0 * nq1, o_tmp2, faceids, outarray);
2781}
void v_ReOrientTracePhysMap(const StdRegions::Orientation orient, Array< OneD, int > &idmap, const int nq0, const int nq1) override
void GetTracePhysMap(const int edge, Array< OneD, int > &outarray)
Definition: Expansion.h:209
int GetTraceNumPoints(const int i) const
This function returns the number of quadrature points belonging to the i-th trace.
Definition: StdExpansion.h:281
void Gathr(I n, const T *x, const I *y, T *z)
Gather vector z[i] = x[y[i]].
Definition: Vmath.hpp:507
void Scatr(int n, const T *x, const int *y, T *z)
Scatter vector z[y[i]] = x[i].
Definition: Vmath.hpp:539

References Nektar::StdRegions::eDir1FwdDir2_Dir2FwdDir1, Nektar::StdRegions::eNoOrientation, Vmath::Gathr(), GetGeom3D(), Nektar::StdRegions::StdExpansion::GetTraceNumPoints(), Nektar::LocalRegions::Expansion::GetTraceOrient(), Nektar::LocalRegions::Expansion::GetTracePhysMap(), Nektar::LibUtilities::Interp2D(), Nektar::StdRegions::StdExpansion::m_base, Vmath::Scatr(), and v_ReOrientTracePhysMap().

◆ v_NormVectorIProductWRTBase()

void Nektar::LocalRegions::Expansion3D::v_NormVectorIProductWRTBase ( const Array< OneD, const Array< OneD, NekDouble > > &  Fvec,
Array< OneD, NekDouble > &  outarray 
)
overridevirtual

Reimplemented from Nektar::StdRegions::StdExpansion.

Definition at line 2938 of file Expansion3D.cpp.

2941{
2942 NormVectorIProductWRTBase(Fvec[0], Fvec[1], Fvec[2], outarray);
2943}
void NormVectorIProductWRTBase(const Array< OneD, const NekDouble > &Fx, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:613

References Nektar::StdRegions::StdExpansion::NormVectorIProductWRTBase().

◆ v_ReOrientTracePhysMap()

void Nektar::LocalRegions::Expansion3D::v_ReOrientTracePhysMap ( const StdRegions::Orientation  orient,
Array< OneD, int > &  idmap,
const int  nq0,
const int  nq1 
)
overridevirtual

Reimplemented from Nektar::LocalRegions::Expansion.

Definition at line 2800 of file Expansion3D.cpp.

2803{
2804
2805 if (idmap.size() != nq0 * nq1)
2806 {
2807 idmap = Array<OneD, int>(nq0 * nq1);
2808 }
2809
2810 if (GetNverts() == 3) // Tri face
2811 {
2812 switch (orient)
2813 {
2815 // eseentially straight copy
2816 for (int i = 0; i < nq0 * nq1; ++i)
2817 {
2818 idmap[i] = i;
2819 }
2820 break;
2822 // reverse
2823 for (int j = 0; j < nq1; ++j)
2824 {
2825 for (int i = 0; i < nq0; ++i)
2826 {
2827 idmap[j * nq0 + i] = nq0 - 1 - i + j * nq0;
2828 }
2829 }
2830 break;
2831 default:
2832 ASSERTL0(false,
2833 "Case not supposed to be used in this function");
2834 }
2835 }
2836 else
2837 {
2838 switch (orient)
2839 {
2841 // eseentially straight copy
2842 for (int i = 0; i < nq0 * nq1; ++i)
2843 {
2844 idmap[i] = i;
2845 }
2846 break;
2848 {
2849 // Direction A negative and B positive
2850 for (int j = 0; j < nq1; j++)
2851 {
2852 for (int i = 0; i < nq0; ++i)
2853 {
2854 idmap[j * nq0 + i] = nq0 - 1 - i + j * nq0;
2855 }
2856 }
2857 }
2858 break;
2860 {
2861 // Direction A positive and B negative
2862 for (int j = 0; j < nq1; j++)
2863 {
2864 for (int i = 0; i < nq0; ++i)
2865 {
2866 idmap[j * nq0 + i] = nq0 * (nq1 - 1 - j) + i;
2867 }
2868 }
2869 }
2870 break;
2872 {
2873 // Direction A negative and B negative
2874 for (int j = 0; j < nq1; j++)
2875 {
2876 for (int i = 0; i < nq0; ++i)
2877 {
2878 idmap[j * nq0 + i] = nq0 * nq1 - 1 - j * nq0 - i;
2879 }
2880 }
2881 }
2882 break;
2884 {
2885 // Transposed, Direction A and B positive
2886 for (int i = 0; i < nq0; ++i)
2887 {
2888 for (int j = 0; j < nq1; ++j)
2889 {
2890 idmap[i * nq1 + j] = i + j * nq0;
2891 }
2892 }
2893 }
2894 break;
2896 {
2897 // Transposed, Direction A positive and B negative
2898 for (int i = 0; i < nq0; ++i)
2899 {
2900 for (int j = 0; j < nq1; ++j)
2901 {
2902 idmap[i * nq1 + j] = nq0 - 1 - i + j * nq0;
2903 }
2904 }
2905 }
2906 break;
2908 {
2909 // Transposed, Direction A negative and B positive
2910 for (int i = 0; i < nq0; ++i)
2911 {
2912 for (int j = 0; j < nq1; ++j)
2913 {
2914 idmap[i * nq1 + j] = i + nq0 * (nq1 - 1) - j * nq0;
2915 }
2916 }
2917 }
2918 break;
2920 {
2921 // Transposed, Direction A and B negative
2922 for (int i = 0; i < nq0; ++i)
2923 {
2924 for (int j = 0; j < nq1; ++j)
2925 {
2926 idmap[i * nq1 + j] = nq0 * nq1 - 1 - i - j * nq0;
2927 }
2928 }
2929 }
2930 break;
2931 default:
2932 ASSERTL0(false, "Unknow orientation");
2933 break;
2934 }
2935 }
2936}

References ASSERTL0, Nektar::StdRegions::eDir1BwdDir1_Dir2BwdDir2, Nektar::StdRegions::eDir1BwdDir1_Dir2FwdDir2, Nektar::StdRegions::eDir1BwdDir2_Dir2BwdDir1, Nektar::StdRegions::eDir1BwdDir2_Dir2FwdDir1, Nektar::StdRegions::eDir1FwdDir1_Dir2BwdDir2, Nektar::StdRegions::eDir1FwdDir1_Dir2FwdDir2, Nektar::StdRegions::eDir1FwdDir2_Dir2BwdDir1, Nektar::StdRegions::eDir1FwdDir2_Dir2FwdDir1, and Nektar::StdRegions::StdExpansion::GetNverts().

Referenced by v_GetTracePhysVals().

◆ v_TraceNormLen()

void Nektar::LocalRegions::Expansion3D::v_TraceNormLen ( const int  traceid,
NekDouble h,
NekDouble p 
)
overrideprotectedvirtual

Reimplemented from Nektar::LocalRegions::Expansion.

Definition at line 2988 of file Expansion3D.cpp.

2989{
2991
2992 int nverts = geom->GetFace(traceid)->GetNumVerts();
2993
2994 SpatialDomains::PointGeom tn1, tn2, normal;
2995 tn1.Sub(*(geom->GetFace(traceid)->GetVertex(1)),
2996 *(geom->GetFace(traceid)->GetVertex(0)));
2997
2998 tn2.Sub(*(geom->GetFace(traceid)->GetVertex(nverts - 1)),
2999 *(geom->GetFace(traceid)->GetVertex(0)));
3000
3001 normal.Mult(tn1, tn2);
3002
3003 // normalise normal
3004 NekDouble mag = normal.dot(normal);
3005 mag = 1.0 / sqrt(mag);
3006 normal.UpdatePosition(normal.x() * mag, normal.y() * mag, normal.z() * mag);
3007
3008 SpatialDomains::PointGeom Dx;
3009 h = 0.0;
3010 p = 0.0;
3011 for (int i = 0; i < nverts; ++i)
3012 {
3013 // vertices on edges
3014 int edgid = geom->GetEdgeNormalToFaceVert(traceid, i);
3015
3016 // vector along noramal edge to each vertex
3017 Dx.Sub(*(geom->GetEdge(edgid)->GetVertex(0)),
3018 *(geom->GetEdge(edgid)->GetVertex(1)));
3019
3020 // calculate perpendicular distance of normal length
3021 // from first vertex
3022 h += fabs(normal.dot(Dx));
3023 }
3024
3025 h /= static_cast<NekDouble>(nverts);
3026
3027 // find normal basis direction
3028 int dir0 = geom->GetDir(traceid, 0);
3029 int dir1 = geom->GetDir(traceid, 1);
3030 int dirn;
3031 for (dirn = 0; dirn < 3; ++dirn)
3032 {
3033 if ((dirn != dir0) && (dirn != dir1))
3034 {
3035 break;
3036 }
3037 }
3038 p = (NekDouble)(GetBasisNumModes(dirn) - 1);
3039}
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:294

References Nektar::SpatialDomains::PointGeom::dot(), Nektar::StdRegions::StdExpansion::GetBasisNumModes(), Nektar::LocalRegions::Expansion::GetGeom(), Nektar::SpatialDomains::PointGeom::Mult(), CellMLToNektar.cellml_metadata::p, tinysimd::sqrt(), Nektar::SpatialDomains::PointGeom::Sub(), Nektar::SpatialDomains::PointGeom::UpdatePosition(), Nektar::NekPoint< data_type >::x(), Nektar::NekPoint< data_type >::y(), and Nektar::NekPoint< data_type >::z().

Member Data Documentation

◆ m_faceNormals

std::map<int, NormalVector> Nektar::LocalRegions::Expansion3D::m_faceNormals
protected

Definition at line 117 of file Expansion3D.h.

◆ m_requireNeg

std::vector<bool> Nektar::LocalRegions::Expansion3D::m_requireNeg
private

Definition at line 168 of file Expansion3D.h.

Referenced by v_AddFaceNormBoundaryInt().