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

This class is the abstractio n of a global discontinuous two- dimensional spectral/hp element expansion which approximates the solution of a set of partial differential equations. More...

#include <DisContField.h>

Inheritance diagram for Nektar::MultiRegions::DisContField:
[legend]

Public Member Functions

 DisContField ()
 Default constructor. More...
 
 DisContField (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &graph, const std::string &variable, const bool SetUpJustDG=true, const bool DeclareCoeffPhysArrays=true, const Collections::ImplementationType ImpType=Collections::eNoImpType)
 Constructs a 1D discontinuous field based on a mesh and boundary conditions. More...
 
 DisContField (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &graph1D, const SpatialDomains::CompositeMap &domain, const SpatialDomains::BoundaryConditions &Allbcs, const std::string &variable, bool SetToOneSpaceDimensions=false, const Collections::ImplementationType ImpType=Collections::eNoImpType)
 Constructor for a DisContField from a List of subdomains New Constructor for arterial network. More...
 
 DisContField (const DisContField &In, const bool DeclareCoeffPhysArrays=true)
 Constructs a 1D discontinuous field based on an existing field. More...
 
 DisContField (const DisContField &In, const SpatialDomains::MeshGraphSharedPtr &graph, const std::string &variable, const bool SetUpJustDG=false, const bool DeclareCoeffPhysArrays=true)
 
 DisContField (const ExpList &In)
 Constructs a 1D discontinuous field based on an existing field. (needed in order to use ContField( const ExpList &In) constructor. More...
 
virtual ~DisContField ()
 Destructor. More...
 
GlobalLinSysSharedPtr GetGlobalBndLinSys (const GlobalLinSysKey &mkey)
 For a given key, returns the associated global linear system. More...
 
bool SameTypeOfBoundaryConditions (const DisContField &In)
 Check to see if expansion has the same BCs as In. More...
 
std::vector< bool > & GetNegatedFluxNormal (void)
 
NekDouble L2_DGDeriv (const int dir, const Array< OneD, const NekDouble > &soln)
 Calculate the \( L^2 \) error of the \( Q_{\rm dir} \) derivative using the consistent DG evaluation of \( Q_{\rm dir} \). More...
 
void EvaluateHDGPostProcessing (Array< OneD, NekDouble > &outarray)
 Evaluate HDG post-processing to increase polynomial order of solution. More...
 
- Public Member Functions inherited from Nektar::MultiRegions::ExpList
 ExpList (const ExpansionType Type=eNoType)
 The default constructor using a type. More...
 
 ExpList (const ExpList &in, const bool DeclareCoeffPhysArrays=true)
 The copy constructor. More...
 
 ExpList (const ExpList &in, const std::vector< unsigned int > &eIDs, const bool DeclareCoeffPhysArrays=true, const Collections::ImplementationType ImpType=Collections::eNoImpType)
 Constructor copying only elements defined in eIds. More...
 
 ExpList (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &graph, const bool DeclareCoeffPhysArrays=true, const std::string &var="DefaultVar", const Collections::ImplementationType ImpType=Collections::eNoImpType)
 Generate an ExpList from a meshgraph graph and session file. More...
 
 ExpList (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::ExpansionInfoMap &expansions, const bool DeclareCoeffPhysArrays=true, const Collections::ImplementationType ImpType=Collections::eNoImpType)
 Sets up a list of local expansions based on an expansion Map. More...
 
 ExpList (const SpatialDomains::PointGeomSharedPtr &geom)
 Specialised constructors for 0D Expansions Wrapper around LocalRegion::PointExp - used in PrePacing.cpp. More...
 
 ExpList (const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, const ExpListSharedPtr > &bndConstraint, const Array< OneD, const SpatialDomains ::BoundaryConditionShPtr > &bndCond, const LocalRegions::ExpansionVector &locexp, const SpatialDomains::MeshGraphSharedPtr &graph, const bool DeclareCoeffPhysArrays=true, const std::string variable="DefaultVar", const Collections::ImplementationType ImpType=Collections::eNoImpType)
 Generate expansions for the trace space expansions used in DisContField. More...
 
 ExpList (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::CompositeMap &domain, const SpatialDomains::MeshGraphSharedPtr &graph, const bool DeclareCoeffPhysArrays=true, const std::string variable="DefaultVar", bool SetToOneSpaceDimension=false, const LibUtilities::CommSharedPtr comm=LibUtilities::CommSharedPtr(), const Collections::ImplementationType ImpType=Collections::eNoImpType)
 Constructor based on domain information only for 1D & 2D boundary conditions. More...
 
virtual ~ExpList ()
 The default destructor. More...
 
int GetNcoeffs (void) const
 Returns the total number of local degrees of freedom \(N_{\mathrm{eof}}=\sum_{e=1}^{{N_{\mathrm{el}}}}N^{e}_m\). More...
 
int GetNcoeffs (const int eid) const
 Returns the total number of local degrees of freedom for element eid. More...
 
ExpansionType GetExpType (void)
 Returns the type of the expansion. More...
 
void SetExpType (ExpansionType Type)
 Returns the type of the expansion. More...
 
int EvalBasisNumModesMax (void) const
 Evaulates the maximum number of modes in the elemental basis order over all elements. More...
 
const Array< OneD, int > EvalBasisNumModesMaxPerExp (void) const
 Returns the vector of the number of modes in the elemental basis order over all elements. More...
 
int GetTotPoints (void) const
 Returns the total number of quadrature points m_npoints \(=Q_{\mathrm{tot}}\). More...
 
int GetTotPoints (const int eid) const
 Returns the total number of quadrature points for eid's element \(=Q_{\mathrm{tot}}\). More...
 
int GetNpoints (void) const
 Returns the total number of quadrature points m_npoints \(=Q_{\mathrm{tot}}\). More...
 
int Get1DScaledTotPoints (const NekDouble scale) const
 Returns the total number of qudature points scaled by the factor scale on each 1D direction. More...
 
void SetWaveSpace (const bool wavespace)
 Sets the wave space to the one of the possible configuration true or false. More...
 
void SetModifiedBasis (const bool modbasis)
 Set Modified Basis for the stability analysis. More...
 
void SetPhys (int i, NekDouble val)
 Set the i th value of m_phys to value val. More...
 
bool GetWaveSpace (void) const
 This function returns the third direction expansion condition, which can be in wave space (coefficient) or not It is stored in the variable m_WaveSpace. More...
 
void SetPhys (const Array< OneD, const NekDouble > &inarray)
 Fills the array m_phys. More...
 
void SetPhysArray (Array< OneD, NekDouble > &inarray)
 Sets the array m_phys. More...
 
void SetPhysState (const bool physState)
 This function manually sets whether the array of physical values \(\boldsymbol{u}_l\) (implemented as m_phys) is filled or not. More...
 
bool GetPhysState (void) const
 This function indicates whether the array of physical values \(\boldsymbol{u}_l\) (implemented as m_phys) is filled or not. More...
 
void MultiplyByQuadratureMetric (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 multiply the metric jacobi and quadrature weights More...
 
void DivideByQuadratureMetric (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 Divided by the metric jacobi and quadrature weights. More...
 
void IProductWRTBase_IterPerExp (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 This function calculates the inner product of a function \(f(\boldsymbol{x})\) with respect to all local expansion modes \(\phi_n^e(\boldsymbol{x})\). More...
 
void IProductWRTBase (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
void IProductWRTDerivBase (const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 This function calculates the inner product of a function \(f(\boldsymbol{x})\) with respect to the derivative (in direction. More...
 
void IProductWRTDirectionalDerivBase (const Array< OneD, const NekDouble > &direction, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 Directional derivative along a given direction. More...
 
void IProductWRTDerivBase (const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, NekDouble > &outarray)
 This function calculates the inner product of a function \(f(\boldsymbol{x})\) with respect to the derivative of all local expansion modes \(\phi_n^e(\boldsymbol{x})\). More...
 
void FwdTrans_IterPerExp (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 This function elementally evaluates the forward transformation of a function \(u(\boldsymbol{x})\) onto the global spectral/hp expansion. More...
 
void FwdTrans (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
void ExponentialFilter (Array< OneD, NekDouble > &array, const NekDouble alpha, const NekDouble exponent, const NekDouble cutoff)
 
void MultiplyByElmtInvMass (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 This function elementally mulplies the coefficient space of Sin my the elemental inverse of the mass matrix. More...
 
void MultiplyByInvMassMatrix (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
void MultiplyByMassMatrix (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
void SmoothField (Array< OneD, NekDouble > &field)
 Smooth a field across elements. More...
 
void HelmSolve (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::ConstFactorMap &factors, const StdRegions::VarCoeffMap &varcoeff=StdRegions::NullVarCoeffMap, const MultiRegions::VarFactorsMap &varfactors=MultiRegions::NullVarFactorsMap, const Array< OneD, const NekDouble > &dirForcing=NullNekDouble1DArray, const bool PhysSpaceForcing=true)
 Solve helmholtz problem. More...
 
void LinearAdvectionDiffusionReactionSolve (const Array< OneD, Array< OneD, NekDouble > > &velocity, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const NekDouble lambda, const Array< OneD, const NekDouble > &dirForcing=NullNekDouble1DArray)
 Solve Advection Diffusion Reaction. More...
 
void LinearAdvectionReactionSolve (const Array< OneD, Array< OneD, NekDouble > > &velocity, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const NekDouble lambda, const Array< OneD, const NekDouble > &dirForcing=NullNekDouble1DArray)
 Solve Advection Diffusion Reaction. More...
 
void FwdTrans_BndConstrained (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
void BwdTrans_IterPerExp (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 This function elementally evaluates the backward transformation of the global spectral/hp element expansion. More...
 
void BwdTrans (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
void GetCoords (Array< OneD, NekDouble > &coord_0, Array< OneD, NekDouble > &coord_1=NullNekDouble1DArray, Array< OneD, NekDouble > &coord_2=NullNekDouble1DArray)
 This function calculates the coordinates of all the elemental quadrature points \(\boldsymbol{x}_i\). More...
 
void HomogeneousFwdTrans (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool Shuff=true, bool UnShuff=true)
 
void HomogeneousBwdTrans (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool Shuff=true, bool UnShuff=true)
 
void DealiasedProd (const Array< OneD, NekDouble > &inarray1, const Array< OneD, NekDouble > &inarray2, Array< OneD, NekDouble > &outarray)
 
void DealiasedDotProd (const Array< OneD, Array< OneD, NekDouble > > &inarray1, const Array< OneD, Array< OneD, NekDouble > > &inarray2, Array< OneD, Array< OneD, NekDouble > > &outarray)
 
void GetBCValues (Array< OneD, NekDouble > &BndVals, const Array< OneD, NekDouble > &TotField, int BndID)
 
void NormVectorIProductWRTBase (Array< OneD, const NekDouble > &V1, Array< OneD, const NekDouble > &V2, Array< OneD, NekDouble > &outarray, int BndID)
 
void NormVectorIProductWRTBase (Array< OneD, Array< OneD, NekDouble > > &V, Array< OneD, NekDouble > &outarray)
 
void ApplyGeomInfo ()
 Apply geometry information to each expansion. More...
 
void Reset ()
 Reset geometry information and reset matrices. More...
 
void WriteTecplotHeader (std::ostream &outfile, std::string var="")
 
void WriteTecplotZone (std::ostream &outfile, int expansion=-1)
 
void WriteTecplotField (std::ostream &outfile, int expansion=-1)
 
void WriteTecplotConnectivity (std::ostream &outfile, int expansion=-1)
 
void WriteVtkHeader (std::ostream &outfile)
 
void WriteVtkFooter (std::ostream &outfile)
 
void WriteVtkPieceHeader (std::ostream &outfile, int expansion, int istrip=0)
 
void WriteVtkPieceFooter (std::ostream &outfile, int expansion)
 
void WriteVtkPieceData (std::ostream &outfile, int expansion, std::string var="v")
 
int GetCoordim (int eid)
 This function returns the dimension of the coordinates of the element eid. More...
 
void SetCoeff (int i, NekDouble val)
 Set the i th coefficiient in m_coeffs to value val. More...
 
void SetCoeffs (int i, NekDouble val)
 Set the i th coefficiient in m_coeffs to value val. More...
 
void SetCoeffsArray (Array< OneD, NekDouble > &inarray)
 Set the m_coeffs array to inarray. More...
 
int GetShapeDimension ()
 This function returns the dimension of the shape of the element eid. More...
 
const Array< OneD, const NekDouble > & GetCoeffs () const
 This function returns (a reference to) the array \(\boldsymbol{\hat{u}}_l\) (implemented as m_coeffs) containing all local expansion coefficients. More...
 
void ImposeDirichletConditions (Array< OneD, NekDouble > &outarray)
 Impose Dirichlet Boundary Conditions onto Array. More...
 
void FillBndCondFromField (void)
 Fill Bnd Condition expansion from the values stored in expansion. More...
 
void FillBndCondFromField (const int nreg)
 Fill Bnd Condition expansion in nreg from the values stored in expansion. More...
 
void LocalToGlobal (bool useComm=true)
 Gathers the global coefficients \(\boldsymbol{\hat{u}}_g\) from the local coefficients \(\boldsymbol{\hat{u}}_l\). More...
 
void LocalToGlobal (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool useComm=true)
 
void GlobalToLocal (void)
 Scatters from the global coefficients \(\boldsymbol{\hat{u}}_g\) to the local coefficients \(\boldsymbol{\hat{u}}_l\). More...
 
void GlobalToLocal (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
NekDouble GetCoeff (int i)
 Get the i th value (coefficient) of m_coeffs. More...
 
NekDouble GetCoeffs (int i)
 Get the i th value (coefficient) of m_coeffs. More...
 
const Array< OneD, const NekDouble > & GetPhys () const
 This function returns (a reference to) the array \(\boldsymbol{u}_l\) (implemented as m_phys) containing the function \(u^{\delta}(\boldsymbol{x})\) evaluated at the quadrature points. More...
 
NekDouble Linf (const Array< OneD, const NekDouble > &inarray, const Array< OneD, const NekDouble > &soln=NullNekDouble1DArray)
 This function calculates the \(L_\infty\) error of the global spectral/hp element approximation. More...
 
NekDouble L2 (const Array< OneD, const NekDouble > &inarray, const Array< OneD, const NekDouble > &soln=NullNekDouble1DArray)
 This function calculates the \(L_2\) error with respect to soln of the global spectral/hp element approximation. More...
 
NekDouble H1 (const Array< OneD, const NekDouble > &inarray, const Array< OneD, const NekDouble > &soln=NullNekDouble1DArray)
 Calculates the \(H^1\) error of the global spectral/hp element approximation. More...
 
NekDouble Integral ()
 
NekDouble Integral (const Array< OneD, const NekDouble > &inarray)
 
NekDouble VectorFlux (const Array< OneD, Array< OneD, NekDouble > > &inarray)
 
Array< OneD, const NekDoubleHomogeneousEnergy (void)
 This function calculates the energy associated with each one of the modesof a 3D homogeneous nD expansion. More...
 
void SetHomo1DSpecVanVisc (Array< OneD, NekDouble > visc)
 This function sets the Spectral Vanishing Viscosity in homogeneous1D expansion. More...
 
Array< OneD, const unsigned int > GetZIDs (void)
 This function returns a vector containing the wave numbers in z-direction associated with the 3D homogenous expansion. Required if a parellelisation is applied in the Fourier direction. More...
 
LibUtilities::TranspositionSharedPtr GetTransposition (void)
 This function returns the transposition class associated with the homogeneous expansion. More...
 
NekDouble GetHomoLen (void)
 This function returns the Width of homogeneous direction associated with the homogeneous expansion. More...
 
void SetHomoLen (const NekDouble lhom)
 This function sets the Width of homogeneous direction associated with the homogeneous expansion. More...
 
Array< OneD, const unsigned int > GetYIDs (void)
 This function returns a vector containing the wave numbers in y-direction associated with the 3D homogenous expansion. Required if a parellelisation is applied in the Fourier direction. More...
 
void PhysInterp1DScaled (const NekDouble scale, const Array< OneD, NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 This function interpolates the physical space points in inarray to outarray using the same points defined in the expansion but where the number of points are rescaled by 1DScale. More...
 
void PhysGalerkinProjection1DScaled (const NekDouble scale, const Array< OneD, NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 This function Galerkin projects the physical space points in inarray to outarray where inarray is assumed to be defined in the expansion but where the number of points are rescaled by 1DScale. More...
 
int GetExpSize (void)
 This function returns the number of elements in the expansion. More...
 
int GetNumElmts (void)
 This function returns the number of elements in the expansion which may be different for a homogeoenous extended expansionp. More...
 
const std::shared_ptr< LocalRegions::ExpansionVectorGetExp () const
 This function returns the vector of elements in the expansion. More...
 
LocalRegions::ExpansionSharedPtrGetExp (int n) const
 This function returns (a shared pointer to) the local elemental expansion of the \(n^{\mathrm{th}}\) element. More...
 
LocalRegions::ExpansionSharedPtrGetExp (const Array< OneD, const NekDouble > &gloCoord)
 This function returns (a shared pointer to) the local elemental expansion containing the arbitrary point given by gloCoord. More...
 
int GetExpIndex (const Array< OneD, const NekDouble > &gloCoord, NekDouble tol=0.0, bool returnNearestElmt=false, int cachedId=-1, NekDouble maxDistance=1e6)
 This function returns the index of the local elemental expansion containing the arbitrary point given by gloCoord, within a distance tolerance of tol. More...
 
int GetExpIndex (const Array< OneD, const NekDouble > &gloCoords, Array< OneD, NekDouble > &locCoords, NekDouble tol=0.0, bool returnNearestElmt=false, int cachedId=-1, NekDouble maxDistance=1e6)
 
NekDouble PhysEvaluate (const Array< OneD, const NekDouble > &coords, const Array< OneD, const NekDouble > &phys)
 
int GetCoeff_Offset (int n) const
 Get the start offset position for a global list of m_coeffs correspoinding to element n. More...
 
int GetPhys_Offset (int n) const
 Get the start offset position for a global list of m_phys correspoinding to element n. More...
 
Array< OneD, NekDouble > & UpdateCoeffs ()
 This function returns (a reference to) the array \(\boldsymbol{\hat{u}}_l\) (implemented as m_coeffs) containing all local expansion coefficients. More...
 
Array< OneD, NekDouble > & UpdatePhys ()
 This function returns (a reference to) the array \(\boldsymbol{u}_l\) (implemented as m_phys) containing the function \(u^{\delta}(\boldsymbol{x})\) evaluated at the quadrature points. More...
 
void PhysDeriv (Direction edir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d)
 
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)
 This function discretely evaluates the derivative of a function \(f(\boldsymbol{x})\) on the domain consisting of all elements of the expansion. More...
 
void PhysDeriv (const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d)
 
void CurlCurl (Array< OneD, Array< OneD, NekDouble > > &Vel, Array< OneD, Array< OneD, NekDouble > > &Q)
 
void PhysDirectionalDeriv (const Array< OneD, const NekDouble > &direction, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
void GetMovingFrames (const SpatialDomains::GeomMMF MMFdir, const Array< OneD, const NekDouble > &CircCentre, Array< OneD, Array< OneD, NekDouble > > &outarray)
 
const Array< OneD, const std::shared_ptr< ExpList > > & GetBndCondExpansions ()
 
const Array< OneD, const NekDouble > & GetBndCondBwdWeight ()
 Get the weight value for boundary conditions. More...
 
void SetBndCondBwdWeight (const int index, const NekDouble value)
 Set the weight value for boundary conditions. More...
 
std::shared_ptr< ExpList > & UpdateBndCondExpansion (int i)
 
void Upwind (const Array< OneD, const NekDouble > &Vn, const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &Upwind)
 
void Upwind (const Array< OneD, const Array< OneD, NekDouble > > &Vec, const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &Upwind)
 
std::shared_ptr< ExpList > & GetTrace ()
 
std::shared_ptr< AssemblyMapDG > & GetTraceMap (void)
 
const Array< OneD, const int > & GetTraceBndMap (void)
 
void GetNormals (Array< OneD, Array< OneD, NekDouble > > &normals)
 
void GetElmtNormalLength (Array< OneD, NekDouble > &lengthsFwd, Array< OneD, NekDouble > &lengthsBwd)
 Get the length of elements in boundary normal direction. More...
 
void GetBwdWeight (Array< OneD, NekDouble > &weightAver, Array< OneD, NekDouble > &weightJump)
 Get the weight value for boundary conditions for boundary average and jump calculations. More...
 
void AddTraceIntegral (const Array< OneD, const NekDouble > &Fx, const Array< OneD, const NekDouble > &Fy, Array< OneD, NekDouble > &outarray)
 
void AddTraceIntegral (const Array< OneD, const NekDouble > &Fn, Array< OneD, NekDouble > &outarray)
 
void AddFwdBwdTraceIntegral (const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &outarray)
 
void GetFwdBwdTracePhys (Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd)
 
void GetFwdBwdTracePhys (const Array< OneD, const NekDouble > &field, Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd, bool FillBnd=true, bool PutFwdInBwdOnBCs=false, bool DoExchange=true)
 
void FillBwdWithBoundCond (const Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd, bool PutFwdInBwdOnBCs=false)
 
void AddTraceQuadPhysToField (const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &field)
 Add Fwd and Bwd value to field, a reverse procedure of GetFwdBwdTracePhys. More...
 
void AddTraceQuadPhysToOffDiag (const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &field)
 
void GetLocTraceFromTracePts (const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &locTraceFwd, Array< OneD, NekDouble > &locTraceBwd)
 
void FillBwdWithBwdWeight (Array< OneD, NekDouble > &weightave, Array< OneD, NekDouble > &weightjmp)
 Fill Bwd with boundary conditions. More...
 
void PeriodicBwdCopy (const Array< OneD, const NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd)
 Copy and fill the Periodic boundaries. More...
 
const std::vector< bool > & GetLeftAdjacentFaces (void) const
 
void ExtractTracePhys (Array< OneD, NekDouble > &outarray)
 
void ExtractTracePhys (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
const Array< OneD, const SpatialDomains::BoundaryConditionShPtr > & GetBndConditions ()
 
Array< OneD, SpatialDomains::BoundaryConditionShPtr > & UpdateBndConditions ()
 
void EvaluateBoundaryConditions (const NekDouble time=0.0, const std::string varName="", const NekDouble=NekConstants::kNekUnsetDouble, const NekDouble=NekConstants::kNekUnsetDouble)
 
void GeneralMatrixOp (const GlobalMatrixKey &gkey, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 This function calculates the result of the multiplication of a matrix of type specified by mkey with a vector given by inarray. More...
 
void GeneralMatrixOp_IterPerExp (const GlobalMatrixKey &gkey, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
void SetUpPhysNormals ()
 
void GetBoundaryToElmtMap (Array< OneD, int > &ElmtID, Array< OneD, int > &EdgeID)
 
void GetBndElmtExpansion (int i, std::shared_ptr< ExpList > &result, const bool DeclareCoeffPhysArrays=true)
 
void ExtractElmtToBndPhys (int i, const Array< OneD, NekDouble > &elmt, Array< OneD, NekDouble > &boundary)
 
void ExtractPhysToBndElmt (int i, const Array< OneD, const NekDouble > &phys, Array< OneD, NekDouble > &bndElmt)
 
void ExtractPhysToBnd (int i, const Array< OneD, const NekDouble > &phys, Array< OneD, NekDouble > &bnd)
 
void GetBoundaryNormals (int i, Array< OneD, Array< OneD, NekDouble > > &normals)
 
void GeneralGetFieldDefinitions (std::vector< LibUtilities::FieldDefinitionsSharedPtr > &fielddef, int NumHomoDir=0, Array< OneD, LibUtilities::BasisSharedPtr > &HomoBasis=LibUtilities::NullBasisSharedPtr1DArray, std::vector< NekDouble > &HomoLen=LibUtilities::NullNekDoubleVector, bool homoStrips=false, std::vector< unsigned int > &HomoSIDs=LibUtilities::NullUnsignedIntVector, std::vector< unsigned int > &HomoZIDs=LibUtilities::NullUnsignedIntVector, std::vector< unsigned int > &HomoYIDs=LibUtilities::NullUnsignedIntVector)
 
std::map< int, RobinBCInfoSharedPtrGetRobinBCInfo ()
 
void GetPeriodicEntities (PeriodicMap &periodicVerts, PeriodicMap &periodicEdges, PeriodicMap &periodicFaces=NullPeriodicMap)
 
std::vector< LibUtilities::FieldDefinitionsSharedPtrGetFieldDefinitions ()
 
void GetFieldDefinitions (std::vector< LibUtilities::FieldDefinitionsSharedPtr > &fielddef)
 
void AppendFieldData (LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector< NekDouble > &fielddata)
 Append the element data listed in elements fielddef->m_ElementIDs onto fielddata. More...
 
void AppendFieldData (LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector< NekDouble > &fielddata, Array< OneD, NekDouble > &coeffs)
 Append the data in coeffs listed in elements fielddef->m_ElementIDs onto fielddata. More...
 
void ExtractElmtDataToCoeffs (LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector< NekDouble > &fielddata, std::string &field, Array< OneD, NekDouble > &coeffs)
 Extract the data in fielddata into the coeffs using the basic ExpList Elemental expansions rather than planes in homogeneous case. More...
 
void ExtractCoeffsToCoeffs (const std::shared_ptr< ExpList > &fromExpList, const Array< OneD, const NekDouble > &fromCoeffs, Array< OneD, NekDouble > &toCoeffs)
 Extract the data from fromField using fromExpList the coeffs using the basic ExpList Elemental expansions rather than planes in homogeneous case. More...
 
void ExtractDataToCoeffs (LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector< NekDouble > &fielddata, std::string &field, Array< OneD, NekDouble > &coeffs)
 Extract the data in fielddata into the coeffs. More...
 
void GenerateElementVector (const int ElementID, const NekDouble scalar1, const NekDouble scalar2, Array< OneD, NekDouble > &outarray)
 Generate vector v such that v[i] = scalar1 if i is in the element < ElementID. Otherwise, v[i] = scalar2. More...
 
std::shared_ptr< ExpListGetSharedThisPtr ()
 Returns a shared pointer to the current object. More...
 
std::shared_ptr< LibUtilities::SessionReaderGetSession () const
 Returns the session object. More...
 
std::shared_ptr< LibUtilities::CommGetComm ()
 Returns the comm object. More...
 
SpatialDomains::MeshGraphSharedPtr GetGraph ()
 
LibUtilities::BasisSharedPtr GetHomogeneousBasis (void)
 
std::shared_ptr< ExpList > & GetPlane (int n)
 
void CreateCollections (Collections::ImplementationType ImpType=Collections::eNoImpType)
 Construct collections of elements containing a single element type and polynomial order from the list of expansions. More...
 
void ClearGlobalLinSysManager (void)
 
const Array< OneD, const std::pair< int, int > > & GetCoeffsToElmt () const
 Get m_coeffs to elemental value map. More...
 
void AddTraceJacToElmtJac (const Array< OneD, const DNekMatSharedPtr > &FwdMat, const Array< OneD, const DNekMatSharedPtr > &BwdMat, Array< OneD, DNekMatSharedPtr > &fieldMat)
 inverse process of v_GetFwdBwdTracePhys. Given Trace integration of Fwd and Bwd Jacobian, with dimension NtotalTrace*TraceCoef*TracePhys. return Elemental Jacobian matrix with dimension NtotalElement*ElementCoef*ElementPhys. More...
 
void GetMatIpwrtDeriveBase (const Array< OneD, const Array< OneD, NekDouble > > &inarray, const int nDirctn, Array< OneD, DNekMatSharedPtr > &mtxPerVar)
 
void GetMatIpwrtDeriveBase (const TensorOfArray3D< NekDouble > &inarray, Array< OneD, DNekMatSharedPtr > &mtxPerVar)
 
void GetDiagMatIpwrtBase (const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, DNekMatSharedPtr > &mtxPerVar)
 
void AddTraceIntegralToOffDiag (const Array< OneD, const NekDouble > &FwdFlux, const Array< OneD, const NekDouble > &BwdFlux, Array< OneD, NekDouble > &outarray)
 
void AddRightIPTPhysDerivBase (const int dir, const Array< OneD, const DNekMatSharedPtr > ElmtJacQuad, Array< OneD, DNekMatSharedPtr > ElmtJacCoef)
 
void AddRightIPTBaseMatrix (const Array< OneD, const DNekMatSharedPtr > ElmtJacQuad, Array< OneD, DNekMatSharedPtr > ElmtJacCoef)
 
const LocTraceToTraceMapSharedPtrGetLocTraceToTraceMap () const
 

Public Attributes

Array< OneD, int > m_BCtoElmMap
 
Array< OneD, int > m_BCtoTraceMap
 

Protected Member Functions

void GenerateBoundaryConditionExpansion (const SpatialDomains::MeshGraphSharedPtr &graph1D, const SpatialDomains::BoundaryConditions &bcs, const std::string variable, const bool DeclareCoeffPhysArrays=true)
 Discretises the boundary conditions. More...
 
void FindPeriodicTraces (const SpatialDomains::BoundaryConditions &bcs, const std::string variable)
 Generate a associative map of periodic vertices in a mesh. More...
 
virtual ExpListSharedPtrv_GetTrace ()
 
virtual AssemblyMapDGSharedPtrv_GetTraceMap (void)
 
virtual const LocTraceToTraceMapSharedPtrv_GetLocTraceToTraceMap (void) const
 
virtual void v_AddTraceIntegral (const Array< OneD, const NekDouble > &Fn, Array< OneD, NekDouble > &outarray)
 Add trace contributions into elemental coefficient spaces. More...
 
virtual void v_AddFwdBwdTraceIntegral (const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &outarray)
 Add trace contributions into elemental coefficient spaces. More...
 
virtual void v_AddTraceQuadPhysToField (const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &field)
 
virtual void v_ExtractTracePhys (Array< OneD, NekDouble > &outarray)
 
virtual void v_ExtractTracePhys (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 This method extracts the trace (verts in 1D) from the field inarray and puts the values in outarray. More...
 
virtual void v_GetLocTraceFromTracePts (const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &locTraceFwd, Array< OneD, NekDouble > &locTraceBwd)
 
void GenerateFieldBnd1D (SpatialDomains::BoundaryConditions &bcs, const std::string variable)
 
virtual std::map< int, RobinBCInfoSharedPtrv_GetRobinBCInfo ()
 
virtual const Array< OneD, const MultiRegions::ExpListSharedPtr > & v_GetBndCondExpansions ()
 
virtual const Array< OneD, const SpatialDomains::BoundaryConditionShPtr > & v_GetBndConditions ()
 
virtual MultiRegions::ExpListSharedPtrv_UpdateBndCondExpansion (int i)
 
virtual Array< OneD, SpatialDomains::BoundaryConditionShPtr > & v_UpdateBndConditions ()
 
virtual void v_GetBoundaryToElmtMap (Array< OneD, int > &ElmtID, Array< OneD, int > &TraceID)
 
virtual void v_GetBndElmtExpansion (int i, std::shared_ptr< ExpList > &result, const bool DeclareCoeffPhysArrays)
 
virtual void v_Reset ()
 Reset this field, so that geometry information can be updated. More...
 
virtual void v_EvaluateBoundaryConditions (const NekDouble time=0.0, const std::string varName="", const NekDouble x2_in=NekConstants::kNekUnsetDouble, const NekDouble x3_in=NekConstants::kNekUnsetDouble)
 Evaluate all boundary conditions at a given time.. More...
 
virtual void v_HelmSolve (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::ConstFactorMap &factors, const StdRegions::VarCoeffMap &varcoeff, const MultiRegions::VarFactorsMap &varfactors, const Array< OneD, const NekDouble > &dirForcing, const bool PhysSpaceForcing)
 Solve the Helmholtz equation. More...
 
virtual void v_PeriodicBwdCopy (const Array< OneD, const NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd)
 
virtual void v_FillBwdWithBwdWeight (Array< OneD, NekDouble > &weightave, Array< OneD, NekDouble > &weightjmp)
 Fill the weight with m_bndCondBndWeight. More...
 
virtual void v_GetFwdBwdTracePhys (Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd)
 
virtual void v_GetFwdBwdTracePhys (const Array< OneD, const NekDouble > &field, Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd, bool FillBnd=true, bool PutFwdInBwdOnBCs=false, bool DoExchange=true)
 This method extracts the "forward" and "backward" trace data from the array field and puts the data into output vectors Fwd and Bwd. More...
 
virtual void v_FillBwdWithBoundCond (const Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd, bool PutFwdInBwdOnBCs)
 
virtual const Array< OneD, const NekDouble > & v_GetBndCondBwdWeight ()
 
virtual void v_SetBndCondBwdWeight (const int index, const NekDouble value)
 
void SetUpDG (const std::string="DefaultVar")
 Set up all DG member variables and maps. More...
 
bool IsLeftAdjacentTrace (const int n, const int e)
 
virtual void v_GetPeriodicEntities (PeriodicMap &periodicVerts, PeriodicMap &periodicEdges, PeriodicMap &periodicFaces)
 Obtain a copy of the periodic edges and vertices for this field. More...
 
- Protected Member Functions inherited from Nektar::MultiRegions::ExpList
std::shared_ptr< DNekMatGenGlobalMatrixFull (const GlobalLinSysKey &mkey, const std::shared_ptr< AssemblyMapCG > &locToGloMap)
 
const DNekScalBlkMatSharedPtr GenBlockMatrix (const GlobalMatrixKey &gkey)
 This function assembles the block diagonal matrix of local matrices of the type mtype. More...
 
const DNekScalBlkMatSharedPtrGetBlockMatrix (const GlobalMatrixKey &gkey)
 
void MultiplyByBlockMatrix (const GlobalMatrixKey &gkey, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
std::shared_ptr< GlobalMatrixGenGlobalMatrix (const GlobalMatrixKey &mkey, const std::shared_ptr< AssemblyMapCG > &locToGloMap)
 Generates a global matrix from the given key and map. More...
 
void GlobalEigenSystem (const std::shared_ptr< DNekMat > &Gmat, Array< OneD, NekDouble > &EigValsReal, Array< OneD, NekDouble > &EigValsImag, Array< OneD, NekDouble > &EigVecs=NullNekDouble1DArray)
 
std::shared_ptr< GlobalLinSysGenGlobalLinSys (const GlobalLinSysKey &mkey, const std::shared_ptr< AssemblyMapCG > &locToGloMap)
 This operation constructs the global linear system of type mkey. More...
 
std::shared_ptr< GlobalLinSysGenGlobalBndLinSys (const GlobalLinSysKey &mkey, const AssemblyMapSharedPtr &locToGloMap)
 Generate a GlobalLinSys from information provided by the key "mkey" and the mapping provided in LocToGloBaseMap. More...
 
virtual int v_GetNumElmts (void)
 
virtual void v_Upwind (const Array< OneD, const Array< OneD, NekDouble > > &Vec, const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &Upwind)
 
virtual void v_Upwind (const Array< OneD, const NekDouble > &Vn, const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &Upwind)
 
virtual const Array< OneD, const int > & v_GetTraceBndMap ()
 
virtual void v_GetNormals (Array< OneD, Array< OneD, NekDouble > > &normals)
 Populate normals with the normals of all expansions. More...
 
virtual void v_AddTraceIntegral (const Array< OneD, const NekDouble > &Fx, const Array< OneD, const NekDouble > &Fy, Array< OneD, NekDouble > &outarray)
 
virtual void v_AddTraceQuadPhysToOffDiag (const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &field)
 
virtual const std::vector< bool > & v_GetLeftAdjacentFaces (void) const
 
virtual void v_MultiplyByInvMassMatrix (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
virtual void v_LinearAdvectionDiffusionReactionSolve (const Array< OneD, Array< OneD, NekDouble > > &velocity, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const NekDouble lambda, const Array< OneD, const NekDouble > &dirForcing=NullNekDouble1DArray)
 
virtual void v_LinearAdvectionReactionSolve (const Array< OneD, Array< OneD, NekDouble > > &velocity, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const NekDouble lambda, const Array< OneD, const NekDouble > &dirForcing=NullNekDouble1DArray)
 
virtual void v_ImposeDirichletConditions (Array< OneD, NekDouble > &outarray)
 
virtual void v_FillBndCondFromField ()
 
virtual void v_FillBndCondFromField (const int nreg)
 
virtual void v_LocalToGlobal (bool UseComm)
 
virtual void v_LocalToGlobal (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool UseComm)
 
virtual void v_GlobalToLocal (void)
 
virtual void v_GlobalToLocal (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
virtual void v_BwdTrans (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
virtual void v_BwdTrans_IterPerExp (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
virtual void v_FwdTrans (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
virtual void v_FwdTrans_IterPerExp (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
virtual void v_FwdTrans_BndConstrained (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
virtual void v_SmoothField (Array< OneD, NekDouble > &field)
 
virtual void v_IProductWRTBase (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
virtual void v_IProductWRTBase_IterPerExp (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
virtual void v_GeneralMatrixOp (const GlobalMatrixKey &gkey, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
virtual void v_GetCoords (Array< OneD, NekDouble > &coord_0, Array< OneD, NekDouble > &coord_1, Array< OneD, NekDouble > &coord_2=NullNekDouble1DArray)
 
virtual void v_PhysDeriv (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2)
 
virtual void v_PhysDeriv (const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d)
 
virtual void v_PhysDeriv (Direction edir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d)
 
virtual void v_CurlCurl (Array< OneD, Array< OneD, NekDouble > > &Vel, Array< OneD, Array< OneD, NekDouble > > &Q)
 
virtual void v_PhysDirectionalDeriv (const Array< OneD, const NekDouble > &direction, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
virtual void v_GetMovingFrames (const SpatialDomains::GeomMMF MMFdir, const Array< OneD, const NekDouble > &CircCentre, Array< OneD, Array< OneD, NekDouble > > &outarray)
 
virtual void v_HomogeneousFwdTrans (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool Shuff=true, bool UnShuff=true)
 
virtual void v_HomogeneousBwdTrans (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool Shuff=true, bool UnShuff=true)
 
virtual void v_DealiasedProd (const Array< OneD, NekDouble > &inarray1, const Array< OneD, NekDouble > &inarray2, Array< OneD, NekDouble > &outarray)
 
virtual void v_DealiasedDotProd (const Array< OneD, Array< OneD, NekDouble > > &inarray1, const Array< OneD, Array< OneD, NekDouble > > &inarray2, Array< OneD, Array< OneD, NekDouble > > &outarray)
 
virtual void v_GetBCValues (Array< OneD, NekDouble > &BndVals, const Array< OneD, NekDouble > &TotField, int BndID)
 
virtual void v_NormVectorIProductWRTBase (Array< OneD, const NekDouble > &V1, Array< OneD, const NekDouble > &V2, Array< OneD, NekDouble > &outarray, int BndID)
 
virtual void v_NormVectorIProductWRTBase (Array< OneD, Array< OneD, NekDouble > > &V, Array< OneD, NekDouble > &outarray)
 
virtual void v_SetUpPhysNormals ()
 
virtual void v_ExtractElmtToBndPhys (const int i, const Array< OneD, NekDouble > &elmt, Array< OneD, NekDouble > &boundary)
 
virtual void v_ExtractPhysToBndElmt (const int i, const Array< OneD, const NekDouble > &phys, Array< OneD, NekDouble > &bndElmt)
 
virtual void v_ExtractPhysToBnd (const int i, const Array< OneD, const NekDouble > &phys, Array< OneD, NekDouble > &bnd)
 
virtual void v_GetBoundaryNormals (int i, Array< OneD, Array< OneD, NekDouble > > &normals)
 
virtual std::vector< LibUtilities::FieldDefinitionsSharedPtrv_GetFieldDefinitions (void)
 
virtual void v_GetFieldDefinitions (std::vector< LibUtilities::FieldDefinitionsSharedPtr > &fielddef)
 
virtual void v_AppendFieldData (LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector< NekDouble > &fielddata)
 
virtual void v_AppendFieldData (LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector< NekDouble > &fielddata, Array< OneD, NekDouble > &coeffs)
 
virtual void v_ExtractDataToCoeffs (LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector< NekDouble > &fielddata, std::string &field, Array< OneD, NekDouble > &coeffs)
 Extract data from raw field data into expansion list. More...
 
virtual void v_ExtractCoeffsToCoeffs (const std::shared_ptr< ExpList > &fromExpList, const Array< OneD, const NekDouble > &fromCoeffs, Array< OneD, NekDouble > &toCoeffs)
 
virtual void v_WriteTecplotHeader (std::ostream &outfile, std::string var="")
 
virtual void v_WriteTecplotZone (std::ostream &outfile, int expansion)
 
virtual void v_WriteTecplotField (std::ostream &outfile, int expansion)
 
virtual void v_WriteTecplotConnectivity (std::ostream &outfile, int expansion)
 
virtual void v_WriteVtkPieceHeader (std::ostream &outfile, int expansion, int istrip)
 
virtual void v_WriteVtkPieceData (std::ostream &outfile, int expansion, std::string var)
 
virtual NekDouble v_L2 (const Array< OneD, const NekDouble > &phys, const Array< OneD, const NekDouble > &soln=NullNekDouble1DArray)
 
virtual NekDouble v_Integral (const Array< OneD, const NekDouble > &inarray)
 
virtual NekDouble v_VectorFlux (const Array< OneD, Array< OneD, NekDouble > > &inarray)
 
virtual Array< OneD, const NekDoublev_HomogeneousEnergy (void)
 
virtual LibUtilities::TranspositionSharedPtr v_GetTransposition (void)
 
virtual NekDouble v_GetHomoLen (void)
 
virtual void v_SetHomoLen (const NekDouble lhom)
 
virtual Array< OneD, const unsigned int > v_GetZIDs (void)
 
virtual Array< OneD, const unsigned int > v_GetYIDs (void)
 
virtual void v_PhysInterp1DScaled (const NekDouble scale, const Array< OneD, NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
virtual void v_PhysGalerkinProjection1DScaled (const NekDouble scale, const Array< OneD, NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
virtual void v_ClearGlobalLinSysManager (void)
 
void ExtractFileBCs (const std::string &fileName, LibUtilities::CommSharedPtr comm, const std::string &varName, const std::shared_ptr< ExpList > locExpList)
 

Protected Attributes

int m_numDirBndCondExpansions
 The number of boundary segments on which Dirichlet boundary conditions are imposed. More...
 
Array< OneD, MultiRegions::ExpListSharedPtrm_bndCondExpansions
 An object which contains the discretised boundary conditions. More...
 
Array< OneD, NekDoublem_bndCondBndWeight
 
Array< OneD, SpatialDomains::BoundaryConditionShPtrm_bndConditions
 An array which contains the information about the boundary condition on the different boundary regions. More...
 
GlobalLinSysMapShPtr m_globalBndMat
 Global boundary matrix. More...
 
ExpListSharedPtr m_trace
 Trace space storage for points between elements. More...
 
AssemblyMapDGSharedPtr m_traceMap
 Local to global DG mapping for trace space. More...
 
std::set< int > m_boundaryTraces
 A set storing the global IDs of any boundary Verts. More...
 
PeriodicMap m_periodicVerts
 A map which identifies groups of periodic vertices. More...
 
PeriodicMap m_periodicEdges
 A map which identifies pairs of periodic edges. More...
 
PeriodicMap m_periodicFaces
 A map which identifies pairs of periodic faces. More...
 
std::vector< int > m_periodicFwdCopy
 A vector indicating degress of freedom which need to be copied from forwards to backwards space in case of a periodic boundary condition. More...
 
std::vector< int > m_periodicBwdCopy
 
std::vector< bool > m_leftAdjacentTraces
 
LocTraceToTraceMapSharedPtr m_locTraceToTraceMap
 
- Protected Attributes inherited from Nektar::MultiRegions::ExpList
ExpansionType m_expType
 Exapnsion type. More...
 
LibUtilities::CommSharedPtr m_comm
 Communicator. More...
 
LibUtilities::SessionReaderSharedPtr m_session
 Session. More...
 
SpatialDomains::MeshGraphSharedPtr m_graph
 Mesh associated with this expansion list. More...
 
int m_ncoeffs
 The total number of local degrees of freedom. m_ncoeffs \(=N_{\mathrm{eof}}=\sum_{e=1}^{{N_{\mathrm{el}}}}N^{e}_l\). More...
 
int m_npoints
 
Array< OneD, NekDoublem_coeffs
 Concatenation of all local expansion coefficients. More...
 
Array< OneD, NekDoublem_phys
 The global expansion evaluated at the quadrature points. More...
 
bool m_physState
 The state of the array m_phys. More...
 
std::shared_ptr< LocalRegions::ExpansionVectorm_exp
 The list of local expansions. More...
 
Collections::CollectionVector m_collections
 
std::vector< bool > m_collectionsDoInit
 Vector of bools to act as an initialise on first call flag. More...
 
std::vector< int > m_coll_coeff_offset
 Offset of elemental data into the array m_coeffs. More...
 
std::vector< int > m_coll_phys_offset
 Offset of elemental data into the array m_phys. More...
 
Array< OneD, int > m_coeff_offset
 Offset of elemental data into the array m_coeffs. More...
 
Array< OneD, int > m_phys_offset
 Offset of elemental data into the array m_phys. More...
 
Array< OneD, std::pair< int, int > > m_coeffsToElmt
 m_coeffs to elemental value map More...
 
BlockMatrixMapShPtr m_blockMat
 
bool m_WaveSpace
 
std::unordered_map< int, int > m_elmtToExpId
 Mapping from geometry ID of element to index inside m_exp. More...
 

Private Member Functions

SpatialDomains::BoundaryConditionsSharedPtr GetDomainBCs (const SpatialDomains::CompositeMap &domain, const SpatialDomains::BoundaryConditions &Allbcs, const std::string &variable)
 
virtual void v_AddTraceIntegralToOffDiag (const Array< OneD, const NekDouble > &FwdFlux, const Array< OneD, const NekDouble > &BwdFlux, Array< OneD, NekDouble > &outarray)
 

Private Attributes

std::vector< bool > m_negatedFluxNormal
 

Additional Inherited Members

- Static Protected Member Functions inherited from Nektar::MultiRegions::ExpList
static SpatialDomains::BoundaryConditionShPtr GetBoundaryCondition (const SpatialDomains::BoundaryConditionCollection &collection, unsigned int index, const std::string &variable)
 

Detailed Description

This class is the abstractio n of a global discontinuous two- dimensional spectral/hp element expansion which approximates the solution of a set of partial differential equations.

This class augments the list of local expansions inherited from ExpList with boundary conditions. Inter-element boundaries are handled using an discontinuous Galerkin scheme.

Definition at line 55 of file DisContField.h.

Constructor & Destructor Documentation

◆ DisContField() [1/6]

Nektar::MultiRegions::DisContField::DisContField ( )

Default constructor.

Constructs an empty expansion list with no boundary conditions.

Definition at line 65 of file DisContField.cpp.

65  :
66  ExpList(),
70  {
71  }
Array< OneD, SpatialDomains::BoundaryConditionShPtr > m_bndConditions
An array which contains the information about the boundary condition on the different boundary region...
Definition: DisContField.h:148
Array< OneD, MultiRegions::ExpListSharedPtr > m_bndCondExpansions
An object which contains the discretised boundary conditions.
Definition: DisContField.h:142
ExpListSharedPtr m_trace
Trace space storage for points between elements.
Definition: DisContField.h:154
ExpList(const ExpansionType Type=eNoType)
The default constructor using a type.
Definition: ExpList.cpp:104
static ExpListSharedPtr NullExpListSharedPtr
Definition: ExpList.h:1792

◆ DisContField() [2/6]

Nektar::MultiRegions::DisContField::DisContField ( const LibUtilities::SessionReaderSharedPtr pSession,
const SpatialDomains::MeshGraphSharedPtr graph,
const std::string &  variable,
const bool  SetUpJustDG = true,
const bool  DeclareCoeffPhysArrays = true,
const Collections::ImplementationType  ImpType = Collections::eNoImpType 
)

Constructs a 1D discontinuous field based on a mesh and boundary conditions.

An expansion list for the boundary expansions is generated first for the field. These are subsequently evaluated for time zero. The trace map is then constructed.

Parameters
graph1DA mesh, containing information about the domain and the spectral/hp element expansions.
bcsInformation about the enforced boundary conditions.
variableThe session variable associated with the boundary conditions to enforce.
solnTypeType of global system to use.

Definition at line 86 of file DisContField.cpp.

93  : ExpList(pSession,graph,DeclareCoeffPhysArrays, variable, ImpType),
95  {
96  if (variable.compare("DefaultVar") != 0)
97  {
98  SpatialDomains::BoundaryConditions bcs(m_session, graph);
99 
100  GenerateBoundaryConditionExpansion(graph,bcs,variable,
101  DeclareCoeffPhysArrays);
102  if(DeclareCoeffPhysArrays)
103  {
104  EvaluateBoundaryConditions(0.0, variable);
105  }
106  ApplyGeomInfo();
107  FindPeriodicTraces(bcs,variable);
108  }
109 
110  if(SetUpJustDG)
111  {
112  SetUpDG(variable);
113  }
114  else
115  {
116  int i,cnt;
117  Array<OneD, int> ElmtID, TraceID;
118  GetBoundaryToElmtMap(ElmtID, TraceID);
119 
120  for(cnt = i = 0; i < m_bndCondExpansions.size(); ++i)
121  {
123  locExpList = m_bndCondExpansions[i];
124 
125  for(int e = 0; e < locExpList->GetExpSize(); ++e)
126  {
127  (*m_exp)[ElmtID[cnt+e]]->SetTraceExp
128  (TraceID[cnt+e], locExpList->GetExp(e));
129  locExpList->GetExp(e)->SetAdjacentElementExp
130  (TraceID[cnt+e], (*m_exp)[ElmtID[cnt+e]]);
131  }
132  cnt += m_bndCondExpansions[i]->GetExpSize();
133  }
134 
135  if(m_session->DefinesSolverInfo("PROJECTION"))
136  {
137  std::string ProjectStr =
138  m_session->GetSolverInfo("PROJECTION");
139  if((ProjectStr == "MixedCGDG") ||
140  (ProjectStr == "Mixed_CG_Discontinuous"))
141  {
142  SetUpDG();
143  }
144  else
145  {
147  }
148  }
149  else
150  {
152  }
153  }
154  }
void FindPeriodicTraces(const SpatialDomains::BoundaryConditions &bcs, const std::string variable)
Generate a associative map of periodic vertices in a mesh.
void SetUpDG(const std::string="DefaultVar")
Set up all DG member variables and maps.
void GenerateBoundaryConditionExpansion(const SpatialDomains::MeshGraphSharedPtr &graph1D, const SpatialDomains::BoundaryConditions &bcs, const std::string variable, const bool DeclareCoeffPhysArrays=true)
Discretises the boundary conditions.
void GetBoundaryToElmtMap(Array< OneD, int > &ElmtID, Array< OneD, int > &EdgeID)
Definition: ExpList.h:2709
std::shared_ptr< LocalRegions::ExpansionVector > m_exp
The list of local expansions.
Definition: ExpList.h:1290
void EvaluateBoundaryConditions(const NekDouble time=0.0, const std::string varName="", const NekDouble=NekConstants::kNekUnsetDouble, const NekDouble=NekConstants::kNekUnsetDouble)
Definition: ExpList.h:2657
void ApplyGeomInfo()
Apply geometry information to each expansion.
Definition: ExpList.cpp:2567
LibUtilities::SessionReaderSharedPtr m_session
Session.
Definition: ExpList.h:1223
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.

References Nektar::MultiRegions::ExpList::ApplyGeomInfo(), Nektar::MultiRegions::ExpList::EvaluateBoundaryConditions(), FindPeriodicTraces(), GenerateBoundaryConditionExpansion(), Nektar::MultiRegions::ExpList::GetBoundaryToElmtMap(), m_bndCondExpansions, Nektar::MultiRegions::ExpList::m_exp, Nektar::MultiRegions::ExpList::m_session, SetUpDG(), and Nektar::MultiRegions::ExpList::SetUpPhysNormals().

◆ DisContField() [3/6]

Nektar::MultiRegions::DisContField::DisContField ( const LibUtilities::SessionReaderSharedPtr pSession,
const SpatialDomains::MeshGraphSharedPtr graph1D,
const SpatialDomains::CompositeMap domain,
const SpatialDomains::BoundaryConditions Allbcs,
const std::string &  variable,
bool  SetToOneSpaceDimension = false,
const Collections::ImplementationType  ImpType = Collections::eNoImpType 
)

Constructor for a DisContField from a List of subdomains New Constructor for arterial network.

Constructor for use in multidomain computations where a domain list can be passed instead of graph1D

Parameters
domainSubdomain specified in the inputfile from which the DisContField is set up

Definition at line 584 of file DisContField.cpp.

591  :
592  ExpList(pSession,domain,graph1D,true,variable,
593  SetToOneSpaceDimension, pSession->GetComm(),ImpType)
594  {
595  if (variable.compare("DefaultVar") != 0)
596  {
598  DomBCs = GetDomainBCs(domain,Allbcs,variable);
599 
601  EvaluateBoundaryConditions(0.0, variable);
602  ApplyGeomInfo();
603  FindPeriodicTraces(*DomBCs,variable);
604  }
605 
606  SetUpDG(variable);
607  }
SpatialDomains::BoundaryConditionsSharedPtr GetDomainBCs(const SpatialDomains::CompositeMap &domain, const SpatialDomains::BoundaryConditions &Allbcs, const std::string &variable)
SpatialDomains::MeshGraphSharedPtr m_graph
Mesh associated with this expansion list.
Definition: ExpList.h:1226
std::shared_ptr< BoundaryConditions > BoundaryConditionsSharedPtr
Definition: Conditions.h:289

References Nektar::MultiRegions::ExpList::ApplyGeomInfo(), Nektar::MultiRegions::ExpList::EvaluateBoundaryConditions(), FindPeriodicTraces(), GenerateBoundaryConditionExpansion(), GetDomainBCs(), Nektar::MultiRegions::ExpList::m_graph, and SetUpDG().

◆ DisContField() [4/6]

Nektar::MultiRegions::DisContField::DisContField ( const DisContField In,
const bool  DeclareCoeffPhysArrays = true 
)

Constructs a 1D discontinuous field based on an existing field.

Constructs a field as a copy of an existing field.

Parameters
InExisting DisContField object to copy.

Definition at line 613 of file DisContField.cpp.

614  :
615  ExpList(In,DeclareCoeffPhysArrays),
616  m_bndCondExpansions(In.m_bndCondExpansions),
617  m_bndConditions(In.m_bndConditions),
618  m_globalBndMat(In.m_globalBndMat),
619  m_traceMap(In.m_traceMap),
620  m_boundaryTraces(In.m_boundaryTraces),
621  m_periodicVerts(In.m_periodicVerts),
622  m_periodicFwdCopy(In.m_periodicFwdCopy),
623  m_periodicBwdCopy(In.m_periodicBwdCopy),
624  m_leftAdjacentTraces(In.m_leftAdjacentTraces),
625  m_locTraceToTraceMap (In.m_locTraceToTraceMap)
626  {
627  if (In.m_trace)
628  {
630  (*In.m_trace,DeclareCoeffPhysArrays);
631  }
632  }
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
std::vector< int > m_periodicBwdCopy
Definition: DisContField.h:185
std::set< int > m_boundaryTraces
A set storing the global IDs of any boundary Verts.
Definition: DisContField.h:162
LocTraceToTraceMapSharedPtr m_locTraceToTraceMap
Definition: DisContField.h:197
AssemblyMapDGSharedPtr m_traceMap
Local to global DG mapping for trace space.
Definition: DisContField.h:157
GlobalLinSysMapShPtr m_globalBndMat
Global boundary matrix.
Definition: DisContField.h:151
PeriodicMap m_periodicVerts
A map which identifies groups of periodic vertices.
Definition: DisContField.h:167
std::vector< int > m_periodicFwdCopy
A vector indicating degress of freedom which need to be copied from forwards to backwards space in ca...
Definition: DisContField.h:184
std::vector< bool > m_leftAdjacentTraces
Definition: DisContField.h:191

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), and m_trace.

◆ DisContField() [5/6]

Nektar::MultiRegions::DisContField::DisContField ( const DisContField In,
const SpatialDomains::MeshGraphSharedPtr graph,
const std::string &  variable,
const bool  SetUpJustDG = false,
const bool  DeclareCoeffPhysArrays = true 
)

Definition at line 639 of file DisContField.cpp.

644  :
645  ExpList(In,DeclareCoeffPhysArrays)
646  {
647 
649 
650  // Set up boundary conditions for this variable.
651  // Do not set up BCs if default variable
652  if(variable.compare("DefaultVar") != 0)
653  {
654  SpatialDomains::BoundaryConditions bcs(m_session, graph);
655  GenerateBoundaryConditionExpansion(graph, bcs, variable);
656 
657  if (DeclareCoeffPhysArrays)
658  {
659  EvaluateBoundaryConditions(0.0, variable);
660  }
661 
663  {
664  // Find periodic edges for this variable.
665  FindPeriodicTraces(bcs, variable);
666 
667  if(SetUpJustDG)
668  {
669  SetUpDG(variable);
670  }
671  else
672  {
673  // set elmt edges to point to robin bc edges if required
674  int i, cnt = 0;
675  Array<OneD, int> ElmtID,TraceID;
676  GetBoundaryToElmtMap(ElmtID,TraceID);
677 
678  for (i = 0; i < m_bndCondExpansions.size(); ++i)
679  {
681 
682  int e;
683  locExpList = m_bndCondExpansions[i];
684 
685  for(e = 0; e < locExpList->GetExpSize(); ++e)
686  {
687  (*m_exp)[ElmtID[cnt+e]]->SetTraceExp
688  (TraceID[cnt+e], locExpList->GetExp(e));
689  locExpList->GetExp(e)->SetAdjacentElementExp
690  (TraceID[cnt+e], (*m_exp)[ElmtID[cnt+e]]);
691  }
692 
693  cnt += m_bndCondExpansions[i]->GetExpSize();
694  }
695 
696 
697  if (m_session->DefinesSolverInfo("PROJECTION"))
698  {
699  std::string ProjectStr =
700  m_session->GetSolverInfo("PROJECTION");
701 
702  if ((ProjectStr == "MixedCGDG") ||
703  (ProjectStr == "Mixed_CG_Discontinuous"))
704  {
705  SetUpDG();
706  }
707  else
708  {
710  }
711  }
712  else
713  {
715  }
716  }
717  }
718  else
719  {
720  m_globalBndMat = In.m_globalBndMat;
721  m_trace = In.m_trace;
722  m_traceMap = In.m_traceMap;
723  m_locTraceToTraceMap = In.m_locTraceToTraceMap;
724  m_periodicVerts = In.m_periodicVerts;
725  m_periodicEdges = In.m_periodicEdges;
726  m_periodicFaces = In.m_periodicFaces;
727  m_periodicFwdCopy = In.m_periodicFwdCopy;
728  m_periodicBwdCopy = In.m_periodicBwdCopy;
729  m_boundaryTraces = In.m_boundaryTraces;
730  m_leftAdjacentTraces = In.m_leftAdjacentTraces;
731 
732  if (SetUpJustDG == false)
733  {
734  // set elmt edges to point to robin bc edges if required
735  int i, cnt = 0;
736  Array<OneD, int> ElmtID, TraceID;
737  GetBoundaryToElmtMap(ElmtID, TraceID);
738 
739  for (i = 0; i < m_bndCondExpansions.size(); ++i)
740  {
742 
743  int e;
744  locExpList = m_bndCondExpansions[i];
745 
746  for (e = 0; e < locExpList->GetExpSize(); ++e)
747  {
748  (*m_exp)[ElmtID[cnt+e]]->SetTraceExp
749  (TraceID[cnt+e], locExpList->GetExp(e));
750  locExpList->GetExp(e)->SetAdjacentElementExp
751  (TraceID[cnt+e], (*m_exp)[ElmtID[cnt+e]]);
752  }
753  cnt += m_bndCondExpansions[i]->GetExpSize();
754  }
755 
757  }
758  }
759  }
760  }
PeriodicMap m_periodicEdges
A map which identifies pairs of periodic edges.
Definition: DisContField.h:172
PeriodicMap m_periodicFaces
A map which identifies pairs of periodic faces.
Definition: DisContField.h:177
bool SameTypeOfBoundaryConditions(const DisContField &In)
Check to see if expansion has the same BCs as In.

References Nektar::MultiRegions::ExpList::EvaluateBoundaryConditions(), FindPeriodicTraces(), GenerateBoundaryConditionExpansion(), Nektar::MultiRegions::ExpList::GetBoundaryToElmtMap(), m_bndCondExpansions, m_boundaryTraces, Nektar::MultiRegions::ExpList::m_exp, m_globalBndMat, m_leftAdjacentTraces, m_locTraceToTraceMap, m_periodicBwdCopy, m_periodicEdges, m_periodicFaces, m_periodicFwdCopy, m_periodicVerts, Nektar::MultiRegions::ExpList::m_session, m_trace, m_traceMap, Nektar::MultiRegions::NullExpListSharedPtr, SameTypeOfBoundaryConditions(), SetUpDG(), and Nektar::MultiRegions::ExpList::SetUpPhysNormals().

◆ DisContField() [6/6]

Nektar::MultiRegions::DisContField::DisContField ( const ExpList In)

Constructs a 1D discontinuous field based on an existing field. (needed in order to use ContField( const ExpList &In) constructor.

Constructs a field as a copy of an existing explist field.

Parameters
InExisting ExpList object to copy.

Definition at line 766 of file DisContField.cpp.

766  :
767  ExpList(In)
768  {
769  }

◆ ~DisContField()

Nektar::MultiRegions::DisContField::~DisContField ( )
virtual

Destructor.

Definition at line 774 of file DisContField.cpp.

775  {
776  }

Member Function Documentation

◆ EvaluateHDGPostProcessing()

void Nektar::MultiRegions::DisContField::EvaluateHDGPostProcessing ( Array< OneD, NekDouble > &  outarray)

Evaluate HDG post-processing to increase polynomial order of solution.

This function takes the solution (assumed to be one order lower) in physical space, and postprocesses at the current polynomial order by solving the system:

\[ \begin{aligned} (\nabla w, \nabla u^*) &= (\nabla w, u), \\ \langle \nabla u^*, 1 \rangle &= \langle \nabla u, 1 \rangle \end{aligned} \]

where \( u \) corresponds with the current solution as stored inside m_coeffs.

Parameters
outarrayThe resulting field \( u^* \).

Definition at line 4192 of file DisContField.cpp.

4194  {
4195  int i,cnt,e,ncoeff_trace;
4196  Array<OneD, NekDouble> force, out_tmp, qrhs, qrhs1;
4197  Array<OneD, Array< OneD, LocalRegions::ExpansionSharedPtr> >
4198  &elmtToTrace = m_traceMap->GetElmtToTrace();
4199 
4200  StdRegions::Orientation edgedir;
4201 
4202  int nq_elmt, nm_elmt;
4203  int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs();
4204  Array<OneD, NekDouble> loc_lambda(LocBndCoeffs), trace_lambda;
4205  Array<OneD, NekDouble> tmp_coeffs;
4206  m_traceMap->GlobalToLocalBnd(m_trace->GetCoeffs(),loc_lambda);
4207 
4208  trace_lambda = loc_lambda;
4209 
4210  int dim = (m_expType == e2D)? 2:3;
4211 
4212  int num_points[] = {0,0,0};
4213  int num_modes [] = {0,0,0};
4214 
4215  // Calculate Q using standard DG formulation.
4216  for(i = cnt = 0; i < GetExpSize(); ++i)
4217  {
4218  nq_elmt = (*m_exp)[i]->GetTotPoints();
4219  nm_elmt = (*m_exp)[i]->GetNcoeffs();
4220  qrhs = Array<OneD, NekDouble>(nq_elmt);
4221  qrhs1 = Array<OneD, NekDouble>(nq_elmt);
4222  force = Array<OneD, NekDouble>(2*nm_elmt);
4223  out_tmp = force + nm_elmt;
4225 
4226  for(int j= 0; j < dim; ++j)
4227  {
4228  num_points[j] = (*m_exp)[i]->GetBasis(j)->GetNumPoints();
4229  num_modes[j] = (*m_exp)[i]->GetBasis(j)->GetNumModes();
4230  }
4231 
4232  // Probably a better way of setting up lambda than this. Note
4233  // cannot use PutCoeffsInToElmts since lambda space is mapped
4234  // during the solve.
4235  int nTraces = (*m_exp)[i]->GetNtraces();
4236  Array<OneD, Array<OneD, NekDouble> > traceCoeffs(nTraces);
4237 
4238  for(e = 0; e < (*m_exp)[i]->GetNtraces(); ++e)
4239  {
4240  edgedir = (*m_exp)[i]->GetTraceOrient(e);
4241  ncoeff_trace = elmtToTrace[i][e]->GetNcoeffs();
4242  traceCoeffs[e] = Array<OneD, NekDouble>(ncoeff_trace);
4243  Vmath::Vcopy(ncoeff_trace, trace_lambda, 1, traceCoeffs[e], 1);
4244  if(dim == 2)
4245  {
4246  elmtToTrace[i][e]->SetCoeffsToOrientation(edgedir,
4247  traceCoeffs[e], traceCoeffs[e]);
4248  }
4249  else
4250  {
4251  (*m_exp)[i]->as<LocalRegions::Expansion3D>()->
4252  SetFaceToGeomOrientation(e,traceCoeffs[e]);
4253 
4254  }
4255  trace_lambda = trace_lambda + ncoeff_trace;
4256  }
4257 
4258  //creating orthogonal expansion (checking if we have quads or triangles)
4259  LibUtilities::ShapeType shape = (*m_exp)[i]->DetShapeType();
4260  switch(shape)
4261  {
4263  {
4264  const LibUtilities::PointsKey PkeyQ1(num_points[0],
4266  const LibUtilities::PointsKey PkeyQ2(num_points[1],
4268  LibUtilities::BasisKey BkeyQ1(LibUtilities::eOrtho_A, num_modes[0],
4269  PkeyQ1);
4270  LibUtilities::BasisKey BkeyQ2(LibUtilities::eOrtho_A, num_modes[1],
4271  PkeyQ2);
4273  std::dynamic_pointer_cast<SpatialDomains::QuadGeom>
4274  ((*m_exp)[i]->GetGeom());
4276  (BkeyQ1, BkeyQ2, qGeom);
4277  }
4278  break;
4280  {
4281  const LibUtilities::PointsKey PkeyT1(num_points[0],
4283  const LibUtilities::PointsKey PkeyT2(num_points[1],
4285  LibUtilities::BasisKey BkeyT1(LibUtilities::eOrtho_A, num_modes[0],
4286  PkeyT1);
4287  LibUtilities::BasisKey BkeyT2(LibUtilities::eOrtho_B, num_modes[1],
4288  PkeyT2);
4289  SpatialDomains::TriGeomSharedPtr tGeom = std::dynamic_pointer_cast
4290  <SpatialDomains::TriGeom>((*m_exp)[i]->GetGeom());
4292  (BkeyT1, BkeyT2, tGeom);
4293  }
4294  break;
4296  {
4297  const LibUtilities::PointsKey PkeyH1(num_points[0],
4299  const LibUtilities::PointsKey PkeyH2(num_points[1],
4301  const LibUtilities::PointsKey PkeyH3(num_points[2],
4303  LibUtilities::BasisKey BkeyH1(LibUtilities::eOrtho_A,
4304  num_modes[0], PkeyH1);
4305  LibUtilities::BasisKey BkeyH2(LibUtilities::eOrtho_A,
4306  num_modes[1], PkeyH2);
4307  LibUtilities::BasisKey BkeyH3(LibUtilities::eOrtho_A,
4308  num_modes[2], PkeyH3);
4310  std::dynamic_pointer_cast<SpatialDomains::HexGeom>
4311  ((*m_exp)[i]->GetGeom());
4313  (BkeyH1, BkeyH2, BkeyH3, hGeom);
4314  }
4315  break;
4317  {
4318  const LibUtilities::PointsKey PkeyT1(num_points[0],
4320  const LibUtilities::PointsKey PkeyT2(num_points[1],
4322  const LibUtilities::PointsKey PkeyT3(num_points[2],
4324  LibUtilities::BasisKey BkeyT1(LibUtilities::eOrtho_A,
4325  num_modes[0], PkeyT1);
4326  LibUtilities::BasisKey BkeyT2(LibUtilities::eOrtho_B,
4327  num_modes[1], PkeyT2);
4328  LibUtilities::BasisKey BkeyT3(LibUtilities::eOrtho_C,
4329  num_modes[2], PkeyT3);
4331  std::dynamic_pointer_cast<SpatialDomains::TetGeom>
4332  ((*m_exp)[i]->GetGeom());
4334  (BkeyT1, BkeyT2, BkeyT3, tGeom);
4335  }
4336  break;
4337  default:
4339  "Wrong shape type, HDG postprocessing is not "
4340  "implemented");
4341  };
4342 
4343 
4344  //DGDeriv
4345  // (d/dx w, d/dx q_0)
4346  (*m_exp)[i]->DGDeriv(
4347  0,tmp_coeffs = m_coeffs + m_coeff_offset[i],
4348  elmtToTrace[i], traceCoeffs, out_tmp);
4349  (*m_exp)[i]->BwdTrans(out_tmp,qrhs);
4350  ppExp->IProductWRTDerivBase(0,qrhs,force);
4351 
4352  // + (d/dy w, d/dy q_1)
4353  (*m_exp)[i]->DGDeriv(
4354  1,tmp_coeffs = m_coeffs + m_coeff_offset[i],
4355  elmtToTrace[i], traceCoeffs, out_tmp);
4356 
4357  (*m_exp)[i]->BwdTrans(out_tmp,qrhs);
4358  ppExp->IProductWRTDerivBase(1,qrhs,out_tmp);
4359 
4360  Vmath::Vadd(nm_elmt,force,1,out_tmp,1,force,1);
4361 
4362  // determine force[0] = (1,u)
4363  (*m_exp)[i]->BwdTrans(
4364  tmp_coeffs = m_coeffs + m_coeff_offset[i],qrhs);
4365  force[0] = (*m_exp)[i]->Integral(qrhs);
4366 
4367  // multiply by inverse Laplacian matrix
4368  // get matrix inverse
4369  LocalRegions::MatrixKey lapkey(StdRegions::eInvLaplacianWithUnityMean,
4370  ppExp->DetShapeType(), *ppExp);
4371  DNekScalMatSharedPtr lapsys = ppExp->GetLocMatrix(lapkey);
4372 
4373  NekVector<NekDouble> in (nm_elmt,force,eWrapper);
4374  NekVector<NekDouble> out(nm_elmt);
4375 
4376  out = (*lapsys)*in;
4377 
4378  // Transforming back to modified basis
4379  Array<OneD, NekDouble> work(nq_elmt);
4380  ppExp->BwdTrans(out.GetPtr(), work);
4381  (*m_exp)[i]->FwdTrans(work, tmp_coeffs = outarray + m_coeff_offset[i]);
4382  }
4383  }
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition: ErrorUtil.hpp:209
Array< OneD, NekDouble > m_coeffs
Concatenation of all local expansion coefficients.
Definition: ExpList.h:1252
int GetExpSize(void)
This function returns the number of elements in the expansion.
Definition: ExpList.h:2401
Array< OneD, int > m_coeff_offset
Offset of elemental data into the array m_coeffs.
Definition: ExpList.h:1304
ExpansionType m_expType
Exapnsion type.
Definition: ExpList.h:1212
@ eGaussRadauMAlpha2Beta0
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:59
@ eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:51
@ eGaussRadauMAlpha1Beta0
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:58
@ eOrtho_A
Principle Orthogonal Functions .
Definition: BasisType.h:45
@ eOrtho_C
Principle Orthogonal Functions .
Definition: BasisType.h:47
@ eOrtho_B
Principle Orthogonal Functions .
Definition: BasisType.h:46
std::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
std::shared_ptr< QuadGeom > QuadGeomSharedPtr
Definition: HexGeom.h:46
std::shared_ptr< TetGeom > TetGeomSharedPtr
Definition: TetGeom.h:82
std::shared_ptr< HexGeom > HexGeomSharedPtr
Definition: HexGeom.h:84
std::shared_ptr< TriGeom > TriGeomSharedPtr
Definition: TriGeom.h:58
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:322
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1199

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::MultiRegions::e2D, Nektar::ErrorUtil::efatal, Nektar::LibUtilities::eGaussLobattoLegendre, Nektar::LibUtilities::eGaussRadauMAlpha1Beta0, Nektar::LibUtilities::eGaussRadauMAlpha2Beta0, Nektar::LibUtilities::eHexahedron, Nektar::StdRegions::eInvLaplacianWithUnityMean, Nektar::LibUtilities::eOrtho_A, Nektar::LibUtilities::eOrtho_B, Nektar::LibUtilities::eOrtho_C, Nektar::LibUtilities::eQuadrilateral, Nektar::LibUtilities::eTetrahedron, Nektar::LibUtilities::eTriangle, Nektar::eWrapper, Nektar::MultiRegions::ExpList::GetExpSize(), Nektar::NekVector< DataType >::GetPtr(), Nektar::MultiRegions::ExpList::m_coeff_offset, Nektar::MultiRegions::ExpList::m_coeffs, Nektar::MultiRegions::ExpList::m_exp, Nektar::MultiRegions::ExpList::m_expType, m_trace, m_traceMap, NEKERROR, Vmath::Vadd(), and Vmath::Vcopy().

◆ FindPeriodicTraces()

void Nektar::MultiRegions::DisContField::FindPeriodicTraces ( const SpatialDomains::BoundaryConditions bcs,
const std::string  variable 
)
protected

Generate a associative map of periodic vertices in a mesh.

Determine the periodic faces, edges and vertices for the given graph.

Parameters
bcsInformation about the boundary conditions.
variableSpecifies the field.

Definition at line 852 of file DisContField.cpp.

855  {
857  = bcs.GetBoundaryRegions();
859  = bcs.GetBoundaryConditions();
860 
862  m_session->GetComm()->GetRowComm();
863 
864 
865  switch(m_expType)
866  {
867  case e1D:
868  {
869  int i, region1ID, region2ID;
870 
872  map<int,int> BregionToVertMap;
873 
874  // Construct list of all periodic Region and their
875  // global vertex on this process.
876  for (auto &it : bregions)
877  {
878  locBCond = GetBoundaryCondition(bconditions,
879  it.first, variable);
880 
881  if (locBCond->GetBoundaryConditionType()
883  {
884  continue;
885  }
886  int id = it.second->begin()->second->
887  m_geomVec[0]->GetGlobalID();
888 
889  BregionToVertMap[it.first] = id;
890  }
891 
892  set<int> islocal;
893 
894  int n = vComm->GetSize();
895  int p = vComm->GetRank();
896 
897  Array<OneD, int> nregions(n, 0);
898  nregions[p] = BregionToVertMap.size();
899  vComm->AllReduce(nregions, LibUtilities::ReduceSum);
900 
901  int totRegions = Vmath::Vsum(n, nregions, 1);
902 
903  Array<OneD, int> regOffset(n, 0);
904 
905  for (i = 1; i < n; ++i)
906  {
907  regOffset[i] = regOffset[i-1] + nregions[i-1];
908  }
909 
910  Array<OneD, int> bregmap(totRegions, 0);
911  Array<OneD, int> bregid (totRegions, 0);
912 
913  i = regOffset[p];
914  for (auto &iit : BregionToVertMap)
915  {
916  bregid [i ] = iit.first;
917  bregmap[i++] = iit.second;
918  islocal.insert(iit.first);
919  }
920 
921  vComm->AllReduce(bregmap, LibUtilities::ReduceSum);
922  vComm->AllReduce(bregid, LibUtilities::ReduceSum);
923 
924  for (int i = 0; i < totRegions; ++i)
925  {
926  BregionToVertMap[bregid[i]] = bregmap[i];
927  }
928 
929  // Construct list of all periodic pairs local to this process.
930  for (auto &it : bregions)
931  {
932  locBCond = GetBoundaryCondition(bconditions,
933  it.first, variable);
934 
935  if (locBCond->GetBoundaryConditionType()
937  {
938  continue;
939  }
940 
941  // Identify periodic boundary region IDs.
942  region1ID = it.first;
943  region2ID = std::static_pointer_cast<
944  SpatialDomains::PeriodicBoundaryCondition>(
945  locBCond)->m_connectedBoundaryRegion;
946 
947  ASSERTL0(BregionToVertMap.count(region1ID) != 0,
948  "Cannot determine vertex of region1ID");
949 
950  ASSERTL0(BregionToVertMap.count(region2ID) != 0,
951  "Cannot determine vertex of region2ID");
952 
953  PeriodicEntity ent(BregionToVertMap[region2ID],
955  islocal.count(region2ID) != 0);
956 
957  m_periodicVerts[BregionToVertMap[region1ID]].push_back(ent);
958  }
959  }
960  break;
961  case e2D:
962  {
963  int region1ID, region2ID, i, j, k, cnt;
965 
967  m_graph->GetCompositeOrdering();
969  m_graph->GetBndRegionOrdering();
971  m_graph->GetComposites();
972 
973  // Unique collection of pairs of periodic composites
974  // (i.e. if composites 1 and 2 are periodic then this
975  // map will contain either the pair (1,2) or (2,1) but
976  // not both).
977  map<int,int> perComps;
978  map<int,vector<int>> allVerts;
979  set<int> locVerts;
980  map<int, pair<int, StdRegions::Orientation>> allEdges;
981 
982  // Set up a set of all local verts and edges.
983  for(i = 0; i < (*m_exp).size(); ++i)
984  {
985  for(j = 0; j < (*m_exp)[i]->GetNverts(); ++j)
986  {
987  int id = (*m_exp)[i]->GetGeom()->GetVid(j);
988  locVerts.insert(id);
989  }
990  }
991 
992  // Construct list of all periodic pairs local to this process.
993  for (auto &it : bregions)
994  {
995  locBCond = GetBoundaryCondition(
996  bconditions, it.first, variable);
997 
998  if (locBCond->GetBoundaryConditionType()
1000  {
1001  continue;
1002  }
1003 
1004  // Identify periodic boundary region IDs.
1005  region1ID = it.first;
1006  region2ID = std::static_pointer_cast<
1007  SpatialDomains::PeriodicBoundaryCondition>(
1008  locBCond)->m_connectedBoundaryRegion;
1009 
1010  // From this identify composites. Note that in
1011  // serial this will be an empty map.
1012  int cId1, cId2;
1013  if (vComm->GetSize() == 1)
1014  {
1015  cId1 = it.second->begin()->first;
1016  cId2 = bregions.find(region2ID)->second->begin()->first;
1017  }
1018  else
1019  {
1020  cId1 = bndRegOrder.find(region1ID)->second[0];
1021  cId2 = bndRegOrder.find(region2ID)->second[0];
1022  }
1023 
1024  ASSERTL0(it.second->size() == 1,
1025  "Boundary region "+boost::lexical_cast<string>(
1026  region1ID)+" should only contain 1 composite.");
1027 
1028  vector<unsigned int> tmpOrder;
1029 
1030  // Construct set containing all periodic edgesd on
1031  // this process
1033  c = it.second->begin()->second;
1034 
1035  for (i = 0; i < c->m_geomVec.size(); ++i)
1036  {
1038  std::dynamic_pointer_cast<
1039  SpatialDomains::SegGeom>(c->m_geomVec[i]);
1040  ASSERTL0(segGeom, "Unable to cast to shared ptr");
1041 
1043  m_graph->GetElementsFromEdge(segGeom);
1044  ASSERTL0(elmt->size() == 1,
1045  "The periodic boundaries belong to "
1046  "more than one element of the mesh");
1047 
1049  std::dynamic_pointer_cast
1050  <SpatialDomains::Geometry2D>(elmt->at(0).first);
1051 
1052  allEdges[c->m_geomVec[i]->GetGlobalID()] =
1053  make_pair(elmt->at(0).second,
1054  geom->GetEorient(elmt->at(0).second));
1055 
1056  // In serial mesh partitioning will not have occurred so
1057  // need to fill composite ordering map manually.
1058  if (vComm->GetSize() == 1)
1059  {
1060  tmpOrder.push_back(c->m_geomVec[i]->GetGlobalID());
1061  }
1062 
1063  vector<int> vertList(2);
1064  vertList[0] = segGeom->GetVid(0);
1065  vertList[1] = segGeom->GetVid(1);
1066  allVerts[c->m_geomVec[i]->GetGlobalID()] = vertList;
1067  }
1068 
1069  if (vComm->GetSize() == 1)
1070  {
1071  compOrder[it.second->begin()->first] = tmpOrder;
1072  }
1073 
1074 
1075  // See if we already have either region1 or
1076  // region2 stored in perComps map.
1077  if (perComps.count(cId1) == 0)
1078  {
1079  if (perComps.count(cId2) == 0)
1080  {
1081  perComps[cId1] = cId2;
1082  }
1083  else
1084  {
1085  std::stringstream ss;
1086  ss << "Boundary region " << cId2 << " should be "
1087  << "periodic with " << perComps[cId2] << " but "
1088  << "found " << cId1 << " instead!";
1089  ASSERTL0(perComps[cId2] == cId1, ss.str());
1090  }
1091  }
1092  else
1093  {
1094  std::stringstream ss;
1095  ss << "Boundary region " << cId1 << " should be "
1096  << "periodic with " << perComps[cId1] << " but "
1097  << "found " << cId2 << " instead!";
1098  ASSERTL0(perComps[cId1] == cId1, ss.str());
1099  }
1100  }
1101 
1102 
1103  // Process local edge list to obtain relative edge orientations.
1104  int n = vComm->GetSize();
1105  int p = vComm->GetRank();
1106  int totEdges;
1107  Array<OneD, int> edgecounts(n,0);
1108  Array<OneD, int> edgeoffset(n,0);
1109  Array<OneD, int> vertoffset(n,0);
1110 
1111  edgecounts[p] = allEdges.size();
1112  vComm->AllReduce(edgecounts, LibUtilities::ReduceSum);
1113 
1114  edgeoffset[0] = 0;
1115  for (i = 1; i < n; ++i)
1116  {
1117  edgeoffset[i] = edgeoffset[i-1] + edgecounts[i-1];
1118  }
1119 
1120  totEdges = Vmath::Vsum(n, edgecounts, 1);
1121  Array<OneD, int> edgeIds (totEdges, 0);
1122  Array<OneD, int> edgeIdx (totEdges, 0);
1123  Array<OneD, int> edgeOrient(totEdges, 0);
1124  Array<OneD, int> edgeVerts (totEdges, 0);
1125 
1126  auto sIt = allEdges.begin();
1127 
1128  for (i = 0; sIt != allEdges.end(); ++sIt)
1129  {
1130  edgeIds [edgeoffset[p] + i ] = sIt->first;
1131  edgeIdx [edgeoffset[p] + i ] = sIt->second.first;
1132  edgeOrient[edgeoffset[p] + i ] = sIt->second.second;
1133  edgeVerts [edgeoffset[p] + i++] = allVerts[sIt->first].size();
1134  }
1135 
1136  vComm->AllReduce(edgeIds, LibUtilities::ReduceSum);
1137  vComm->AllReduce(edgeIdx, LibUtilities::ReduceSum);
1138  vComm->AllReduce(edgeOrient, LibUtilities::ReduceSum);
1139  vComm->AllReduce(edgeVerts, LibUtilities::ReduceSum);
1140 
1141  // Calculate number of vertices on each processor.
1142  Array<OneD, int> procVerts(n,0);
1143  int nTotVerts;
1144 
1145  // Note if there are no periodic edges at all calling Vsum will
1146  // cause a segfault.
1147  if (totEdges > 0)
1148  {
1149  nTotVerts = Vmath::Vsum(totEdges, edgeVerts, 1);
1150  }
1151  else
1152  {
1153  nTotVerts = 0;
1154  }
1155 
1156  for (i = 0; i < n; ++i)
1157  {
1158  if (edgecounts[i] > 0)
1159  {
1160  procVerts[i] = Vmath::Vsum(
1161  edgecounts[i], edgeVerts + edgeoffset[i], 1);
1162  }
1163  else
1164  {
1165  procVerts[i] = 0;
1166  }
1167  }
1168  vertoffset[0] = 0;
1169 
1170  for (i = 1; i < n; ++i)
1171  {
1172  vertoffset[i] = vertoffset[i-1] + procVerts[i-1];
1173  }
1174 
1175  Array<OneD, int> traceIds(nTotVerts, 0);
1176  for (i = 0, sIt = allEdges.begin(); sIt != allEdges.end(); ++sIt)
1177  {
1178  for (j = 0; j < allVerts[sIt->first].size(); ++j)
1179  {
1180  traceIds[vertoffset[p] + i++] = allVerts[sIt->first][j];
1181  }
1182  }
1183 
1184  vComm->AllReduce(traceIds, LibUtilities::ReduceSum);
1185 
1186  // For simplicity's sake create a map of edge id -> orientation.
1187  map<int, StdRegions::Orientation> orientMap;
1188  map<int, vector<int> > vertMap;
1189 
1190  for (cnt = i = 0; i < totEdges; ++i)
1191  {
1192  ASSERTL0(orientMap.count(edgeIds[i]) == 0,
1193  "Already found edge in orientation map!");
1194 
1195  // Work out relative orientations. To avoid having
1196  // to exchange vertex locations, we figure out if
1197  // the edges are backwards or forwards orientated
1198  // with respect to a forwards orientation that is
1199  // CCW. Since our local geometries are
1200  // forwards-orientated with respect to the
1201  // Cartesian axes, we need to invert the
1202  // orientation for the top and left edges of a
1203  // quad and the left edge of a triangle.
1205  (StdRegions::Orientation)edgeOrient[i];
1206 
1207  if (edgeIdx[i] > 1)
1208  {
1209  o = o == StdRegions::eForwards ?
1211  }
1212 
1213  orientMap[edgeIds[i]] = o;
1214 
1215  vector<int> verts(edgeVerts[i]);
1216 
1217  for (j = 0; j < edgeVerts[i]; ++j)
1218  {
1219  verts[j] = traceIds[cnt++];
1220  }
1221  vertMap[edgeIds[i]] = verts;
1222  }
1223 
1224  // Go through list of composites and figure out which
1225  // edges are parallel from original ordering in
1226  // session file. This includes composites which are
1227  // not necessarily on this process.
1228  map<int,int> allCompPairs;
1229 
1230  // Store temporary map of periodic vertices which will hold all
1231  // periodic vertices on the entire mesh so that doubly periodic
1232  // vertices can be counted properly across partitions. Local
1233  // vertices are copied into m_periodicVerts at the end of the
1234  // function.
1235  PeriodicMap periodicVerts;
1236 
1237  for (auto &cIt : perComps)
1238  {
1240  const int id1 = cIt.first;
1241  const int id2 = cIt.second;
1242  std::string id1s = boost::lexical_cast<string>(id1);
1243  std::string id2s = boost::lexical_cast<string>(id2);
1244 
1245  if (compMap.count(id1) > 0)
1246  {
1247  c[0] = compMap[id1];
1248  }
1249 
1250  if (compMap.count(id2) > 0)
1251  {
1252  c[1] = compMap[id2];
1253  }
1254 
1255  ASSERTL0(c[0] || c[1],
1256  "Both composites not found on this process!");
1257 
1258  // Loop over composite ordering to construct list of all
1259  // periodic edges regardless of whether they are on this
1260  // process.
1261  map<int,int> compPairs;
1262 
1263  ASSERTL0(compOrder.count(id1) > 0,
1264  "Unable to find composite "+id1s+" in order map.");
1265  ASSERTL0(compOrder.count(id2) > 0,
1266  "Unable to find composite "+id2s+" in order map.");
1267  ASSERTL0(compOrder[id1].size() == compOrder[id2].size(),
1268  "Periodic composites "+id1s+" and "+id2s+
1269  " should have the same number of elements.");
1270  ASSERTL0(compOrder[id1].size() > 0,
1271  "Periodic composites "+id1s+" and "+id2s+
1272  " are empty!");
1273 
1274  // TODO: Add more checks.
1275  for (i = 0; i < compOrder[id1].size(); ++i)
1276  {
1277  int eId1 = compOrder[id1][i];
1278  int eId2 = compOrder[id2][i];
1279 
1280  ASSERTL0(compPairs.count(eId1) == 0,
1281  "Already paired.");
1282 
1283  if (compPairs.count(eId2) != 0)
1284  {
1285  ASSERTL0(compPairs[eId2] == eId1, "Pairing incorrect");
1286  }
1287  compPairs[eId1] = eId2;
1288  }
1289 
1290  // Construct set of all edges that we have locally on this
1291  // processor.
1292  set<int> locEdges;
1293  for (i = 0; i < 2; ++i)
1294  {
1295  if (!c[i])
1296  {
1297  continue;
1298  }
1299 
1300  if (c[i]->m_geomVec.size() > 0)
1301  {
1302  for (j = 0; j < c[i]->m_geomVec.size(); ++j)
1303  {
1304  locEdges.insert(c[i]->m_geomVec[j]->GetGlobalID());
1305  }
1306  }
1307  }
1308 
1309  // Loop over all edges in the geometry composite.
1310  for (auto &pIt : compPairs)
1311  {
1312  int ids [2] = {pIt.first, pIt.second};
1313  bool local[2] = {locEdges.count(pIt.first) > 0,
1314  locEdges.count(pIt.second) > 0};
1315 
1316  ASSERTL0(orientMap.count(ids[0]) > 0 &&
1317  orientMap.count(ids[1]) > 0,
1318  "Unable to find edge in orientation map.");
1319 
1320  allCompPairs[pIt.first ] = pIt.second;
1321  allCompPairs[pIt.second] = pIt.first;
1322 
1323  for (i = 0; i < 2; ++i)
1324  {
1325  if (!local[i])
1326  {
1327  continue;
1328  }
1329 
1330  int other = (i+1) % 2;
1331 
1333  orientMap[ids[i]] == orientMap[ids[other]] ?
1336 
1337  PeriodicEntity ent(ids [other], o,
1338  local[other]);
1339  m_periodicEdges[ids[i]].push_back(ent);
1340  }
1341 
1342  for (i = 0; i < 2; ++i)
1343  {
1344  int other = (i+1) % 2;
1345 
1347  orientMap[ids[i]] == orientMap[ids[other]] ?
1350 
1351  // Determine periodic vertices.
1352  vector<int> perVerts1 = vertMap[ids[i]];
1353  vector<int> perVerts2 = vertMap[ids[other]];
1354 
1355  map<int, pair<int, bool> > tmpMap;
1356  if (o == StdRegions::eForwards)
1357  {
1358  tmpMap[perVerts1[0]] = make_pair(perVerts2[0],
1359  locVerts.count(perVerts2[0]) > 0);
1360  tmpMap[perVerts1[1]] = make_pair(perVerts2[1],
1361  locVerts.count(perVerts2[1]) > 0);
1362  }
1363  else
1364  {
1365  tmpMap[perVerts1[0]] = make_pair(perVerts2[1],
1366  locVerts.count(perVerts2[1]) > 0);
1367  tmpMap[perVerts1[1]] = make_pair(
1368  perVerts2[0], locVerts.count(perVerts2[0]) > 0);
1369  }
1370 
1371  for (auto &mIt : tmpMap)
1372  {
1373  // See if this vertex has been recorded already.
1374  PeriodicEntity ent2(mIt.second.first,
1376  mIt.second.second);
1377  auto perIt = periodicVerts.find(mIt.first);
1378 
1379  if (perIt == periodicVerts.end())
1380  {
1381  periodicVerts[mIt.first].push_back(ent2);
1382  perIt = periodicVerts.find(mIt.first);
1383  }
1384  else
1385  {
1386  bool doAdd = true;
1387  for (j = 0; j < perIt->second.size(); ++j)
1388  {
1389  if (perIt->second[j].id ==
1390  mIt.second.first)
1391  {
1392  doAdd = false;
1393  break;
1394  }
1395  }
1396 
1397  if (doAdd)
1398  {
1399  perIt->second.push_back(ent2);
1400  }
1401  }
1402  }
1403  }
1404  }
1405  }
1406 
1407  // Distribute the size of the periodic boundary to all
1408  // processes. We assume that all processes who own periodic
1409  // faces have the same local copy of allCompPairs
1410  int NPairs = allCompPairs.size();
1411  vComm->AllReduce(NPairs, LibUtilities::ReduceMax);
1412 
1413  // Check that the previous assertion regarding allCompPairs
1414  // is correct
1415  ASSERTL0(
1416  allCompPairs.size() == NPairs || allCompPairs.size() == 0,
1417  "Local copy of allCompPairs not the same for all ranks"
1418  );
1419 
1420  // Allocate local vectors that will contain the content of
1421  // allCompPairs if process owns faces on periodic boundary
1422  Array<OneD, int> first(NPairs, -1);
1423  Array<OneD, int> second(NPairs, -1);
1424  cnt = 0;
1425  for(const auto &it : allCompPairs)
1426  {
1427  first[cnt] = it.first;
1428  second[cnt++] = it.second;
1429  }
1430 
1431  // Distribute the content in first and second to all processes
1432  vComm->AllReduce(first, LibUtilities::ReduceMax);
1433  vComm->AllReduce(second, LibUtilities::ReduceMax);
1434 
1435  // Check that the MPI Allreduce routine worked
1436  ASSERTL0(std::count(first.begin(), first.end(), -1) == 0,
1437  "Distribution of allCompPairs failed");
1438  ASSERTL0(std::count(second.begin(), second.end(), -1) == 0,
1439  "Distribution of allCompParis failed")
1440 
1441  // Put content back in allCompPairs
1442  allCompPairs.clear();
1443  for(cnt = 0; cnt < NPairs; ++cnt)
1444  {
1445  allCompPairs[first[cnt]] = second[cnt];
1446  }
1447 
1448  // Search for periodic vertices and edges which are not in
1449  // a periodic composite but lie in this process. First, loop
1450  // over all information we have from other processors.
1451  for (cnt = i = 0; i < totEdges; ++i)
1452  {
1453  int edgeId = edgeIds[i];
1454 
1455  ASSERTL0(allCompPairs.count(edgeId) > 0,
1456  "Unable to find matching periodic edge.");
1457 
1458  int perEdgeId = allCompPairs[edgeId];
1459 
1460  for (j = 0; j < edgeVerts[i]; ++j, ++cnt)
1461  {
1462  int vId = traceIds[cnt];
1463 
1464  auto perId = periodicVerts.find(vId);
1465 
1466  if (perId == periodicVerts.end())
1467  {
1468  // This vertex is not included in the
1469  // map. Figure out which vertex it is
1470  // supposed to be periodic with. perEdgeId
1471  // is the edge ID which is periodic with
1472  // edgeId. The logic is much the same as
1473  // the loop above.
1474  int perVertexId =
1475  orientMap[edgeId] == orientMap[perEdgeId] ?
1476  vertMap[perEdgeId][(j+1)%2]:
1477  vertMap[perEdgeId][j];
1478 
1479  PeriodicEntity ent(perVertexId,
1481  locVerts.count(perVertexId) > 0);
1482 
1483  periodicVerts[vId].push_back(ent);
1484  }
1485  }
1486  }
1487 
1488  // Loop over all periodic vertices to complete connectivity
1489  // information.
1490  for (auto &perIt : periodicVerts)
1491  {
1492  // Loop over associated vertices.
1493  for (i = 0; i < perIt.second.size(); ++i)
1494  {
1495  auto perIt2 = periodicVerts.find(perIt.second[i].id);
1496  ASSERTL0(perIt2 != periodicVerts.end(),
1497  "Couldn't find periodic vertex.");
1498 
1499  for (j = 0; j < perIt2->second.size(); ++j)
1500  {
1501  if (perIt2->second[j].id == perIt.first)
1502  {
1503  continue;
1504  }
1505 
1506  bool doAdd = true;
1507  for (k = 0; k < perIt.second.size(); ++k)
1508  {
1509  if (perIt2->second[j].id == perIt.second[k].id)
1510  {
1511  doAdd = false;
1512  break;
1513  }
1514  }
1515 
1516  if (doAdd)
1517  {
1518  perIt.second.push_back(perIt2->second[j]);
1519  }
1520  }
1521  }
1522  }
1523 
1524  // Do one final loop over periodic vertices to remove non-local
1525  // vertices from map.
1526  for (auto &perIt : periodicVerts)
1527  {
1528  if (locVerts.count(perIt.first) > 0)
1529  {
1530  m_periodicVerts.insert(perIt);
1531  }
1532  }
1533  }
1534  break;
1535  case e3D:
1536  {
1538  m_graph->GetCompositeOrdering();
1539  SpatialDomains::BndRegionOrdering bndRegOrder =
1540  m_graph->GetBndRegionOrdering();
1542  m_graph->GetComposites();
1543 
1544  // perComps: Stores a unique collection of pairs of periodic
1545  // composites (i.e. if composites 1 and 2 are periodic then this map
1546  // will contain either the pair (1,2) or (2,1) but not both).
1547  //
1548  // The four maps allVerts, allCoord, allEdges and allOrient map a
1549  // periodic face to a vector containing the vertex ids of the face;
1550  // their coordinates; the edge ids of the face; and their
1551  // orientation within that face respectively.
1552  //
1553  // Finally the three sets locVerts, locEdges and locFaces store any
1554  // vertices, edges and faces that belong to a periodic composite and
1555  // lie on this process.
1556  map<int,RotPeriodicInfo> rotComp;
1557  map<int,int> perComps;
1558  map<int,vector<int> > allVerts;
1559  map<int,SpatialDomains::PointGeomVector> allCoord;
1560  map<int,vector<int> > allEdges;
1561  map<int,vector<StdRegions::Orientation> > allOrient;
1562  set<int> locVerts;
1563  set<int> locEdges;
1564  set<int> locFaces;
1565 
1566  int region1ID, region2ID, i, j, k, cnt;
1568 
1569  // Set up a set of all local verts and edges.
1570  for(i = 0; i < (*m_exp).size(); ++i)
1571  {
1572  for(j = 0; j < (*m_exp)[i]->GetNverts(); ++j)
1573  {
1574  int id = (*m_exp)[i]->GetGeom()->GetVid(j);
1575  locVerts.insert(id);
1576  }
1577 
1578  for(j = 0; j < (*m_exp)[i]->GetGeom()->GetNumEdges(); ++j)
1579  {
1580  int id = (*m_exp)[i]->GetGeom()->GetEid(j);
1581  locEdges.insert(id);
1582  }
1583  }
1584 
1585  // Begin by populating the perComps map. We loop over all periodic
1586  // boundary conditions and determine the composite associated with
1587  // it, then fill out the all* maps.
1588  for (auto &it : bregions)
1589  {
1590 
1591  locBCond = GetBoundaryCondition
1592  (bconditions, it.first, variable);
1593 
1594  if (locBCond->GetBoundaryConditionType()
1596  {
1597  continue;
1598  }
1599 
1600  // Identify periodic boundary region IDs.
1601  region1ID = it.first;
1602  region2ID = std::static_pointer_cast<
1603  SpatialDomains::PeriodicBoundaryCondition>
1604  (locBCond)->m_connectedBoundaryRegion;
1605 
1606  // Check the region only contains a single composite.
1607  ASSERTL0(it.second->size() == 1,
1608  "Boundary region "+boost::lexical_cast<string>
1609  (region1ID)+" should only contain 1 composite.");
1610 
1611  // From this identify composites by looking at the original
1612  // boundary region ordering. Note that in serial the mesh
1613  // partitioner is not run, so this map will be empty and
1614  // therefore needs to be populated by using the corresponding
1615  // boundary region.
1616  int cId1, cId2;
1617  if (vComm->GetSize() == 1)
1618  {
1619  cId1 = it.second->begin()->first;
1620  cId2 = bregions.find(region2ID)->second->begin()->first;
1621  }
1622  else
1623  {
1624  cId1 = bndRegOrder.find(region1ID)->second[0];
1625  cId2 = bndRegOrder.find(region2ID)->second[0];
1626  }
1627 
1628  // check to see if boundary is rotationally aligned
1629  if(boost::icontains(locBCond->GetUserDefined(),"Rotated"))
1630  {
1631  vector<string> tmpstr;
1632 
1633  boost::split(tmpstr,locBCond->GetUserDefined(),
1634  boost::is_any_of(":"));
1635 
1636  if(boost::iequals(tmpstr[0],"Rotated"))
1637  {
1638  ASSERTL1(tmpstr.size() > 2,
1639  "Expected Rotated user defined string to "
1640  "contain direction and rotation angle "
1641  "and optionally a tolerance, "
1642  "i.e. Rotated:dir:PI/2:1e-6");
1643 
1644 
1645  ASSERTL1((tmpstr[1] == "x")||(tmpstr[1] == "y")
1646  ||(tmpstr[1] == "z"), "Rotated Dir is "
1647  "not specified as x,y or z");
1648 
1649  RotPeriodicInfo RotInfo;
1650  RotInfo.m_dir = (tmpstr[1] == "x")? 0:
1651  (tmpstr[1] == "y")? 1:2;
1652 
1653  LibUtilities::Interpreter strEval;
1654  int ExprId = strEval.DefineFunction("", tmpstr[2]);
1655  RotInfo.m_angle = strEval.Evaluate(ExprId);
1656 
1657  if(tmpstr.size() == 4)
1658  {
1659  try {
1660  RotInfo.m_tol = boost::lexical_cast
1661  <NekDouble>(tmpstr[3]);
1662  }
1663  catch (...) {
1665  "failed to cast tolerance input "
1666  "to a double value in Rotated"
1667  "boundary information");
1668  }
1669  }
1670  else
1671  {
1672  RotInfo.m_tol = 1e-8;
1673  }
1674  rotComp[cId1] = RotInfo;
1675  }
1676  }
1677 
1678  SpatialDomains::CompositeSharedPtr c = it.second->begin()->second;
1679 
1680  vector<unsigned int> tmpOrder;
1681 
1682  // store the rotation info of this
1683 
1684  // From the composite, we now construct the allVerts, allEdges
1685  // and allCoord map so that they can be transferred across
1686  // processors. We also populate the locFaces set to store a
1687  // record of all faces local to this process.
1688  for (i = 0; i < c->m_geomVec.size(); ++i)
1689  {
1691  std::dynamic_pointer_cast<
1692  SpatialDomains::Geometry2D>(c->m_geomVec[i]);
1693  ASSERTL1(faceGeom, "Unable to cast to shared ptr");
1694 
1695  // Get geometry ID of this face and store in locFaces.
1696  int faceId = c->m_geomVec[i]->GetGlobalID();
1697  locFaces.insert(faceId);
1698 
1699  // In serial, mesh partitioning will not have occurred so
1700  // need to fill composite ordering map manually.
1701  if (vComm->GetSize() == 1)
1702  {
1703  tmpOrder.push_back(c->m_geomVec[i]->GetGlobalID());
1704  }
1705 
1706  // Loop over vertices and edges of the face to populate
1707  // allVerts, allEdges and allCoord maps.
1708  vector<int> vertList, edgeList;
1710  vector<StdRegions::Orientation> orientVec;
1711  for (j = 0; j < faceGeom->GetNumVerts(); ++j)
1712  {
1713  vertList .push_back(faceGeom->GetVid (j));
1714  edgeList .push_back(faceGeom->GetEid (j));
1715  coordVec .push_back(faceGeom->GetVertex(j));
1716  orientVec.push_back(faceGeom->GetEorient(j));
1717  }
1718 
1719  allVerts [faceId] = vertList;
1720  allEdges [faceId] = edgeList;
1721  allCoord [faceId] = coordVec;
1722  allOrient[faceId] = orientVec;
1723  }
1724 
1725  // In serial, record the composite ordering in compOrder for
1726  // later in the routine.
1727  if (vComm->GetSize() == 1)
1728  {
1729  compOrder[it.second->begin()->first] = tmpOrder;
1730  }
1731 
1732  // See if we already have either region1 or region2 stored in
1733  // perComps map already and do a sanity check to ensure regions
1734  // are mutually periodic.
1735  if (perComps.count(cId1) == 0)
1736  {
1737  if (perComps.count(cId2) == 0)
1738  {
1739  perComps[cId1] = cId2;
1740  }
1741  else
1742  {
1743  std::stringstream ss;
1744  ss << "Boundary region " << cId2 << " should be "
1745  << "periodic with " << perComps[cId2] << " but "
1746  << "found " << cId1 << " instead!";
1747  ASSERTL0(perComps[cId2] == cId1, ss.str());
1748  }
1749  }
1750  else
1751  {
1752  std::stringstream ss;
1753  ss << "Boundary region " << cId1 << " should be "
1754  << "periodic with " << perComps[cId1] << " but "
1755  << "found " << cId2 << " instead!";
1756  ASSERTL0(perComps[cId1] == cId1, ss.str());
1757  }
1758  }
1759 
1760  // The next routines process local face lists to exchange vertices,
1761  // edges and faces.
1762  int n = vComm->GetSize();
1763  int p = vComm->GetRank();
1764  int totFaces;
1765  Array<OneD, int> facecounts(n,0);
1766  Array<OneD, int> vertcounts(n,0);
1767  Array<OneD, int> faceoffset(n,0);
1768  Array<OneD, int> vertoffset(n,0);
1769 
1770  Array<OneD, int> rotcounts(n,0);
1771  Array<OneD, int> rotoffset(n,0);
1772 
1773  rotcounts[p] = rotComp.size();
1774  vComm->AllReduce(rotcounts, LibUtilities::ReduceSum);
1775  int totrot = Vmath::Vsum(n,rotcounts,1);
1776 
1777  if(totrot)
1778  {
1779  for (i = 1; i < n ; ++i)
1780  {
1781  rotoffset[i] = rotoffset[i-1] + rotcounts[i-1];
1782  }
1783 
1784  Array<OneD, int> compid(totrot,0);
1785  Array<OneD, int> rotdir(totrot,0);
1786  Array<OneD, NekDouble> rotangle(totrot,0.0);
1787  Array<OneD, NekDouble> rottol(totrot,0.0);
1788 
1789  // fill in rotational informaiton
1790  auto rIt = rotComp.begin();
1791 
1792  for(i = 0; rIt != rotComp.end(); ++rIt)
1793  {
1794  compid [rotoffset[p] + i ] = rIt->first;
1795  rotdir [rotoffset[p] + i ] = rIt->second.m_dir;
1796  rotangle[rotoffset[p] + i ] = rIt->second.m_angle;
1797  rottol [rotoffset[p] + i++] = rIt->second.m_tol;
1798  }
1799 
1800  vComm->AllReduce(compid, LibUtilities::ReduceSum);
1801  vComm->AllReduce(rotdir, LibUtilities::ReduceSum);
1802  vComm->AllReduce(rotangle, LibUtilities::ReduceSum);
1803  vComm->AllReduce(rottol, LibUtilities::ReduceSum);
1804 
1805  // Fill in full rotational composite list
1806  for(i =0; i < totrot; ++i)
1807  {
1808  RotPeriodicInfo rinfo(rotdir[i],rotangle[i], rottol[i]);
1809 
1810  rotComp[compid[i]] = rinfo;
1811  }
1812  }
1813 
1814  // First exchange the number of faces on each process.
1815  facecounts[p] = locFaces.size();
1816  vComm->AllReduce(facecounts, LibUtilities::ReduceSum);
1817 
1818  // Set up an offset map to allow us to distribute face IDs to all
1819  // processors.
1820  faceoffset[0] = 0;
1821  for (i = 1; i < n; ++i)
1822  {
1823  faceoffset[i] = faceoffset[i-1] + facecounts[i-1];
1824  }
1825 
1826  // Calculate total number of faces.
1827  totFaces = Vmath::Vsum(n, facecounts, 1);
1828 
1829  // faceIds holds face IDs for each periodic face. faceVerts holds
1830  // the number of vertices in this face.
1831  Array<OneD, int> faceIds (totFaces, 0);
1832  Array<OneD, int> faceVerts(totFaces, 0);
1833 
1834  // Process p writes IDs of its faces into position faceoffset[p] of
1835  // faceIds which allows us to perform an AllReduce to distribute
1836  // information amongst processors.
1837  auto sIt = locFaces.begin();
1838  for (i = 0; sIt != locFaces.end(); ++sIt)
1839  {
1840  faceIds [faceoffset[p] + i ] = *sIt;
1841  faceVerts[faceoffset[p] + i++] = allVerts[*sIt].size();
1842  }
1843 
1844  vComm->AllReduce(faceIds, LibUtilities::ReduceSum);
1845  vComm->AllReduce(faceVerts, LibUtilities::ReduceSum);
1846 
1847  // procVerts holds number of vertices (and also edges since each
1848  // face is 2D) on each process.
1849  Array<OneD, int> procVerts(n,0);
1850  int nTotVerts;
1851 
1852  // Note if there are no periodic faces at all calling Vsum will
1853  // cause a segfault.
1854  if (totFaces > 0)
1855  {
1856  // Calculate number of vertices on each processor.
1857  nTotVerts = Vmath::Vsum(totFaces, faceVerts, 1);
1858  }
1859  else
1860  {
1861  nTotVerts = 0;
1862  }
1863 
1864  for (i = 0; i < n; ++i)
1865  {
1866  if (facecounts[i] > 0)
1867  {
1868  procVerts[i] = Vmath::Vsum
1869  (facecounts[i], faceVerts + faceoffset[i], 1);
1870  }
1871  else
1872  {
1873  procVerts[i] = 0;
1874  }
1875  }
1876 
1877  // vertoffset is defined in the same manner as edgeoffset
1878  // beforehand.
1879  vertoffset[0] = 0;
1880  for (i = 1; i < n; ++i)
1881  {
1882  vertoffset[i] = vertoffset[i-1] + procVerts[i-1];
1883  }
1884 
1885  // At this point we exchange all vertex IDs, edge IDs and vertex
1886  // coordinates for each face. The coordinates are necessary because
1887  // we need to calculate relative face orientations between periodic
1888  // faces to determined edge and vertex connectivity.
1889  Array<OneD, int> vertIds(nTotVerts, 0);
1890  Array<OneD, int> edgeIds(nTotVerts, 0);
1891  Array<OneD, int> edgeOrt(nTotVerts, 0);
1892  Array<OneD, NekDouble> vertX (nTotVerts, 0.0);
1893  Array<OneD, NekDouble> vertY (nTotVerts, 0.0);
1894  Array<OneD, NekDouble> vertZ (nTotVerts, 0.0);
1895 
1896  for (cnt = 0, sIt = locFaces.begin();
1897  sIt != locFaces.end(); ++sIt)
1898  {
1899  for (j = 0; j < allVerts[*sIt].size(); ++j)
1900  {
1901  int vertId = allVerts[*sIt][j];
1902  vertIds[vertoffset[p] + cnt ] = vertId;
1903  vertX [vertoffset[p] + cnt ] = (*allCoord[*sIt][j])(0);
1904  vertY [vertoffset[p] + cnt ] = (*allCoord[*sIt][j])(1);
1905  vertZ [vertoffset[p] + cnt ] = (*allCoord[*sIt][j])(2);
1906  edgeIds[vertoffset[p] + cnt ] = allEdges [*sIt][j];
1907  edgeOrt[vertoffset[p] + cnt++] = allOrient[*sIt][j];
1908  }
1909  }
1910 
1911  vComm->AllReduce(vertIds, LibUtilities::ReduceSum);
1912  vComm->AllReduce(vertX, LibUtilities::ReduceSum);
1913  vComm->AllReduce(vertY, LibUtilities::ReduceSum);
1914  vComm->AllReduce(vertZ, LibUtilities::ReduceSum);
1915  vComm->AllReduce(edgeIds, LibUtilities::ReduceSum);
1916  vComm->AllReduce(edgeOrt, LibUtilities::ReduceSum);
1917 
1918  // Finally now we have all of this information, we construct maps
1919  // which make accessing the information easier. These are
1920  // conceptually the same as all* maps at the beginning of the
1921  // routine, but now hold information for all periodic vertices.
1922  map<int, vector<int> > vertMap;
1923  map<int, vector<int> > edgeMap;
1924  map<int, SpatialDomains::PointGeomVector> coordMap;
1925 
1926  // These final two maps are required for determining the relative
1927  // orientation of periodic edges. vCoMap associates vertex IDs with
1928  // their coordinates, and eIdMap maps an edge ID to the two vertices
1929  // which construct it.
1930  map<int, SpatialDomains::PointGeomSharedPtr> vCoMap;
1931  map<int, pair<int, int> > eIdMap;
1932 
1933  for (cnt = i = 0; i < totFaces; ++i)
1934  {
1935  vector<int> edges(faceVerts[i]);
1936  vector<int> verts(faceVerts[i]);
1937  SpatialDomains::PointGeomVector coord(faceVerts[i]);
1938 
1939  // Keep track of cnt to enable correct edge vertices to be
1940  // inserted into eIdMap.
1941  int tmp = cnt;
1942  for (j = 0; j < faceVerts[i]; ++j, ++cnt)
1943  {
1944  edges[j] = edgeIds[cnt];
1945  verts[j] = vertIds[cnt];
1947  ::AllocateSharedPtr(3, verts[j], vertX[cnt],
1948  vertY[cnt], vertZ[cnt]);
1949  vCoMap[vertIds[cnt]] = coord[j];
1950 
1951  // Try to insert edge into the eIdMap to avoid re-inserting.
1952  auto testIns = eIdMap.insert
1953  ( make_pair( edgeIds[cnt], make_pair(vertIds[tmp+j],
1954  vertIds[tmp+((j+1) % faceVerts[i])])));
1955 
1956  if (testIns.second == false)
1957  {
1958  continue;
1959  }
1960 
1961  // If the edge is reversed with respect to the face, then
1962  // swap the edges so that we have the original ordering of
1963  // the edge in the 3D element. This is necessary to properly
1964  // determine edge orientation. Note that the logic relies on
1965  // the fact that all edge forward directions are CCW
1966  // orientated: we use a tensor product ordering for 2D
1967  // elements so need to reverse this for edge IDs 2 and 3.
1968  StdRegions::Orientation edgeOrient =
1969  static_cast<StdRegions::Orientation>(edgeOrt[cnt]);
1970  if (j > 1)
1971  {
1972  edgeOrient = edgeOrient == StdRegions::eForwards ?
1974  }
1975 
1976  if (edgeOrient == StdRegions::eBackwards)
1977  {
1978  swap(testIns.first->second.first,
1979  testIns.first->second.second);
1980  }
1981  }
1982 
1983  vertMap [faceIds[i]] = verts;
1984  edgeMap [faceIds[i]] = edges;
1985  coordMap[faceIds[i]] = coord;
1986  }
1987 
1988  // Go through list of composites and figure out which edges are
1989  // parallel from original ordering in session file. This includes
1990  // composites which are not necessarily on this process.
1991 
1992  // Store temporary map of periodic vertices which will hold all
1993  // periodic vertices on the entire mesh so that doubly periodic
1994  // vertices/edges can be counted properly across partitions. Local
1995  // vertices/edges are copied into m_periodicVerts and
1996  // m_periodicEdges at the end of the function.
1997  PeriodicMap periodicVerts, periodicEdges;
1998 
1999  // Construct two maps which determine how vertices and edges of
2000  // faces connect given a specific face orientation. The key of the
2001  // map is the number of vertices in the face, used to determine
2002  // difference between tris and quads.
2003  map<int, map<StdRegions::Orientation, vector<int> > > vmap;
2004  map<int, map<StdRegions::Orientation, vector<int> > > emap;
2005 
2006  map<StdRegions::Orientation, vector<int> > quadVertMap;
2007  quadVertMap[StdRegions::eDir1FwdDir1_Dir2FwdDir2] = {0,1,2,3};
2008  quadVertMap[StdRegions::eDir1FwdDir1_Dir2BwdDir2] = {3,2,1,0};
2009  quadVertMap[StdRegions::eDir1BwdDir1_Dir2FwdDir2] = {1,0,3,2};
2010  quadVertMap[StdRegions::eDir1BwdDir1_Dir2BwdDir2] = {2,3,0,1};
2011  quadVertMap[StdRegions::eDir1FwdDir2_Dir2FwdDir1] = {0,3,2,1};
2012  quadVertMap[StdRegions::eDir1FwdDir2_Dir2BwdDir1] = {1,2,3,0};
2013  quadVertMap[StdRegions::eDir1BwdDir2_Dir2FwdDir1] = {3,0,1,2};
2014  quadVertMap[StdRegions::eDir1BwdDir2_Dir2BwdDir1] = {2,1,0,3};
2015 
2016  map<StdRegions::Orientation, vector<int> > quadEdgeMap;
2017  quadEdgeMap[StdRegions::eDir1FwdDir1_Dir2FwdDir2] = {0,1,2,3};
2018  quadEdgeMap[StdRegions::eDir1FwdDir1_Dir2BwdDir2] = {2,1,0,3};
2019  quadEdgeMap[StdRegions::eDir1BwdDir1_Dir2FwdDir2] = {0,3,2,1};
2020  quadEdgeMap[StdRegions::eDir1BwdDir1_Dir2BwdDir2] = {2,3,0,1};
2021  quadEdgeMap[StdRegions::eDir1FwdDir2_Dir2FwdDir1] = {3,2,1,0};
2022  quadEdgeMap[StdRegions::eDir1FwdDir2_Dir2BwdDir1] = {1,2,3,0};
2023  quadEdgeMap[StdRegions::eDir1BwdDir2_Dir2FwdDir1] = {3,0,1,2};
2024  quadEdgeMap[StdRegions::eDir1BwdDir2_Dir2BwdDir1] = {1,0,3,2};
2025 
2026  map<StdRegions::Orientation, vector<int> > triVertMap;
2027  triVertMap[StdRegions::eDir1FwdDir1_Dir2FwdDir2] = {0,1,2};
2028  triVertMap[StdRegions::eDir1BwdDir1_Dir2FwdDir2] = {1,0,2};
2029 
2030  map<StdRegions::Orientation, vector<int> > triEdgeMap;
2031  triEdgeMap[StdRegions::eDir1FwdDir1_Dir2FwdDir2] = {0,1,2};
2032  triEdgeMap[StdRegions::eDir1BwdDir1_Dir2FwdDir2] = {0,2,1};
2033 
2034  vmap[3] = triVertMap;
2035  vmap[4] = quadVertMap;
2036  emap[3] = triEdgeMap;
2037  emap[4] = quadEdgeMap;
2038 
2039  map<int,int> allCompPairs;
2040 
2041  // Collect composite id's of each periodic face for use if rotation is required
2042  map<int,int> fIdToCompId;
2043 
2044  // Finally we have enough information to populate the periodic
2045  // vertex, edge and face maps. Begin by looping over all pairs of
2046  // periodic composites to determine pairs of periodic faces.
2047  for (auto &cIt : perComps)
2048  {
2050  const int id1 = cIt.first;
2051  const int id2 = cIt.second;
2052  std::string id1s = boost::lexical_cast<string>(id1);
2053  std::string id2s = boost::lexical_cast<string>(id2);
2054 
2055  if (compMap.count(id1) > 0)
2056  {
2057  c[0] = compMap[id1];
2058  }
2059 
2060  if (compMap.count(id2) > 0)
2061  {
2062  c[1] = compMap[id2];
2063  }
2064 
2065  ASSERTL0(c[0] || c[1],
2066  "Neither composite not found on this process!");
2067 
2068  // Loop over composite ordering to construct list of all
2069  // periodic faces, regardless of whether they are on this
2070  // process.
2071  map<int,int> compPairs;
2072 
2073 
2074  ASSERTL0(compOrder.count(id1) > 0,
2075  "Unable to find composite "+id1s+" in order map.");
2076  ASSERTL0(compOrder.count(id2) > 0,
2077  "Unable to find composite "+id2s+" in order map.");
2078  ASSERTL0(compOrder[id1].size() == compOrder[id2].size(),
2079  "Periodic composites "+id1s+" and "+id2s+
2080  " should have the same number of elements.");
2081  ASSERTL0(compOrder[id1].size() > 0,
2082  "Periodic composites "+id1s+" and "+id2s+
2083  " are empty!");
2084 
2085  // Look up composite ordering to determine pairs.
2086  for (i = 0; i < compOrder[id1].size(); ++i)
2087  {
2088  int eId1 = compOrder[id1][i];
2089  int eId2 = compOrder[id2][i];
2090 
2091  ASSERTL0(compPairs.count(eId1) == 0,
2092  "Already paired.");
2093 
2094  // Sanity check that the faces are mutually periodic.
2095  if (compPairs.count(eId2) != 0)
2096  {
2097  ASSERTL0(compPairs[eId2] == eId1, "Pairing incorrect");
2098  }
2099  compPairs[eId1] = eId2;
2100 
2101  // store a map of face ids to composite ids
2102  fIdToCompId[eId1] = id1;
2103  fIdToCompId[eId2] = id2;
2104  }
2105 
2106  // Now that we have all pairs of periodic faces, loop over the
2107  // ones local on this process and populate face/edge/vertex
2108  // maps.
2109  for (auto &pIt : compPairs)
2110  {
2111  int ids [2] = {pIt.first, pIt.second};
2112  bool local[2] = {locFaces.count(pIt.first) > 0,
2113  locFaces.count(pIt.second) > 0};
2114 
2115  ASSERTL0(coordMap.count(ids[0]) > 0 &&
2116  coordMap.count(ids[1]) > 0,
2117  "Unable to find face in coordinate map");
2118 
2119  allCompPairs[pIt.first ] = pIt.second;
2120  allCompPairs[pIt.second] = pIt.first;
2121 
2122  // Loop up coordinates of the faces, check they have the
2123  // same number of vertices.
2125  = { coordMap[ids[0]], coordMap[ids[1]] };
2126 
2127  ASSERTL0(tmpVec[0].size() == tmpVec[1].size(),
2128  "Two periodic faces have different number "
2129  "of vertices!");
2130 
2131  // o will store relative orientation of faces. Note that in
2132  // some transpose cases (Dir1FwdDir2_Dir2BwdDir1 and
2133  // Dir1BwdDir1_Dir2FwdDir1) it seems orientation will be
2134  // different going from face1->face2 instead of face2->face1
2135  // (check this).
2137  bool rotbnd = false;
2138  int dir = 0;
2139  NekDouble angle = 0.0;
2140  NekDouble sign = 0.0;
2141  NekDouble tol = 1e-8;
2142 
2143  // check to see if perioid boundary is rotated
2144  if(rotComp.count(fIdToCompId[pIt.first]))
2145  {
2146  rotbnd = true;
2147  dir = rotComp[fIdToCompId[pIt.first]].m_dir;
2148  angle = rotComp[fIdToCompId[pIt.first]].m_angle;
2149  tol = rotComp[fIdToCompId[pIt.first]].m_tol;
2150  }
2151 
2152  // Record periodic faces.
2153  for (i = 0; i < 2; ++i)
2154  {
2155  if (!local[i])
2156  {
2157  continue;
2158  }
2159 
2160  // Reference to the other face.
2161  int other = (i+1) % 2;
2162 
2163  // angle is set up for i = 0 to i = 1
2164  sign = (i == 0)? 1.0:-1.0;
2165 
2166  // Calculate relative face orientation.
2167  if (tmpVec[0].size() == 3)
2168  {
2170  (tmpVec[i], tmpVec[other],
2171  rotbnd, dir, sign*angle, tol);
2172  }
2173  else
2174  {
2176  tmpVec[i], tmpVec[other],
2177  rotbnd,dir,sign*angle,tol);
2178  }
2179 
2180  // Record face ID, orientation and whether other face is
2181  // local.
2182  PeriodicEntity ent(ids [other], o,
2183  local[other]);
2184  m_periodicFaces[ids[i]].push_back(ent);
2185  }
2186 
2187  int nFaceVerts = vertMap[ids[0]].size();
2188 
2189  // Determine periodic vertices.
2190  for (i = 0; i < 2; ++i)
2191  {
2192  int other = (i+1) % 2;
2193 
2194  // angle is set up for i = 0 to i = 1
2195  sign = (i == 0)? 1.0:-1.0;
2196 
2197  // Calculate relative face orientation.
2198  if (tmpVec[0].size() == 3)
2199  {
2201  tmpVec[i], tmpVec[other], rotbnd, dir,
2202  sign*angle, tol);
2203  }
2204  else
2205  {
2207  tmpVec[i], tmpVec[other], rotbnd, dir,
2208  sign*angle, tol);
2209  }
2210 
2211  if (nFaceVerts == 3)
2212  {
2215  "Unsupported face orientation for face "+
2216  boost::lexical_cast<string>(ids[i]));
2217  }
2218 
2219  // Look up vertices for this face.
2220  vector<int> per1 = vertMap[ids[i]];
2221  vector<int> per2 = vertMap[ids[other]];
2222 
2223  // tmpMap will hold the pairs of vertices which are
2224  // periodic.
2225  map<int, pair<int, bool> > tmpMap;
2226 
2227  // Use vmap to determine which vertices connect given
2228  // the orientation o.
2229  for (j = 0; j < nFaceVerts; ++j)
2230  {
2231  int v = vmap[nFaceVerts][o][j];
2232  tmpMap[per1[j]] = make_pair
2233  (per2[v], locVerts.count(per2[v]) > 0);
2234  }
2235 
2236  // Now loop over tmpMap to associate periodic vertices.
2237  for (auto &mIt : tmpMap)
2238  {
2239  PeriodicEntity ent2(mIt.second.first,
2241  mIt.second.second);
2242 
2243  // See if this vertex has been recorded already.
2244  auto perIt = periodicVerts.find(mIt.first);
2245 
2246  if (perIt == periodicVerts.end())
2247  {
2248  // Vertex is new - just record this entity as
2249  // usual.
2250  periodicVerts[mIt.first].push_back(ent2);
2251  perIt = periodicVerts.find(mIt.first);
2252  }
2253  else
2254  {
2255  // Vertex is known - loop over the vertices
2256  // inside the record and potentially add vertex
2257  // mIt.second to the list.
2258  for (k = 0; k < perIt->second.size(); ++k)
2259  {
2260  if (perIt->second[k].id == mIt.second.first)
2261  {
2262  break;
2263  }
2264  }
2265 
2266  if (k == perIt->second.size())
2267  {
2268  perIt->second.push_back(ent2);
2269  }
2270  }
2271  }
2272  }
2273 
2274  // Determine periodic edges. Logic is the same as above,
2275  // and perhaps should be condensed to avoid replication.
2276  for (i = 0; i < 2; ++i)
2277  {
2278  int other = (i+1) % 2;
2279 
2280  // angle is set up for i = 0 to i = 1
2281  sign = (i == 0)? 1.0:-1.0;
2282 
2283  if (tmpVec[0].size() == 3)
2284  {
2286  tmpVec[i], tmpVec[other], rotbnd, dir,
2287  sign*angle, tol);
2288  }
2289  else
2290  {
2292  tmpVec[i], tmpVec[other], rotbnd, dir,
2293  sign*angle, tol);
2294  }
2295 
2296  vector<int> per1 = edgeMap[ids[i]];
2297  vector<int> per2 = edgeMap[ids[other]];
2298 
2299  map<int, pair<int, bool> > tmpMap;
2300 
2301  for (j = 0; j < nFaceVerts; ++j)
2302  {
2303  int e = emap[nFaceVerts][o][j];
2304  tmpMap[per1[j]] = make_pair
2305  (per2[e], locEdges.count(per2[e]) > 0);
2306  }
2307 
2308  for (auto &mIt : tmpMap)
2309  {
2310  // Note we assume orientation of edges is forwards -
2311  // this may be reversed later.
2312  PeriodicEntity ent2(mIt.second.first,
2314  mIt.second.second);
2315  auto perIt = periodicEdges.find(mIt.first);
2316 
2317  if (perIt == periodicEdges.end())
2318  {
2319  periodicEdges[mIt.first].push_back(ent2);
2320  perIt = periodicEdges.find(mIt.first);
2321  }
2322  else
2323  {
2324  for (k = 0; k < perIt->second.size(); ++k)
2325  {
2326  if (perIt->second[k].id == mIt.second.first)
2327  {
2328  break;
2329  }
2330  }
2331 
2332  if (k == perIt->second.size())
2333  {
2334  perIt->second.push_back(ent2);
2335  }
2336  }
2337  }
2338  }
2339  }
2340  }
2341 
2342  // Make sure that the nubmer of face pairs and the
2343  // face Id to composite Id map match in size
2344  ASSERTL1(allCompPairs.size() == fIdToCompId.size(),
2345  "At this point the size of allCompPairs "
2346  "should have been the same as fIdToCompId");
2347 
2348  // Distribute the size of the periodic boundary to all
2349  // processes. We assume that all processes who own periodic
2350  // faces have the same local copy of allCompPairs
2351  int NPairs = allCompPairs.size();
2352  vComm->AllReduce(NPairs, LibUtilities::ReduceMax);
2353 
2354  // Check that the previous assertion regarding allCompPairs
2355  // is correct
2356  ASSERTL0(
2357  allCompPairs.size() == NPairs || allCompPairs.size() == 0,
2358  "Local copy of allCompPairs not the same for all ranks"
2359  );
2360 
2361  // Allocate local vectors that will contain the content of
2362  // allCompPairs if process owns faces on periodic boundary
2363  Array<OneD, int> first(NPairs, -1);
2364  Array<OneD, int> second(NPairs, -1);
2365  cnt = 0;
2366  for(const auto &it : allCompPairs)
2367  {
2368  first[cnt] = it.first;
2369  second[cnt++] = it.second;
2370  }
2371 
2372  // Distribute the content in first and second to all processes
2373  vComm->AllReduce(first, LibUtilities::ReduceMax);
2374  vComm->AllReduce(second, LibUtilities::ReduceMax);
2375 
2376  // Check that the MPI Allreduce routine worked
2377  ASSERTL0(std::count(first.begin(), first.end(), -1) == 0,
2378  "Distribution of allCompPairs failed");
2379  ASSERTL0(std::count(second.begin(), second.end(), -1) == 0,
2380  "Distribution of allCompPairs failed")
2381 
2382  // Put content back in allCompPairs
2383  allCompPairs.clear();
2384  for(cnt = 0; cnt < NPairs; ++cnt)
2385  {
2386  allCompPairs[first[cnt]] = second[cnt];
2387  }
2388 
2389  // Store face ID to composite ID map for rotational boundaries
2390  if(rotComp.size())
2391  {
2392  // Set values to -1
2393  std::fill(first.begin(), first.end(), -1);
2394  std::fill(second.begin(), second.end(), -1);
2395 
2396  cnt = 0;
2397  for (const auto &pIt : fIdToCompId)
2398  {
2399  first [cnt ] = pIt.first;
2400  second[cnt++] = pIt.second;
2401  }
2402 
2403  vComm->AllReduce(first, LibUtilities::ReduceMax);
2404  vComm->AllReduce(second, LibUtilities::ReduceMax);
2405 
2406  // Check that the MPI Allreduce routine worked
2407  ASSERTL0(std::count(first.begin(), first.end(), -1) == 0,
2408  "Distribution of fIdToCompId failed");
2409  ASSERTL0(std::count(second.begin(), second.end(), -1) == 0,
2410  "Distribution of fIdToCompId failed")
2411 
2412  fIdToCompId.clear();
2413  for(cnt = 0; cnt < NPairs; ++cnt)
2414  {
2415  fIdToCompId[first[cnt]] = second[cnt];
2416  }
2417  }
2418 
2419  // also will need an edge id to composite id at end of routine
2420  map<int,int> eIdToCompId;
2421 
2422  // Search for periodic vertices and edges which are not
2423  // in a periodic composite but lie in this process. First,
2424  // loop over all information we have from other
2425  // processors.
2426  for (cnt = i = 0; i < totFaces; ++i)
2427  {
2428  bool rotbnd = false;
2429  int dir = 0;
2430  NekDouble angle = 0.0;
2431  NekDouble tol = 1e-8;
2432 
2433  int faceId = faceIds[i];
2434 
2435  ASSERTL0(allCompPairs.count(faceId) > 0,
2436  "Unable to find matching periodic face.");
2437 
2438  int perFaceId = allCompPairs[faceId];
2439 
2440  // check to see if periodic boundary is rotated
2441  ASSERTL1((rotComp.size() == 0) ||
2442  fIdToCompId.count(faceId) > 0,"Face " +
2443  boost::lexical_cast<string>(faceId) +
2444  " not found in fIdtoCompId map");
2445  if(rotComp.count(fIdToCompId[faceId]))
2446  {
2447  rotbnd = true;
2448  dir = rotComp[fIdToCompId[faceId]].m_dir;
2449  angle = rotComp[fIdToCompId[faceId]].m_angle;
2450  tol = rotComp[fIdToCompId[faceId]].m_tol;
2451  }
2452 
2453  for (j = 0; j < faceVerts[i]; ++j, ++cnt)
2454  {
2455  int vId = vertIds[cnt];
2456 
2457  auto perId = periodicVerts.find(vId);
2458 
2459  if (perId == periodicVerts.end())
2460  {
2461 
2462  // This vertex is not included in the
2463  // map. Figure out which vertex it is supposed
2464  // to be periodic with. perFaceId is the face
2465  // ID which is periodic with faceId. The logic
2466  // is much the same as the loop above.
2468  = { coordMap[faceId], coordMap[perFaceId] };
2469 
2470  int nFaceVerts = tmpVec[0].size();
2471  StdRegions::Orientation o = nFaceVerts == 3 ?
2473  tmpVec[0], tmpVec[1], rotbnd, dir, angle, tol):
2474  SpatialDomains::QuadGeom::GetFaceOrientation(
2475  tmpVec[0], tmpVec[1], rotbnd, dir, angle, tol);
2476 
2477  // Use vmap to determine which vertex of the other face
2478  // should be periodic with this one.
2479  int perVertexId = vertMap[perFaceId][vmap[nFaceVerts][o][j]];
2480 
2481 
2482  PeriodicEntity ent(perVertexId,
2484  locVerts.count(perVertexId) > 0);
2485 
2486  periodicVerts[vId].push_back(ent);
2487  }
2488 
2489  int eId = edgeIds[cnt];
2490 
2491  perId = periodicEdges.find(eId);
2492 
2493  // this map is required at very end to determine rotation of edges.
2494  if(rotbnd)
2495  {
2496  eIdToCompId[eId] = fIdToCompId[faceId];
2497  }
2498 
2499  if (perId == periodicEdges.end())
2500  {
2501  // This edge is not included in the map. Figure
2502  // out which edge it is supposed to be periodic
2503  // with. perFaceId is the face ID which is
2504  // periodic with faceId. The logic is much the
2505  // same as the loop above.
2507  = { coordMap[faceId], coordMap[perFaceId] };
2508 
2509  int nFaceEdges = tmpVec[0].size();
2510  StdRegions::Orientation o = nFaceEdges == 3 ?
2512  tmpVec[0], tmpVec[1], rotbnd, dir, angle, tol):
2513  SpatialDomains::QuadGeom::GetFaceOrientation(
2514  tmpVec[0], tmpVec[1], rotbnd, dir, angle, tol);
2515 
2516  // Use emap to determine which edge of the other
2517  // face should be periodic with this one.
2518  int perEdgeId = edgeMap[perFaceId][emap[nFaceEdges][o][j]];
2519 
2520  PeriodicEntity ent(perEdgeId,
2522  locEdges.count(perEdgeId) > 0);
2523 
2524  periodicEdges[eId].push_back(ent);
2525 
2526 
2527  // this map is required at very end to
2528  // determine rotation of edges.
2529  if(rotbnd)
2530  {
2531  eIdToCompId[perEdgeId] = fIdToCompId[perFaceId];
2532  }
2533  }
2534  }
2535  }
2536 
2537  // Finally, we must loop over the periodicVerts and periodicEdges
2538  // map to complete connectivity information.
2539  for (auto &perIt : periodicVerts)
2540  {
2541  // For each vertex that is periodic with this one...
2542  for (i = 0; i < perIt.second.size(); ++i)
2543  {
2544  // Find the vertex in the periodicVerts map...
2545  auto perIt2 = periodicVerts.find(perIt.second[i].id);
2546  ASSERTL0(perIt2 != periodicVerts.end(),
2547  "Couldn't find periodic vertex.");
2548 
2549  // Now search through this vertex's list and make sure that
2550  // we have a record of any vertices which aren't in the
2551  // original list.
2552  for (j = 0; j < perIt2->second.size(); ++j)
2553  {
2554  if (perIt2->second[j].id == perIt.first)
2555  {
2556  continue;
2557  }
2558 
2559  for (k = 0; k < perIt.second.size(); ++k)
2560  {
2561  if (perIt2->second[j].id == perIt.second[k].id)
2562  {
2563  break;
2564  }
2565  }
2566 
2567  if (k == perIt.second.size())
2568  {
2569  perIt.second.push_back(perIt2->second[j]);
2570  }
2571  }
2572  }
2573  }
2574 
2575  for (auto &perIt : periodicEdges)
2576  {
2577  for (i = 0; i < perIt.second.size(); ++i)
2578  {
2579  auto perIt2 = periodicEdges.find(perIt.second[i].id);
2580  ASSERTL0(perIt2 != periodicEdges.end(),
2581  "Couldn't find periodic edge.");
2582 
2583  for (j = 0; j < perIt2->second.size(); ++j)
2584  {
2585  if (perIt2->second[j].id == perIt.first)
2586  {
2587  continue;
2588  }
2589 
2590  for (k = 0; k < perIt.second.size(); ++k)
2591  {
2592  if (perIt2->second[j].id == perIt.second[k].id)
2593  {
2594  break;
2595  }
2596  }
2597 
2598  if (k == perIt.second.size())
2599  {
2600  perIt.second.push_back(perIt2->second[j]);
2601  }
2602  }
2603  }
2604  }
2605 
2606  // Loop over periodic edges to determine relative edge orientations.
2607  for (auto &perIt : periodicEdges)
2608  {
2609  bool rotbnd = false;
2610  int dir = 0;
2611  NekDouble angle = 0.0;
2612  NekDouble tol = 1e-8;
2613 
2614 
2615  // Find edge coordinates
2616  auto eIt = eIdMap.find(perIt.first);
2617  SpatialDomains::PointGeom v[2] = {
2618  *vCoMap[eIt->second.first],
2619  *vCoMap[eIt->second.second]
2620  };
2621 
2622  // check to see if perioid boundary is rotated
2623  if(rotComp.count(eIdToCompId[perIt.first]))
2624  {
2625  rotbnd = true;
2626  dir = rotComp[eIdToCompId[perIt.first]].m_dir;
2627  angle = rotComp[eIdToCompId[perIt.first]].m_angle;
2628  tol = rotComp[eIdToCompId[perIt.first]].m_tol;
2629  }
2630 
2631  // Loop over each edge, and construct a vector that takes us
2632  // from one vertex to another. Use this to figure out which
2633  // vertex maps to which.
2634  for (i = 0; i < perIt.second.size(); ++i)
2635  {
2636  eIt = eIdMap.find(perIt.second[i].id);
2637 
2638  SpatialDomains::PointGeom w[2] = {
2639  *vCoMap[eIt->second.first],
2640  *vCoMap[eIt->second.second]
2641  };
2642 
2643  int vMap[2] = {-1,-1};
2644  if(rotbnd)
2645  {
2646 
2647  SpatialDomains::PointGeom r;
2648 
2649  r.Rotate(v[0],dir,angle);
2650 
2651  if(r.dist(w[0])< tol)
2652  {
2653  vMap[0] = 0;
2654  }
2655  else
2656  {
2657  r.Rotate(v[1],dir,angle);
2658  if(r.dist(w[0]) < tol)
2659  {
2660  vMap[0] = 1;
2661  }
2662  else
2663  {
2665  "Unable to align rotationally "
2666  "periodic edge vertex");
2667  }
2668  }
2669  }
2670  else // translation test
2671  {
2672  NekDouble cx = 0.5*(w[0](0)-v[0](0)+w[1](0)-v[1](0));
2673  NekDouble cy = 0.5*(w[0](1)-v[0](1)+w[1](1)-v[1](1));
2674  NekDouble cz = 0.5*(w[0](2)-v[0](2)+w[1](2)-v[1](2));
2675 
2676  for (j = 0; j < 2; ++j)
2677  {
2678  NekDouble x = v[j](0);
2679  NekDouble y = v[j](1);
2680  NekDouble z = v[j](2);
2681  for (k = 0; k < 2; ++k)
2682  {
2683  NekDouble x1 = w[k](0)-cx;
2684  NekDouble y1 = w[k](1)-cy;
2685  NekDouble z1 = w[k](2)-cz;
2686 
2687  if (sqrt((x1-x)*(x1-x)+(y1-y)*(y1-y)+(z1-z)*(z1-z))
2688  < 1e-8)
2689  {
2690  vMap[k] = j;
2691  break;
2692  }
2693  }
2694  }
2695 
2696  // Sanity check the map.
2697  ASSERTL0(vMap[0] >= 0 && vMap[1] >= 0,
2698  "Unable to align periodic edge vertex.");
2699  ASSERTL0((vMap[0] == 0 || vMap[0] == 1) &&
2700  (vMap[1] == 0 || vMap[1] == 1) &&
2701  (vMap[0] != vMap[1]),
2702  "Unable to align periodic edge vertex.");
2703  }
2704 
2705  // If 0 -> 0 then edges are aligned already; otherwise
2706  // reverse the orientation.
2707  if (vMap[0] != 0)
2708  {
2709  perIt.second[i].orient = StdRegions::eBackwards;
2710  }
2711  }
2712  }
2713 
2714  // Do one final loop over periodic vertices/edges to remove
2715  // non-local vertices/edges from map.
2716  for (auto &perIt : periodicVerts)
2717  {
2718  if (locVerts.count(perIt.first) > 0)
2719  {
2720  m_periodicVerts.insert(perIt);
2721  }
2722  }
2723 
2724  for (auto &perIt : periodicEdges)
2725  {
2726  if (locEdges.count(perIt.first) > 0)
2727  {
2728  m_periodicEdges.insert(perIt);
2729  }
2730  }
2731  }
2732  break;
2733  default:
2734  ASSERTL1(false,"Not setup for this expansion");
2735  break;
2736  }
2737  }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:250
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:15
static SpatialDomains::BoundaryConditionShPtr GetBoundaryCondition(const SpatialDomains::BoundaryConditionCollection &collection, unsigned int index, const std::string &variable)
Definition: ExpList.cpp:4983
static StdRegions::Orientation GetFaceOrientation(const QuadGeom &face1, const QuadGeom &face2, bool doRot=false, int dir=0, NekDouble angle=0.0, NekDouble tol=1e-8)
Get the orientation of face1.
Definition: QuadGeom.cpp:138
static StdRegions::Orientation GetFaceOrientation(const TriGeom &face1, const TriGeom &face2, bool doRot, int dir, NekDouble angle, NekDouble tol)
Definition: TriGeom.cpp:118
std::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:54
std::map< int, std::vector< PeriodicEntity > > PeriodicMap
std::map< int, BoundaryRegionShPtr > BoundaryRegionCollection
Definition: Conditions.h:217
std::map< int, std::vector< unsigned int > > CompositeOrdering
Definition: MeshGraph.h:109
std::shared_ptr< BoundaryConditionBase > BoundaryConditionShPtr
Definition: Conditions.h:219
std::shared_ptr< std::vector< std::pair< GeometrySharedPtr, int > > > GeometryLinkSharedPtr
Definition: MeshGraph.h:170
std::map< int, std::vector< unsigned int > > BndRegionOrdering
Definition: MeshGraph.h:110
std::vector< PointGeomSharedPtr > PointGeomVector
Definition: Geometry2D.h:64
std::map< int, CompositeSharedPtr > CompositeMap
Definition: MeshGraph.h:138
std::map< int, BoundaryConditionMapShPtr > BoundaryConditionCollection
Definition: Conditions.h:226
std::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:62
std::shared_ptr< Geometry2D > Geometry2DSharedPtr
Definition: Geometry.h:65
std::shared_ptr< Composite > CompositeSharedPtr
Definition: MeshGraph.h:137
double NekDouble
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.cpp:846
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:267

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, ASSERTL1, Nektar::LibUtilities::Interpreter::DefineFunction(), Nektar::SpatialDomains::PointGeom::dist(), Nektar::MultiRegions::e1D, Nektar::MultiRegions::e2D, Nektar::MultiRegions::e3D, Nektar::StdRegions::eBackwards, 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, Nektar::ErrorUtil::efatal, Nektar::StdRegions::eForwards, Nektar::StdRegions::eNoOrientation, Nektar::SpatialDomains::ePeriodic, Nektar::LibUtilities::Interpreter::Evaluate(), Nektar::MultiRegions::ExpList::GetBoundaryCondition(), Nektar::SpatialDomains::BoundaryConditions::GetBoundaryConditions(), Nektar::SpatialDomains::BoundaryConditions::GetBoundaryRegions(), Nektar::SpatialDomains::QuadGeom::GetFaceOrientation(), Nektar::SpatialDomains::TriGeom::GetFaceOrientation(), Nektar::MultiRegions::RotPeriodicInfo::m_angle, Nektar::MultiRegions::RotPeriodicInfo::m_dir, Nektar::MultiRegions::ExpList::m_expType, Nektar::MultiRegions::ExpList::m_graph, m_periodicEdges, m_periodicFaces, m_periodicVerts, Nektar::MultiRegions::ExpList::m_session, Nektar::MultiRegions::RotPeriodicInfo::m_tol, NEKERROR, CellMLToNektar.cellml_metadata::p, Nektar::LibUtilities::ReduceMax, Nektar::LibUtilities::ReduceSum, Nektar::SpatialDomains::PointGeom::Rotate(), sign, tinysimd::sqrt(), and Vmath::Vsum().

Referenced by DisContField().

◆ GenerateBoundaryConditionExpansion()

void Nektar::MultiRegions::DisContField::GenerateBoundaryConditionExpansion ( const SpatialDomains::MeshGraphSharedPtr graph,
const SpatialDomains::BoundaryConditions bcs,
const std::string  variable,
const bool  DeclareCoeffPhysArrays = true 
)
protected

Discretises the boundary conditions.

This function discretises the boundary conditions by setting up a list of one-dimensions lower boundary expansions.

According to their boundary region, the separate boundary expansions are bundled together in an object of the class

Parameters
graphA mesh, containing information about the domain and the spectral/hp element expansions.
bcsInformation about the enforced boundary conditions.
variableThe session variable associated with the boundary conditions to enforce.
DeclareCoeffPhysArraysbool to identify if array space should be setup. Default is true.

Definition at line 796 of file DisContField.cpp.

801  {
802  int cnt = 0;
806  = bcs.GetBoundaryRegions();
808  = bcs.GetBoundaryConditions();
809 
811  Array<OneD, MultiRegions::ExpListSharedPtr>(bregions.size());
813  Array<OneD, SpatialDomains::BoundaryConditionShPtr>(bregions.size());
814 
815  m_bndCondBndWeight = Array<OneD, NekDouble> {bregions.size(), 0.0};
816 
817  // count the number of non-periodic boundary points
818  for (auto &it : bregions)
819  {
820  bc = GetBoundaryCondition(bconditions, it.first, variable);
821 
823  ::AllocateSharedPtr(m_session, *(it.second), graph,
824  DeclareCoeffPhysArrays, variable,
825  false, bc->GetComm());
826 
827  m_bndCondExpansions[cnt] = locExpList;
828  m_bndConditions[cnt] = bc;
829 
830  std::string type = m_bndConditions[cnt]->GetUserDefined();
831 
832  // Set up normals on non-Dirichlet boundary conditions. Second
833  // two conditions ideally should be in local solver setup (when
834  // made into factory)
835  if(bc->GetBoundaryConditionType() != SpatialDomains::eDirichlet
836  || boost::iequals(type,"I") || boost::iequals(type,"CalcBC"))
837  {
839  }
840  cnt ++;
841  }
842  }
Array< OneD, NekDouble > m_bndCondBndWeight
Definition: DisContField.h:144

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::SpatialDomains::eDirichlet, Nektar::MultiRegions::ExpList::GetBoundaryCondition(), Nektar::SpatialDomains::BoundaryConditions::GetBoundaryConditions(), Nektar::SpatialDomains::BoundaryConditions::GetBoundaryRegions(), m_bndCondBndWeight, m_bndCondExpansions, m_bndConditions, Nektar::MultiRegions::ExpList::m_session, and Nektar::MultiRegions::ExpList::SetUpPhysNormals().

Referenced by DisContField().

◆ GenerateFieldBnd1D()

void Nektar::MultiRegions::DisContField::GenerateFieldBnd1D ( SpatialDomains::BoundaryConditions bcs,
const std::string  variable 
)
protected

◆ GetDomainBCs()

SpatialDomains::BoundaryConditionsSharedPtr Nektar::MultiRegions::DisContField::GetDomainBCs ( const SpatialDomains::CompositeMap domain,
const SpatialDomains::BoundaryConditions Allbcs,
const std::string &  variable 
)
private

Definition at line 470 of file DisContField.cpp.

474  {
476 
478 
479  map<int,int> GeometryToRegionsMap;
480 
482  = Allbcs.GetBoundaryRegions();
484  = Allbcs.GetBoundaryConditions();
485 
486  // Set up a map of all boundary regions
487  for(auto &it : bregions)
488  {
489  for (auto &bregionIt : *it.second)
490  {
491  // can assume that all regions only contain one point in 1D
492  // Really do not need loop above
493  int id = bregionIt.second->m_geomVec[0]->GetGlobalID();
494  GeometryToRegionsMap[id] = it.first;
495  }
496  }
497 
498  map<int,SpatialDomains::GeometrySharedPtr> EndOfDomain;
499 
500  // Now find out which points in domain have only one vertex
501  for(auto &domIt : domain)
502  {
503  SpatialDomains::CompositeSharedPtr geomvector = domIt.second;
504  for(int i = 0; i < geomvector->m_geomVec.size(); ++i)
505  {
506  for(int j = 0; j < 2; ++j)
507  {
508  int vid = geomvector->m_geomVec[i]->GetVid(j);
509  if(EndOfDomain.count(vid) == 0)
510  {
511  EndOfDomain[vid] = geomvector->m_geomVec[i]->GetVertex(j);
512  }
513  else
514  {
515  EndOfDomain.erase(vid);
516  }
517  }
518  }
519  }
520  ASSERTL1(EndOfDomain.size() == 2,"Did not find two ends of domain");
521 
522  int numNewBc = 1;
523  for(auto &regIt : EndOfDomain)
524  {
525  if(GeometryToRegionsMap.count(regIt.first) != 0)
526  {
527  // Set up boundary condition up
528  auto iter = GeometryToRegionsMap.find(regIt.first);
529  ASSERTL1(iter != GeometryToRegionsMap.end(),
530  "Failied to find GeometryToRegionMap");
531 
532  int regionId = iter->second;
533  auto bregionsIter = bregions.find(regionId);
534  ASSERTL1(bregionsIter != bregions.end(),
535  "Failed to find boundary region");
536 
538  bregionsIter->second;
539  returnval->AddBoundaryRegions(regionId, breg);
540 
541  auto bconditionsIter = bconditions.find(regionId);
542  ASSERTL1(bconditionsIter != bconditions.end(),
543  "Failed to find boundary collection");
544 
546  bconditionsIter->second;
547  returnval->AddBoundaryConditions(regionId,bcond);
548  }
549  else // Set up an undefined region.
550  {
552 
553  // Set up Composite (GemetryVector) to contain vertex and put into bRegion
557  gvec->m_geomVec.push_back(regIt.second);
558  (*breg)[regIt.first] = gvec;
559 
560  returnval->AddBoundaryRegions(bregions.size()+numNewBc,breg);
561 
563 
564  // Set up just boundary condition for this variable.
566  (*bCondition)[variable] = notDefinedCondition;
567 
568  returnval->AddBoundaryConditions(bregions.size()+numNewBc,bCondition);
569  ++numNewBc;
570 
571  }
572  }
573 
574  return returnval;
575  }
std::shared_ptr< BoundaryRegion > BoundaryRegionShPtr
Definition: Conditions.h:215
std::shared_ptr< BoundaryConditionMap > BoundaryConditionMapShPtr
Definition: Conditions.h:225

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL1, Nektar::SpatialDomains::BoundaryConditions::GetBoundaryConditions(), Nektar::SpatialDomains::BoundaryConditions::GetBoundaryRegions(), and Nektar::MultiRegions::ExpList::m_session.

Referenced by DisContField().

◆ GetGlobalBndLinSys()

GlobalLinSysSharedPtr Nektar::MultiRegions::DisContField::GetGlobalBndLinSys ( const GlobalLinSysKey mkey)

For a given key, returns the associated global linear system.

Definition at line 2742 of file DisContField.cpp.

2744  {
2745  ASSERTL0(mkey.GetMatrixType() == StdRegions::eHybridDGHelmBndLam,
2746  "Routine currently only tested for HybridDGHelmholtz");
2747 
2748  ASSERTL1(mkey.GetGlobalSysSolnType() != eDirectFullMatrix,
2749  "Full matrix global systems are not supported for HDG "
2750  "expansions");
2751 
2752  ASSERTL1(mkey.GetGlobalSysSolnType()
2753  ==m_traceMap->GetGlobalSysSolnType(),
2754  "The local to global map is not set up for the requested "
2755  "solution type");
2756 
2757  GlobalLinSysSharedPtr glo_matrix;
2758  auto matrixIter = m_globalBndMat->find(mkey);
2759 
2760  if (matrixIter == m_globalBndMat->end())
2761  {
2762  glo_matrix = GenGlobalBndLinSys(mkey,m_traceMap);
2763  (*m_globalBndMat)[mkey] = glo_matrix;
2764  }
2765  else
2766  {
2767  glo_matrix = matrixIter->second;
2768  }
2769 
2770  return glo_matrix;
2771  }
std::shared_ptr< GlobalLinSys > GenGlobalBndLinSys(const GlobalLinSysKey &mkey, const AssemblyMapSharedPtr &locToGloMap)
Generate a GlobalLinSys from information provided by the key "mkey" and the mapping provided in LocTo...
Definition: ExpList.cpp:2329
std::shared_ptr< GlobalLinSys > GlobalLinSysSharedPtr
Pointer to a GlobalLinSys object.
Definition: GlobalLinSys.h:50

References ASSERTL0, ASSERTL1, Nektar::MultiRegions::eDirectFullMatrix, Nektar::StdRegions::eHybridDGHelmBndLam, Nektar::MultiRegions::ExpList::GenGlobalBndLinSys(), Nektar::MultiRegions::GlobalLinSysKey::GetGlobalSysSolnType(), Nektar::MultiRegions::GlobalMatrixKey::GetMatrixType(), m_globalBndMat, and m_traceMap.

Referenced by v_HelmSolve().

◆ GetNegatedFluxNormal()

vector< bool > & Nektar::MultiRegions::DisContField::GetNegatedFluxNormal ( void  )

Definition at line 2809 of file DisContField.cpp.

2810  {
2811  if(m_negatedFluxNormal.size() == 0)
2812  {
2813  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
2814  &elmtToTrace = m_traceMap->GetElmtToTrace();
2815 
2816  m_negatedFluxNormal.resize(2*GetExpSize());
2817 
2818  for(int i = 0; i < GetExpSize(); ++i)
2819  {
2820 
2821  for(int v = 0; v < 2; ++v)
2822  {
2824  elmtToTrace[i][v]->as<LocalRegions::Expansion0D>();
2825 
2826  if(vertExp->GetLeftAdjacentElementExp()->GetGeom()
2827  ->GetGlobalID() !=
2828  (*m_exp)[i]->GetGeom()->GetGlobalID())
2829  {
2830  m_negatedFluxNormal[2*i+v] = true;
2831  }
2832  else
2833  {
2834  m_negatedFluxNormal[2*i+v] = false;
2835  }
2836  }
2837  }
2838 
2839  }
2840 
2841  return m_negatedFluxNormal;
2842  }
std::vector< bool > m_negatedFluxNormal
Definition: DisContField.h:373
std::shared_ptr< Expansion0D > Expansion0DSharedPtr
Definition: Expansion0D.h:48

References Nektar::MultiRegions::ExpList::GetExpSize(), m_negatedFluxNormal, and m_traceMap.

Referenced by v_AddTraceIntegral().

◆ IsLeftAdjacentTrace()

bool Nektar::MultiRegions::DisContField::IsLeftAdjacentTrace ( const int  n,
const int  e 
)
protected

Definition at line 428 of file DisContField.cpp.

429  {
431  m_traceMap->GetElmtToTrace()[n][e];
432 
433  PeriodicMap periodicTraces = (m_expType == e1D)? m_periodicVerts:
435 
436  bool fwd = true;
437  if (traceEl->GetLeftAdjacentElementTrace () == -1 ||
438  traceEl->GetRightAdjacentElementTrace() == -1)
439  {
440  // Boundary edge (1 connected element). Do nothing in
441  // serial.
442  auto it = m_boundaryTraces.find(traceEl->GetElmtId());
443 
444  // If the edge does not have a boundary condition set on
445  // it, then assume it is a partition edge or periodic.
446  if (it == m_boundaryTraces.end())
447  {
448  fwd = true;
449  }
450  }
451  else if (traceEl->GetLeftAdjacentElementTrace () != -1 &&
452  traceEl->GetRightAdjacentElementTrace() != -1)
453  {
454  // Non-boundary edge (2 connected elements).
455  fwd = (traceEl->GetLeftAdjacentElementExp().get())
456  == (*m_exp)[n].get();
457  }
458  else
459  {
460  ASSERTL2(false, "Unconnected trace element!");
461  }
462 
463  return fwd;
464  }
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed to...
Definition: ErrorUtil.hpp:274

References ASSERTL2, Nektar::MultiRegions::e1D, Nektar::MultiRegions::e2D, m_boundaryTraces, Nektar::MultiRegions::ExpList::m_exp, Nektar::MultiRegions::ExpList::m_expType, m_periodicEdges, m_periodicFaces, m_periodicVerts, and m_traceMap.

Referenced by SetUpDG().

◆ L2_DGDeriv()

NekDouble Nektar::MultiRegions::DisContField::L2_DGDeriv ( const int  dir,
const Array< OneD, const NekDouble > &  soln 
)

Calculate the \( L^2 \) error of the \( Q_{\rm dir} \) derivative using the consistent DG evaluation of \( Q_{\rm dir} \).

The solution provided is of the primative variation at the quadrature points and the derivative is compared to the discrete derivative at these points, which is likely to be undesirable unless using a much higher number of quadrature points than the polynomial order used to evaluate \( Q_{\rm dir} \).

Definition at line 4117 of file DisContField.cpp.

4120  {
4121  int i,e,ncoeff_edge;
4122  Array<OneD, const NekDouble> tmp_coeffs;
4123  Array<OneD, NekDouble> out_d(m_ncoeffs), out_tmp;
4124 
4125  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
4126  &elmtToTrace = m_traceMap->GetElmtToTrace();
4127 
4128  StdRegions::Orientation edgedir;
4129 
4130  int cnt;
4131  int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs();
4132  Array<OneD, NekDouble> loc_lambda(LocBndCoeffs), edge_lambda;
4133 
4134 
4135  m_traceMap->GlobalToLocalBnd(m_trace->GetCoeffs(),loc_lambda);
4136 
4137  edge_lambda = loc_lambda;
4138 
4139  // Calculate Q using standard DG formulation.
4140  for(i = cnt = 0; i < GetExpSize(); ++i)
4141  {
4142  // Probably a better way of setting up lambda than this.
4143  // Note cannot use PutCoeffsInToElmts since lambda space
4144  // is mapped during the solve.
4145  int nEdges = (*m_exp)[i]->GetGeom()->GetNumEdges();
4146  Array<OneD, Array<OneD, NekDouble> > edgeCoeffs(nEdges);
4147 
4148  for(e = 0; e < nEdges; ++e)
4149  {
4150  edgedir = (*m_exp)[i]->GetTraceOrient(e);
4151  ncoeff_edge = elmtToTrace[i][e]->GetNcoeffs();
4152  edgeCoeffs[e] = Array<OneD, NekDouble>(ncoeff_edge);
4153  Vmath::Vcopy(ncoeff_edge, edge_lambda, 1, edgeCoeffs[e], 1);
4154  elmtToTrace[i][e]->SetCoeffsToOrientation(
4155  edgedir, edgeCoeffs[e], edgeCoeffs[e]);
4156  edge_lambda = edge_lambda + ncoeff_edge;
4157  }
4158 
4159  (*m_exp)[i]->DGDeriv(dir,
4160  tmp_coeffs=m_coeffs+m_coeff_offset[i],
4161  elmtToTrace[i],
4162  edgeCoeffs,
4163  out_tmp = out_d+cnt);
4164  cnt += (*m_exp)[i]->GetNcoeffs();
4165  }
4166 
4167  BwdTrans(out_d,m_phys);
4168  Vmath::Vsub(m_npoints,m_phys,1,soln,1,m_phys,1);
4169  return L2(m_phys);
4170  }
int m_ncoeffs
The total number of local degrees of freedom. m_ncoeffs .
Definition: ExpList.h:1230
NekDouble L2(const Array< OneD, const NekDouble > &inarray, const Array< OneD, const NekDouble > &soln=NullNekDouble1DArray)
This function calculates the error with respect to soln of the global spectral/hp element approximat...
Definition: ExpList.h:596
void BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.h:2015
Array< OneD, NekDouble > m_phys
The global expansion evaluated at the quadrature points.
Definition: ExpList.h:1269
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:372

References Nektar::MultiRegions::ExpList::BwdTrans(), Nektar::MultiRegions::ExpList::GetExpSize(), Nektar::MultiRegions::ExpList::L2(), Nektar::MultiRegions::ExpList::m_coeff_offset, Nektar::MultiRegions::ExpList::m_coeffs, Nektar::MultiRegions::ExpList::m_ncoeffs, Nektar::MultiRegions::ExpList::m_npoints, Nektar::MultiRegions::ExpList::m_phys, m_trace, m_traceMap, Vmath::Vcopy(), and Vmath::Vsub().

◆ SameTypeOfBoundaryConditions()

bool Nektar::MultiRegions::DisContField::SameTypeOfBoundaryConditions ( const DisContField In)

Check to see if expansion has the same BCs as In.

For each boundary region, checks that the types and number of boundary expansions in that region match.

Parameters
InField to compare with.
Returns
True if boundary conditions match.

Definition at line 2779 of file DisContField.cpp.

2781  {
2782  int i;
2783  bool returnval = true;
2784 
2785  for(i = 0; i < m_bndConditions.size(); ++i)
2786  {
2787  // check to see if boundary condition type is the same
2788  // and there are the same number of boundary
2789  // conditions in the boundary definition.
2790  if((m_bndConditions[i]->GetBoundaryConditionType()
2791  != In.m_bndConditions[i]->GetBoundaryConditionType())||
2792  (m_bndCondExpansions[i]->GetExpSize()
2793  != In.m_bndCondExpansions[i]->GetExpSize()))
2794  {
2795  returnval = false;
2796  break;
2797  }
2798  }
2799 
2800  // Compare with all other processes. Return true only if all
2801  // processes report having the same boundary conditions.
2802  int vSame = (returnval?1:0);
2803  m_comm->GetRowComm()->AllReduce(vSame, LibUtilities::ReduceMin);
2804 
2805  return (vSame == 1);
2806  }
LibUtilities::CommSharedPtr m_comm
Communicator.
Definition: ExpList.h:1220

References m_bndCondExpansions, m_bndConditions, Nektar::MultiRegions::ExpList::m_comm, and Nektar::LibUtilities::ReduceMin.

Referenced by Nektar::MultiRegions::ContField::ContField(), and DisContField().

◆ SetUpDG()

void Nektar::MultiRegions::DisContField::SetUpDG ( const std::string  variable = "DefaultVar")
protected

Set up all DG member variables and maps.

Definition at line 159 of file DisContField.cpp.

160  {
161  // Check for multiple calls
163  {
164  return;
165  }
166 
167  // Set up matrix map
170 
171  // Set up trace space
174  *m_exp,m_graph);
175 
176  PeriodicMap periodicTraces = (m_expType == e1D)? m_periodicVerts:
178 
182  periodicTraces, variable);
183 
184  if (m_session->DefinesCmdLineArgument("verbose"))
185  {
186  m_traceMap->PrintStats(std::cout, variable);
187  }
188 
189  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
190  &elmtToTrace = m_traceMap->GetElmtToTrace();
191 
192  // Scatter trace to elements. For each element, we find
193  // the trace associated to each elemental trace. The element
194  // then retains a pointer to the elemental trace, to
195  // ensure uniqueness of normals when retrieving from two
196  // adjoining elements which do not lie in a plane.
197 
198  for (int i = 0; i < m_exp->size(); ++i)
199  {
200  for (int j = 0; j < (*m_exp)[i]->GetNtraces(); ++j)
201  {
202  LocalRegions::ExpansionSharedPtr exp = elmtToTrace[i][j];;
203 
204  exp->SetAdjacentElementExp(j, (*m_exp)[i]);
205  (*m_exp)[i]->SetTraceExp (j, exp);
206  }
207  }
208 
209  // Set up physical normals
211 
212  int cnt, n;
213 
214  // Identify boundary trace
215  for(cnt = 0, n = 0; n < m_bndCondExpansions.size(); ++n)
216  {
217  if (m_bndConditions[n]->GetBoundaryConditionType() !=
219  {
220  for(int v = 0; v < m_bndCondExpansions[n]->GetExpSize(); ++v)
221  {
222  m_boundaryTraces.insert
223  (m_traceMap->GetBndCondIDToGlobalTraceID(cnt+v));
224  }
225  cnt += m_bndCondExpansions[n]->GetExpSize();
226  }
227  }
228 
229  // Set up information for periodic boundary conditions.
230  std::unordered_map<int,pair<int,int> > perTraceToExpMap;
231  for (cnt = n = 0; n < m_exp->size(); ++n)
232  {
233  for (int v = 0; v < (*m_exp)[n]->GetNtraces(); ++v, ++cnt)
234  {
235  auto it = periodicTraces.find
236  ((*m_exp)[n]->GetGeom()->GetTid(v));
237 
238  if (it != periodicTraces.end())
239  {
240  perTraceToExpMap[it->first] = make_pair(n,v);
241  }
242  }
243  }
244 
245  // Set up left-adjacent tracelist.
246  m_leftAdjacentTraces.resize(cnt);
247 
248  // count size of trace
249  for (cnt = n = 0; n < m_exp->size(); ++n)
250  {
251  for (int v = 0; v < (*m_exp)[n]->GetNtraces(); ++v, ++cnt)
252  {
254  }
255  }
256 
257 
258  // Set up mapping to copy Fwd of periodic bcs to Bwd of other edge.
259  for (cnt = n = 0; n < m_exp->size(); ++n)
260  {
261  for (int v = 0; v < (*m_exp)[n]->GetNtraces(); ++v, ++cnt)
262  {
263  int GeomId = (*m_exp)[n]->GetGeom()->GetTid(v);
264 
265  // Check to see if this trace is periodic.
266  auto it = periodicTraces.find(GeomId);
267 
268  if (it != periodicTraces.end())
269  {
270  const PeriodicEntity &ent = it->second[0];
271  auto it2 = perTraceToExpMap.find(ent.id);
272 
273  if (it2 == perTraceToExpMap.end())
274  {
275  if (m_session->GetComm()->GetRowComm()->GetSize()
276  > 1 && !ent.isLocal)
277  {
278  continue;
279  }
280  else
281  {
282  ASSERTL1(false, "Periodic trace not found!");
283  }
284  }
285 
287  "Periodic trace in non-forward space?");
288 
289  int offset = m_trace->GetPhys_Offset
290  ((m_traceMap->GetElmtToTrace())
291  [n][v]->GetElmtId());
292 
293  int offset2 = m_trace->GetPhys_Offset
294  ((m_traceMap->GetElmtToTrace())
295  [it2->second.first]
296  [it2->second.second]->GetElmtId());
297 
298 
299  switch(m_expType)
300  {
301  case e1D:
302  {
303  m_periodicFwdCopy.push_back(offset);
304  m_periodicBwdCopy.push_back(offset2);
305  }
306  break;
307  case e2D:
308  {
309  // Calculate relative orientations between trace to
310  // calculate copying map.
311  int nquad = elmtToTrace[n][v]->GetNumPoints(0);
312 
313  vector<int> tmpBwd(nquad);
314  vector<int> tmpFwd(nquad);
315 
316  if (ent.orient == StdRegions::eForwards)
317  {
318  for (int i = 0; i < nquad; ++i)
319  {
320  tmpBwd[i] = offset2 + i;
321  tmpFwd[i] = offset + i;
322  }
323  }
324  else
325  {
326  for (int i = 0; i < nquad; ++i)
327  {
328  tmpBwd[i] = offset2 + i;
329  tmpFwd[i] = offset + nquad - i - 1;
330  }
331  }
332 
333  for (int i = 0; i < nquad; ++i)
334  {
335  m_periodicFwdCopy.push_back(tmpFwd[i]);
336  m_periodicBwdCopy.push_back(tmpBwd[i]);
337  }
338  }
339  break;
340  case e3D:
341  {
342  // Calculate relative orientations between faces to
343  // calculate copying map.
344  int nquad1 = elmtToTrace[n][v]->GetNumPoints(0);
345  int nquad2 = elmtToTrace[n][v]->GetNumPoints(1);
346 
347  vector<int> tmpBwd(nquad1*nquad2);
348  vector<int> tmpFwd(nquad1*nquad2);
349 
350  if (ent.orient == StdRegions::eDir1FwdDir2_Dir2FwdDir1 ||
351  ent.orient == StdRegions::eDir1BwdDir2_Dir2FwdDir1 ||
352  ent.orient == StdRegions::eDir1FwdDir2_Dir2BwdDir1 ||
354  {
355  for (int i = 0; i < nquad2; ++i)
356  {
357  for (int j = 0; j < nquad1; ++j)
358  {
359  tmpBwd[i*nquad1+j] = offset2 + i*nquad1+j;
360  tmpFwd[i*nquad1+j] = offset + j*nquad2+i;
361  }
362  }
363  }
364  else
365  {
366  for (int i = 0; i < nquad2; ++i)
367  {
368  for (int j = 0; j < nquad1; ++j)
369  {
370  tmpBwd[i*nquad1+j] = offset2 + i*nquad1+j;
371  tmpFwd[i*nquad1+j] = offset + i*nquad1+j;
372  }
373  }
374  }
375 
376  if (ent.orient==StdRegions::eDir1BwdDir1_Dir2FwdDir2 ||
380  {
381  // Reverse x direction
382  for (int i = 0; i < nquad2; ++i)
383  {
384  for (int j = 0; j < nquad1/2; ++j)
385  {
386  swap(tmpFwd[i*nquad1 + j],
387  tmpFwd[i*nquad1 + nquad1-j-1]);
388  }
389  }
390  }
391 
392  if (ent.orient==StdRegions::eDir1FwdDir1_Dir2BwdDir2 ||
396  {
397  // Reverse y direction
398  for (int j = 0; j < nquad1; ++j)
399  {
400  for (int i = 0; i < nquad2/2; ++i)
401  {
402  swap(tmpFwd[i*nquad1 + j],
403  tmpFwd[(nquad2-i-1)*nquad1 + j]);
404  }
405  }
406  }
407 
408  for (int i = 0; i < nquad1*nquad2; ++i)
409  {
410  m_periodicFwdCopy.push_back(tmpFwd[i]);
411  m_periodicBwdCopy.push_back(tmpBwd[i]);
412  }
413  }
414  break;
415  default:
416  ASSERTL1(false,"not set up");
417  }
418  }
419  }
420  }
421 
423  AllocateSharedPtr(*this, m_trace, elmtToTrace,
425  }
bool IsLeftAdjacentTrace(const int n, const int e)

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL1, Nektar::MultiRegions::e1D, Nektar::MultiRegions::e2D, Nektar::MultiRegions::e3D, Nektar::StdRegions::eDir1BwdDir1_Dir2BwdDir2, Nektar::StdRegions::eDir1BwdDir1_Dir2FwdDir2, Nektar::StdRegions::eDir1BwdDir2_Dir2BwdDir1, Nektar::StdRegions::eDir1BwdDir2_Dir2FwdDir1, Nektar::StdRegions::eDir1FwdDir1_Dir2BwdDir2, Nektar::StdRegions::eDir1FwdDir2_Dir2BwdDir1, Nektar::StdRegions::eDir1FwdDir2_Dir2FwdDir1, Nektar::StdRegions::eForwards, Nektar::SpatialDomains::ePeriodic, Nektar::MultiRegions::PeriodicEntity::id, IsLeftAdjacentTrace(), Nektar::MultiRegions::PeriodicEntity::isLocal, m_bndCondExpansions, m_bndConditions, m_boundaryTraces, Nektar::MultiRegions::ExpList::m_exp, Nektar::MultiRegions::ExpList::m_expType, m_globalBndMat, Nektar::MultiRegions::ExpList::m_graph, m_leftAdjacentTraces, m_locTraceToTraceMap, m_periodicBwdCopy, m_periodicEdges, m_periodicFaces, m_periodicFwdCopy, m_periodicVerts, Nektar::MultiRegions::ExpList::m_session, m_trace, m_traceMap, Nektar::MultiRegions::NullExpListSharedPtr, Nektar::MultiRegions::PeriodicEntity::orient, and Nektar::MultiRegions::ExpList::SetUpPhysNormals().

Referenced by DisContField(), and v_GetTrace().

◆ v_AddFwdBwdTraceIntegral()

void Nektar::MultiRegions::DisContField::v_AddFwdBwdTraceIntegral ( const Array< OneD, const NekDouble > &  Fwd,
const Array< OneD, const NekDouble > &  Bwd,
Array< OneD, NekDouble > &  outarray 
)
protectedvirtual

Add trace contributions into elemental coefficient spaces.

Given some quantity \( \vec{q} \), calculate the elemental integral

\[ \int_{\Omega^e} \vec{q}, \mathrm{d}S \]

and adds this to the coefficient space provided by outarray. The value of q is determined from the routine IsLeftAdjacentTrace() which if true we use Fwd else we use Bwd

Parameters
FwdThe trace quantities associated with left (fwd) adjancent elmt.
BwdThe trace quantities associated with right (bwd) adjacent elet.
outarrayResulting coefficient space.

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 3362 of file DisContField.cpp.

3366  {
3367 
3368  ASSERTL0(m_expType != e1D, "This method is not setup or "
3369  "tested for 1D expansion");
3370 
3371  Array<OneD, NekDouble> Coeffs(m_trace->GetNcoeffs());
3372 
3373  m_trace->IProductWRTBase(Fwd,Coeffs);
3374  m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(0,Coeffs,outarray);
3375  m_trace->IProductWRTBase(Bwd,Coeffs);
3376  m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(1,Coeffs,outarray);
3377  }

References ASSERTL0, Nektar::MultiRegions::e1D, Nektar::MultiRegions::ExpList::m_expType, m_locTraceToTraceMap, and m_trace.

◆ v_AddTraceIntegral()

void Nektar::MultiRegions::DisContField::v_AddTraceIntegral ( const Array< OneD, const NekDouble > &  Fn,
Array< OneD, NekDouble > &  outarray 
)
protectedvirtual

Add trace contributions into elemental coefficient spaces.

Given some quantity \( \vec{Fn} \), which conatins this routine calculates the integral

\[ \int_{\Omega^e} \vec{Fn}, \mathrm{d}S \]

and adds this to the coefficient space provided by outarray.

See also
Expansion3D::AddFaceNormBoundaryInt
Parameters
FnThe trace quantities.
outarrayResulting 3D coefficient space.

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 3167 of file DisContField.cpp.

3170  {
3171  if(m_expType == e1D)
3172  {
3173  int n,offset, t_offset;
3174 
3175  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
3176  &elmtToTrace = m_traceMap->GetElmtToTrace();
3177 
3178  vector<bool> negatedFluxNormal = GetNegatedFluxNormal();
3179 
3180  for (n = 0; n < GetExpSize(); ++n)
3181  {
3182  // Number of coefficients on each element
3183  int e_ncoeffs = (*m_exp)[n]->GetNcoeffs();
3184 
3185  offset = GetCoeff_Offset(n);
3186 
3187  // Implementation for every points except Gauss points
3188  if ((*m_exp)[n]->GetBasis(0)->GetBasisType() !=
3190  {
3191  t_offset = GetTrace()->GetCoeff_Offset
3192  (elmtToTrace[n][0]->GetElmtId());
3193  if(negatedFluxNormal[2*n])
3194  {
3195  outarray[offset] -= Fn[t_offset];
3196  }
3197  else
3198  {
3199  outarray[offset] += Fn[t_offset];
3200  }
3201 
3202  t_offset = GetTrace()->GetCoeff_Offset
3203  (elmtToTrace[n][1]->GetElmtId());
3204 
3205  if(negatedFluxNormal[2*n+1])
3206  {
3207  outarray[offset+(*m_exp)[n]->GetVertexMap(1)] -=
3208  Fn[t_offset];
3209  }
3210  else
3211  {
3212  outarray[offset+(*m_exp)[n]->GetVertexMap(1)] +=
3213  Fn[t_offset];
3214  }
3215 
3216  }
3217  else
3218  {
3219  int j;
3220  static DNekMatSharedPtr m_Ixm, m_Ixp;
3221  static int sav_ncoeffs = 0;
3222  if(!m_Ixm.get() || e_ncoeffs != sav_ncoeffs)
3223  {
3225  const LibUtilities::PointsKey
3226  BS_p(e_ncoeffs,LibUtilities::eGaussGaussLegendre);
3227  const LibUtilities::BasisKey
3228  BS_k(LibUtilities::eGauss_Lagrange,e_ncoeffs,BS_p);
3229 
3230  BASE = LibUtilities::BasisManager()[BS_k];
3231 
3232  Array<OneD, NekDouble> coords(1, 0.0);
3233 
3234  coords[0] = -1.0;
3235  m_Ixm = BASE->GetI(coords);
3236 
3237  coords[0] = 1.0;
3238  m_Ixp = BASE->GetI(coords);
3239 
3240  sav_ncoeffs = e_ncoeffs;
3241  }
3242 
3243  t_offset = GetTrace()->GetCoeff_Offset
3244  (elmtToTrace[n][0]->GetElmtId());
3245 
3246  if(negatedFluxNormal[2*n])
3247  {
3248  for (j = 0; j < e_ncoeffs; j++)
3249  {
3250  outarray[offset + j] -=
3251  (m_Ixm->GetPtr())[j] * Fn[t_offset];
3252  }
3253  }
3254  else
3255  {
3256  for (j = 0; j < e_ncoeffs; j++)
3257  {
3258  outarray[offset + j] +=
3259  (m_Ixm->GetPtr())[j] * Fn[t_offset];
3260  }
3261  }
3262 
3263  t_offset = GetTrace()->GetCoeff_Offset
3264  (elmtToTrace[n][1]->GetElmtId());
3265 
3266  if (negatedFluxNormal[2*n+1])
3267  {
3268  for (j = 0; j < e_ncoeffs; j++)
3269  {
3270  outarray[offset + j] -=
3271  (m_Ixp->GetPtr())[j] * Fn[t_offset];
3272  }
3273  }
3274  else
3275  {
3276  for (j = 0; j < e_ncoeffs; j++)
3277  {
3278  outarray[offset + j] +=
3279  (m_Ixp->GetPtr())[j] * Fn[t_offset];
3280  }
3281  }
3282  }
3283  }
3284  }
3285  else // other dimensions
3286  {
3287  // Basis definition on each element
3288  LibUtilities::BasisSharedPtr basis = (*m_exp)[0]->GetBasis(0);
3289  if((m_expType != e1D)&&
3290  (basis->GetBasisType() != LibUtilities::eGauss_Lagrange))
3291  {
3292  Array<OneD, NekDouble> Fcoeffs(m_trace->GetNcoeffs());
3293  m_trace->IProductWRTBase(Fn, Fcoeffs);
3294 
3295  m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(Fcoeffs,
3296  outarray);
3297  }
3298  else
3299  {
3300  int e, n, offset, t_offset;
3301  Array<OneD, NekDouble> e_outarray;
3302  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
3303  &elmtToTrace = m_traceMap->GetElmtToTrace();
3304 
3305  if(m_expType == e2D)
3306  {
3307  for(n = 0; n < GetExpSize(); ++n)
3308  {
3309  offset = GetCoeff_Offset(n);
3310  for(e = 0; e < (*m_exp)[n]->GetNtraces(); ++e)
3311  {
3312  t_offset = GetTrace()->GetPhys_Offset
3313  (elmtToTrace[n][e]->GetElmtId());
3314  (*m_exp)[n]->AddEdgeNormBoundaryInt
3315  (e, elmtToTrace[n][e],
3316  Fn+t_offset,
3317  e_outarray = outarray+offset);
3318  }
3319  }
3320  }
3321  else if (m_expType == e3D)
3322  {
3323  for(n = 0; n < GetExpSize(); ++n)
3324  {
3325  offset = GetCoeff_Offset(n);
3326  for(e = 0; e < (*m_exp)[n]->GetNtraces(); ++e)
3327  {
3328  t_offset = GetTrace()->GetPhys_Offset
3329  (elmtToTrace[n][e]->GetElmtId());
3330  (*m_exp)[n]->AddFaceNormBoundaryInt
3331  (e, elmtToTrace[n][e],
3332  Fn+t_offset,
3333  e_outarray = outarray+offset);
3334  }
3335  }
3336 
3337  }
3338  }
3339  }
3340  }
std::vector< bool > & GetNegatedFluxNormal(void)
std::shared_ptr< ExpList > & GetTrace()
Definition: ExpList.h:2525
int GetCoeff_Offset(int n) const
Get the start offset position for a global list of m_coeffs correspoinding to element n.
Definition: ExpList.h:2431
BasisManagerT & BasisManager(void)
std::shared_ptr< Basis > BasisSharedPtr
@ eGaussGaussLegendre
1D Gauss-Gauss-Legendre quadrature points
Definition: PointsType.h:48
@ eGauss_Lagrange
Lagrange Polynomials using the Gauss points .
Definition: BasisType.h:55
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:69

References Nektar::LibUtilities::BasisManager(), Nektar::MultiRegions::e1D, Nektar::MultiRegions::e2D, Nektar::MultiRegions::e3D, Nektar::LibUtilities::eGauss_Lagrange, Nektar::LibUtilities::eGaussGaussLegendre, Nektar::MultiRegions::ExpList::GetCoeff_Offset(), Nektar::MultiRegions::ExpList::GetExpSize(), GetNegatedFluxNormal(), Nektar::MultiRegions::ExpList::GetTrace(), Nektar::MultiRegions::ExpList::m_exp, Nektar::MultiRegions::ExpList::m_expType, m_locTraceToTraceMap, m_trace, and m_traceMap.

◆ v_AddTraceIntegralToOffDiag()

void Nektar::MultiRegions::DisContField::v_AddTraceIntegralToOffDiag ( const Array< OneD, const NekDouble > &  FwdFlux,
const Array< OneD, const NekDouble > &  BwdFlux,
Array< OneD, NekDouble > &  outarray 
)
privatevirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 4428 of file DisContField.cpp.

4432  {
4433  Array<OneD, NekDouble> FCoeffs(m_trace->GetNcoeffs());
4434 
4435  m_trace->IProductWRTBase(FwdFlux,FCoeffs);
4436  m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(1, FCoeffs,
4437  outarray);
4438  m_trace->IProductWRTBase(BwdFlux,FCoeffs);
4439  m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(0, FCoeffs,
4440  outarray);
4441  }

References m_locTraceToTraceMap, and m_trace.

◆ v_AddTraceQuadPhysToField()

void Nektar::MultiRegions::DisContField::v_AddTraceQuadPhysToField ( const Array< OneD, const NekDouble > &  Fwd,
const Array< OneD, const NekDouble > &  Bwd,
Array< OneD, NekDouble > &  field 
)
protectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 3054 of file DisContField.cpp.

3058  {
3059  // Basis definition on each element
3060  LibUtilities::BasisSharedPtr basis = (*m_exp)[0]->GetBasis(0);
3061  if (basis->GetBasisType() != LibUtilities::eGauss_Lagrange)
3062  {
3063  Array<OneD, NekDouble> edgevals(m_locTraceToTraceMap->
3064  GetNLocTracePts(), 0.0);
3065 
3066  Array<OneD, NekDouble> invals = edgevals +
3067  m_locTraceToTraceMap->GetNFwdLocTracePts();
3068  m_locTraceToTraceMap->RightIPTWLocEdgesToTraceInterpMat(
3069  1, Bwd, invals);
3070 
3071  m_locTraceToTraceMap->RightIPTWLocEdgesToTraceInterpMat(
3072  0, Fwd, edgevals);
3073 
3074  m_locTraceToTraceMap->AddLocTracesToField(edgevals, field);
3075  }
3076  else
3077  {
3079  "v_AddTraceQuadPhysToField not coded for eGauss_Lagrange");
3080  }
3081  }

References Nektar::ErrorUtil::efatal, Nektar::LibUtilities::eGauss_Lagrange, m_locTraceToTraceMap, and NEKERROR.

◆ v_EvaluateBoundaryConditions()

void Nektar::MultiRegions::DisContField::v_EvaluateBoundaryConditions ( const NekDouble  time = 0.0,
const std::string  varName = "",
const NekDouble  x2_in = NekConstants::kNekUnsetDouble,
const NekDouble  x3_in = NekConstants::kNekUnsetDouble 
)
protectedvirtual

Evaluate all boundary conditions at a given time..

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 3555 of file DisContField.cpp.

3560  {
3561  int i;
3562  int npoints;
3563 
3564  MultiRegions::ExpListSharedPtr locExpList;
3565 
3566  for (i = 0; i < m_bndCondExpansions.size(); ++i)
3567  {
3568  if (time == 0.0 || m_bndConditions[i]->IsTimeDependent())
3569  {
3570  m_bndCondBndWeight[i] = 1.0;
3571  locExpList = m_bndCondExpansions[i];
3572 
3573  npoints = locExpList->GetNpoints();
3574  Array<OneD, NekDouble> x0(npoints, 0.0);
3575  Array<OneD, NekDouble> x1(npoints, 0.0);
3576  Array<OneD, NekDouble> x2(npoints, 0.0);
3577 
3578  locExpList->GetCoords(x0, x1, x2);
3579 
3580  if (x2_in != NekConstants::kNekUnsetDouble && x3_in !=
3582  {
3583  Vmath::Fill(npoints,x2_in,x1,1);
3584  Vmath::Fill(npoints,x3_in,x2,1);
3585  }
3586  else if(x2_in != NekConstants::kNekUnsetDouble)
3587  {
3588  Vmath::Fill(npoints,x2_in,x2,1);
3589  }
3590 
3591  // treat 1D expansions separately since we only
3592  // require an evaluation at a point rather than
3593  // any projections or inner products that are not
3594  // available in a PointExp
3595  if(m_expType == e1D)
3596  {
3597  if (m_bndConditions[i]->GetBoundaryConditionType() ==
3599  {
3600  m_bndCondExpansions[i]->SetCoeff
3601  (0,(std::static_pointer_cast<SpatialDomains
3602  ::DirichletBoundaryCondition>
3603  (m_bndConditions[i])
3604  ->m_dirichletCondition).Evaluate
3605  (x0[0],x1[0],x2[0],time));
3606  m_bndCondExpansions[i]->SetPhys
3607  (0,m_bndCondExpansions[i]->GetCoeff(0));
3608  }
3609  else if (m_bndConditions[i]->GetBoundaryConditionType()
3611  {
3612  m_bndCondExpansions[i]->SetCoeff
3613  (0,(std::static_pointer_cast<SpatialDomains
3614  ::NeumannBoundaryCondition>
3615  (m_bndConditions[i])
3616  ->m_neumannCondition).Evaluate
3617  (x0[0],x1[0],x2[0],time));
3618  }
3619  else if (m_bndConditions[i]->GetBoundaryConditionType()
3621  {
3622  m_bndCondExpansions[i]->SetCoeff
3623  (0,(std::static_pointer_cast<SpatialDomains
3624  ::RobinBoundaryCondition>
3625  (m_bndConditions[i])
3626  ->m_robinFunction).Evaluate
3627  (x0[0],x1[0],x2[0],time));
3628 
3629  }
3630  else if (m_bndConditions[i]->GetBoundaryConditionType()
3632  {
3633  continue;
3634  }
3635  else if (m_bndConditions[i]->GetBoundaryConditionType()
3637  {
3638  }
3639  else
3640  {
3642  "This type of BC not implemented yet");
3643  }
3644  }
3645  else // 2D and 3D versions
3646  {
3647  if (m_bndConditions[i]->GetBoundaryConditionType()
3649  {
3651  = std::static_pointer_cast<
3652  SpatialDomains::DirichletBoundaryCondition>
3653  (m_bndConditions[i]);
3654 
3655  Array<OneD, NekDouble> valuesFile(npoints, 1.0),
3656  valuesExp(npoints, 1.0);
3657 
3658  string filebcs = bcPtr->m_filename;
3659  string exprbcs = bcPtr->m_expr;
3660 
3661  if (filebcs != "")
3662  {
3664  (filebcs, bcPtr->GetComm(), varName, locExpList);
3665  valuesFile = locExpList->GetPhys();
3666  }
3667 
3668  if (exprbcs != "")
3669  {
3670  LibUtilities::Equation condition =
3671  std::static_pointer_cast<
3672  SpatialDomains::DirichletBoundaryCondition>
3673  (m_bndConditions[i])->m_dirichletCondition;
3674 
3675  condition.Evaluate(x0, x1, x2, time, valuesExp);
3676  }
3677 
3678  Vmath::Vmul(npoints, valuesExp, 1, valuesFile, 1,
3679  locExpList->UpdatePhys(), 1);
3680 
3681  locExpList->FwdTrans_BndConstrained
3682  (locExpList->GetPhys(),
3683  locExpList->UpdateCoeffs());
3684  }
3685  else if (m_bndConditions[i]->GetBoundaryConditionType()
3687  {
3689  bcPtr = std::static_pointer_cast<
3690  SpatialDomains::NeumannBoundaryCondition>
3691  (m_bndConditions[i]);
3692  string filebcs = bcPtr->m_filename;
3693  if (filebcs != "")
3694  {
3696  (filebcs, bcPtr->GetComm(), varName, locExpList);
3697  }
3698  else
3699  {
3700  LibUtilities::Equation condition =
3701  std::static_pointer_cast<
3702  SpatialDomains::NeumannBoundaryCondition>
3703  (m_bndConditions[i])->m_neumannCondition;
3704  condition.Evaluate(x0, x1, x2, time,
3705  locExpList->UpdatePhys());
3706  }
3707 
3708  locExpList->IProductWRTBase(locExpList->GetPhys(),
3709  locExpList->UpdateCoeffs());
3710  }
3711  else if (m_bndConditions[i]->GetBoundaryConditionType()
3713  {
3715  bcPtr = std::static_pointer_cast<
3716  SpatialDomains::RobinBoundaryCondition>
3717  (m_bndConditions[i]);
3718 
3719  string filebcs = bcPtr->m_filename;
3720 
3721  if (filebcs != "")
3722  {
3724  (filebcs, bcPtr->GetComm(), varName,
3725  locExpList);
3726  }
3727  else
3728  {
3729  LibUtilities::Equation condition =
3730  std::static_pointer_cast<
3731  SpatialDomains::RobinBoundaryCondition>
3732  (m_bndConditions[i])->m_robinFunction;
3733  condition.Evaluate(x0, x1, x2, time,
3734  locExpList->UpdatePhys());
3735  }
3736 
3737  locExpList->IProductWRTBase
3738  (locExpList->GetPhys(),
3739  locExpList->UpdateCoeffs());
3740  }
3741  else if (m_bndConditions[i]->GetBoundaryConditionType()
3743  {
3744  continue;
3745  }
3746  else
3747  {
3749  "This type of BC not implemented yet");
3750  }
3751  }
3752  }
3753  }
3754  }
void ExtractFileBCs(const std::string &fileName, LibUtilities::CommSharedPtr comm, const std::string &varName, const std::shared_ptr< ExpList > locExpList)
Definition: ExpList.cpp:3207
NekDouble GetCoeff(int i)
Get the i th value (coefficient) of m_coeffs.
Definition: ExpList.h:2370
static const NekDouble kNekUnsetDouble
std::shared_ptr< DirichletBoundaryCondition > DirichletBCShPtr
Definition: Conditions.h:220
std::shared_ptr< NeumannBoundaryCondition > NeumannBCShPtr
Definition: Conditions.h:221
std::shared_ptr< RobinBoundaryCondition > RobinBCShPtr
Definition: Conditions.h:222
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.cpp:192
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:45

References Nektar::MultiRegions::e1D, Nektar::SpatialDomains::eDirichlet, Nektar::ErrorUtil::efatal, Nektar::SpatialDomains::eNeumann, Nektar::SpatialDomains::eNotDefined, Nektar::SpatialDomains::ePeriodic, Nektar::SpatialDomains::eRobin, Nektar::LibUtilities::Equation::Evaluate(), Nektar::MultiRegions::ExpList::ExtractFileBCs(), Vmath::Fill(), Nektar::MultiRegions::ExpList::GetCoeff(), Nektar::NekConstants::kNekUnsetDouble, m_bndCondBndWeight, m_bndCondExpansions, m_bndConditions, Nektar::MultiRegions::ExpList::m_expType, Nektar::SpatialDomains::NeumannBoundaryCondition::m_filename, Nektar::SpatialDomains::RobinBoundaryCondition::m_filename, NEKERROR, and Vmath::Vmul().

◆ v_ExtractTracePhys() [1/2]

void Nektar::MultiRegions::DisContField::v_ExtractTracePhys ( Array< OneD, NekDouble > &  outarray)
protectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 3083 of file DisContField.cpp.

3084  {
3085  ASSERTL1(m_physState == true,"local physical space is not true ");
3086  v_ExtractTracePhys(m_phys, outarray);
3087  }
virtual void v_ExtractTracePhys(Array< OneD, NekDouble > &outarray)
bool m_physState
The state of the array m_phys.
Definition: ExpList.h:1278

References ASSERTL1, Nektar::MultiRegions::ExpList::m_phys, and Nektar::MultiRegions::ExpList::m_physState.

◆ v_ExtractTracePhys() [2/2]

void Nektar::MultiRegions::DisContField::v_ExtractTracePhys ( const Array< OneD, const NekDouble > &  inarray,
Array< OneD, NekDouble > &  outarray 
)
protectedvirtual

This method extracts the trace (verts in 1D) from the field inarray and puts the values in outarray.

It assumes the field is C0 continuous so that it can overwrite the edge data when visited by the two adjacent elements.

Parameters
inarrayAn array containing the 1D data from which we wish to extract the edge data.
outarrayThe resulting edge information.

This will not work for non-boundary expansions

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 3103 of file DisContField.cpp.

3106  {
3107  LibUtilities::BasisSharedPtr basis = (*m_exp)[0]->GetBasis(0);
3108  if( (basis->GetBasisType() != LibUtilities::eGauss_Lagrange))
3109  {
3110  Vmath::Zero(outarray.size(), outarray, 1);
3111 
3112  Array<OneD, NekDouble> tracevals(
3113  m_locTraceToTraceMap->GetNFwdLocTracePts());
3114  m_locTraceToTraceMap->FwdLocTracesFromField(inarray,tracevals);
3115  m_locTraceToTraceMap->InterpLocTracesToTrace(0,tracevals,outarray);
3116  m_traceMap->GetAssemblyCommDG()->
3117  PerformExchange(outarray, outarray);
3118  }
3119  else
3120  {
3121 
3122  // Loop over elemente and collect forward expansion
3123  int nexp = GetExpSize();
3124  int n,p,offset,phys_offset;
3125  Array<OneD,NekDouble> t_tmp;
3126 
3127  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
3128  &elmtToTrace = m_traceMap->GetElmtToTrace();
3129 
3130  ASSERTL1(outarray.size() >= m_trace->GetNpoints(),
3131  "input array is of insufficient length");
3132 
3133  for (n = 0; n < nexp; ++n)
3134  {
3135  phys_offset = GetPhys_Offset(n);
3136 
3137  for (p = 0; p < (*m_exp)[n]->GetNtraces(); ++p)
3138  {
3139  offset = m_trace->GetPhys_Offset
3140  (elmtToTrace[n][p]->GetElmtId());
3141  (*m_exp)[n]->GetTracePhysVals(p,elmtToTrace[n][p],
3142  inarray + phys_offset,
3143  t_tmp = outarray +offset);
3144  }
3145  }
3146  }
3147  }
int GetPhys_Offset(int n) const
Get the start offset position for a global list of m_phys correspoinding to element n.
Definition: ExpList.h:2439
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:436

References ASSERTL1, Nektar::LibUtilities::eGauss_Lagrange, Nektar::MultiRegions::ExpList::GetExpSize(), Nektar::MultiRegions::ExpList::GetPhys_Offset(), m_locTraceToTraceMap, m_trace, m_traceMap, CellMLToNektar.cellml_metadata::p, and Vmath::Zero().

◆ v_FillBwdWithBoundCond()

void Nektar::MultiRegions::DisContField::v_FillBwdWithBoundCond ( const Array< OneD, NekDouble > &  Fwd,
Array< OneD, NekDouble > &  Bwd,
bool  PutFwdInBwdOnBCs 
)
protectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 2942 of file DisContField.cpp.

2946  {
2947  // Fill boundary conditions into missing elements
2948  if (PutFwdInBwdOnBCs) // just set Bwd value to be Fwd value on BCs
2949  {
2950  // Fill boundary conditions into missing elements
2951  for (int n = 0, cnt = 0; n < m_bndCondExpansions.size(); ++n)
2952  {
2953  if (m_bndConditions[n]->GetBoundaryConditionType() ==
2955  {
2956  auto ne = m_bndCondExpansions[n]->GetExpSize();
2957  for (int e = 0; e < ne; ++e)
2958  {
2959  auto npts = m_bndCondExpansions[n]->
2960  GetExp(e)->GetTotPoints();
2961  auto id2 = m_trace->GetPhys_Offset(m_traceMap->
2962  GetBndCondIDToGlobalTraceID(cnt+e));
2963  Vmath::Vcopy(npts, &Fwd[id2], 1, &Bwd[id2], 1);
2964  }
2965 
2966  cnt += ne;
2967  }
2968  else if (m_bndConditions[n]->GetBoundaryConditionType() ==
2970  m_bndConditions[n]->GetBoundaryConditionType() ==
2972  {
2973  auto ne = m_bndCondExpansions[n]->GetExpSize();
2974  for (int e = 0; e < ne; ++e)
2975  {
2976  auto npts = m_bndCondExpansions[n]->
2977  GetExp(e)->GetTotPoints();
2978  auto id2 = m_trace->GetPhys_Offset(m_traceMap->
2979  GetBndCondIDToGlobalTraceID(cnt+e));
2980  Vmath::Vcopy(npts, &Fwd[id2], 1, &Bwd[id2], 1);
2981  }
2982  cnt += ne;
2983  }
2984  else if (m_bndConditions[n]->GetBoundaryConditionType() !=
2986  {
2987  ASSERTL0(false,
2988  "Method not set up for this "
2989  "boundary condition.");
2990  }
2991  }
2992  }
2993  else
2994  {
2995  for (int n = 0, cnt = 0; n < m_bndCondExpansions.size(); ++n)
2996  {
2997  if (m_bndConditions[n]->GetBoundaryConditionType() ==
2999  {
3000  auto ne = m_bndCondExpansions[n]->GetExpSize();
3001  for (int e = 0; e < ne; ++e)
3002  {
3003  auto npts = m_bndCondExpansions[n]->GetExp(e)
3004  ->GetTotPoints();
3005  auto id1 = m_bndCondExpansions[n]
3006  ->GetPhys_Offset(e);
3007  auto id2 = m_trace->GetPhys_Offset(m_traceMap->
3008  GetBndCondIDToGlobalTraceID(cnt+e));
3009  Vmath::Vcopy(npts,
3010  &(m_bndCondExpansions[n]->GetPhys())
3011  [id1], 1, &Bwd[id2], 1);
3012  }
3013  cnt += ne;
3014  }
3015  else if (m_bndConditions[n]->GetBoundaryConditionType()
3017  m_bndConditions[n]->GetBoundaryConditionType()
3019  {
3020  auto ne = m_bndCondExpansions[n]->GetExpSize();
3021  for(int e = 0; e < ne; ++e)
3022  {
3023  auto npts = m_bndCondExpansions[n]->GetExp(e)
3024  ->GetTotPoints();
3025  auto id1 = m_bndCondExpansions[n]->
3026  GetPhys_Offset(e);
3028  GetPhys())[id1] == 0.0,
3029  "Method not set up for non-zero "
3030  "Neumann boundary condition");
3031  auto id2 = m_trace->GetPhys_Offset(
3032  m_traceMap->GetBndCondIDToGlobalTraceID(cnt+e));
3033  Vmath::Vcopy(npts, &Fwd[id2], 1, &Bwd[id2], 1);
3034  }
3035 
3036  cnt += ne;
3037  }
3038  else if (m_bndConditions[n]->GetBoundaryConditionType()
3040  {
3041  // Do nothing
3042  }
3043  else if (m_bndConditions[n]->GetBoundaryConditionType() !=
3045  {
3047  "Method not set up for this boundary "
3048  "condition.");
3049  }
3050  }
3051  }
3052  }
const Array< OneD, const NekDouble > & GetPhys() const
This function returns (a reference to) the array (implemented as m_phys) containing the function ev...
Definition: ExpList.h:2392
const std::shared_ptr< LocalRegions::ExpansionVector > GetExp() const
This function returns the vector of elements in the expansion.
Definition: ExpList.h:2422

References ASSERTL0, Nektar::SpatialDomains::eDirichlet, Nektar::ErrorUtil::efatal, Nektar::SpatialDomains::eNeumann, Nektar::SpatialDomains::eNotDefined, Nektar::SpatialDomains::ePeriodic, Nektar::SpatialDomains::eRobin, Nektar::MultiRegions::ExpList::GetExp(), Nektar::MultiRegions::ExpList::GetPhys(), Nektar::MultiRegions::ExpList::GetPhys_Offset(), m_bndCondExpansions, m_bndConditions, m_trace, m_traceMap, NEKERROR, and Vmath::Vcopy().

Referenced by v_GetFwdBwdTracePhys().

◆ v_FillBwdWithBwdWeight()

void Nektar::MultiRegions::DisContField::v_FillBwdWithBwdWeight ( Array< OneD, NekDouble > &  weightave,
Array< OneD, NekDouble > &  weightjmp 
)
protectedvirtual

Fill the weight with m_bndCondBndWeight.

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 3759 of file DisContField.cpp.

3762  {
3763  int cnt;
3764  int npts;
3765  int e = 0;
3766 
3767  // Fill boundary conditions into missing elements
3768  int id2 = 0;
3769 
3770  for (int n = cnt = 0; n < m_bndCondExpansions.size(); ++n)
3771  {
3772 
3773  if (m_bndConditions[n]->GetBoundaryConditionType() ==
3775  {
3776  for (e = 0; e < m_bndCondExpansions[n]->GetExpSize(); ++e)
3777  {
3778  npts = m_bndCondExpansions[n]->
3779  GetExp(e)->GetTotPoints();
3780  id2 = m_trace->GetPhys_Offset(m_traceMap->
3781  GetBndCondIDToGlobalTraceID(cnt+e));
3783  &weightave[id2], 1);
3784  Vmath::Fill(npts, 0.0, &weightjmp[id2], 1);
3785 
3786  }
3787 
3788  cnt += e;
3789  }
3790  else if (m_bndConditions[n]->GetBoundaryConditionType() ==
3792  m_bndConditions[n]->GetBoundaryConditionType() ==
3794  {
3795  for (e = 0; e < m_bndCondExpansions[n]->GetExpSize(); ++e)
3796  {
3797  npts = m_bndCondExpansions[n]->
3798  GetExp(e)->GetTotPoints();
3799  id2 = m_trace->GetPhys_Offset(m_traceMap->
3800  GetBndCondIDToGlobalTraceID(cnt+e));
3801  Vmath::Fill(npts,
3802  m_bndCondBndWeight[n],
3803  &weightave[id2], 1);
3804  Vmath::Fill(npts, 0.0, &weightjmp[id2], 1);
3805  }
3806 
3807  cnt += e;
3808  }
3809  else if (m_bndConditions[n]->GetBoundaryConditionType() !=
3811  {
3813  "Method not set up for this boundary condition.");
3814  }
3815  }
3816  }

References Nektar::SpatialDomains::eDirichlet, Nektar::ErrorUtil::efatal, Nektar::SpatialDomains::eNeumann, Nektar::SpatialDomains::ePeriodic, Nektar::SpatialDomains::eRobin, Vmath::Fill(), Nektar::MultiRegions::ExpList::GetExp(), m_bndCondBndWeight, m_bndCondExpansions, m_bndConditions, m_trace, m_traceMap, and NEKERROR.

◆ v_GetBndCondBwdWeight()

const Array< OneD, const NekDouble > & Nektar::MultiRegions::DisContField::v_GetBndCondBwdWeight ( )
inlineprotectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 399 of file DisContField.h.

400  {
401  return m_bndCondBndWeight;
402  }

References m_bndCondBndWeight.

◆ v_GetBndCondExpansions()

virtual const Array<OneD,const MultiRegions::ExpListSharedPtr>& Nektar::MultiRegions::DisContField::v_GetBndCondExpansions ( void  )
inlineprotectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 276 of file DisContField.h.

277  {
278  return m_bndCondExpansions;
279  }

References m_bndCondExpansions.

◆ v_GetBndConditions()

virtual const Array<OneD,const SpatialDomains::BoundaryConditionShPtr>& Nektar::MultiRegions::DisContField::v_GetBndConditions ( void  )
inlineprotectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Reimplemented in Nektar::MultiRegions::ContField.

Definition at line 282 of file DisContField.h.

283  {
284  return m_bndConditions;
285  }

References m_bndConditions.

◆ v_GetBndElmtExpansion()

void Nektar::MultiRegions::DisContField::v_GetBndElmtExpansion ( int  i,
std::shared_ptr< ExpList > &  result,
const bool  DeclareCoeffPhysArrays 
)
protectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 3971 of file DisContField.cpp.

3974  {
3975  int n, cnt, nq;
3976  int offsetOld, offsetNew;
3977  std::vector<unsigned int> eIDs;
3978 
3979  Array<OneD, int> ElmtID,TraceID;
3980  GetBoundaryToElmtMap(ElmtID,TraceID);
3981 
3982  // Skip other boundary regions
3983  for (cnt = n = 0; n < i; ++n)
3984  {
3985  cnt += m_bndCondExpansions[n]->GetExpSize();
3986  }
3987 
3988  // Populate eIDs with information from BoundaryToElmtMap
3989  for (n = 0; n < m_bndCondExpansions[i]->GetExpSize(); ++n)
3990  {
3991  eIDs.push_back(ElmtID[cnt+n]);
3992  }
3993 
3994  // Create expansion list
3995  result =
3997  (*this, eIDs, DeclareCoeffPhysArrays);
3998 
3999  // Copy phys and coeffs to new explist
4000  if( DeclareCoeffPhysArrays)
4001  {
4002  Array<OneD, NekDouble> tmp1, tmp2;
4003  for (n = 0; n < result->GetExpSize(); ++n)
4004  {
4005  nq = GetExp(ElmtID[cnt+n])->GetTotPoints();
4006  offsetOld = GetPhys_Offset(ElmtID[cnt+n]);
4007  offsetNew = result->GetPhys_Offset(n);
4008  Vmath::Vcopy(nq, tmp1 = GetPhys()+ offsetOld, 1,
4009  tmp2 = result->UpdatePhys()+ offsetNew, 1);
4010 
4011  nq = GetExp(ElmtID[cnt+n])->GetNcoeffs();
4012  offsetOld = GetCoeff_Offset(ElmtID[cnt+n]);
4013  offsetNew = result->GetCoeff_Offset(n);
4014  Vmath::Vcopy(nq, tmp1 = GetCoeffs()+ offsetOld, 1,
4015  tmp2 = result->UpdateCoeffs()+ offsetNew, 1);
4016  }
4017  }
4018  }
const Array< OneD, const NekDouble > & GetCoeffs() const
This function returns (a reference to) the array (implemented as m_coeffs) containing all local expa...
Definition: ExpList.h:2293

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::MultiRegions::ExpList::GetBoundaryToElmtMap(), Nektar::MultiRegions::ExpList::GetCoeff_Offset(), Nektar::MultiRegions::ExpList::GetCoeffs(), Nektar::MultiRegions::ExpList::GetExp(), Nektar::MultiRegions::ExpList::GetPhys(), Nektar::MultiRegions::ExpList::GetPhys_Offset(), m_bndCondExpansions, and Vmath::Vcopy().

◆ v_GetBoundaryToElmtMap()

void Nektar::MultiRegions::DisContField::v_GetBoundaryToElmtMap ( Array< OneD, int > &  ElmtID,
Array< OneD, int > &  TraceID 
)
protectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 3821 of file DisContField.cpp.

3824  {
3825 
3826  if (m_BCtoElmMap.size() == 0)
3827  {
3828  switch(m_expType)
3829  {
3830  case e1D:
3831  {
3832  map<int, int> VertGID;
3833  int i,n,id;
3834  int bid,cnt,Vid;
3835  int nbcs = m_bndConditions.size();
3836 
3837  // make sure arrays are of sufficient length
3838  m_BCtoElmMap = Array<OneD, int>(nbcs,-1);
3839  m_BCtoTraceMap = Array<OneD, int>(nbcs);
3840 
3841  // setup map of all global ids along boundary
3842  cnt = 0;
3843  for (n = 0; n < m_bndCondExpansions.size(); ++n)
3844  {
3845  Vid = m_bndCondExpansions[n]->GetExp(0)->
3846  GetGeom()->GetVertex(0)->GetVid();
3847  VertGID[Vid] = cnt++;
3848  }
3849 
3850  // Loop over elements and find verts that match;
3852  for(cnt = n = 0; n < GetExpSize(); ++n)
3853  {
3854  exp = (*m_exp)[n];
3855  for(i = 0; i < exp->GetNverts(); ++i)
3856  {
3857  id = exp->GetGeom()->GetVid(i);
3858 
3859  if (VertGID.count(id) > 0)
3860  {
3861  bid = VertGID.find(id)->second;
3862  ASSERTL1(m_BCtoElmMap[bid] == -1,
3863  "Edge already set");
3864  m_BCtoElmMap [bid] = n;
3865  m_BCtoTraceMap[bid] = i;
3866  cnt ++;
3867  }
3868  }
3869  }
3870  ASSERTL1(cnt == nbcs,
3871  "Failed to visit all boundary condtiions");
3872  }
3873  break;
3874  case e2D:
3875  {
3876  map<int, int> globalIdMap;
3877  int i,n;
3878  int cnt;
3879  int nbcs = 0;
3880 
3881  // Populate global ID map (takes global geometry
3882  // ID to local expansion list ID).
3883  for (i = 0; i < GetExpSize(); ++i)
3884  {
3885  globalIdMap[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
3886  }
3887 
3888  // Determine number of boundary condition expansions.
3889  for(i = 0; i < m_bndConditions.size(); ++i)
3890  {
3891  nbcs += m_bndCondExpansions[i]->GetExpSize();
3892  }
3893 
3894  // Initialize arrays
3895  m_BCtoElmMap = Array<OneD, int>(nbcs);
3896  m_BCtoTraceMap = Array<OneD, int>(nbcs);
3897 
3899  cnt = 0;
3900  for (n = 0; n < m_bndCondExpansions.size(); ++n)
3901  {
3902  for (i = 0; i < m_bndCondExpansions[n]->GetExpSize();
3903  ++i, ++cnt)
3904  {
3905  exp1d = m_bndCondExpansions[n]->GetExp(i)->
3906  as<LocalRegions::Expansion1D>();
3907 
3908  // Use edge to element map from MeshGraph.
3910  m_graph->GetElementsFromEdge(exp1d->GetGeom1D());
3911 
3912  m_BCtoElmMap[cnt] = globalIdMap[
3913  (*tmp)[0].first->GetGlobalID()];
3914  m_BCtoTraceMap[cnt] = (*tmp)[0].second;
3915  }
3916  }
3917  }
3918  break;
3919  case e3D:
3920  {
3921  map<int,int> globalIdMap;
3922  int i, n;
3923  int cnt;
3924  int nbcs = 0;
3925 
3926  // Populate global ID map (takes global geometry ID to local
3927  // expansion list ID).
3929  for (i = 0; i < GetExpSize(); ++i)
3930  {
3931  globalIdMap[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
3932  }
3933 
3934  // Determine number of boundary condition expansions.
3935  for(i = 0; i < m_bndConditions.size(); ++i)
3936  {
3937  nbcs += m_bndCondExpansions[i]->GetExpSize();
3938  }
3939 
3940  // Initialize arrays
3941  m_BCtoElmMap = Array<OneD, int>(nbcs);
3942  m_BCtoTraceMap = Array<OneD, int>(nbcs);
3943 
3945  for(cnt = n = 0; n < m_bndCondExpansions.size(); ++n)
3946  {
3947  for(i = 0; i < m_bndCondExpansions[n]->GetExpSize(); ++i, ++cnt)
3948  {
3949  exp2d = m_bndCondExpansions[n]->GetExp(i)->
3950  as<LocalRegions::Expansion2D>();
3951 
3953  m_graph->GetElementsFromFace(exp2d->GetGeom2D());
3954  m_BCtoElmMap[cnt] = globalIdMap
3955  [tmp->at(0).first->GetGlobalID()];
3956  m_BCtoTraceMap[cnt] = tmp->at(0).second;
3957  }
3958  }
3959  }
3960  break;
3961  default:
3962  ASSERTL1(false,"Not setup for this expansion");
3963  break;
3964  }
3965  }
3966 
3967  ElmtID = m_BCtoElmMap;
3968  TraceID = m_BCtoTraceMap;
3969  }
Array< OneD, int > m_BCtoTraceMap
Definition: DisContField.h:59
std::shared_ptr< Expansion2D > Expansion2DSharedPtr
Definition: Expansion1D.h:47
std::shared_ptr< Expansion1D > Expansion1DSharedPtr
Definition: Expansion1D.h:51
std::shared_ptr< Expansion3D > Expansion3DSharedPtr
Definition: Expansion2D.h:49

References ASSERTL1, Nektar::MultiRegions::e1D, Nektar::MultiRegions::e2D, Nektar::MultiRegions::e3D, Nektar::MultiRegions::ExpList::GetExpSize(), m_BCtoElmMap, m_BCtoTraceMap, m_bndCondExpansions, m_bndConditions, Nektar::MultiRegions::ExpList::m_expType, and Nektar::MultiRegions::ExpList::m_graph.

◆ v_GetFwdBwdTracePhys() [1/2]

void Nektar::MultiRegions::DisContField::v_GetFwdBwdTracePhys ( Array< OneD, NekDouble > &  Fwd,
Array< OneD, NekDouble > &  Bwd 
)
inlineprotectedvirtual

Generate the forward or backward state for each trace point.

Parameters
FwdForward state.
BwdBackward state.

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 392 of file DisContField.h.

394  {
395  v_GetFwdBwdTracePhys(m_phys,Fwd,Bwd);
396  }
virtual void v_GetFwdBwdTracePhys(Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd)
Definition: DisContField.h:392

References Nektar::MultiRegions::ExpList::m_phys.

◆ v_GetFwdBwdTracePhys() [2/2]

void Nektar::MultiRegions::DisContField::v_GetFwdBwdTracePhys ( const Array< OneD, const NekDouble > &  field,
Array< OneD, NekDouble > &  Fwd,
Array< OneD, NekDouble > &  Bwd,
bool  FillBnd = true,
bool  PutFwdInBwdOnBCs = false,
bool  DoExchange = true 
)
protectedvirtual

This method extracts the "forward" and "backward" trace data from the array field and puts the data into output vectors Fwd and Bwd.

We first define the convention which defines "forwards" and "backwards". First an association is made between the vertex/edge/face of each element and its corresponding vertex/edge/face in the trace space using the mapping m_traceMap. The element can either be left-adjacent or right-adjacent to this trace face (see Expansion2D::GetLeftAdjacentElementExp). Boundary faces are always left-adjacent since left-adjacency is populated first.

If the element is left-adjacent we extract the trace data from field into the forward trace space Fwd; otherwise, we place it in the backwards trace space Bwd. In this way, we form a unique set of trace normals since these are always extracted from left-adjacent elements.

Parameters
fieldis a NekDouble array which contains the fielddata from which we wish to extract the backward and forward orientated trace/face arrays.
Returns
Updates a NekDouble array Fwd and Bwd

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 2871 of file DisContField.cpp.

2878  {
2879  // Is this zeroing necessary?
2880  // Zero forward/backward vectors.
2881  Vmath::Zero(Fwd.size(), Fwd, 1);
2882  Vmath::Zero(Bwd.size(), Bwd, 1);
2883 
2884  // Basis definition on each element
2885  LibUtilities::BasisSharedPtr basis = (*m_exp)[0]->GetBasis(0);
2886  if((basis->GetBasisType() != LibUtilities::eGauss_Lagrange))
2887  {
2888  // blocked routine
2889  Array<OneD, NekDouble> tracevals(m_locTraceToTraceMap->
2890  GetNLocTracePts());
2891 
2892  m_locTraceToTraceMap->LocTracesFromField(field, tracevals);
2893  m_locTraceToTraceMap->InterpLocTracesToTrace(0, tracevals, Fwd);
2894 
2895  Array<OneD, NekDouble> invals = tracevals + m_locTraceToTraceMap->
2896  GetNFwdLocTracePts();
2897  m_locTraceToTraceMap->InterpLocTracesToTrace(1, invals, Bwd);
2898  }
2899  else
2900  {
2901  // Loop over elements and collect forward expansion
2902  auto nexp = GetExpSize();
2903  Array<OneD,NekDouble> e_tmp;
2905 
2906  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
2907  &elmtToTrace = m_traceMap->GetElmtToTrace();
2908 
2909  for (int n = 0, cnt = 0; n < nexp; ++n)
2910  {
2911  exp = (*m_exp)[n];
2912  auto phys_offset = GetPhys_Offset(n);
2913 
2914  for(int e = 0; e < exp->GetNtraces(); ++e, ++cnt)
2915  {
2916  auto offset = m_trace->GetPhys_Offset(
2917  elmtToTrace[n][e]->GetElmtId());
2918 
2919  e_tmp = (m_leftAdjacentTraces[cnt])? Fwd + offset:
2920  Bwd + offset;
2921 
2922  exp->GetTracePhysVals(e, elmtToTrace[n][e],
2923  field + phys_offset, e_tmp);
2924  }
2925  }
2926  }
2927 
2929 
2930  if (FillBnd)
2931  {
2932  v_FillBwdWithBoundCond(Fwd, Bwd, PutFwdInBwdOnBCs);
2933  }
2934 
2935  if(DoExchange)
2936  {
2937  // Do parallel exchange for forwards/backwards spaces.
2938  m_traceMap->GetAssemblyCommDG()->PerformExchange(Fwd, Bwd);
2939  }
2940  }
virtual void v_PeriodicBwdCopy(const Array< OneD, const NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd)
Definition: DisContField.h:411
virtual void v_FillBwdWithBoundCond(const Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd, bool PutFwdInBwdOnBCs)

References Nektar::LibUtilities::eGauss_Lagrange, Nektar::MultiRegions::ExpList::GetExpSize(), Nektar::MultiRegions::ExpList::GetPhys_Offset(), m_leftAdjacentTraces, m_locTraceToTraceMap, m_trace, m_traceMap, v_FillBwdWithBoundCond(), v_PeriodicBwdCopy(), and Vmath::Zero().

◆ v_GetLocTraceFromTracePts()

void Nektar::MultiRegions::DisContField::v_GetLocTraceFromTracePts ( const Array< OneD, const NekDouble > &  Fwd,
const Array< OneD, const NekDouble > &  Bwd,
Array< OneD, NekDouble > &  locTraceFwd,
Array< OneD, NekDouble > &  locTraceBwd 
)
protectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 4385 of file DisContField.cpp.

4390  {
4391  if (NullNekDouble1DArray != locTraceBwd)
4392  {
4393  switch(m_expType)
4394  {
4395  case e2D:
4396  m_locTraceToTraceMap->RightIPTWLocEdgesToTraceInterpMat(
4397  1, Bwd, locTraceBwd);
4398  break;
4399  case e3D:
4400  m_locTraceToTraceMap->RightIPTWLocFacesToTraceInterpMat(
4401  1, Bwd, locTraceBwd);
4402  break;
4403  default:
4405  "GetLocTraceFromTracePts not defined");
4406  }
4407  }
4408 
4409  if (NullNekDouble1DArray != locTraceFwd)
4410  {
4411  switch(m_expType)
4412  {
4413  case e2D:
4414  m_locTraceToTraceMap->RightIPTWLocEdgesToTraceInterpMat(
4415  0, Fwd, locTraceFwd);
4416  break;
4417  case e3D:
4418  m_locTraceToTraceMap->RightIPTWLocFacesToTraceInterpMat(
4419  0, Fwd, locTraceFwd);
4420  break;
4421  default:
4423  "GetLocTraceFromTracePts not defined");
4424  }
4425  }
4426  }
static Array< OneD, NekDouble > NullNekDouble1DArray

References Nektar::MultiRegions::e2D, Nektar::MultiRegions::e3D, Nektar::ErrorUtil::efatal, Nektar::MultiRegions::ExpList::m_expType, m_locTraceToTraceMap, NEKERROR, and Nektar::NullNekDouble1DArray.

◆ v_GetLocTraceToTraceMap()

virtual const LocTraceToTraceMapSharedPtr& Nektar::MultiRegions::DisContField::v_GetLocTraceToTraceMap ( void  ) const
inlineprotectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 228 of file DisContField.h.

229  {
230  return m_locTraceToTraceMap;
231  }

References m_locTraceToTraceMap.

◆ v_GetPeriodicEntities()

virtual void Nektar::MultiRegions::DisContField::v_GetPeriodicEntities ( PeriodicMap periodicVerts,
PeriodicMap periodicEdges,
PeriodicMap periodicFaces 
)
inlineprotectedvirtual

Obtain a copy of the periodic edges and vertices for this field.

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 361 of file DisContField.h.

365  {
366  periodicVerts = m_periodicVerts;
367  periodicEdges = m_periodicEdges;
368  periodicFaces = m_periodicFaces;
369  }

References m_periodicEdges, m_periodicFaces, and m_periodicVerts.

◆ v_GetRobinBCInfo()

map< int, RobinBCInfoSharedPtr > Nektar::MultiRegions::DisContField::v_GetRobinBCInfo ( void  )
protectedvirtual

Search through the edge expansions and identify which ones have Robin/Mixed type boundary conditions. If find a Robin boundary then store the edge id of the boundary condition and the array of points of the physical space boundary condition which are hold the boundary condition primitive variable coefficient at the quatrature points

Returns
std map containing the robin boundary condition info using a key of the element id

There is a next member to allow for more than one Robin boundary condition per element

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 4048 of file DisContField.cpp.

4049  {
4050  int i,cnt;
4051  map<int, RobinBCInfoSharedPtr> returnval;
4052  Array<OneD, int> ElmtID,TraceID;
4053  GetBoundaryToElmtMap(ElmtID,TraceID);
4054 
4055  for(cnt = i = 0; i < m_bndCondExpansions.size(); ++i)
4056  {
4057  MultiRegions::ExpListSharedPtr locExpList;
4058 
4059  if(m_bndConditions[i]->GetBoundaryConditionType() ==
4061  {
4062  int e,elmtid;
4063  Array<OneD, NekDouble> Array_tmp;
4064 
4065  locExpList = m_bndCondExpansions[i];
4066 
4067  int npoints = locExpList->GetNpoints();
4068  Array<OneD, NekDouble> x0(npoints, 0.0);
4069  Array<OneD, NekDouble> x1(npoints, 0.0);
4070  Array<OneD, NekDouble> x2(npoints, 0.0);
4071  Array<OneD, NekDouble> coeffphys(npoints);
4072 
4073  locExpList->GetCoords(x0, x1, x2);
4074 
4075  LibUtilities::Equation coeffeqn =
4076  std::static_pointer_cast<
4077  SpatialDomains::RobinBoundaryCondition>
4078  (m_bndConditions[i])->m_robinPrimitiveCoeff;
4079 
4080  // evalaute coefficient
4081  coeffeqn.Evaluate(x0, x1, x2, 0.0, coeffphys);
4082 
4083  for(e = 0; e < locExpList->GetExpSize(); ++e)
4084  {
4085  RobinBCInfoSharedPtr rInfo =
4088  TraceID[cnt+e],
4089  Array_tmp = coeffphys +
4090  locExpList->GetPhys_Offset(e));
4091 
4092  elmtid = ElmtID[cnt+e];
4093  // make link list if necessary
4094  if(returnval.count(elmtid) != 0)
4095  {
4096  rInfo->next = returnval.find(elmtid)->second;
4097  }
4098  returnval[elmtid] = rInfo;
4099  }
4100  }
4101  cnt += m_bndCondExpansions[i]->GetExpSize();
4102  }
4103 
4104  return returnval;
4105  }
std::shared_ptr< RobinBCInfo > RobinBCInfoSharedPtr

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::SpatialDomains::eRobin, Nektar::LibUtilities::Equation::Evaluate(), Nektar::MultiRegions::ExpList::GetBoundaryToElmtMap(), m_bndCondExpansions, and m_bndConditions.

◆ v_GetTrace()

virtual ExpListSharedPtr& Nektar::MultiRegions::DisContField::v_GetTrace ( )
inlineprotectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 212 of file DisContField.h.

213  {
215  {
216  SetUpDG();
217  }
218 
219  return m_trace;
220  }

References m_trace, Nektar::MultiRegions::NullExpListSharedPtr, and SetUpDG().

◆ v_GetTraceMap()

virtual AssemblyMapDGSharedPtr& Nektar::MultiRegions::DisContField::v_GetTraceMap ( void  )
inlineprotectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 222 of file DisContField.h.

223  {
224  return m_traceMap;
225  }

References m_traceMap.

◆ v_HelmSolve()

void Nektar::MultiRegions::DisContField::v_HelmSolve ( const Array< OneD, const NekDouble > &  inarray,
Array< OneD, NekDouble > &  outarray,
const StdRegions::ConstFactorMap factors,
const StdRegions::VarCoeffMap varcoeff,
const MultiRegions::VarFactorsMap varfactors,
const Array< OneD, const NekDouble > &  dirForcing,
const bool  PhysSpaceForcing 
)
protectedvirtual

Solve the Helmholtz equation.

Reimplemented from Nektar::MultiRegions::ExpList.

Reimplemented in Nektar::MultiRegions::ContField.

Definition at line 3379 of file DisContField.cpp.

3387  {
3388  boost::ignore_unused(varfactors,dirForcing);
3389  int i,n,cnt,nbndry;
3390  int nexp = GetExpSize();
3391 
3392  Array<OneD,NekDouble> f(m_ncoeffs);
3393  DNekVec F(m_ncoeffs,f,eWrapper);
3394  Array<OneD,NekDouble> e_f, e_l;
3395 
3396  //----------------------------------
3397  // Setup RHS Inner product if required
3398  //----------------------------------
3399  if(PhysSpaceForcing)
3400  {
3401  IProductWRTBase(inarray,f);
3402  Vmath::Neg(m_ncoeffs,f,1);
3403  }
3404  else
3405  {
3406  Vmath::Smul(m_ncoeffs,-1.0,inarray,1,f,1);
3407  }
3408 
3409  //----------------------------------
3410  // Solve continuous Boundary System
3411  //----------------------------------
3412  int GloBndDofs = m_traceMap->GetNumGlobalBndCoeffs();
3413  int NumDirBCs = m_traceMap->GetNumLocalDirBndCoeffs();
3414  int e_ncoeffs;
3415 
3416  GlobalMatrixKey HDGLamToUKey(StdRegions::eHybridDGLamToU,
3417  NullAssemblyMapSharedPtr,factors,varcoeff);
3418  const DNekScalBlkMatSharedPtr &HDGLamToU =
3419  GetBlockMatrix(HDGLamToUKey);
3420 
3421  // Retrieve number of local trace space coefficients N_{\lambda},
3422  // and set up local elemental trace solution \lambda^e.
3423  int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs();
3424  Array<OneD, NekDouble> bndrhs(LocBndCoeffs,0.0);
3425  Array<OneD, NekDouble> loclambda(LocBndCoeffs,0.0);
3426  DNekVec LocLambda(LocBndCoeffs,loclambda,eWrapper);
3427 
3428  //----------------------------------
3429  // Evaluate Trace Forcing
3430  // Kirby et al, 2010, P23, Step 5.
3431  //----------------------------------
3432  // Determing <u_lam,f> terms using HDGLamToU matrix
3433  for (cnt = n = 0; n < nexp; ++n)
3434  {
3435  nbndry = (*m_exp)[n]->NumDGBndryCoeffs();
3436 
3437  e_ncoeffs = (*m_exp)[n]->GetNcoeffs();
3438  e_f = f + m_coeff_offset[n];
3439  e_l = bndrhs + cnt;
3440 
3441  // use outarray as tmp space
3442  DNekVec Floc (nbndry, e_l, eWrapper);
3443  DNekVec ElmtFce (e_ncoeffs, e_f, eWrapper);
3444  Floc = Transpose(*(HDGLamToU->GetBlock(n,n)))*ElmtFce;
3445 
3446  cnt += nbndry;
3447  }
3448 
3449  Array<OneD, const int> bndCondMap =
3450  m_traceMap->GetBndCondCoeffsToLocalTraceMap();
3451  Array<OneD, const NekDouble> Sign =
3452  m_traceMap->GetLocalToGlobalBndSign();
3453 
3454  // Copy Dirichlet boundary conditions and weak forcing
3455  // into trace space
3456  int locid;
3457  cnt = 0;
3458  for(i = 0; i < m_bndCondExpansions.size(); ++i)
3459  {
3460  Array<OneD, const NekDouble> bndcoeffs =
3461  m_bndCondExpansions[i]->GetCoeffs();
3462 
3463  if(m_bndConditions[i]->GetBoundaryConditionType() ==
3465  {
3466  for(int j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); ++j)
3467  {
3468  locid = bndCondMap[cnt + j];
3469  loclambda[locid] = Sign[locid]*bndcoeffs[j];
3470  }
3471  }
3472  else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3474  m_bndConditions[i]->GetBoundaryConditionType() ==
3476  {
3477  //Add weak boundary condition to trace forcing
3478  for(int j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); ++j)
3479  {
3480  locid = bndCondMap[cnt + j];
3481  bndrhs[locid] += Sign[locid]*bndcoeffs[j];
3482  }
3483  }
3484 
3485  cnt += (m_bndCondExpansions[i])->GetNcoeffs();
3486  }
3487 
3488  //----------------------------------
3489  // Solve trace problem: \Lambda = K^{-1} F
3490  // K is the HybridDGHelmBndLam matrix.
3491  //----------------------------------
3492  if(GloBndDofs - NumDirBCs > 0)
3493  {
3494  GlobalLinSysKey key(StdRegions::eHybridDGHelmBndLam,
3495  m_traceMap,factors,varcoeff);
3497  LinSys->Solve(bndrhs,loclambda,m_traceMap);
3498 
3499  // For consistency with previous version put global
3500  // solution into m_trace->m_coeffs
3501  m_traceMap->LocalToGlobal(loclambda,m_trace->UpdateCoeffs());
3502  }
3503 
3504  //----------------------------------
3505  // Internal element solves
3506  //----------------------------------
3507  GlobalMatrixKey invHDGhelmkey(StdRegions::eInvHybridDGHelmholtz,
3509  factors, varcoeff);
3510 
3511  const DNekScalBlkMatSharedPtr& InvHDGHelm =
3512  GetBlockMatrix(invHDGhelmkey);
3513  DNekVec out(m_ncoeffs,outarray,eWrapper);
3514  Vmath::Zero(m_ncoeffs,outarray,1);
3515 
3516  // out = u_f + u_lam = (*InvHDGHelm)*f + (LamtoU)*Lam
3517  out = (*InvHDGHelm)*F + (*HDGLamToU)*LocLambda;
3518  }
GlobalLinSysSharedPtr GetGlobalBndLinSys(const GlobalLinSysKey &mkey)
For a given key, returns the associated global linear system.
void IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.h:1956
int GetNcoeffs(void) const
Returns the total number of local degrees of freedom .
Definition: ExpList.h:1800
const DNekScalBlkMatSharedPtr & GetBlockMatrix(const GlobalMatrixKey &gkey)
Definition: ExpList.cpp:1946
static AssemblyMapSharedPtr NullAssemblyMapSharedPtr
Definition: AssemblyMap.h:54
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:73
NekVector< NekDouble > DNekVec
Definition: NekTypeDefs.hpp:48
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:461
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:225

References Nektar::SpatialDomains::eDirichlet, Nektar::StdRegions::eHybridDGHelmBndLam, Nektar::StdRegions::eHybridDGLamToU, Nektar::StdRegions::eInvHybridDGHelmholtz, Nektar::SpatialDomains::eNeumann, Nektar::SpatialDomains::eRobin, Nektar::eWrapper, Nektar::MultiRegions::ExpList::GetBlockMatrix(), Nektar::MultiRegions::ExpList::GetExpSize(), GetGlobalBndLinSys(), Nektar::MultiRegions::ExpList::GetNcoeffs(), Nektar::MultiRegions::ExpList::IProductWRTBase(), m_bndCondExpansions, m_bndConditions, Nektar::MultiRegions::ExpList::m_coeff_offset, Nektar::MultiRegions::ExpList::m_ncoeffs, m_trace, m_traceMap, Vmath::Neg(), Nektar::MultiRegions::NullAssemblyMapSharedPtr, Vmath::Smul(), Nektar::Transpose(), and Vmath::Zero().

◆ v_PeriodicBwdCopy()

void Nektar::MultiRegions::DisContField::v_PeriodicBwdCopy ( const Array< OneD, const NekDouble > &  Fwd,
Array< OneD, NekDouble > &  Bwd 
)
inlineprotectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 411 of file DisContField.h.

414  {
415  for (int n = 0; n < m_periodicFwdCopy.size(); ++n)
416  {
417  Bwd[m_periodicBwdCopy[n]] = Fwd[m_periodicFwdCopy[n]];
418  }
419  }

References m_periodicBwdCopy, and m_periodicFwdCopy.

Referenced by v_GetFwdBwdTracePhys().

◆ v_Reset()

void Nektar::MultiRegions::DisContField::v_Reset ( )
protectedvirtual

Reset this field, so that geometry information can be updated.

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 4023 of file DisContField.cpp.

4024  {
4025  ExpList::v_Reset();
4026 
4027  // Reset boundary condition expansions.
4028  for (int n = 0; n < m_bndCondExpansions.size(); ++n)
4029  {
4030  m_bndCondExpansions[n]->Reset();
4031  }
4032  }
virtual void v_Reset()
Reset geometry information, metrics, matrix managers and geometry information.
Definition: ExpList.cpp:2580

References m_bndCondExpansions, and Nektar::MultiRegions::ExpList::v_Reset().

◆ v_SetBndCondBwdWeight()

void Nektar::MultiRegions::DisContField::v_SetBndCondBwdWeight ( const int  index,
const NekDouble  value 
)
inlineprotectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 404 of file DisContField.h.

407  {
408  m_bndCondBndWeight[index] = value;
409  }

References m_bndCondBndWeight.

◆ v_UpdateBndCondExpansion()

virtual MultiRegions::ExpListSharedPtr& Nektar::MultiRegions::DisContField::v_UpdateBndCondExpansion ( int  i)
inlineprotectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 288 of file DisContField.h.

289  {
290  return m_bndCondExpansions[i];
291  }

References m_bndCondExpansions.

◆ v_UpdateBndConditions()

virtual Array<OneD, SpatialDomains::BoundaryConditionShPtr>& Nektar::MultiRegions::DisContField::v_UpdateBndConditions ( )
inlineprotectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 294 of file DisContField.h.

295  {
296  return m_bndConditions;
297  }

References m_bndConditions.

Member Data Documentation

◆ m_BCtoElmMap

Array<OneD, int> Nektar::MultiRegions::DisContField::m_BCtoElmMap

Definition at line 58 of file DisContField.h.

Referenced by v_GetBoundaryToElmtMap().

◆ m_BCtoTraceMap

Array<OneD, int> Nektar::MultiRegions::DisContField::m_BCtoTraceMap

Definition at line 59 of file DisContField.h.

Referenced by v_GetBoundaryToElmtMap().

◆ m_bndCondBndWeight

Array<OneD,NekDouble> Nektar::MultiRegions::DisContField::m_bndCondBndWeight
protected

◆ m_bndCondExpansions

Array<OneD,MultiRegions::ExpListSharedPtr> Nektar::MultiRegions::DisContField::m_bndCondExpansions
protected

◆ m_bndConditions

Array<OneD,SpatialDomains::BoundaryConditionShPtr> Nektar::MultiRegions::DisContField::m_bndConditions
protected

◆ m_boundaryTraces

std::set<int> Nektar::MultiRegions::DisContField::m_boundaryTraces
protected

A set storing the global IDs of any boundary Verts.

Definition at line 162 of file DisContField.h.

Referenced by DisContField(), IsLeftAdjacentTrace(), and SetUpDG().

◆ m_globalBndMat

GlobalLinSysMapShPtr Nektar::MultiRegions::DisContField::m_globalBndMat
protected

Global boundary matrix.

Definition at line 151 of file DisContField.h.

Referenced by DisContField(), GetGlobalBndLinSys(), and SetUpDG().

◆ m_leftAdjacentTraces

std::vector<bool> Nektar::MultiRegions::DisContField::m_leftAdjacentTraces
protected

Definition at line 191 of file DisContField.h.

Referenced by DisContField(), SetUpDG(), and v_GetFwdBwdTracePhys().

◆ m_locTraceToTraceMap

LocTraceToTraceMapSharedPtr Nektar::MultiRegions::DisContField::m_locTraceToTraceMap
protected

Map of local trace (the points at the edge,face of the element) to the trace space discretisation

Definition at line 197 of file DisContField.h.

Referenced by DisContField(), SetUpDG(), v_AddFwdBwdTraceIntegral(), v_AddTraceIntegral(), v_AddTraceIntegralToOffDiag(), v_AddTraceQuadPhysToField(), v_ExtractTracePhys(), v_GetFwdBwdTracePhys(), v_GetLocTraceFromTracePts(), and v_GetLocTraceToTraceMap().

◆ m_negatedFluxNormal

std::vector<bool> Nektar::MultiRegions::DisContField::m_negatedFluxNormal
private

Definition at line 373 of file DisContField.h.

Referenced by GetNegatedFluxNormal().

◆ m_numDirBndCondExpansions

int Nektar::MultiRegions::DisContField::m_numDirBndCondExpansions
protected

The number of boundary segments on which Dirichlet boundary conditions are imposed.

Definition at line 129 of file DisContField.h.

◆ m_periodicBwdCopy

std::vector<int> Nektar::MultiRegions::DisContField::m_periodicBwdCopy
protected

Definition at line 185 of file DisContField.h.

Referenced by DisContField(), SetUpDG(), and v_PeriodicBwdCopy().

◆ m_periodicEdges

PeriodicMap Nektar::MultiRegions::DisContField::m_periodicEdges
protected

A map which identifies pairs of periodic edges.

Definition at line 172 of file DisContField.h.

Referenced by Nektar::MultiRegions::ContField::ContField(), DisContField(), FindPeriodicTraces(), IsLeftAdjacentTrace(), SetUpDG(), and v_GetPeriodicEntities().

◆ m_periodicFaces

PeriodicMap Nektar::MultiRegions::DisContField::m_periodicFaces
protected

A map which identifies pairs of periodic faces.

Definition at line 177 of file DisContField.h.

Referenced by Nektar::MultiRegions::ContField::ContField(), DisContField(), FindPeriodicTraces(), IsLeftAdjacentTrace(), SetUpDG(), and v_GetPeriodicEntities().

◆ m_periodicFwdCopy

std::vector<int> Nektar::MultiRegions::DisContField::m_periodicFwdCopy
protected

A vector indicating degress of freedom which need to be copied from forwards to backwards space in case of a periodic boundary condition.

Definition at line 184 of file DisContField.h.

Referenced by DisContField(), SetUpDG(), and v_PeriodicBwdCopy().

◆ m_periodicVerts

PeriodicMap Nektar::MultiRegions::DisContField::m_periodicVerts
protected

A map which identifies groups of periodic vertices.

Definition at line 167 of file DisContField.h.

Referenced by Nektar::MultiRegions::ContField::ContField(), DisContField(), FindPeriodicTraces(), IsLeftAdjacentTrace(), SetUpDG(), and v_GetPeriodicEntities().

◆ m_trace

ExpListSharedPtr Nektar::MultiRegions::DisContField::m_trace
protected

◆ m_traceMap

AssemblyMapDGSharedPtr Nektar::MultiRegions::DisContField::m_traceMap
protected