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

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).
 
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.
 
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}}\).
 
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 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 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)
 
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_GetTracePhysMap (const int edge, Array< OneD, int > &outarray)
 
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_AddRobinTraceContribution (const int traceid, const Array< OneD, const NekDouble > &primCoeffs, const Array< OneD, NekDouble > &incoeffs, Array< OneD, NekDouble > &coeffs)
 
- 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.
 
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 NekDouble v_StdPhysEvaluate (const Array< OneD, const NekDouble > &Lcoord, const Array< OneD, const NekDouble > &physvals)
 
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.
 
template<int DIR>
NekDouble BaryEvaluateBasis (const NekDouble &coord, const int &mode)
 
template<int DIR, bool DERIV = false, bool DERIV2 = false>
NekDouble BaryEvaluate (const NekDouble &coord, const NekDouble *physvals)
 Helper function to pass an unused value by reference into BaryEvaluate.
 
template<int DIR, bool DERIV = false, bool DERIV2 = false>
NekDouble BaryEvaluate (const NekDouble &coord, const NekDouble *physvals, NekDouble &deriv)
 
- 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.
 
NekDouble v_PhysEvaluateInterp (const Array< OneD, DNekMatSharedPtr > &I, const Array< OneD, const NekDouble > &physvals) 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.
 
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
 
void v_PhysInterp (std::shared_ptr< StdExpansion > fromExp, const Array< OneD, const NekDouble > &fromData, Array< OneD, NekDouble > &toData) 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::Geometrym_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)
 
- 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::Geometry3D pGeom)
inline

Definition at line 59 of file Expansion3D.h.

61 {
62 }
std::vector< bool > m_requireNeg
Expansion(SpatialDomains::Geometry *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.
int GetBasisNumModes(const int dir) const
This function returns the number of expansion modes in the dir direction.
std::shared_ptr< IndexMapValues > IndexMapValuesSharedPtr

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

Referenced by AddNormTraceInt(), and 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)
ExpansionSharedPtr GetTraceExp(const int traceid)
Definition Expansion.h:420
DNekScalMatSharedPtr GetLocMatrix(const LocalRegions::MatrixKey &mkey)
Definition Expansion.cpp:84
Array< OneD, NekDouble > GetMFDiv(const int dir, const StdRegions::VarCoeffMap &varcoeffs)
const NormalVector & GetTraceNormal(const int id)
Array< OneD, NekDouble > GetMF(const int dir, const int shapedim, const StdRegions::VarCoeffMap &varcoeffs)
int GetNcoeffs(void) const
This function returns the total number of coefficients used in the expansion.
std::shared_ptr< Expansion > ExpansionSharedPtr
Definition Expansion.h:66
static ConstFactorMap NullConstFactorMap
std::map< StdRegions::VarCoeffType, VarCoeffEntry > VarCoeffMap
NekMatrix< NekMatrix< NekDouble, StandardMatrixTag >, ScaledMatrixTag > DNekScalMat
NekVector< NekDouble > DNekVec
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.

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 advection matrix
644 // Check for varcoeffs not required;
645 // assume advection velocity is always time-dependent
646 MatrixKey advkey(mkey, StdRegions::eLinearAdvection);
647 DNekScalMat &AdvMat = *GetLocMatrix(advkey);
648
649 int rows = MassMat.GetRows();
650 int cols = MassMat.GetColumns();
651
652 DNekMatSharedPtr adr =
654
655 NekDouble one = 1.0;
656 (*adr) = -lambda * MassMat + AdvMat;
657
659
660 // Clear memory for time-dependent matrices
661 DropLocMatrix(advkey);
662 if (!massVarcoeffs.empty())
663 {
664 DropLocMatrix(masskey);
665 }
666 }
667 break;
669 {
670 NekDouble lambda = mkey.GetConstFactor(StdRegions::eFactorLambda);
671
672 // Construct mass matrix
673 // Check for mass-specific varcoeffs to avoid unncessary
674 // re-computation of the elemental matrix every time step
676 if (mkey.HasVarCoeff(StdRegions::eVarCoeffMass))
677 {
678 massVarcoeffs[StdRegions::eVarCoeffMass] =
679 mkey.GetVarCoeff(StdRegions::eVarCoeffMass);
680 }
681 MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this,
682 mkey.GetConstFactors(), massVarcoeffs);
683 DNekScalMat &MassMat = *GetLocMatrix(masskey);
684
685 // Construct laplacian matrix (Check for varcoeffs)
686 // Take all varcoeffs if one or more are detected
687 // TODO We might want to have a map
688 // from MatrixType to Vector of Varcoeffs and vice-versa
690 if ((mkey.HasVarCoeff(StdRegions::eVarCoeffLaplacian)) ||
691 (mkey.HasVarCoeff(StdRegions::eVarCoeffD00)) ||
692 (mkey.HasVarCoeff(StdRegions::eVarCoeffD01)) ||
693 (mkey.HasVarCoeff(StdRegions::eVarCoeffD10)) ||
694 (mkey.HasVarCoeff(StdRegions::eVarCoeffD02)) ||
695 (mkey.HasVarCoeff(StdRegions::eVarCoeffD20)) ||
696 (mkey.HasVarCoeff(StdRegions::eVarCoeffD11)) ||
697 (mkey.HasVarCoeff(StdRegions::eVarCoeffD12)) ||
698 (mkey.HasVarCoeff(StdRegions::eVarCoeffD21)) ||
699 (mkey.HasVarCoeff(StdRegions::eVarCoeffD22)))
700 {
701 lapVarcoeffs = mkey.GetVarCoeffs();
702 }
703 MatrixKey lapkey(StdRegions::eLaplacian, mkey.GetShapeType(), *this,
704 mkey.GetConstFactors(), lapVarcoeffs);
705 DNekScalMat &LapMat = *GetLocMatrix(lapkey);
706
707 // Construct advection matrix
708 // Check for varcoeffs not required;
709 // assume advection velocity is always time-dependent
710 MatrixKey advkey(mkey, StdRegions::eLinearAdvection);
711 DNekScalMat &AdvMat = *GetLocMatrix(advkey);
712
713 int rows = LapMat.GetRows();
714 int cols = LapMat.GetColumns();
715
716 DNekMatSharedPtr adr =
718
719 NekDouble one = 1.0;
720 (*adr) = LapMat - lambda * MassMat + AdvMat;
721
723
724 // Clear memory for time-dependent matrices
725 DropLocMatrix(advkey);
726 if (!massVarcoeffs.empty())
727 {
728 DropLocMatrix(masskey);
729 }
730 if (!lapVarcoeffs.empty())
731 {
732 DropLocMatrix(lapkey);
733 }
734 }
735 break;
737 {
738 // Copied mostly from ADR solve to have fine-grain control
739 // over updating only advection matrix, relevant for performance!
740 NekDouble lambda = mkey.GetConstFactor(StdRegions::eFactorLambda);
741
742 // Construct mass matrix (Check for varcoeffs)
744 if (mkey.HasVarCoeff(StdRegions::eVarCoeffMass))
745 {
746 massVarcoeffs[StdRegions::eVarCoeffMass] =
747 mkey.GetVarCoeff(StdRegions::eVarCoeffMass);
748 }
749 MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this,
750 mkey.GetConstFactors(), massVarcoeffs);
751 DNekScalMat &MassMat = *GetLocMatrix(masskey);
752
753 // Construct laplacian matrix (Check for varcoeffs)
755 if ((mkey.HasVarCoeff(StdRegions::eVarCoeffLaplacian)) ||
756 (mkey.HasVarCoeff(StdRegions::eVarCoeffD00)) ||
757 (mkey.HasVarCoeff(StdRegions::eVarCoeffD01)) ||
758 (mkey.HasVarCoeff(StdRegions::eVarCoeffD10)) ||
759 (mkey.HasVarCoeff(StdRegions::eVarCoeffD02)) ||
760 (mkey.HasVarCoeff(StdRegions::eVarCoeffD20)) ||
761 (mkey.HasVarCoeff(StdRegions::eVarCoeffD11)) ||
762 (mkey.HasVarCoeff(StdRegions::eVarCoeffD12)) ||
763 (mkey.HasVarCoeff(StdRegions::eVarCoeffD21)) ||
764 (mkey.HasVarCoeff(StdRegions::eVarCoeffD22)))
765 {
766 lapVarcoeffs = mkey.GetVarCoeffs();
767 }
768 MatrixKey lapkey(StdRegions::eLaplacian, mkey.GetShapeType(), *this,
769 mkey.GetConstFactors(), lapVarcoeffs);
770 DNekScalMat &LapMat = *GetLocMatrix(lapkey);
771
772 // Construct advection matrix
773 // (assume advection velocity defined and non-zero)
774 // Could check L2(AdvectionVelocity) or HasVarCoeff
775 MatrixKey advkey(mkey, StdRegions::eLinearAdvection);
776 DNekScalMat &AdvMat = *GetLocMatrix(advkey);
777
778 // Generate a local copy of traceMat
779 MatrixKey gjpkey(StdRegions::eNormDerivOnTrace, mkey.GetShapeType(),
780 *this, mkey.GetConstFactors());
781 DNekScalMat &NDTraceMat = *GetLocMatrix(gjpkey);
782
783 NekDouble gjpfactor = mkey.GetConstFactor(StdRegions::eFactorGJP);
784 ASSERTL1(mkey.ConstFactorExists(StdRegions::eFactorGJP),
785 "Need to specify eFactorGJP to construct "
786 "a LinearAdvectionDiffusionReactionGJP matrix");
787
788 int rows = LapMat.GetRows();
789 int cols = LapMat.GetColumns();
790
791 DNekMatSharedPtr adr =
793
794 NekDouble one = 1.0;
795 (*adr) =
796 LapMat - lambda * MassMat + AdvMat + gjpfactor * NDTraceMat;
797
799
800 // Clear memory
801 DropLocMatrix(advkey);
802 DropLocMatrix(masskey);
803 DropLocMatrix(lapkey);
804 DropLocMatrix(gjpkey);
805 }
806 break;
808 {
809 NekDouble one = 1.0;
811
813 }
814 break;
816 {
817 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
818 {
819 NekDouble one = 1.0;
820 DNekMatSharedPtr mat = GenMatrix(mkey);
821 returnval =
823 }
824 else
825 {
826 NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
827 DNekMatSharedPtr mat = GetStdMatrix(mkey);
828 returnval =
830 }
831 }
832 break;
834 {
835 NekDouble one = 1.0;
836
838 *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
839 DNekMatSharedPtr mat = GenMatrix(hkey);
840
841 mat->Invert();
843 }
844 break;
846 {
847 NekDouble one = 1.0;
848 MatrixKey helmkey(StdRegions::eHelmholtz, mkey.GetShapeType(),
849 *this, mkey.GetConstFactors(),
850 mkey.GetVarCoeffs());
851 DNekScalBlkMatSharedPtr helmStatCond =
852 GetLocStaticCondMatrix(helmkey);
853 DNekScalMatSharedPtr A = helmStatCond->GetBlock(0, 0);
855
857 }
858 break;
860 {
861 NekDouble one = 1.0;
862 MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this);
863 DNekScalBlkMatSharedPtr massStatCond =
864 GetLocStaticCondMatrix(masskey);
865 DNekScalMatSharedPtr A = massStatCond->GetBlock(0, 0);
867
869 }
870 break;
872 {
873 NekDouble one = 1.0;
874 MatrixKey helmkey(StdRegions::eHelmholtz, mkey.GetShapeType(),
875 *this, mkey.GetConstFactors(),
876 mkey.GetVarCoeffs());
877 DNekScalBlkMatSharedPtr helmStatCond =
878 GetLocStaticCondMatrix(helmkey);
879
880 DNekScalMatSharedPtr A = helmStatCond->GetBlock(0, 0);
881
883
885 BuildTransformationMatrix(A, mkey.GetMatrixType());
886
888
889 // Free memory,
890 // need to build keys to delete Mass and Laplacian
891 MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this);
892 MatrixKey lapkey(StdRegions::eLaplacian, mkey.GetShapeType(), *this,
893 mkey.GetConstFactors(), mkey.GetVarCoeffs());
895 DropLocMatrix(helmkey);
896 DropLocMatrix(lapkey);
897 DropLocMatrix(masskey);
898 }
899 break;
901 {
902 NekDouble one = 1.0;
903 MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this);
905 DNekScalMatSharedPtr A = StatCond->GetBlock(0, 0);
906
909 BuildTransformationMatrix(A, mkey.GetMatrixType());
910
912 }
913 break;
914 default:
915 {
916 NekDouble one = 1.0;
917 DNekMatSharedPtr mat = GenMatrix(mkey);
918
920 }
921 break;
922 }
923
924 return returnval;
925}
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed to...
DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey) override
void DropLocMatrix(const LocalRegions::MatrixKey &mkey)
Definition Expansion.cpp:90
DNekMatSharedPtr BuildVertexMatrix(const DNekScalMatSharedPtr &r_bnd)
SpatialDomains::GeomFactorsSharedPtr m_metricinfo
Definition Expansion.h:280
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)
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
const LibUtilities::PointsKeyVector GetPointsKeys() const
DNekScalBlkMatSharedPtr GetLocStaticCondMatrix(const LocalRegions::MatrixKey &mkey)
DNekMatSharedPtr GenMatrix(const StdMatrixKey &mkey)
std::vector< PointsKey > PointsKeyVector
Definition Points.h:231
@ eNoGeomType
No type defined.
@ eDeformed
Geometry is curved or has non-constant factors.
static VarCoeffMap NullVarCoeffMap
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
NekMatrix< NekDouble, StandardMatrixTag > DNekMat
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
std::shared_ptr< DNekMat > DNekMatSharedPtr
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::eLinearAdvectionReaction, 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 2545 of file Expansion3D.cpp.

2546{
2547 int n, j;
2548 int nEdgeCoeffs;
2549 int nBndCoeffs = NumBndryCoeffs();
2550
2551 Array<OneD, unsigned int> bmap(nBndCoeffs);
2552 GetBoundaryMap(bmap);
2553
2554 // Map from full system to statically condensed system (i.e reverse
2555 // GetBoundaryMap)
2556 map<int, int> invmap;
2557 for (j = 0; j < nBndCoeffs; ++j)
2558 {
2559 invmap[bmap[j]] = j;
2560 }
2561
2562 // Number of interior edge coefficients
2563 nEdgeCoeffs = GetEdgeNcoeffs(eid) - 2;
2564
2565 const SpatialDomains::Geometry3D *geom = GetGeom3D();
2566
2567 Array<OneD, unsigned int> edgemaparray(nEdgeCoeffs);
2568 StdRegions::Orientation eOrient = geom->GetEorient(eid);
2569 Array<OneD, unsigned int> maparray = Array<OneD, unsigned int>(nEdgeCoeffs);
2570 Array<OneD, int> signarray = Array<OneD, int>(nEdgeCoeffs, 1);
2571
2572 // maparray is the location of the edge within the matrix
2573 GetEdgeInteriorToElementMap(eid, maparray, signarray, eOrient);
2574
2575 for (n = 0; n < nEdgeCoeffs; ++n)
2576 {
2577 edgemaparray[n] = invmap[maparray[n]];
2578 }
2579
2580 return edgemaparray;
2581}
SpatialDomains::Geometry3D * GetGeom3D() const
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)

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

Referenced by v_BuildInverseTransformationMatrix(), and v_BuildTransformationMatrix().

◆ GetGeom3D()

SpatialDomains::Geometry3D * Nektar::LocalRegions::Expansion3D::GetGeom3D ( ) const
inline

Definition at line 176 of file Expansion3D.h.

177{
178 return static_cast<SpatialDomains::Geometry3D *>(m_geom);
179}
SpatialDomains::Geometry * m_geom
Definition Expansion.h:279

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

2701{
2702 int n, j;
2703 int nEdgeCoeffs;
2704 int nFaceCoeffs;
2705
2706 int nBndCoeffs = NumBndryCoeffs();
2707
2708 Array<OneD, unsigned int> bmap(nBndCoeffs);
2709 GetBoundaryMap(bmap);
2710
2711 // Map from full system to statically condensed system (i.e reverse
2712 // GetBoundaryMap)
2713 map<int, int> reversemap;
2714 for (j = 0; j < bmap.size(); ++j)
2715 {
2716 reversemap[bmap[j]] = j;
2717 }
2718
2719 int nverts = GetNverts();
2720 vmap = Array<OneD, unsigned int>(nverts);
2721 for (n = 0; n < nverts; ++n)
2722 {
2723 int id = GetVertexMap(n);
2724 vmap[n] = reversemap[id]; // not sure what should be true here.
2725 }
2726
2727 int nedges = GetNedges();
2728 emap = Array<OneD, Array<OneD, unsigned int>>(nedges);
2729
2730 for (int eid = 0; eid < nedges; ++eid)
2731 {
2732 // Number of interior edge coefficients
2733 nEdgeCoeffs = GetEdgeNcoeffs(eid) - 2;
2734
2735 Array<OneD, unsigned int> edgemaparray(nEdgeCoeffs);
2736 Array<OneD, unsigned int> maparray =
2737 Array<OneD, unsigned int>(nEdgeCoeffs);
2738 Array<OneD, int> signarray = Array<OneD, int>(nEdgeCoeffs, 1);
2739
2740 // maparray is the location of the edge within the matrix
2741 GetEdgeInteriorToElementMap(eid, maparray, signarray,
2743
2744 for (n = 0; n < nEdgeCoeffs; ++n)
2745 {
2746 edgemaparray[n] = reversemap[maparray[n]];
2747 }
2748 emap[eid] = edgemaparray;
2749 }
2750
2751 int nfaces = GetNtraces();
2752 fmap = Array<OneD, Array<OneD, unsigned int>>(nfaces);
2753
2754 for (int fid = 0; fid < nfaces; ++fid)
2755 {
2756 // Number of interior face coefficients
2757 nFaceCoeffs = GetTraceIntNcoeffs(fid);
2758
2759 Array<OneD, unsigned int> facemaparray(nFaceCoeffs);
2760 Array<OneD, unsigned int> maparray =
2761 Array<OneD, unsigned int>(nFaceCoeffs);
2762 Array<OneD, int> signarray = Array<OneD, int>(nFaceCoeffs, 1);
2763
2764 // maparray is the location of the face within the matrix
2765 GetTraceInteriorToElementMap(fid, maparray, signarray,
2767
2768 for (n = 0; n < nFaceCoeffs; ++n)
2769 {
2770 facemaparray[n] = reversemap[maparray[n]];
2771 }
2772
2773 fmap[fid] = facemaparray;
2774 }
2775}
int GetNedges() const
return the number of edges in 3D expansion
int GetTraceIntNcoeffs(const int i) const
int GetVertexMap(const int localVertexId, bool useCoeffPacking=false)
int GetNverts() const
This function returns the number of vertices of the expansion domain.
void GetTraceInteriorToElementMap(const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, const Orientation traceOrient=eForwards)

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

3006{
3007 int nquad_f = FaceExp_f->GetNumPoints(0) * FaceExp_f->GetNumPoints(1);
3008 int coordim = GetCoordim();
3009
3010 int nquad0 = m_base[0]->GetNumPoints();
3011 int nquad1 = m_base[1]->GetNumPoints();
3012 int nquad2 = m_base[2]->GetNumPoints();
3013 int nqtot = nquad0 * nquad1 * nquad2;
3014
3015 StdRegions::VarCoeffType MMFCoeffs[15] = {
3024
3025 StdRegions::VarCoeffMap::const_iterator MFdir;
3026
3027 Array<OneD, NekDouble> nFacecdotMF(nqtot, 0.0);
3028 Array<OneD, NekDouble> tmp(nqtot);
3029 Array<OneD, NekDouble> tmp_f(nquad_f);
3030 for (int k = 0; k < coordim; k++)
3031 {
3032 MFdir = varcoeffs.find(MMFCoeffs[dir * 5 + k]);
3033 tmp = MFdir->second.GetValue();
3034
3035 GetPhysFaceVarCoeffsFromElement(face, FaceExp_f, tmp, tmp_f);
3036
3037 Vmath::Vvtvp(nquad_f, &tmp_f[0], 1, &normals[k][0], 1, &nFacecdotMF[0],
3038 1, &nFacecdotMF[0], 1);
3039 }
3040
3041 return nFacecdotMF;
3042}
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)
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 2583 of file Expansion3D.cpp.

2585{
2586 int n, j;
2587 int nFaceCoeffs;
2588
2589 int nBndCoeffs = NumBndryCoeffs();
2590
2591 Array<OneD, unsigned int> bmap(nBndCoeffs);
2592 GetBoundaryMap(bmap);
2593
2594 // Map from full system to statically condensed system (i.e reverse
2595 // GetBoundaryMap)
2596 map<int, int> reversemap;
2597 for (j = 0; j < bmap.size(); ++j)
2598 {
2599 reversemap[bmap[j]] = j;
2600 }
2601
2602 // Number of interior face coefficients
2603 nFaceCoeffs = GetTraceIntNcoeffs(fid);
2604
2606 Array<OneD, unsigned int> maparray = Array<OneD, unsigned int>(nFaceCoeffs);
2607 Array<OneD, int> signarray = Array<OneD, int>(nFaceCoeffs, 1);
2608
2609 if (faceOrient == StdRegions::eNoOrientation)
2610 {
2611 fOrient = GetTraceOrient(fid);
2612 }
2613 else
2614 {
2615 fOrient = faceOrient;
2616 }
2617
2618 // maparray is the location of the face within the matrix
2619 GetTraceInteriorToElementMap(fid, maparray, signarray, fOrient);
2620
2621 Array<OneD, unsigned int> facemaparray;
2622 int locP1, locP2;
2623 GetTraceNumModes(fid, locP1, locP2, fOrient);
2624
2625 if (P1 == -1)
2626 {
2627 P1 = locP1;
2628 }
2629 else
2630 {
2631 ASSERTL1(P1 <= locP1, "Expect value of passed P1 to "
2632 "be lower or equal to face num modes");
2633 }
2634
2635 if (P2 == -1)
2636 {
2637 P2 = locP2;
2638 }
2639 else
2640 {
2641 ASSERTL1(P2 <= locP2, "Expect value of passed P2 to "
2642 "be lower or equal to face num modes");
2643 }
2644
2645 switch (GetGeom3D()->GetFace(fid)->GetShapeType())
2646 {
2648 {
2649 if (((P1 - 3) > 0) && ((P2 - 3) > 0))
2650 {
2651 facemaparray = Array<OneD, unsigned int>(
2653 P2 - 3));
2654 int cnt = 0;
2655 int cnt1 = 0;
2656 for (n = 0; n < P1 - 3; ++n)
2657 {
2658 for (int m = 0; m < P2 - 3 - n; ++m, ++cnt)
2659 {
2660 facemaparray[cnt] = reversemap[maparray[cnt1 + m]];
2661 }
2662 cnt1 += locP2 - 3 - n;
2663 }
2664 }
2665 }
2666 break;
2668 {
2669 if (((P1 - 2) > 0) && ((P2 - 2) > 0))
2670 {
2671 facemaparray = Array<OneD, unsigned int>(
2673 P2 - 2));
2674 int cnt = 0;
2675 int cnt1 = 0;
2676 for (n = 0; n < P2 - 2; ++n)
2677 {
2678 for (int m = 0; m < P1 - 2; ++m, ++cnt)
2679 {
2680 facemaparray[cnt] = reversemap[maparray[cnt1 + m]];
2681 }
2682 cnt1 += locP1 - 2;
2683 }
2684 }
2685 }
2686 break;
2687 default:
2688 {
2689 ASSERTL0(false, "Invalid shape type.");
2690 }
2691 break;
2692 }
2693
2694 return facemaparray;
2695}
#define ASSERTL0(condition, msg)
void GetTraceNumModes(const int tid, int &numModes0, int &numModes1, const Orientation traceOrient=eDir1FwdDir1_Dir2FwdDir2)
int getNumberOfCoefficients(int Na, int Nb)
int getNumberOfCoefficients(int Na, int Nb)

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

1708{
1709 int i, j;
1710
1711 /*
1712 * Coming into this routine, the velocity V will have been
1713 * multiplied by the trace normals to give the input vector Vn. By
1714 * convention, these normals are inwards facing for elements which
1715 * have FaceExp as their right-adjacent face. This conditional
1716 * statement therefore determines whether the normals must be
1717 * negated, since the integral being performed here requires an
1718 * outwards facing normal.
1719 */
1720 if (m_requireNeg.size() == 0)
1721 {
1722 m_requireNeg.resize(GetNtraces());
1723
1724 for (i = 0; i < GetNtraces(); ++i)
1725 {
1726 m_requireNeg[i] = false;
1727
1728 ExpansionSharedPtr faceExp = m_traceExp[i].lock();
1729
1730 if (faceExp->GetRightAdjacentElementExp())
1731 {
1732 if (faceExp->GetRightAdjacentElementExp()
1733 ->GetGeom()
1734 ->GetGlobalID() == GetGeom()->GetGlobalID())
1735 {
1736 m_requireNeg[i] = true;
1737 }
1738 }
1739 }
1740 }
1741
1742 IndexMapKey ikey(eFaceToElement, DetShapeType(), GetBasisNumModes(0),
1744 GetTraceOrient(face));
1746
1747 int order_e = (*map).size(); // Order of the element
1748 int n_coeffs = FaceExp->GetNcoeffs();
1749
1750 Array<OneD, NekDouble> faceCoeffs(n_coeffs);
1751
1752 if (n_coeffs != order_e) // Going to orthogonal space
1753 {
1754 FaceExp->FwdTrans(Fn, faceCoeffs);
1755
1756 int NumModesElementMax = FaceExp->GetBasis(0)->GetNumModes();
1757 int NumModesElementMin = m_base[0]->GetNumModes();
1758
1759 FaceExp->ReduceOrderCoeffs(NumModesElementMin, faceCoeffs, faceCoeffs);
1760
1761 StdRegions::StdMatrixKey masskey(StdRegions::eMass,
1762 FaceExp->DetShapeType(), *FaceExp);
1763 FaceExp->MassMatrixOp(faceCoeffs, faceCoeffs, masskey);
1764
1765 // Reorder coefficients for the lower degree face.
1766 int offset1 = 0, offset2 = 0;
1767
1768 if (FaceExp->DetShapeType() == LibUtilities::eQuadrilateral)
1769 {
1770 for (i = 0; i < NumModesElementMin; ++i)
1771 {
1772 for (j = 0; j < NumModesElementMin; ++j)
1773 {
1774 faceCoeffs[offset1 + j] = faceCoeffs[offset2 + j];
1775 }
1776 offset1 += NumModesElementMin;
1777 offset2 += NumModesElementMax;
1778 }
1779
1780 // Extract lower degree modes. TODO: Check this is correct.
1781 for (i = NumModesElementMin; i < NumModesElementMax; ++i)
1782 {
1783 for (j = NumModesElementMin; j < NumModesElementMax; ++j)
1784 {
1785 faceCoeffs[i * NumModesElementMax + j] = 0.0;
1786 }
1787 }
1788 }
1789
1790 if (FaceExp->DetShapeType() == LibUtilities::eTriangle)
1791 {
1792
1793 // Reorder coefficients for the lower degree face.
1794 int offset1 = 0, offset2 = 0;
1795
1796 for (i = 0; i < NumModesElementMin; ++i)
1797 {
1798 for (j = 0; j < NumModesElementMin - i; ++j)
1799 {
1800 faceCoeffs[offset1 + j] = faceCoeffs[offset2 + j];
1801 }
1802 offset1 += NumModesElementMin - i;
1803 offset2 += NumModesElementMax - i;
1804 }
1805 }
1806 }
1807 else
1808 {
1809 FaceExp->IProductWRTBase(Fn, faceCoeffs);
1810 }
1811
1812 if (m_requireNeg[face])
1813 {
1814 for (i = 0; i < order_e; ++i)
1815 {
1816 outarray[(*map)[i].index] -= (*map)[i].sign * faceCoeffs[i];
1817 }
1818 }
1819 else
1820 {
1821 for (i = 0; i < order_e; ++i)
1822 {
1823 outarray[(*map)[i].index] += (*map)[i].sign * faceCoeffs[i];
1824 }
1825 }
1826}
SpatialDomains::Geometry * GetGeom() const
std::map< int, ExpansionWeakPtr > m_traceExp
Definition Expansion.h:278
int GetGlobalID(void) const
Get the ID of this object.
Definition Geometry.h:322

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::SpatialDomains::Geometry::GetGlobalID(), 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 1862 of file Expansion3D.cpp.

1865{
1867 "Not set up for non boundary-interior expansions");
1868 ASSERTL1(inoutmat->GetRows() == inoutmat->GetColumns(),
1869 "Assuming that input matrix was square");
1870
1871 int i, j;
1872 int id1, id2;
1873 ExpansionSharedPtr faceExp = m_traceExp[face].lock();
1874 int order_f = faceExp->GetNcoeffs();
1875
1876 Array<OneD, unsigned int> map;
1877 Array<OneD, int> sign;
1878
1879 StdRegions::VarCoeffMap varcoeffs;
1880 varcoeffs[StdRegions::eVarCoeffMass] = primCoeffs;
1881
1882 LibUtilities::ShapeType shapeType = faceExp->DetShapeType();
1883
1884 LocalRegions::MatrixKey mkey(StdRegions::eMass, shapeType, *faceExp,
1886
1887 DNekScalMat &facemat = *faceExp->GetLocMatrix(mkey);
1888
1889 // Now need to identify a map which takes the local face
1890 // mass matrix to the matrix stored in inoutmat;
1891 // This can currently be deduced from the size of the matrix
1892
1893 // - if inoutmat.m_rows() == v_NCoeffs() it is a full
1894 // matrix system
1895
1896 // - if inoutmat.m_rows() == v_GetNverts() it is a vertex space
1897 // preconditioner.
1898
1899 // - if inoutmat.m_rows() == v_NumBndCoeffs() it is a
1900 // boundary CG system
1901
1902 // - if inoutmat.m_rows() == v_NumDGBndCoeffs() it is a
1903 // trace DG system; still needs implementing.
1904 int rows = inoutmat->GetRows();
1905
1906 if (rows == GetNcoeffs())
1907 {
1908 GetTraceToElementMap(face, map, sign, GetTraceOrient(face));
1909 }
1910 else if (rows == GetNverts())
1911 {
1912 int nfvert = faceExp->GetNverts();
1913
1914 // Need to find where linear vertices are in facemat
1915 Array<OneD, unsigned int> linmap;
1916 Array<OneD, int> linsign;
1917
1918 // Use a linear expansion to get correct mapping
1919 GetLinStdExp()->GetTraceToElementMap(face, linmap, linsign,
1920 GetTraceOrient(face));
1921
1922 // zero out sign map to remove all other modes
1923 sign = Array<OneD, int>(order_f, 0);
1924 map = Array<OneD, unsigned int>(order_f, (unsigned int)0);
1925
1926 int fmap;
1927 // Reset sign map to only have contribution from vertices
1928 for (i = 0; i < nfvert; ++i)
1929 {
1930 fmap = faceExp->GetVertexMap(i, true);
1931 sign[fmap] = 1;
1932
1933 // need to reset map
1934 map[fmap] = linmap[i];
1935 }
1936 }
1937 else if (rows == NumBndryCoeffs())
1938 {
1939 int nbndry = NumBndryCoeffs();
1940 Array<OneD, unsigned int> bmap(nbndry);
1941
1942 GetTraceToElementMap(face, map, sign, GetTraceOrient(face));
1943 GetBoundaryMap(bmap);
1944
1945 for (i = 0; i < order_f; ++i)
1946 {
1947 for (j = 0; j < nbndry; ++j)
1948 {
1949 if (map[i] == bmap[j])
1950 {
1951 map[i] = j;
1952 break;
1953 }
1954 }
1955 ASSERTL1(j != nbndry, "Did not find number in map");
1956 }
1957 }
1958 else if (rows == NumDGBndryCoeffs())
1959 {
1960 // possibly this should be a separate method
1961 int cnt = 0;
1962 map = Array<OneD, unsigned int>(order_f);
1963 sign = Array<OneD, int>(order_f, 1);
1964
1965 IndexMapKey ikey1(eFaceToElement, DetShapeType(), GetBasisNumModes(0),
1967 GetTraceOrient(face));
1969 IndexMapKey ikey2(eFaceToElement, DetShapeType(), GetBasisNumModes(0),
1973
1974 ASSERTL1((*map1).size() == (*map2).size(),
1975 "There is an error with the GetTraceToElementMap");
1976
1977 for (i = 0; i < face; ++i)
1978 {
1979 cnt += GetTraceNcoeffs(i);
1980 }
1981
1982 for (i = 0; i < (*map1).size(); ++i)
1983 {
1984 int idx = -1;
1985
1986 for (j = 0; j < (*map2).size(); ++j)
1987 {
1988 if ((*map1)[i].index == (*map2)[j].index)
1989 {
1990 idx = j;
1991 break;
1992 }
1993 }
1994
1995 ASSERTL2(idx >= 0, "Index not found");
1996 map[i] = idx + cnt;
1997 sign[i] = (*map2)[idx].sign;
1998 }
1999 }
2000 else
2001 {
2002 ASSERTL0(false, "Could not identify matrix type from dimension");
2003 }
2004
2005 for (i = 0; i < order_f; ++i)
2006 {
2007 id1 = map[i];
2008 for (j = 0; j < order_f; ++j)
2009 {
2010 id2 = map[j];
2011 (*inoutmat)(id1, id2) += facemat(i, j) * sign[i] * sign[j];
2012 }
2013 }
2014}
std::shared_ptr< StdExpansion > GetLinStdExp(void) const

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

2430{
2431 int i, j, n, eid = 0, fid = 0;
2432 int nCoeffs = NumBndryCoeffs();
2433 NekDouble MatrixValue;
2434 DNekScalMat &R = (*transformationmatrix);
2435
2436 // Define storage for vertex transpose matrix and zero all entries
2437 MatrixStorage storage = eFULL;
2438
2439 // Inverse transformation matrix
2440 DNekMatSharedPtr inversetransformationmatrix =
2441 MemoryManager<DNekMat>::AllocateSharedPtr(nCoeffs, nCoeffs, 0.0,
2442 storage);
2443 DNekMat &InvR = (*inversetransformationmatrix);
2444
2445 int nVerts = GetNverts();
2446 int nEdges = GetNedges();
2447 int nFaces = GetNtraces();
2448
2449 int nedgemodes = 0;
2450 int nfacemodes = 0;
2451 int nedgemodestotal = 0;
2452 int nfacemodestotal = 0;
2453
2454 for (eid = 0; eid < nEdges; ++eid)
2455 {
2456 nedgemodes = GetEdgeNcoeffs(eid) - 2;
2457 nedgemodestotal += nedgemodes;
2458 }
2459
2460 for (fid = 0; fid < nFaces; ++fid)
2461 {
2462 nfacemodes = GetTraceIntNcoeffs(fid);
2463 nfacemodestotal += nfacemodes;
2464 }
2465
2466 Array<OneD, unsigned int> edgemodearray(nedgemodestotal);
2467 Array<OneD, unsigned int> facemodearray(nfacemodestotal);
2468
2469 int offset = 0;
2470
2471 // Create array of edge modes
2472 for (eid = 0; eid < nEdges; ++eid)
2473 {
2474 Array<OneD, unsigned int> edgearray = GetEdgeInverseBoundaryMap(eid);
2475 nedgemodes = GetEdgeNcoeffs(eid) - 2;
2476
2477 // Only copy if there are edge modes
2478 if (nedgemodes)
2479 {
2480 Vmath::Vcopy(nedgemodes, &edgearray[0], 1, &edgemodearray[offset],
2481 1);
2482 }
2483
2484 offset += nedgemodes;
2485 }
2486
2487 offset = 0;
2488
2489 // Create array of face modes
2490 for (fid = 0; fid < nFaces; ++fid)
2491 {
2492 Array<OneD, unsigned int> facearray = GetTraceInverseBoundaryMap(fid);
2493 nfacemodes = GetTraceIntNcoeffs(fid);
2494
2495 // Only copy if there are face modes
2496 if (nfacemodes)
2497 {
2498 Vmath::Vcopy(nfacemodes, &facearray[0], 1, &facemodearray[offset],
2499 1);
2500 }
2501
2502 offset += nfacemodes;
2503 }
2504
2505 // Vertex-edge/face
2506 for (i = 0; i < nVerts; ++i)
2507 {
2508 for (j = 0; j < nedgemodestotal; ++j)
2509 {
2510 InvR.SetValue(GetVertexMap(i), edgemodearray[j],
2511 -R(GetVertexMap(i), edgemodearray[j]));
2512 }
2513 for (j = 0; j < nfacemodestotal; ++j)
2514 {
2515 InvR.SetValue(GetVertexMap(i), facemodearray[j],
2516 -R(GetVertexMap(i), facemodearray[j]));
2517 for (n = 0; n < nedgemodestotal; ++n)
2518 {
2519 MatrixValue = InvR.GetValue(GetVertexMap(i), facemodearray[j]) +
2520 R(GetVertexMap(i), edgemodearray[n]) *
2521 R(edgemodearray[n], facemodearray[j]);
2522 InvR.SetValue(GetVertexMap(i), facemodearray[j], MatrixValue);
2523 }
2524 }
2525 }
2526
2527 // Edge-face contributions
2528 for (i = 0; i < nedgemodestotal; ++i)
2529 {
2530 for (j = 0; j < nfacemodestotal; ++j)
2531 {
2532 InvR.SetValue(edgemodearray[i], facemodearray[j],
2533 -R(edgemodearray[i], facemodearray[j]));
2534 }
2535 }
2536
2537 for (i = 0; i < nCoeffs; ++i)
2538 {
2539 InvR.SetValue(i, i, 1.0);
2540 }
2541
2542 return inversetransformationmatrix;
2543}
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 2046 of file Expansion3D.cpp.

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

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::ErrorUtil::efatal, Nektar::eFULL, Nektar::StdRegions::ePreconR, Nektar::StdRegions::ePreconRMass, Nektar::SpatialDomains::Geometry::GetEdgeFaceMap(), GetEdgeInverseBoundaryMap(), Nektar::StdRegions::StdExpansion3D::GetEdgeNcoeffs(), GetGeom3D(), Nektar::StdRegions::StdExpansion3D::GetNedges(), Nektar::StdRegions::StdExpansion::GetNverts(), Nektar::StdRegions::StdExpansion::GetTraceIntNcoeffs(), GetTraceInverseBoundaryMap(), Nektar::SpatialDomains::Geometry::GetVertexEdgeMap(), Nektar::SpatialDomains::Geometry::GetVertexFaceMap(), 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 2016 of file Expansion3D.cpp.

2018{
2019 MatrixStorage storage = eFULL;
2020 DNekMatSharedPtr vertexmatrix;
2021
2022 int nVerts, vid1, vid2, vMap1, vMap2;
2023 NekDouble VertexValue;
2024
2025 nVerts = GetNverts();
2026
2027 vertexmatrix =
2028 MemoryManager<DNekMat>::AllocateSharedPtr(nVerts, nVerts, 0.0, storage);
2029 DNekMat &VertexMat = (*vertexmatrix);
2030
2031 for (vid1 = 0; vid1 < nVerts; ++vid1)
2032 {
2033 vMap1 = GetVertexMap(vid1, true);
2034
2035 for (vid2 = 0; vid2 < nVerts; ++vid2)
2036 {
2037 vMap2 = GetVertexMap(vid2, true);
2038 VertexValue = (*r_bnd)(vMap1, vMap2);
2039 VertexMat.SetValue(vid1, vid2, VertexValue);
2040 }
2041 }
2042
2043 return vertexmatrix;
2044}

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

1838{
1839 int ncoeffs = GetNcoeffs();
1843
1845 DNekScalMat &Dmat = *GetLocMatrix(DerivType[dir]);
1846
1847 Array<OneD, NekDouble> coeffs = incoeffs;
1848 DNekVec Coeffs(ncoeffs, coeffs, eWrapper);
1849
1850 Coeffs = Transpose(Dmat) * Coeffs;
1851 Vmath::Neg(ncoeffs, coeffs, 1);
1852
1853 // Add the boundary integral including the relevant part of
1854 // the normal
1855 AddNormTraceInt(dir, FaceExp, faceCoeffs, coeffs);
1856
1857 DNekVec Out_d(ncoeffs, out_d, eWrapper);
1858
1859 Out_d = InvMass * Coeffs;
1860}
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 936 of file Expansion3D.cpp.

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

2840{
2841 SpatialDomains::Geometry *faceGeom = m_geom->GetFace(traceid);
2842 if (faceGeom->GetNumVerts() == 3)
2843 {
2845 GetTraceBasisKey(traceid, 0), GetTraceBasisKey(traceid, 1),
2846 m_geom->GetFace(traceid));
2847 }
2848 else
2849 {
2851 GetTraceBasisKey(traceid, 0), GetTraceBasisKey(traceid, 1),
2852 m_geom->GetFace(traceid));
2853 }
2854}
Geometry2D * GetFace(int i) const
Returns face i of this object.
Definition Geometry.h:377

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::SpatialDomains::Geometry::GetFace(), Nektar::SpatialDomains::Geometry::GetNumVerts(), 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 2777 of file Expansion3D.cpp.

2778{
2779 return m_geom->GetForient(face);
2780}
StdRegions::Orientation GetForient(const int i) const
Returns the orientation of face i with respect to the ordering of faces in the standard element.
Definition Geometry.h:395

References Nektar::SpatialDomains::Geometry::GetForient(), and 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 2787 of file Expansion3D.cpp.

2791{
2792 if (orient == StdRegions::eNoOrientation)
2793 {
2794 orient = GetTraceOrient(face);
2795 }
2796
2797 int nq0 = FaceExp->GetNumPoints(0);
2798 int nq1 = FaceExp->GetNumPoints(1);
2799
2800 int nfacepts = GetTraceNumPoints(face);
2801 int dir0 = GetGeom3D()->GetDir(face, 0);
2802 int dir1 = GetGeom3D()->GetDir(face, 1);
2803
2804 Array<OneD, NekDouble> o_tmp(nfacepts);
2805 Array<OneD, NekDouble> o_tmp2(FaceExp->GetTotPoints());
2806 Array<OneD, int> faceids;
2807
2808 // Get local face pts and put into o_tmp
2809 GetTracePhysMap(face, faceids);
2810 // The static cast is necessary because faceids should be
2811 // Array<OneD, size_t> faceids ... or at leaste the same type as
2812 // faceids.size() ...
2813 Vmath::Gathr(static_cast<int>(faceids.size()), inarray, faceids, o_tmp);
2814
2815 int to_id0, to_id1;
2816
2818 {
2819 to_id0 = 0;
2820 to_id1 = 1;
2821 }
2822 else // transpose points key evaluation
2823 {
2824 to_id0 = 1;
2825 to_id1 = 0;
2826 }
2827
2828 // interpolate to points distrbution given in FaceExp
2830 m_base[dir0]->GetPointsKey(), m_base[dir1]->GetPointsKey(),
2831 o_tmp.data(), FaceExp->GetBasis(to_id0)->GetPointsKey(),
2832 FaceExp->GetBasis(to_id1)->GetPointsKey(), o_tmp2.data());
2833
2834 // Reshuffule points as required and put into outarray.
2835 v_ReOrientTracePhysMap(orient, faceids, nq0, nq1);
2836 Vmath::Scatr(nq0 * nq1, o_tmp2, faceids, outarray);
2837}
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 GetDir(const int i, const int j=0) const
Returns the element coordinate direction corresponding to a given face coordinate direction.
Definition Geometry.h:661
int GetTraceNumPoints(const int i) const
This function returns the number of quadrature points belonging to the i-th trace.
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(), Nektar::SpatialDomains::Geometry::GetDir(), 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 2994 of file Expansion3D.cpp.

2997{
2998 NormVectorIProductWRTBase(Fvec[0], Fvec[1], Fvec[2], outarray);
2999}
void NormVectorIProductWRTBase(const Array< OneD, const NekDouble > &Fx, Array< OneD, NekDouble > &outarray)

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

2859{
2860
2861 if (idmap.size() != nq0 * nq1)
2862 {
2863 idmap = Array<OneD, int>(nq0 * nq1);
2864 }
2865
2866 if (GetNverts() == 3) // Tri face
2867 {
2868 switch (orient)
2869 {
2871 // eseentially straight copy
2872 for (int i = 0; i < nq0 * nq1; ++i)
2873 {
2874 idmap[i] = i;
2875 }
2876 break;
2878 // reverse
2879 for (int j = 0; j < nq1; ++j)
2880 {
2881 for (int i = 0; i < nq0; ++i)
2882 {
2883 idmap[j * nq0 + i] = nq0 - 1 - i + j * nq0;
2884 }
2885 }
2886 break;
2887 default:
2888 ASSERTL0(false,
2889 "Case not supposed to be used in this function");
2890 }
2891 }
2892 else
2893 {
2894 switch (orient)
2895 {
2897 // eseentially straight copy
2898 for (int i = 0; i < nq0 * nq1; ++i)
2899 {
2900 idmap[i] = i;
2901 }
2902 break;
2904 {
2905 // Direction A negative and B positive
2906 for (int j = 0; j < nq1; j++)
2907 {
2908 for (int i = 0; i < nq0; ++i)
2909 {
2910 idmap[j * nq0 + i] = nq0 - 1 - i + j * nq0;
2911 }
2912 }
2913 }
2914 break;
2916 {
2917 // Direction A positive and B negative
2918 for (int j = 0; j < nq1; j++)
2919 {
2920 for (int i = 0; i < nq0; ++i)
2921 {
2922 idmap[j * nq0 + i] = nq0 * (nq1 - 1 - j) + i;
2923 }
2924 }
2925 }
2926 break;
2928 {
2929 // Direction A negative and B negative
2930 for (int j = 0; j < nq1; j++)
2931 {
2932 for (int i = 0; i < nq0; ++i)
2933 {
2934 idmap[j * nq0 + i] = nq0 * nq1 - 1 - j * nq0 - i;
2935 }
2936 }
2937 }
2938 break;
2940 {
2941 // Transposed, Direction A and B positive
2942 for (int i = 0; i < nq0; ++i)
2943 {
2944 for (int j = 0; j < nq1; ++j)
2945 {
2946 idmap[i * nq1 + j] = i + j * nq0;
2947 }
2948 }
2949 }
2950 break;
2952 {
2953 // Transposed, Direction A positive and B negative
2954 for (int i = 0; i < nq0; ++i)
2955 {
2956 for (int j = 0; j < nq1; ++j)
2957 {
2958 idmap[i * nq1 + j] = nq0 - 1 - i + j * nq0;
2959 }
2960 }
2961 }
2962 break;
2964 {
2965 // Transposed, Direction A negative and B positive
2966 for (int i = 0; i < nq0; ++i)
2967 {
2968 for (int j = 0; j < nq1; ++j)
2969 {
2970 idmap[i * nq1 + j] = i + nq0 * (nq1 - 1) - j * nq0;
2971 }
2972 }
2973 }
2974 break;
2976 {
2977 // Transposed, Direction A and B negative
2978 for (int i = 0; i < nq0; ++i)
2979 {
2980 for (int j = 0; j < nq1; ++j)
2981 {
2982 idmap[i * nq1 + j] = nq0 * nq1 - 1 - i - j * nq0;
2983 }
2984 }
2985 }
2986 break;
2987 default:
2988 ASSERTL0(false, "Unknow orientation");
2989 break;
2990 }
2991 }
2992}

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

3045{
3046 SpatialDomains::Geometry *geom = GetGeom();
3047
3048 int nverts = geom->GetFace(traceid)->GetNumVerts();
3049
3050 SpatialDomains::PointGeom tn1, tn2, normal;
3051 tn1.Sub(*(geom->GetFace(traceid)->GetVertex(1)),
3052 *(geom->GetFace(traceid)->GetVertex(0)));
3053
3054 tn2.Sub(*(geom->GetFace(traceid)->GetVertex(nverts - 1)),
3055 *(geom->GetFace(traceid)->GetVertex(0)));
3056
3057 normal.Mult(tn1, tn2);
3058
3059 // normalise normal
3060 NekDouble mag = normal.dot(normal);
3061 mag = 1.0 / sqrt(mag);
3062 normal.UpdatePosition(normal.x() * mag, normal.y() * mag, normal.z() * mag);
3063
3064 SpatialDomains::PointGeom Dx;
3065 h = 0.0;
3066 p = 0.0;
3067 for (int i = 0; i < nverts; ++i)
3068 {
3069 // vertices on edges
3070 int edgid = geom->GetEdgeNormalToFaceVert(traceid, i);
3071
3072 // vector along noramal edge to each vertex
3073 Dx.Sub(*(geom->GetEdge(edgid)->GetVertex(0)),
3074 *(geom->GetEdge(edgid)->GetVertex(1)));
3075
3076 // calculate perpendicular distance of normal length
3077 // from first vertex
3078 h += fabs(normal.dot(Dx));
3079 }
3080
3081 h /= static_cast<NekDouble>(nverts);
3082
3083 // find normal basis direction
3084 int dir0 = geom->GetDir(traceid, 0);
3085 int dir1 = geom->GetDir(traceid, 1);
3086 int dirn;
3087 for (dirn = 0; dirn < 3; ++dirn)
3088 {
3089 if ((dirn != dir0) && (dirn != dir1))
3090 {
3091 break;
3092 }
3093 }
3094 p = (NekDouble)(GetBasisNumModes(dirn) - 1);
3095}
scalarT< T > sqrt(scalarT< T > in)
Definition scalar.hpp:290

References Nektar::SpatialDomains::PointGeom::dot(), Nektar::StdRegions::StdExpansion::GetBasisNumModes(), Nektar::SpatialDomains::Geometry::GetDir(), Nektar::SpatialDomains::Geometry::GetEdge(), Nektar::SpatialDomains::Geometry::GetEdgeNormalToFaceVert(), Nektar::SpatialDomains::Geometry::GetFace(), Nektar::LocalRegions::Expansion::GetGeom(), Nektar::SpatialDomains::Geometry::GetNumVerts(), Nektar::SpatialDomains::Geometry::GetVertex(), Nektar::SpatialDomains::PointGeom::Mult(), 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().