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

#include <Expansion3D.h>

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

Public Member Functions

 Expansion3D (SpatialDomains::Geometry3DSharedPtr pGeom)
 
 ~Expansion3D () override=default
 
void SetTraceToGeomOrientation (Array< OneD, NekDouble > &inout)
 Align trace orientation with the geometry orientation. More...
 
void SetFaceToGeomOrientation (const int face, Array< OneD, NekDouble > &inout)
 Align face orientation with the geometry orientation. More...
 
void AddHDGHelmholtzFaceTerms (const NekDouble tau, const int edge, Array< OneD, NekDouble > &facePhys, const StdRegions::VarCoeffMap &dirForcing, Array< OneD, NekDouble > &outarray)
 
void AddNormTraceInt (const int dir, Array< OneD, ExpansionSharedPtr > &FaceExp, Array< OneD, Array< OneD, NekDouble > > &faceCoeffs, Array< OneD, NekDouble > &outarray)
 
void AddNormTraceInt (const int dir, Array< OneD, const NekDouble > &inarray, Array< OneD, ExpansionSharedPtr > &FaceExp, Array< OneD, NekDouble > &outarray, const StdRegions::VarCoeffMap &varcoeffs)
 
void AddFaceBoundaryInt (const int face, ExpansionSharedPtr &FaceExp, Array< OneD, NekDouble > &facePhys, Array< OneD, NekDouble > &outarray, const StdRegions::VarCoeffMap &varcoeffs=StdRegions::NullVarCoeffMap)
 
SpatialDomains::Geometry3DSharedPtr GetGeom3D () const
 
void v_ReOrientTracePhysMap (const StdRegions::Orientation orient, Array< OneD, int > &idmap, const int nq0, const int nq1) override
 
void v_NormVectorIProductWRTBase (const Array< OneD, const Array< OneD, NekDouble > > &Fvec, Array< OneD, NekDouble > &outarray) override
 
Array< OneD, unsigned int > GetEdgeInverseBoundaryMap (int eid)
 
Array< OneD, unsigned int > GetTraceInverseBoundaryMap (int fid, StdRegions::Orientation faceOrient=StdRegions::eNoOrientation, int P1=-1, int P2=-1)
 
void GetInverseBoundaryMaps (Array< OneD, unsigned int > &vmap, Array< OneD, Array< OneD, unsigned int > > &emap, Array< OneD, Array< OneD, unsigned int > > &fmap)
 
DNekScalMatSharedPtr CreateMatrix (const MatrixKey &mkey)
 
- Public Member Functions inherited from Nektar::LocalRegions::Expansion
 Expansion (SpatialDomains::GeometrySharedPtr pGeom)
 
 Expansion (const Expansion &pSrc)
 
 ~Expansion () override
 
void SetTraceExp (const int traceid, ExpansionSharedPtr &f)
 
ExpansionSharedPtr GetTraceExp (const int traceid)
 
DNekScalMatSharedPtr GetLocMatrix (const LocalRegions::MatrixKey &mkey)
 
void DropLocMatrix (const LocalRegions::MatrixKey &mkey)
 
DNekScalMatSharedPtr GetLocMatrix (const StdRegions::MatrixType mtype, const StdRegions::ConstFactorMap &factors=StdRegions::NullConstFactorMap, const StdRegions::VarCoeffMap &varcoeffs=StdRegions::NullVarCoeffMap)
 
SpatialDomains::GeometrySharedPtr GetGeom () const
 
void Reset ()
 
IndexMapValuesSharedPtr CreateIndexMap (const IndexMapKey &ikey)
 
DNekScalBlkMatSharedPtr CreateStaticCondMatrix (const MatrixKey &mkey)
 
const SpatialDomains::GeomFactorsSharedPtrGetMetricInfo () const
 
DNekMatSharedPtr BuildTransformationMatrix (const DNekScalMatSharedPtr &r_bnd, const StdRegions::MatrixType matrixType)
 
DNekMatSharedPtr BuildVertexMatrix (const DNekScalMatSharedPtr &r_bnd)
 
void ExtractDataToCoeffs (const NekDouble *data, const std::vector< unsigned int > &nummodes, const int nmodes_offset, NekDouble *coeffs, std::vector< LibUtilities::BasisType > &fromType)
 
void AddEdgeNormBoundaryInt (const int edge, const std::shared_ptr< Expansion > &EdgeExp, const Array< OneD, const NekDouble > &Fx, const Array< OneD, const NekDouble > &Fy, Array< OneD, NekDouble > &outarray)
 
void AddEdgeNormBoundaryInt (const int edge, const std::shared_ptr< Expansion > &EdgeExp, const Array< OneD, const NekDouble > &Fn, Array< OneD, NekDouble > &outarray)
 
void AddFaceNormBoundaryInt (const int face, const std::shared_ptr< Expansion > &FaceExp, const Array< OneD, const NekDouble > &Fn, Array< OneD, NekDouble > &outarray)
 
void DGDeriv (const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, ExpansionSharedPtr > &EdgeExp, Array< OneD, Array< OneD, NekDouble > > &coeffs, Array< OneD, NekDouble > &outarray)
 
NekDouble VectorFlux (const Array< OneD, Array< OneD, NekDouble > > &vec)
 
void NormalTraceDerivFactors (Array< OneD, Array< OneD, NekDouble > > &factors, Array< OneD, Array< OneD, NekDouble > > &d0factors, Array< OneD, Array< OneD, NekDouble > > &d1factors)
 
IndexMapValuesSharedPtr GetIndexMap (const IndexMapKey &ikey)
 
void AlignVectorToCollapsedDir (const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
 
ExpansionSharedPtr GetLeftAdjacentElementExp () const
 
ExpansionSharedPtr GetRightAdjacentElementExp () const
 
int GetLeftAdjacentElementTrace () const
 
int GetRightAdjacentElementTrace () const
 
void SetAdjacentElementExp (int traceid, ExpansionSharedPtr &e)
 
StdRegions::Orientation GetTraceOrient (int trace)
 
void SetCoeffsToOrientation (StdRegions::Orientation dir, Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
void DivideByQuadratureMetric (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 Divided by the metric jacobi and quadrature weights. More...
 
void GetTraceQFactors (const int trace, Array< OneD, NekDouble > &outarray)
 Extract the metric factors to compute the contravariant fluxes along edge edge and stores them into outarray following the local edge orientation (i.e. anticlockwise convention). More...
 
void GetTracePhysVals (const int trace, const StdRegions::StdExpansionSharedPtr &TraceExp, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, StdRegions::Orientation orient=StdRegions::eNoOrientation)
 
void GetTracePhysMap (const int edge, Array< OneD, int > &outarray)
 
void ReOrientTracePhysMap (const StdRegions::Orientation orient, Array< OneD, int > &idmap, const int nq0, const int nq1)
 
const NormalVectorGetTraceNormal (const int id)
 
void ComputeTraceNormal (const int id)
 
const Array< OneD, const NekDouble > & GetPhysNormals (void)
 
void SetPhysNormals (Array< OneD, const NekDouble > &normal)
 
void SetUpPhysNormals (const int trace)
 
void AddRobinMassMatrix (const int traceid, const Array< OneD, const NekDouble > &primCoeffs, DNekMatSharedPtr &inoutmat)
 
void TraceNormLen (const int traceid, NekDouble &h, NekDouble &p)
 
void AddRobinTraceContribution (const int traceid, const Array< OneD, const NekDouble > &primCoeffs, const Array< OneD, NekDouble > &incoeffs, Array< OneD, NekDouble > &coeffs)
 
const Array< OneD, const NekDouble > & GetElmtBndNormDirElmtLen (const int nbnd) const
 
void StdDerivBaseOnTraceMat (Array< OneD, DNekMatSharedPtr > &DerivMat)
 
- Public Member Functions inherited from Nektar::StdRegions::StdExpansion
 StdExpansion ()
 Default Constructor. More...
 
 StdExpansion (const int numcoeffs, const int numbases, const LibUtilities::BasisKey &Ba=LibUtilities::NullBasisKey, const LibUtilities::BasisKey &Bb=LibUtilities::NullBasisKey, const LibUtilities::BasisKey &Bc=LibUtilities::NullBasisKey)
 Constructor. More...
 
 StdExpansion (const StdExpansion &T)
 Copy Constructor. More...
 
virtual ~StdExpansion ()
 Destructor. More...
 
int GetNumBases () const
 This function returns the number of 1D bases used in the expansion. More...
 
const Array< OneD, const LibUtilities::BasisSharedPtr > & GetBase () const
 This function gets the shared point to basis. More...
 
const LibUtilities::BasisSharedPtrGetBasis (int dir) const
 This function gets the shared point to basis in the dir direction. More...
 
int GetNcoeffs (void) const
 This function returns the total number of coefficients used in the expansion. More...
 
int GetTotPoints () const
 This function returns the total number of quadrature points used in the element. More...
 
LibUtilities::BasisType GetBasisType (const int dir) const
 This function returns the type of basis used in the dir direction. More...
 
int GetBasisNumModes (const int dir) const
 This function returns the number of expansion modes in the dir direction. More...
 
int EvalBasisNumModesMax (void) const
 This function returns the maximum number of expansion modes over all local directions. More...
 
LibUtilities::PointsType GetPointsType (const int dir) const
 This function returns the type of quadrature points used in the dir direction. More...
 
int GetNumPoints (const int dir) const
 This function returns the number of quadrature points in the dir direction. More...
 
const Array< OneD, const NekDouble > & GetPoints (const int dir) const
 This function returns a pointer to the array containing the quadrature points in dir direction. More...
 
int GetNverts () const
 This function returns the number of vertices of the expansion domain. More...
 
int GetTraceNcoeffs (const int i) const
 This function returns the number of expansion coefficients belonging to the i-th trace. More...
 
int GetTraceIntNcoeffs (const int i) const
 
int GetTraceNumPoints (const int i) const
 This function returns the number of quadrature points belonging to the i-th trace. More...
 
const LibUtilities::BasisKey GetTraceBasisKey (const int i, int k=-1) const
 This function returns the basis key belonging to the i-th trace. More...
 
LibUtilities::PointsKey GetTracePointsKey (const int i, int k=-1) const
 This function returns the basis key belonging to the i-th trace. More...
 
int NumBndryCoeffs (void) const
 
int NumDGBndryCoeffs (void) const
 
const LibUtilities::PointsKey GetNodalPointsKey () const
 This function returns the type of expansion Nodal point type if defined. More...
 
int GetNtraces () const
 Returns the number of trace elements connected to this element. More...
 
LibUtilities::ShapeType DetShapeType () const
 This function returns the shape of the expansion domain. More...
 
std::shared_ptr< StdExpansionGetStdExp () const
 
std::shared_ptr< StdExpansionGetLinStdExp (void) const
 
int GetShapeDimension () const
 
bool IsBoundaryInteriorExpansion () const
 
bool IsNodalNonTensorialExp ()
 
void BwdTrans (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 This function performs the Backward transformation from coefficient space to physical space. More...
 
void FwdTrans (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 This function performs the Forward transformation from physical space to coefficient space. More...
 
void FwdTransBndConstrained (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
NekDouble Integral (const Array< OneD, const NekDouble > &inarray)
 This function integrates the specified function over the domain. More...
 
void FillMode (const int mode, Array< OneD, NekDouble > &outarray)
 This function fills the array outarray with the mode-th mode of the expansion. More...
 
void IProductWRTBase (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 this function calculates the inner product of a given function f with the different modes of the expansion More...
 
void IProductWRTBase (const Array< OneD, const NekDouble > &base, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, int coll_check)
 
void IProductWRTDerivBase (const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
void IProductWRTDirectionalDerivBase (const Array< OneD, const NekDouble > &direction, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
int GetElmtId ()
 Get the element id of this expansion when used in a list by returning value of m_elmt_id. More...
 
void SetElmtId (const int id)
 Set the element id of this expansion when used in a list by returning value of m_elmt_id. More...
 
void GetCoords (Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2=NullNekDouble1DArray, Array< OneD, NekDouble > &coords_3=NullNekDouble1DArray)
 this function returns the physical coordinates of the quadrature points of the expansion More...
 
void GetCoord (const Array< OneD, const NekDouble > &Lcoord, Array< OneD, NekDouble > &coord)
 given the coordinates of a point of the element in the local collapsed coordinate system, this function calculates the physical coordinates of the point More...
 
DNekMatSharedPtr GetStdMatrix (const StdMatrixKey &mkey)
 
DNekBlkMatSharedPtr GetStdStaticCondMatrix (const StdMatrixKey &mkey)
 
void NormVectorIProductWRTBase (const Array< OneD, const NekDouble > &Fx, Array< OneD, NekDouble > &outarray)
 
void NormVectorIProductWRTBase (const Array< OneD, const NekDouble > &Fx, const Array< OneD, NekDouble > &Fy, Array< OneD, NekDouble > &outarray)
 
void NormVectorIProductWRTBase (const Array< OneD, const NekDouble > &Fx, const Array< OneD, const NekDouble > &Fy, const Array< OneD, const NekDouble > &Fz, Array< OneD, NekDouble > &outarray)
 
void NormVectorIProductWRTBase (const Array< OneD, const Array< OneD, NekDouble > > &Fvec, Array< OneD, NekDouble > &outarray)
 
DNekScalBlkMatSharedPtr GetLocStaticCondMatrix (const LocalRegions::MatrixKey &mkey)
 
void DropLocStaticCondMatrix (const LocalRegions::MatrixKey &mkey)
 
int CalcNumberOfCoefficients (const std::vector< unsigned int > &nummodes, int &modes_offset)
 
NekDouble StdPhysEvaluate (const Array< OneD, const NekDouble > &Lcoord, const Array< OneD, const NekDouble > &physvals)
 
int GetCoordim ()
 
void GetBoundaryMap (Array< OneD, unsigned int > &outarray)
 
void GetInteriorMap (Array< OneD, unsigned int > &outarray)
 
int GetVertexMap (const int localVertexId, bool useCoeffPacking=false)
 
void GetTraceToElementMap (const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, Orientation traceOrient=eForwards, int P=-1, int Q=-1)
 
void GetTraceCoeffMap (const unsigned int traceid, Array< OneD, unsigned int > &maparray)
 
void GetElmtTraceToTraceMap (const unsigned int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, Orientation traceOrient=eForwards, int P=-1, int Q=-1)
 
void GetTraceInteriorToElementMap (const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, const Orientation traceOrient=eForwards)
 
void GetTraceNumModes (const int tid, int &numModes0, int &numModes1, const Orientation traceOrient=eDir1FwdDir1_Dir2FwdDir2)
 
void MultiplyByQuadratureMetric (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
void MultiplyByStdQuadratureMetric (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
DNekMatSharedPtr CreateGeneralMatrix (const StdMatrixKey &mkey)
 this function generates the mass matrix \(\mathbf{M}[i][j] = \int \phi_i(\mathbf{x}) \phi_j(\mathbf{x}) d\mathbf{x}\) More...
 
void GeneralMatrixOp (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
 
void MassMatrixOp (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
 
void LaplacianMatrixOp (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
 
void ReduceOrderCoeffs (int numMin, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
void SVVLaplacianFilter (Array< OneD, NekDouble > &array, const StdMatrixKey &mkey)
 
void ExponentialFilter (Array< OneD, NekDouble > &array, const NekDouble alpha, const NekDouble exponent, const NekDouble cutoff)
 
void LaplacianMatrixOp (const int k1, const int k2, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
 
void WeakDerivMatrixOp (const int i, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
 
void WeakDirectionalDerivMatrixOp (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
 
void MassLevelCurvatureMatrixOp (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
 
void LinearAdvectionMatrixOp (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
 
void LinearAdvectionDiffusionReactionMatrixOp (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey, bool addDiffusionTerm=true)
 
void HelmholtzMatrixOp (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
 
DNekMatSharedPtr GenMatrix (const StdMatrixKey &mkey)
 
void PhysDeriv (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1=NullNekDouble1DArray, Array< OneD, NekDouble > &out_d2=NullNekDouble1DArray)
 
void PhysDeriv (const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
void PhysDeriv_s (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_ds)
 
void PhysDeriv_n (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_dn)
 
void PhysDirectionalDeriv (const Array< OneD, const NekDouble > &inarray, const Array< OneD, const NekDouble > &direction, Array< OneD, NekDouble > &outarray)
 
void StdPhysDeriv (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1=NullNekDouble1DArray, Array< OneD, NekDouble > &out_d2=NullNekDouble1DArray)
 
void StdPhysDeriv (const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
NekDouble PhysEvaluate (const Array< OneD, const NekDouble > &coords, const Array< OneD, const NekDouble > &physvals)
 This function evaluates the expansion at a single (arbitrary) point of the domain. More...
 
NekDouble PhysEvaluate (const Array< OneD, NekDouble > &coord, const Array< OneD, const NekDouble > &inarray, std::array< NekDouble, 3 > &firstOrderDerivs)
 This function evaluates the first derivative of the expansion at a single (arbitrary) point of the domain. More...
 
NekDouble PhysEvaluate (const Array< OneD, NekDouble > &coord, const Array< OneD, const NekDouble > &inarray, std::array< NekDouble, 3 > &firstOrderDerivs, std::array< NekDouble, 6 > &secondOrderDerivs)
 
NekDouble PhysEvaluate (const Array< OneD, DNekMatSharedPtr > &I, const Array< OneD, const NekDouble > &physvals)
 This function evaluates the expansion at a single (arbitrary) point of the domain. More...
 
NekDouble PhysEvaluateBasis (const Array< OneD, const NekDouble > &coords, int mode)
 This function evaluates the basis function mode mode at a point coords of the domain. More...
 
void LocCoordToLocCollapsed (const Array< OneD, const NekDouble > &xi, Array< OneD, NekDouble > &eta)
 Convert local cartesian coordinate xi into local collapsed coordinates eta. More...
 
void LocCollapsedToLocCoord (const Array< OneD, const NekDouble > &eta, Array< OneD, NekDouble > &xi)
 Convert local collapsed coordinates eta into local cartesian coordinate xi. More...
 
virtual int v_CalcNumberOfCoefficients (const std::vector< unsigned int > &nummodes, int &modes_offset)
 
virtual void v_NormVectorIProductWRTBase (const Array< OneD, const NekDouble > &Fx, Array< OneD, NekDouble > &outarray)
 
virtual void v_NormVectorIProductWRTBase (const Array< OneD, const NekDouble > &Fx, const Array< OneD, const NekDouble > &Fy, Array< OneD, NekDouble > &outarray)
 
virtual void v_NormVectorIProductWRTBase (const Array< OneD, const NekDouble > &Fx, const Array< OneD, const NekDouble > &Fy, const Array< OneD, const NekDouble > &Fz, Array< OneD, NekDouble > &outarray)
 
virtual void v_NormVectorIProductWRTBase (const Array< OneD, const Array< OneD, NekDouble > > &Fvec, Array< OneD, NekDouble > &outarray)
 
virtual DNekScalBlkMatSharedPtr v_GetLocStaticCondMatrix (const LocalRegions::MatrixKey &mkey)
 
virtual void v_DropLocStaticCondMatrix (const LocalRegions::MatrixKey &mkey)
 
NekDouble Linf (const Array< OneD, const NekDouble > &phys, const Array< OneD, const NekDouble > &sol=NullNekDouble1DArray)
 Function to evaluate the discrete \( L_\infty\) error \( |\epsilon|_\infty = \max |u - u_{exact}|\) where \( u_{exact}\) is given by the array sol. More...
 
NekDouble L2 (const Array< OneD, const NekDouble > &phys, const Array< OneD, const NekDouble > &sol=NullNekDouble1DArray)
 Function to evaluate the discrete \( L_2\) error, \( | \epsilon |_{2} = \left [ \int^1_{-1} [u - u_{exact}]^2 dx \right]^{1/2} d\xi_1 \) where \( u_{exact}\) is given by the array sol. More...
 
NekDouble H1 (const Array< OneD, const NekDouble > &phys, const Array< OneD, const NekDouble > &sol=NullNekDouble1DArray)
 Function to evaluate the discrete \( H^1\) error, \( | \epsilon |^1_{2} = \left [ \int^1_{-1} [u - u_{exact}]^2 + \nabla(u - u_{exact})\cdot\nabla(u - u_{exact})\cdot dx \right]^{1/2} d\xi_1 \) where \( u_{exact}\) is given by the array sol. More...
 
const LibUtilities::PointsKeyVector GetPointsKeys () const
 
DNekMatSharedPtr BuildInverseTransformationMatrix (const DNekScalMatSharedPtr &m_transformationmatrix)
 
void PhysInterpToSimplexEquiSpaced (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, int npset=-1)
 This function performs an interpolation from the physical space points provided at input into an array of equispaced points which are not the collapsed coordinate. So for a tetrahedron you will only get a tetrahedral number of values. More...
 
void GetSimplexEquiSpacedConnectivity (Array< OneD, int > &conn, bool standard=true)
 This function provides the connectivity of local simplices (triangles or tets) to connect the equispaced data points provided by PhysInterpToSimplexEquiSpaced. More...
 
void EquiSpacedToCoeffs (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 This function performs a projection/interpolation from the equispaced points sometimes used in post-processing onto the coefficient space. More...
 
template<class T >
std::shared_ptr< T > as ()
 
void IProductWRTBase_SumFac (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
 
void GenStdMatBwdDeriv (const int dir, DNekMatSharedPtr &mat)
 
- Public Member Functions inherited from Nektar::StdRegions::StdExpansion3D
 StdExpansion3D (int numcoeffs, const LibUtilities::BasisKey &Ba, const LibUtilities::BasisKey &Bb, const LibUtilities::BasisKey &Bc)
 
 StdExpansion3D ()=default
 
 StdExpansion3D (const StdExpansion3D &T)=default
 
 ~StdExpansion3D () override=default
 
void PhysTensorDeriv (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray_d1, Array< OneD, NekDouble > &outarray_d2, Array< OneD, NekDouble > &outarray_d3)
 Calculate the 3D derivative in the local tensor/collapsed coordinate at the physical points. More...
 
void BwdTrans_SumFacKernel (const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0, bool doCheckCollDir1, bool doCheckCollDir2)
 
void IProductWRTBase_SumFacKernel (const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0, bool doCheckCollDir1, bool doCheckCollDir2)
 
int GetNedges () const
 return the number of edges in 3D expansion More...
 
int GetEdgeNcoeffs (const int i) const
 This function returns the number of expansion coefficients belonging to the i-th edge. More...
 
void GetEdgeInteriorToElementMap (const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, Orientation traceOrient=eForwards)
 

Protected Member Functions

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

Protected Attributes

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

Private Member Functions

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

Private Attributes

std::vector< bool > m_requireNeg
 

Detailed Description

Definition at line 55 of file Expansion3D.h.

Constructor & Destructor Documentation

◆ Expansion3D()

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

Definition at line 59 of file Expansion3D.h.

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

◆ ~Expansion3D()

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

Member Function Documentation

◆ AddFaceBoundaryInt()

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

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

Definition at line 322 of file Expansion3D.cpp.

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

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

Referenced by AddNormTraceInt().

◆ AddHDGHelmholtzFaceTerms()

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

Definition at line 50 of file Expansion3D.cpp.

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

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

Referenced by v_GenMatrix().

◆ AddNormTraceInt() [1/2]

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

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

@TODO: Document this

Definition at line 232 of file Expansion3D.cpp.

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

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

◆ AddNormTraceInt() [2/2]

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

Definition at line 294 of file Expansion3D.cpp.

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

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

Referenced by v_DGDeriv(), and v_GenMatrix().

◆ CreateMatrix()

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

Definition at line 407 of file Expansion3D.cpp.

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

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL1, ASSERTL2, Nektar::LocalRegions::Expansion::BuildTransformationMatrix(), Nektar::LocalRegions::Expansion::BuildVertexMatrix(), Nektar::StdRegions::StdMatrixKey::ConstFactorExists(), Nektar::StdRegions::StdExpansion::DetShapeType(), Nektar::LocalRegions::Expansion::DropLocMatrix(), Nektar::StdRegions::StdExpansion::DropLocStaticCondMatrix(), 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 2501 of file Expansion3D.cpp.

2502{
2503 int n, j;
2504 int nEdgeCoeffs;
2505 int nBndCoeffs = NumBndryCoeffs();
2506
2507 Array<OneD, unsigned int> bmap(nBndCoeffs);
2508 GetBoundaryMap(bmap);
2509
2510 // Map from full system to statically condensed system (i.e reverse
2511 // GetBoundaryMap)
2512 map<int, int> invmap;
2513 for (j = 0; j < nBndCoeffs; ++j)
2514 {
2515 invmap[bmap[j]] = j;
2516 }
2517
2518 // Number of interior edge coefficients
2519 nEdgeCoeffs = GetEdgeNcoeffs(eid) - 2;
2520
2522
2523 Array<OneD, unsigned int> edgemaparray(nEdgeCoeffs);
2524 StdRegions::Orientation eOrient = geom->GetEorient(eid);
2525 Array<OneD, unsigned int> maparray = Array<OneD, unsigned int>(nEdgeCoeffs);
2526 Array<OneD, int> signarray = Array<OneD, int>(nEdgeCoeffs, 1);
2527
2528 // maparray is the location of the edge within the matrix
2529 GetEdgeInteriorToElementMap(eid, maparray, signarray, eOrient);
2530
2531 for (n = 0; n < nEdgeCoeffs; ++n)
2532 {
2533 edgemaparray[n] = invmap[maparray[n]];
2534 }
2535
2536 return edgemaparray;
2537}
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 2653 of file Expansion3D.cpp.

2657{
2658 int n, j;
2659 int nEdgeCoeffs;
2660 int nFaceCoeffs;
2661
2662 int nBndCoeffs = NumBndryCoeffs();
2663
2664 Array<OneD, unsigned int> bmap(nBndCoeffs);
2665 GetBoundaryMap(bmap);
2666
2667 // Map from full system to statically condensed system (i.e reverse
2668 // GetBoundaryMap)
2669 map<int, int> reversemap;
2670 for (j = 0; j < bmap.size(); ++j)
2671 {
2672 reversemap[bmap[j]] = j;
2673 }
2674
2675 int nverts = GetNverts();
2676 vmap = Array<OneD, unsigned int>(nverts);
2677 for (n = 0; n < nverts; ++n)
2678 {
2679 int id = GetVertexMap(n);
2680 vmap[n] = reversemap[id]; // not sure what should be true here.
2681 }
2682
2683 int nedges = GetNedges();
2684 emap = Array<OneD, Array<OneD, unsigned int>>(nedges);
2685
2686 for (int eid = 0; eid < nedges; ++eid)
2687 {
2688 // Number of interior edge coefficients
2689 nEdgeCoeffs = GetEdgeNcoeffs(eid) - 2;
2690
2691 Array<OneD, unsigned int> edgemaparray(nEdgeCoeffs);
2692 Array<OneD, unsigned int> maparray =
2693 Array<OneD, unsigned int>(nEdgeCoeffs);
2694 Array<OneD, int> signarray = Array<OneD, int>(nEdgeCoeffs, 1);
2695
2696 // maparray is the location of the edge within the matrix
2697 GetEdgeInteriorToElementMap(eid, maparray, signarray,
2699
2700 for (n = 0; n < nEdgeCoeffs; ++n)
2701 {
2702 edgemaparray[n] = reversemap[maparray[n]];
2703 }
2704 emap[eid] = edgemaparray;
2705 }
2706
2707 int nfaces = GetNtraces();
2708 fmap = Array<OneD, Array<OneD, unsigned int>>(nfaces);
2709
2710 for (int fid = 0; fid < nfaces; ++fid)
2711 {
2712 // Number of interior face coefficients
2713 nFaceCoeffs = GetTraceIntNcoeffs(fid);
2714
2715 Array<OneD, unsigned int> facemaparray(nFaceCoeffs);
2716 Array<OneD, unsigned int> maparray =
2717 Array<OneD, unsigned int>(nFaceCoeffs);
2718 Array<OneD, int> signarray = Array<OneD, int>(nFaceCoeffs, 1);
2719
2720 // maparray is the location of the face within the matrix
2721 GetTraceInteriorToElementMap(fid, maparray, signarray,
2723
2724 for (n = 0; n < nFaceCoeffs; ++n)
2725 {
2726 facemaparray[n] = reversemap[maparray[n]];
2727 }
2728
2729 fmap[fid] = facemaparray;
2730 }
2731}
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 2958 of file Expansion3D.cpp.

2962{
2963 int nquad_f = FaceExp_f->GetNumPoints(0) * FaceExp_f->GetNumPoints(1);
2964 int coordim = GetCoordim();
2965
2966 int nquad0 = m_base[0]->GetNumPoints();
2967 int nquad1 = m_base[1]->GetNumPoints();
2968 int nquad2 = m_base[2]->GetNumPoints();
2969 int nqtot = nquad0 * nquad1 * nquad2;
2970
2971 StdRegions::VarCoeffType MMFCoeffs[15] = {
2980
2981 StdRegions::VarCoeffMap::const_iterator MFdir;
2982
2983 Array<OneD, NekDouble> nFacecdotMF(nqtot, 0.0);
2984 Array<OneD, NekDouble> tmp(nqtot);
2985 Array<OneD, NekDouble> tmp_f(nquad_f);
2986 for (int k = 0; k < coordim; k++)
2987 {
2988 MFdir = varcoeffs.find(MMFCoeffs[dir * 5 + k]);
2989 tmp = MFdir->second.GetValue();
2990
2991 GetPhysFaceVarCoeffsFromElement(face, FaceExp_f, tmp, tmp_f);
2992
2993 Vmath::Vvtvp(nquad_f, &tmp_f[0], 1, &normals[k][0], 1, &nFacecdotMF[0],
2994 1, &nFacecdotMF[0], 1);
2995 }
2996
2997 return nFacecdotMF;
2998}
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 2539 of file Expansion3D.cpp.

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

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

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

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

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

1974{
1975 MatrixStorage storage = eFULL;
1976 DNekMatSharedPtr vertexmatrix;
1977
1978 int nVerts, vid1, vid2, vMap1, vMap2;
1979 NekDouble VertexValue;
1980
1981 nVerts = GetNverts();
1982
1983 vertexmatrix =
1984 MemoryManager<DNekMat>::AllocateSharedPtr(nVerts, nVerts, 0.0, storage);
1985 DNekMat &VertexMat = (*vertexmatrix);
1986
1987 for (vid1 = 0; vid1 < nVerts; ++vid1)
1988 {
1989 vMap1 = GetVertexMap(vid1, true);
1990
1991 for (vid2 = 0; vid2 < nVerts; ++vid2)
1992 {
1993 vMap2 = GetVertexMap(vid2, true);
1994 VertexValue = (*r_bnd)(vMap1, vMap2);
1995 VertexMat.SetValue(vid1, vid2, VertexValue);
1996 }
1997 }
1998
1999 return vertexmatrix;
2000}

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

1794{
1795 int ncoeffs = GetNcoeffs();
1799
1801 DNekScalMat &Dmat = *GetLocMatrix(DerivType[dir]);
1802
1803 Array<OneD, NekDouble> coeffs = incoeffs;
1804 DNekVec Coeffs(ncoeffs, coeffs, eWrapper);
1805
1806 Coeffs = Transpose(Dmat) * Coeffs;
1807 Vmath::Neg(ncoeffs, coeffs, 1);
1808
1809 // Add the boundary integral including the relevant part of
1810 // the normal
1811 AddNormTraceInt(dir, FaceExp, faceCoeffs, coeffs);
1812
1813 DNekVec Out_d(ncoeffs, out_d, eWrapper);
1814
1815 Out_d = InvMass * Coeffs;
1816}
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 894 of file Expansion3D.cpp.

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

2796{
2797 SpatialDomains::GeometrySharedPtr faceGeom = m_geom->GetFace(traceid);
2798 if (faceGeom->GetNumVerts() == 3)
2799 {
2801 GetTraceBasisKey(traceid, 0), GetTraceBasisKey(traceid, 1),
2802 m_geom->GetFace(traceid));
2803 }
2804 else
2805 {
2807 GetTraceBasisKey(traceid, 0), GetTraceBasisKey(traceid, 1),
2808 m_geom->GetFace(traceid));
2809 }
2810}
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 2733 of file Expansion3D.cpp.

2734{
2735 return m_geom->GetForient(face);
2736}

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

2747{
2748 if (orient == StdRegions::eNoOrientation)
2749 {
2750 orient = GetTraceOrient(face);
2751 }
2752
2753 int nq0 = FaceExp->GetNumPoints(0);
2754 int nq1 = FaceExp->GetNumPoints(1);
2755
2756 int nfacepts = GetTraceNumPoints(face);
2757 int dir0 = GetGeom3D()->GetDir(face, 0);
2758 int dir1 = GetGeom3D()->GetDir(face, 1);
2759
2760 Array<OneD, NekDouble> o_tmp(nfacepts);
2761 Array<OneD, NekDouble> o_tmp2(FaceExp->GetTotPoints());
2762 Array<OneD, int> faceids;
2763
2764 // Get local face pts and put into o_tmp
2765 GetTracePhysMap(face, faceids);
2766 // The static cast is necessary because faceids should be
2767 // Array<OneD, size_t> faceids ... or at leaste the same type as
2768 // faceids.size() ...
2769 Vmath::Gathr(static_cast<int>(faceids.size()), inarray, faceids, o_tmp);
2770
2771 int to_id0, to_id1;
2772
2774 {
2775 to_id0 = 0;
2776 to_id1 = 1;
2777 }
2778 else // transpose points key evaluation
2779 {
2780 to_id0 = 1;
2781 to_id1 = 0;
2782 }
2783
2784 // interpolate to points distrbution given in FaceExp
2786 m_base[dir0]->GetPointsKey(), m_base[dir1]->GetPointsKey(), o_tmp.get(),
2787 FaceExp->GetBasis(to_id0)->GetPointsKey(),
2788 FaceExp->GetBasis(to_id1)->GetPointsKey(), o_tmp2.get());
2789
2790 // Reshuffule points as required and put into outarray.
2791 v_ReOrientTracePhysMap(orient, faceids, nq0, nq1);
2792 Vmath::Scatr(nq0 * nq1, o_tmp2, faceids, outarray);
2793}
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 2950 of file Expansion3D.cpp.

2953{
2954 NormVectorIProductWRTBase(Fvec[0], Fvec[1], Fvec[2], outarray);
2955}
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 2812 of file Expansion3D.cpp.

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

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

3001{
3003
3004 int nverts = geom->GetFace(traceid)->GetNumVerts();
3005
3006 SpatialDomains::PointGeom tn1, tn2, normal;
3007 tn1.Sub(*(geom->GetFace(traceid)->GetVertex(1)),
3008 *(geom->GetFace(traceid)->GetVertex(0)));
3009
3010 tn2.Sub(*(geom->GetFace(traceid)->GetVertex(nverts - 1)),
3011 *(geom->GetFace(traceid)->GetVertex(0)));
3012
3013 normal.Mult(tn1, tn2);
3014
3015 // normalise normal
3016 NekDouble mag = normal.dot(normal);
3017 mag = 1.0 / sqrt(mag);
3018 normal.UpdatePosition(normal.x() * mag, normal.y() * mag, normal.z() * mag);
3019
3020 SpatialDomains::PointGeom Dx;
3021 h = 0.0;
3022 p = 0.0;
3023 for (int i = 0; i < nverts; ++i)
3024 {
3025 // vertices on edges
3026 int edgid = geom->GetEdgeNormalToFaceVert(traceid, i);
3027
3028 // vector along noramal edge to each vertex
3029 Dx.Sub(*(geom->GetEdge(edgid)->GetVertex(0)),
3030 *(geom->GetEdge(edgid)->GetVertex(1)));
3031
3032 // calculate perpendicular distance of normal length
3033 // from first vertex
3034 h += fabs(normal.dot(Dx));
3035 }
3036
3037 h /= static_cast<NekDouble>(nverts);
3038
3039 // find normal basis direction
3040 int dir0 = geom->GetDir(traceid, 0);
3041 int dir1 = geom->GetDir(traceid, 1);
3042 int dirn;
3043 for (dirn = 0; dirn < 3; ++dirn)
3044 {
3045 if ((dirn != dir0) && (dirn != dir1))
3046 {
3047 break;
3048 }
3049 }
3050 p = (NekDouble)(GetBasisNumModes(dirn) - 1);
3051}
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().