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:403
std::map< StdRegions::VarCoeffType, VarCoeffEntry > VarCoeffMap
Definition: StdRegions.hpp:346
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 (Repeat varcoeff checks)
683 DropLocMatrix(advkey);
684 DropLocMatrix(masskey);
685 DropLocMatrix(lapkey);
686 }
687 break;
689 {
690 // Copied mostly from ADR solve to have fine-grain control
691 // over updating only advection matrix, relevant for performance!
692 NekDouble lambda = mkey.GetConstFactor(StdRegions::eFactorLambda);
693
694 // Construct mass matrix (Check for varcoeffs)
696 if (mkey.HasVarCoeff(StdRegions::eVarCoeffMass))
697 {
698 massVarcoeffs[StdRegions::eVarCoeffMass] =
699 mkey.GetVarCoeff(StdRegions::eVarCoeffMass);
700 }
701 MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this,
702 mkey.GetConstFactors(), massVarcoeffs);
703 DNekScalMat &MassMat = *GetLocMatrix(masskey);
704
705 // Construct laplacian matrix (Check for varcoeffs)
707 if ((mkey.HasVarCoeff(StdRegions::eVarCoeffLaplacian)) ||
708 (mkey.HasVarCoeff(StdRegions::eVarCoeffD00)) ||
709 (mkey.HasVarCoeff(StdRegions::eVarCoeffD01)) ||
710 (mkey.HasVarCoeff(StdRegions::eVarCoeffD10)) ||
711 (mkey.HasVarCoeff(StdRegions::eVarCoeffD02)) ||
712 (mkey.HasVarCoeff(StdRegions::eVarCoeffD20)) ||
713 (mkey.HasVarCoeff(StdRegions::eVarCoeffD11)) ||
714 (mkey.HasVarCoeff(StdRegions::eVarCoeffD12)) ||
715 (mkey.HasVarCoeff(StdRegions::eVarCoeffD21)) ||
716 (mkey.HasVarCoeff(StdRegions::eVarCoeffD22)))
717 {
718 lapVarcoeffs = mkey.GetVarCoeffs();
719 }
720 MatrixKey lapkey(StdRegions::eLaplacian, mkey.GetShapeType(), *this,
721 mkey.GetConstFactors(), lapVarcoeffs);
722 DNekScalMat &LapMat = *GetLocMatrix(lapkey);
723
724 // Construct advection matrix
725 // (assume advection velocity defined and non-zero)
726 // Could check L2(AdvectionVelocity) or HasVarCoeff
727 MatrixKey advkey(mkey, StdRegions::eLinearAdvection);
728 DNekScalMat &AdvMat = *GetLocMatrix(advkey);
729
730 // Generate a local copy of traceMat
731 MatrixKey gjpkey(StdRegions::eNormDerivOnTrace, mkey.GetShapeType(),
732 *this, mkey.GetConstFactors());
733 DNekScalMat &NDTraceMat = *GetLocMatrix(gjpkey);
734
735 NekDouble gjpfactor = mkey.GetConstFactor(StdRegions::eFactorGJP);
736 ASSERTL1(mkey.ConstFactorExists(StdRegions::eFactorGJP),
737 "Need to specify eFactorGJP to construct "
738 "a LinearAdvectionDiffusionReactionGJP matrix");
739
740 int rows = LapMat.GetRows();
741 int cols = LapMat.GetColumns();
742
743 DNekMatSharedPtr adr =
745
746 NekDouble one = 1.0;
747 (*adr) =
748 LapMat - lambda * MassMat + AdvMat + gjpfactor * NDTraceMat;
749
751
752 // Clear memory
753 DropLocMatrix(advkey);
754 DropLocMatrix(masskey);
755 DropLocMatrix(lapkey);
756 DropLocMatrix(gjpkey);
757 }
758 break;
760 {
761 NekDouble one = 1.0;
763
765 }
766 break;
768 {
769 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
770 {
771 NekDouble one = 1.0;
772 DNekMatSharedPtr mat = GenMatrix(mkey);
773 returnval =
775 }
776 else
777 {
778 NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
779 DNekMatSharedPtr mat = GetStdMatrix(mkey);
780 returnval =
782 }
783 }
784 break;
786 {
787 NekDouble one = 1.0;
788
790 *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
791 DNekMatSharedPtr mat = GenMatrix(hkey);
792
793 mat->Invert();
795 }
796 break;
798 {
799 NekDouble one = 1.0;
800 MatrixKey helmkey(StdRegions::eHelmholtz, mkey.GetShapeType(),
801 *this, mkey.GetConstFactors(),
802 mkey.GetVarCoeffs());
803 DNekScalBlkMatSharedPtr helmStatCond =
804 GetLocStaticCondMatrix(helmkey);
805 DNekScalMatSharedPtr A = helmStatCond->GetBlock(0, 0);
807
809 }
810 break;
812 {
813 NekDouble one = 1.0;
814 MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this);
815 DNekScalBlkMatSharedPtr massStatCond =
816 GetLocStaticCondMatrix(masskey);
817 DNekScalMatSharedPtr A = massStatCond->GetBlock(0, 0);
819
821 }
822 break;
824 {
825 NekDouble one = 1.0;
826 MatrixKey helmkey(StdRegions::eHelmholtz, mkey.GetShapeType(),
827 *this, mkey.GetConstFactors(),
828 mkey.GetVarCoeffs());
829 DNekScalBlkMatSharedPtr helmStatCond =
830 GetLocStaticCondMatrix(helmkey);
831 DNekScalMatSharedPtr A = helmStatCond->GetBlock(0, 0);
832
835 BuildTransformationMatrix(A, mkey.GetMatrixType());
836
838 }
839 break;
841 {
842 NekDouble one = 1.0;
843 MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this);
845 DNekScalMatSharedPtr A = StatCond->GetBlock(0, 0);
846
849 BuildTransformationMatrix(A, mkey.GetMatrixType());
850
852 }
853 break;
854 default:
855 {
856 NekDouble one = 1.0;
857 DNekMatSharedPtr mat = GenMatrix(mkey);
858
860 }
861 break;
862 }
863
864 return returnval;
865}
#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:347
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 2483 of file Expansion3D.cpp.

2484{
2485 int n, j;
2486 int nEdgeCoeffs;
2487 int nBndCoeffs = NumBndryCoeffs();
2488
2489 Array<OneD, unsigned int> bmap(nBndCoeffs);
2490 GetBoundaryMap(bmap);
2491
2492 // Map from full system to statically condensed system (i.e reverse
2493 // GetBoundaryMap)
2494 map<int, int> invmap;
2495 for (j = 0; j < nBndCoeffs; ++j)
2496 {
2497 invmap[bmap[j]] = j;
2498 }
2499
2500 // Number of interior edge coefficients
2501 nEdgeCoeffs = GetEdgeNcoeffs(eid) - 2;
2502
2504
2505 Array<OneD, unsigned int> edgemaparray(nEdgeCoeffs);
2506 StdRegions::Orientation eOrient = geom->GetEorient(eid);
2507 Array<OneD, unsigned int> maparray = Array<OneD, unsigned int>(nEdgeCoeffs);
2508 Array<OneD, int> signarray = Array<OneD, int>(nEdgeCoeffs, 1);
2509
2510 // maparray is the location of the edge within the matrix
2511 GetEdgeInteriorToElementMap(eid, maparray, signarray, eOrient);
2512
2513 for (n = 0; n < nEdgeCoeffs; ++n)
2514 {
2515 edgemaparray[n] = invmap[maparray[n]];
2516 }
2517
2518 return edgemaparray;
2519}
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 2635 of file Expansion3D.cpp.

2639{
2640 int n, j;
2641 int nEdgeCoeffs;
2642 int nFaceCoeffs;
2643
2644 int nBndCoeffs = NumBndryCoeffs();
2645
2646 Array<OneD, unsigned int> bmap(nBndCoeffs);
2647 GetBoundaryMap(bmap);
2648
2649 // Map from full system to statically condensed system (i.e reverse
2650 // GetBoundaryMap)
2651 map<int, int> reversemap;
2652 for (j = 0; j < bmap.size(); ++j)
2653 {
2654 reversemap[bmap[j]] = j;
2655 }
2656
2657 int nverts = GetNverts();
2658 vmap = Array<OneD, unsigned int>(nverts);
2659 for (n = 0; n < nverts; ++n)
2660 {
2661 int id = GetVertexMap(n);
2662 vmap[n] = reversemap[id]; // not sure what should be true here.
2663 }
2664
2665 int nedges = GetNedges();
2666 emap = Array<OneD, Array<OneD, unsigned int>>(nedges);
2667
2668 for (int eid = 0; eid < nedges; ++eid)
2669 {
2670 // Number of interior edge coefficients
2671 nEdgeCoeffs = GetEdgeNcoeffs(eid) - 2;
2672
2673 Array<OneD, unsigned int> edgemaparray(nEdgeCoeffs);
2674 Array<OneD, unsigned int> maparray =
2675 Array<OneD, unsigned int>(nEdgeCoeffs);
2676 Array<OneD, int> signarray = Array<OneD, int>(nEdgeCoeffs, 1);
2677
2678 // maparray is the location of the edge within the matrix
2679 GetEdgeInteriorToElementMap(eid, maparray, signarray,
2681
2682 for (n = 0; n < nEdgeCoeffs; ++n)
2683 {
2684 edgemaparray[n] = reversemap[maparray[n]];
2685 }
2686 emap[eid] = edgemaparray;
2687 }
2688
2689 int nfaces = GetNtraces();
2690 fmap = Array<OneD, Array<OneD, unsigned int>>(nfaces);
2691
2692 for (int fid = 0; fid < nfaces; ++fid)
2693 {
2694 // Number of interior face coefficients
2695 nFaceCoeffs = GetTraceIntNcoeffs(fid);
2696
2697 Array<OneD, unsigned int> facemaparray(nFaceCoeffs);
2698 Array<OneD, unsigned int> maparray =
2699 Array<OneD, unsigned int>(nFaceCoeffs);
2700 Array<OneD, int> signarray = Array<OneD, int>(nFaceCoeffs, 1);
2701
2702 // maparray is the location of the face within the matrix
2703 GetTraceInteriorToElementMap(fid, maparray, signarray,
2705
2706 for (n = 0; n < nFaceCoeffs; ++n)
2707 {
2708 facemaparray[n] = reversemap[maparray[n]];
2709 }
2710
2711 fmap[fid] = facemaparray;
2712 }
2713}
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 2940 of file Expansion3D.cpp.

2944{
2945 int nquad_f = FaceExp_f->GetNumPoints(0) * FaceExp_f->GetNumPoints(1);
2946 int coordim = GetCoordim();
2947
2948 int nquad0 = m_base[0]->GetNumPoints();
2949 int nquad1 = m_base[1]->GetNumPoints();
2950 int nquad2 = m_base[2]->GetNumPoints();
2951 int nqtot = nquad0 * nquad1 * nquad2;
2952
2953 StdRegions::VarCoeffType MMFCoeffs[15] = {
2962
2963 StdRegions::VarCoeffMap::const_iterator MFdir;
2964
2965 Array<OneD, NekDouble> nFacecdotMF(nqtot, 0.0);
2966 Array<OneD, NekDouble> tmp(nqtot);
2967 Array<OneD, NekDouble> tmp_f(nquad_f);
2968 for (int k = 0; k < coordim; k++)
2969 {
2970 MFdir = varcoeffs.find(MMFCoeffs[dir * 5 + k]);
2971 tmp = MFdir->second.GetValue();
2972
2973 GetPhysFaceVarCoeffsFromElement(face, FaceExp_f, tmp, tmp_f);
2974
2975 Vmath::Vvtvp(nquad_f, &tmp_f[0], 1, &normals[k][0], 1, &nFacecdotMF[0],
2976 1, &nFacecdotMF[0], 1);
2977 }
2978
2979 return nFacecdotMF;
2980}
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 2521 of file Expansion3D.cpp.

2523{
2524 int n, j;
2525 int nFaceCoeffs;
2526
2527 int nBndCoeffs = NumBndryCoeffs();
2528
2529 Array<OneD, unsigned int> bmap(nBndCoeffs);
2530 GetBoundaryMap(bmap);
2531
2532 // Map from full system to statically condensed system (i.e reverse
2533 // GetBoundaryMap)
2534 map<int, int> reversemap;
2535 for (j = 0; j < bmap.size(); ++j)
2536 {
2537 reversemap[bmap[j]] = j;
2538 }
2539
2540 // Number of interior face coefficients
2541 nFaceCoeffs = GetTraceIntNcoeffs(fid);
2542
2544 Array<OneD, unsigned int> maparray = Array<OneD, unsigned int>(nFaceCoeffs);
2545 Array<OneD, int> signarray = Array<OneD, int>(nFaceCoeffs, 1);
2546
2547 if (faceOrient == StdRegions::eNoOrientation)
2548 {
2549 fOrient = GetTraceOrient(fid);
2550 }
2551 else
2552 {
2553 fOrient = faceOrient;
2554 }
2555
2556 // maparray is the location of the face within the matrix
2557 GetTraceInteriorToElementMap(fid, maparray, signarray, fOrient);
2558
2559 Array<OneD, unsigned int> facemaparray;
2560 int locP1, locP2;
2561 GetTraceNumModes(fid, locP1, locP2, fOrient);
2562
2563 if (P1 == -1)
2564 {
2565 P1 = locP1;
2566 }
2567 else
2568 {
2569 ASSERTL1(P1 <= locP1, "Expect value of passed P1 to "
2570 "be lower or equal to face num modes");
2571 }
2572
2573 if (P2 == -1)
2574 {
2575 P2 = locP2;
2576 }
2577 else
2578 {
2579 ASSERTL1(P2 <= locP2, "Expect value of passed P2 to "
2580 "be lower or equal to face num modes");
2581 }
2582
2583 switch (GetGeom3D()->GetFace(fid)->GetShapeType())
2584 {
2586 {
2587 if (((P1 - 3) > 0) && ((P2 - 3) > 0))
2588 {
2589 facemaparray = Array<OneD, unsigned int>(
2591 P2 - 3));
2592 int cnt = 0;
2593 int cnt1 = 0;
2594 for (n = 0; n < P1 - 3; ++n)
2595 {
2596 for (int m = 0; m < P2 - 3 - n; ++m, ++cnt)
2597 {
2598 facemaparray[cnt] = reversemap[maparray[cnt1 + m]];
2599 }
2600 cnt1 += locP2 - 3 - n;
2601 }
2602 }
2603 }
2604 break;
2606 {
2607 if (((P1 - 2) > 0) && ((P2 - 2) > 0))
2608 {
2609 facemaparray = Array<OneD, unsigned int>(
2611 P2 - 2));
2612 int cnt = 0;
2613 int cnt1 = 0;
2614 for (n = 0; n < P2 - 2; ++n)
2615 {
2616 for (int m = 0; m < P1 - 2; ++m, ++cnt)
2617 {
2618 facemaparray[cnt] = reversemap[maparray[cnt1 + m]];
2619 }
2620 cnt1 += locP1 - 2;
2621 }
2622 }
2623 }
2624 break;
2625 default:
2626 {
2627 ASSERTL0(false, "Invalid shape type.");
2628 }
2629 break;
2630 }
2631
2632 return facemaparray;
2633}
#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 1643 of file Expansion3D.cpp.

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

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

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

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

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

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 1771 of file Expansion3D.cpp.

1776{
1777 int ncoeffs = GetNcoeffs();
1781
1783 DNekScalMat &Dmat = *GetLocMatrix(DerivType[dir]);
1784
1785 Array<OneD, NekDouble> coeffs = incoeffs;
1786 DNekVec Coeffs(ncoeffs, coeffs, eWrapper);
1787
1788 Coeffs = Transpose(Dmat) * Coeffs;
1789 Vmath::Neg(ncoeffs, coeffs, 1);
1790
1791 // Add the boundary integral including the relevant part of
1792 // the normal
1793 AddNormTraceInt(dir, FaceExp, faceCoeffs, coeffs);
1794
1795 DNekVec Out_d(ncoeffs, out_d, eWrapper);
1796
1797 Out_d = InvMass * Coeffs;
1798}
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 876 of file Expansion3D.cpp.

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

2778{
2779 SpatialDomains::GeometrySharedPtr faceGeom = m_geom->GetFace(traceid);
2780 if (faceGeom->GetNumVerts() == 3)
2781 {
2783 GetTraceBasisKey(traceid, 0), GetTraceBasisKey(traceid, 1),
2784 m_geom->GetFace(traceid));
2785 }
2786 else
2787 {
2789 GetTraceBasisKey(traceid, 0), GetTraceBasisKey(traceid, 1),
2790 m_geom->GetFace(traceid));
2791 }
2792}
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 2715 of file Expansion3D.cpp.

2716{
2717 return m_geom->GetForient(face);
2718}

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 2725 of file Expansion3D.cpp.

2729{
2730 if (orient == StdRegions::eNoOrientation)
2731 {
2732 orient = GetTraceOrient(face);
2733 }
2734
2735 int nq0 = FaceExp->GetNumPoints(0);
2736 int nq1 = FaceExp->GetNumPoints(1);
2737
2738 int nfacepts = GetTraceNumPoints(face);
2739 int dir0 = GetGeom3D()->GetDir(face, 0);
2740 int dir1 = GetGeom3D()->GetDir(face, 1);
2741
2742 Array<OneD, NekDouble> o_tmp(nfacepts);
2743 Array<OneD, NekDouble> o_tmp2(FaceExp->GetTotPoints());
2744 Array<OneD, int> faceids;
2745
2746 // Get local face pts and put into o_tmp
2747 GetTracePhysMap(face, faceids);
2748 // The static cast is necessary because faceids should be
2749 // Array<OneD, size_t> faceids ... or at leaste the same type as
2750 // faceids.size() ...
2751 Vmath::Gathr(static_cast<int>(faceids.size()), inarray, faceids, o_tmp);
2752
2753 int to_id0, to_id1;
2754
2756 {
2757 to_id0 = 0;
2758 to_id1 = 1;
2759 }
2760 else // transpose points key evaluation
2761 {
2762 to_id0 = 1;
2763 to_id1 = 0;
2764 }
2765
2766 // interpolate to points distrbution given in FaceExp
2768 m_base[dir0]->GetPointsKey(), m_base[dir1]->GetPointsKey(), o_tmp.get(),
2769 FaceExp->GetBasis(to_id0)->GetPointsKey(),
2770 FaceExp->GetBasis(to_id1)->GetPointsKey(), o_tmp2.get());
2771
2772 // Reshuffule points as required and put into outarray.
2773 v_ReOrientTracePhysMap(orient, faceids, nq0, nq1);
2774 Vmath::Scatr(nq0 * nq1, o_tmp2, faceids, outarray);
2775}
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 2932 of file Expansion3D.cpp.

2935{
2936 NormVectorIProductWRTBase(Fvec[0], Fvec[1], Fvec[2], outarray);
2937}
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 2794 of file Expansion3D.cpp.

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

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 2982 of file Expansion3D.cpp.

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