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, const std::string bcvariable="NotSet")
 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, const LibUtilities::CommSharedPtr &comm, 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...
 
void GetLocTraceToTraceMap (LocTraceToTraceMapSharedPtr &LocTraceToTraceMap)
 
- 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 LibUtilities::CommSharedPtr &comm, 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 LocalRegions::ExpansionVector &locexp, const SpatialDomains::MeshGraphSharedPtr &graph, const bool DeclareCoeffPhysArrays, const std::string variable, const Collections::ImplementationType ImpType=Collections::eNoImpType)
 Generate an trace ExpList from a meshgraph graph and session file. 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 (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 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 FwdTransLocalElmt (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 FwdTransBndConstrained (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
void BwdTrans (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 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 GetCoords (const int eid, Array< OneD, NekDouble > &coord_0, Array< OneD, NekDouble > &coord_1=NullNekDouble1DArray, Array< OneD, NekDouble > &coord_2=NullNekDouble1DArray)
 
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 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 () const
 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
 
std::vector< bool > & GetLeftAdjacentTraces (void)
 

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 std::vector< bool > & v_GetLeftAdjacentTraces (void)
 
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)
 This routine determines if an element is to the "left" of the adjacent trace, which arises from the idea there is a local normal direction between two elements (i.e. on the trace) and one elements would then be the left. More...
 
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_FwdTrans (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
virtual void v_FwdTransLocalElmt (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
virtual void v_FwdTransBndConstrained (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_GetCoords (Array< OneD, NekDouble > &coord_0, Array< OneD, NekDouble > &coord_1, Array< OneD, NekDouble > &coord_2=NullNekDouble1DArray)
 
virtual void v_GetCoords (const int eid, Array< OneD, NekDouble > &xc0, Array< OneD, NekDouble > &xc1, Array< OneD, NekDouble > &xc2)
 
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 ()
 : Set up a normal along the trace elements between two elements at elemental level More...
 
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 64 of file DisContField.cpp.

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

◆ 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,
const std::string  bcvariable = "NotSet" 
)

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 82 of file DisContField.cpp.

88  : ExpList(pSession, graph, DeclareCoeffPhysArrays, variable, ImpType),
90 {
91  std::string bcvar;
92  if (bcvariable == "NotSet")
93  {
94  bcvar = variable;
95  }
96  else
97  {
98  bcvar = bcvariable;
99  }
100 
101  if (bcvar.compare("DefaultVar") != 0)
102  {
103  SpatialDomains::BoundaryConditions bcs(m_session, graph);
104 
105  GenerateBoundaryConditionExpansion(graph, bcs, bcvar,
106  DeclareCoeffPhysArrays);
107  if (DeclareCoeffPhysArrays)
108  {
109  EvaluateBoundaryConditions(0.0, bcvar);
110  }
111  ApplyGeomInfo();
112  FindPeriodicTraces(bcs, bcvar);
113  }
114 
115  if (SetUpJustDG)
116  {
117  SetUpDG(variable);
118  }
119  else
120  {
121  int i, cnt;
122  Array<OneD, int> ElmtID, TraceID;
123  GetBoundaryToElmtMap(ElmtID, TraceID);
124 
125  for (cnt = i = 0; i < m_bndCondExpansions.size(); ++i)
126  {
128  locExpList = m_bndCondExpansions[i];
129 
130  for (int e = 0; e < locExpList->GetExpSize(); ++e)
131  {
132  (*m_exp)[ElmtID[cnt + e]]->SetTraceExp(TraceID[cnt + e],
133  locExpList->GetExp(e));
134  locExpList->GetExp(e)->SetAdjacentElementExp(
135  TraceID[cnt + e], (*m_exp)[ElmtID[cnt + e]]);
136  }
137  cnt += m_bndCondExpansions[i]->GetExpSize();
138  }
139 
140  if (m_session->DefinesSolverInfo("PROJECTION"))
141  {
142  std::string ProjectStr = m_session->GetSolverInfo("PROJECTION");
143  if ((ProjectStr == "MixedCGDG") ||
144  (ProjectStr == "Mixed_CG_Discontinuous"))
145  {
146  SetUpDG();
147  }
148  else
149  {
151  }
152  }
153  else
154  {
156  }
157  }
158 }
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:2450
std::shared_ptr< LocalRegions::ExpansionVector > m_exp
The list of local expansions.
Definition: ExpList.h:1196
void EvaluateBoundaryConditions(const NekDouble time=0.0, const std::string varName="", const NekDouble=NekConstants::kNekUnsetDouble, const NekDouble=NekConstants::kNekUnsetDouble)
Definition: ExpList.h:2437
void ApplyGeomInfo()
Apply geometry information to each expansion.
Definition: ExpList.cpp:2658
LibUtilities::SessionReaderSharedPtr m_session
Session.
Definition: ExpList.h:1129
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,
const LibUtilities::CommSharedPtr comm,
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 605 of file DisContField.cpp.

613  : ExpList(pSession, domain, graph1D, true, variable, SetToOneSpaceDimension,
614  comm, ImpType)
615 {
616  if (variable.compare("DefaultVar") != 0)
617  {
619  GetDomainBCs(domain, Allbcs, variable);
620 
621  GenerateBoundaryConditionExpansion(m_graph, *DomBCs, variable);
622  EvaluateBoundaryConditions(0.0, variable);
623  ApplyGeomInfo();
624  FindPeriodicTraces(*DomBCs, variable);
625  }
626 
627  SetUpDG(variable);
628 }
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:1132
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 634 of file DisContField.cpp.

636  : ExpList(In, DeclareCoeffPhysArrays),
637  m_bndCondExpansions(In.m_bndCondExpansions),
638  m_bndConditions(In.m_bndConditions), m_globalBndMat(In.m_globalBndMat),
639  m_traceMap(In.m_traceMap), m_boundaryTraces(In.m_boundaryTraces),
640  m_periodicVerts(In.m_periodicVerts),
641  m_periodicFwdCopy(In.m_periodicFwdCopy),
642  m_periodicBwdCopy(In.m_periodicBwdCopy),
643  m_leftAdjacentTraces(In.m_leftAdjacentTraces),
644  m_locTraceToTraceMap(In.m_locTraceToTraceMap)
645 {
646  if (In.m_trace)
647  {
649  *In.m_trace, DeclareCoeffPhysArrays);
650  }
651 }
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:186
std::set< int > m_boundaryTraces
A set storing the global IDs of any boundary Verts.
Definition: DisContField.h:163
LocTraceToTraceMapSharedPtr m_locTraceToTraceMap
Definition: DisContField.h:198
AssemblyMapDGSharedPtr m_traceMap
Local to global DG mapping for trace space.
Definition: DisContField.h:158
GlobalLinSysMapShPtr m_globalBndMat
Global boundary matrix.
Definition: DisContField.h:152
PeriodicMap m_periodicVerts
A map which identifies groups of periodic vertices.
Definition: DisContField.h:168
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:185
std::vector< bool > m_leftAdjacentTraces
Definition: DisContField.h:192

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 657 of file DisContField.cpp.

661  : ExpList(In, DeclareCoeffPhysArrays)
662 {
663 
665 
666  // Set up boundary conditions for this variable.
667  // Do not set up BCs if default variable
668  if (variable.compare("DefaultVar") != 0)
669  {
670  SpatialDomains::BoundaryConditions bcs(m_session, graph);
671  GenerateBoundaryConditionExpansion(graph, bcs, variable);
672 
673  if (DeclareCoeffPhysArrays)
674  {
675  EvaluateBoundaryConditions(0.0, variable);
676  }
677 
679  {
680  // Find periodic edges for this variable.
681  FindPeriodicTraces(bcs, variable);
682 
683  if (SetUpJustDG)
684  {
685  SetUpDG(variable);
686  }
687  else
688  {
689  // set elmt edges to point to robin bc edges if required
690  int i, cnt = 0;
691  Array<OneD, int> ElmtID, TraceID;
692  GetBoundaryToElmtMap(ElmtID, TraceID);
693 
694  for (i = 0; i < m_bndCondExpansions.size(); ++i)
695  {
697 
698  int e;
699  locExpList = m_bndCondExpansions[i];
700 
701  for (e = 0; e < locExpList->GetExpSize(); ++e)
702  {
703  (*m_exp)[ElmtID[cnt + e]]->SetTraceExp(
704  TraceID[cnt + e], locExpList->GetExp(e));
705  locExpList->GetExp(e)->SetAdjacentElementExp(
706  TraceID[cnt + e], (*m_exp)[ElmtID[cnt + e]]);
707  }
708 
709  cnt += m_bndCondExpansions[i]->GetExpSize();
710  }
711 
712  if (m_session->DefinesSolverInfo("PROJECTION"))
713  {
714  std::string ProjectStr =
715  m_session->GetSolverInfo("PROJECTION");
716 
717  if ((ProjectStr == "MixedCGDG") ||
718  (ProjectStr == "Mixed_CG_Discontinuous"))
719  {
720  SetUpDG();
721  }
722  else
723  {
725  }
726  }
727  else
728  {
730  }
731  }
732  }
733  else
734  {
735  m_globalBndMat = In.m_globalBndMat;
736  m_trace = In.m_trace;
737  m_traceMap = In.m_traceMap;
738  m_locTraceToTraceMap = In.m_locTraceToTraceMap;
739  m_periodicVerts = In.m_periodicVerts;
740  m_periodicEdges = In.m_periodicEdges;
741  m_periodicFaces = In.m_periodicFaces;
742  m_periodicFwdCopy = In.m_periodicFwdCopy;
743  m_periodicBwdCopy = In.m_periodicBwdCopy;
744  m_boundaryTraces = In.m_boundaryTraces;
745  m_leftAdjacentTraces = In.m_leftAdjacentTraces;
746 
747  if (SetUpJustDG == false)
748  {
749  // set elmt edges to point to robin bc edges if required
750  int i, cnt = 0;
751  Array<OneD, int> ElmtID, TraceID;
752  GetBoundaryToElmtMap(ElmtID, TraceID);
753 
754  for (i = 0; i < m_bndCondExpansions.size(); ++i)
755  {
757 
758  int e;
759  locExpList = m_bndCondExpansions[i];
760 
761  for (e = 0; e < locExpList->GetExpSize(); ++e)
762  {
763  (*m_exp)[ElmtID[cnt + e]]->SetTraceExp(
764  TraceID[cnt + e], locExpList->GetExp(e));
765  locExpList->GetExp(e)->SetAdjacentElementExp(
766  TraceID[cnt + e], (*m_exp)[ElmtID[cnt + e]]);
767  }
768  cnt += m_bndCondExpansions[i]->GetExpSize();
769  }
770 
772  }
773  }
774  }
775 }
PeriodicMap m_periodicEdges
A map which identifies pairs of periodic edges.
Definition: DisContField.h:173
PeriodicMap m_periodicFaces
A map which identifies pairs of periodic faces.
Definition: DisContField.h:178
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 781 of file DisContField.cpp.

781  : ExpList(In)
782 {
783 }

◆ ~DisContField()

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

Destructor.

Definition at line 788 of file DisContField.cpp.

789 {
790 }

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 4116 of file DisContField.cpp.

4117 {
4118  int i, cnt, e, ncoeff_trace;
4119  Array<OneD, NekDouble> force, out_tmp, qrhs, qrhs1;
4120  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
4121  m_traceMap->GetElmtToTrace();
4122 
4123  StdRegions::Orientation edgedir;
4124 
4125  int nq_elmt, nm_elmt;
4126  int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs();
4127  Array<OneD, NekDouble> loc_lambda(LocBndCoeffs), trace_lambda;
4128  Array<OneD, NekDouble> tmp_coeffs;
4129  m_traceMap->GlobalToLocalBnd(m_trace->GetCoeffs(), loc_lambda);
4130 
4131  trace_lambda = loc_lambda;
4132 
4133  int dim = (m_expType == e2D) ? 2 : 3;
4134 
4135  int num_points[] = {0, 0, 0};
4136  int num_modes[] = {0, 0, 0};
4137 
4138  // Calculate Q using standard DG formulation.
4139  for (i = cnt = 0; i < GetExpSize(); ++i)
4140  {
4141  nq_elmt = (*m_exp)[i]->GetTotPoints();
4142  nm_elmt = (*m_exp)[i]->GetNcoeffs();
4143  qrhs = Array<OneD, NekDouble>(nq_elmt);
4144  qrhs1 = Array<OneD, NekDouble>(nq_elmt);
4145  force = Array<OneD, NekDouble>(2 * nm_elmt);
4146  out_tmp = force + nm_elmt;
4148 
4149  for (int j = 0; j < dim; ++j)
4150  {
4151  num_points[j] = (*m_exp)[i]->GetBasis(j)->GetNumPoints();
4152  num_modes[j] = (*m_exp)[i]->GetBasis(j)->GetNumModes();
4153  }
4154 
4155  // Probably a better way of setting up lambda than this. Note
4156  // cannot use PutCoeffsInToElmts since lambda space is mapped
4157  // during the solve.
4158  int nTraces = (*m_exp)[i]->GetNtraces();
4159  Array<OneD, Array<OneD, NekDouble>> traceCoeffs(nTraces);
4160 
4161  for (e = 0; e < (*m_exp)[i]->GetNtraces(); ++e)
4162  {
4163  edgedir = (*m_exp)[i]->GetTraceOrient(e);
4164  ncoeff_trace = elmtToTrace[i][e]->GetNcoeffs();
4165  traceCoeffs[e] = Array<OneD, NekDouble>(ncoeff_trace);
4166  Vmath::Vcopy(ncoeff_trace, trace_lambda, 1, traceCoeffs[e], 1);
4167  if (dim == 2)
4168  {
4169  elmtToTrace[i][e]->SetCoeffsToOrientation(
4170  edgedir, traceCoeffs[e], traceCoeffs[e]);
4171  }
4172  else
4173  {
4174  (*m_exp)[i]
4175  ->as<LocalRegions::Expansion3D>()
4176  ->SetFaceToGeomOrientation(e, traceCoeffs[e]);
4177  }
4178  trace_lambda = trace_lambda + ncoeff_trace;
4179  }
4180 
4181  // creating orthogonal expansion (checking if we have quads or
4182  // triangles)
4183  LibUtilities::ShapeType shape = (*m_exp)[i]->DetShapeType();
4184  switch (shape)
4185  {
4187  {
4188  const LibUtilities::PointsKey PkeyQ1(
4189  num_points[0], LibUtilities::eGaussLobattoLegendre);
4190  const LibUtilities::PointsKey PkeyQ2(
4191  num_points[1], LibUtilities::eGaussLobattoLegendre);
4192  LibUtilities::BasisKey BkeyQ1(LibUtilities::eOrtho_A,
4193  num_modes[0], PkeyQ1);
4194  LibUtilities::BasisKey BkeyQ2(LibUtilities::eOrtho_A,
4195  num_modes[1], PkeyQ2);
4197  std::dynamic_pointer_cast<SpatialDomains::QuadGeom>(
4198  (*m_exp)[i]->GetGeom());
4200  BkeyQ1, BkeyQ2, qGeom);
4201  }
4202  break;
4204  {
4205  const LibUtilities::PointsKey PkeyT1(
4206  num_points[0], LibUtilities::eGaussLobattoLegendre);
4207  const LibUtilities::PointsKey PkeyT2(
4208  num_points[1], LibUtilities::eGaussRadauMAlpha1Beta0);
4209  LibUtilities::BasisKey BkeyT1(LibUtilities::eOrtho_A,
4210  num_modes[0], PkeyT1);
4211  LibUtilities::BasisKey BkeyT2(LibUtilities::eOrtho_B,
4212  num_modes[1], PkeyT2);
4214  std::dynamic_pointer_cast<SpatialDomains::TriGeom>(
4215  (*m_exp)[i]->GetGeom());
4217  BkeyT1, BkeyT2, tGeom);
4218  }
4219  break;
4221  {
4222  const LibUtilities::PointsKey PkeyH1(
4223  num_points[0], LibUtilities::eGaussLobattoLegendre);
4224  const LibUtilities::PointsKey PkeyH2(
4225  num_points[1], LibUtilities::eGaussLobattoLegendre);
4226  const LibUtilities::PointsKey PkeyH3(
4227  num_points[2], LibUtilities::eGaussLobattoLegendre);
4228  LibUtilities::BasisKey BkeyH1(LibUtilities::eOrtho_A,
4229  num_modes[0], PkeyH1);
4230  LibUtilities::BasisKey BkeyH2(LibUtilities::eOrtho_A,
4231  num_modes[1], PkeyH2);
4232  LibUtilities::BasisKey BkeyH3(LibUtilities::eOrtho_A,
4233  num_modes[2], PkeyH3);
4235  std::dynamic_pointer_cast<SpatialDomains::HexGeom>(
4236  (*m_exp)[i]->GetGeom());
4238  BkeyH1, BkeyH2, BkeyH3, hGeom);
4239  }
4240  break;
4242  {
4243  const LibUtilities::PointsKey PkeyT1(
4244  num_points[0], LibUtilities::eGaussLobattoLegendre);
4245  const LibUtilities::PointsKey PkeyT2(
4246  num_points[1], LibUtilities::eGaussRadauMAlpha1Beta0);
4247  const LibUtilities::PointsKey PkeyT3(
4248  num_points[2], LibUtilities::eGaussRadauMAlpha2Beta0);
4249  LibUtilities::BasisKey BkeyT1(LibUtilities::eOrtho_A,
4250  num_modes[0], PkeyT1);
4251  LibUtilities::BasisKey BkeyT2(LibUtilities::eOrtho_B,
4252  num_modes[1], PkeyT2);
4253  LibUtilities::BasisKey BkeyT3(LibUtilities::eOrtho_C,
4254  num_modes[2], PkeyT3);
4256  std::dynamic_pointer_cast<SpatialDomains::TetGeom>(
4257  (*m_exp)[i]->GetGeom());
4259  BkeyT1, BkeyT2, BkeyT3, tGeom);
4260  }
4261  break;
4262  default:
4264  "Wrong shape type, HDG postprocessing is not "
4265  "implemented");
4266  };
4267 
4268  // DGDeriv
4269  // (d/dx w, d/dx q_0)
4270  (*m_exp)[i]->DGDeriv(0, tmp_coeffs = m_coeffs + m_coeff_offset[i],
4271  elmtToTrace[i], traceCoeffs, out_tmp);
4272  (*m_exp)[i]->BwdTrans(out_tmp, qrhs);
4273  ppExp->IProductWRTDerivBase(0, qrhs, force);
4274 
4275  // + (d/dy w, d/dy q_1)
4276  (*m_exp)[i]->DGDeriv(1, tmp_coeffs = m_coeffs + m_coeff_offset[i],
4277  elmtToTrace[i], traceCoeffs, out_tmp);
4278 
4279  (*m_exp)[i]->BwdTrans(out_tmp, qrhs);
4280  ppExp->IProductWRTDerivBase(1, qrhs, out_tmp);
4281 
4282  Vmath::Vadd(nm_elmt, force, 1, out_tmp, 1, force, 1);
4283 
4284  // determine force[0] = (1,u)
4285  (*m_exp)[i]->BwdTrans(tmp_coeffs = m_coeffs + m_coeff_offset[i], qrhs);
4286  force[0] = (*m_exp)[i]->Integral(qrhs);
4287 
4288  // multiply by inverse Laplacian matrix
4289  // get matrix inverse
4290  LocalRegions::MatrixKey lapkey(StdRegions::eInvLaplacianWithUnityMean,
4291  ppExp->DetShapeType(), *ppExp);
4292  DNekScalMatSharedPtr lapsys = ppExp->GetLocMatrix(lapkey);
4293 
4294  NekVector<NekDouble> in(nm_elmt, force, eWrapper);
4295  NekVector<NekDouble> out(nm_elmt);
4296 
4297  out = (*lapsys) * in;
4298 
4299  // Transforming back to modified basis
4300  Array<OneD, NekDouble> work(nq_elmt);
4301  ppExp->BwdTrans(out.GetPtr(), work);
4302  (*m_exp)[i]->FwdTrans(work, tmp_coeffs = outarray + m_coeff_offset[i]);
4303  }
4304 }
#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:1158
int GetExpSize(void)
This function returns the number of elements in the expansion.
Definition: ExpList.h:2204
Array< OneD, int > m_coeff_offset
Offset of elemental data into the array m_coeffs.
Definition: ExpList.h:1210
ExpansionType m_expType
Exapnsion type.
Definition: ExpList.h:1119
@ eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:53
@ eOrtho_A
Principle Orthogonal Functions .
Definition: BasisType.h:44
@ eOrtho_C
Principle Orthogonal Functions .
Definition: BasisType.h:48
@ 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:84
std::shared_ptr< HexGeom > HexGeomSharedPtr
Definition: HexGeom.h:86
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:359
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1255

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::MultiRegions::e2D, Nektar::ErrorUtil::efatal, Nektar::LibUtilities::eGaussLobattoLegendre, 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 862 of file DisContField.cpp.

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

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::ExpList::m_comm, Nektar::MultiRegions::RotPeriodicInfo::m_dir, Nektar::MultiRegions::ExpList::m_expType, Nektar::MultiRegions::ExpList::m_graph, m_periodicEdges, m_periodicFaces, m_periodicVerts, 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 809 of file DisContField.cpp.

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

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 514 of file DisContField.cpp.

518 {
519  boost::ignore_unused(variable);
521 
522  returnval =
524 
525  map<int, int> GeometryToRegionsMap;
526 
528  Allbcs.GetBoundaryRegions();
530  Allbcs.GetBoundaryConditions();
531 
532  // Set up a map of all boundary regions
533  for (auto &it : bregions)
534  {
535  for (auto &bregionIt : *it.second)
536  {
537  // can assume that all regions only contain one point in 1D
538  // Really do not need loop above
539  int id = bregionIt.second->m_geomVec[0]->GetGlobalID();
540  GeometryToRegionsMap[id] = it.first;
541  }
542  }
543 
544  map<int, SpatialDomains::GeometrySharedPtr> EndOfDomain;
545 
546  // Now find out which points in domain have only one vertex
547  for (auto &domIt : domain)
548  {
549  SpatialDomains::CompositeSharedPtr geomvector = domIt.second;
550  for (int i = 0; i < geomvector->m_geomVec.size(); ++i)
551  {
552  for (int j = 0; j < 2; ++j)
553  {
554  int vid = geomvector->m_geomVec[i]->GetVid(j);
555  if (EndOfDomain.count(vid) == 0)
556  {
557  EndOfDomain[vid] = geomvector->m_geomVec[i]->GetVertex(j);
558  }
559  else
560  {
561  EndOfDomain.erase(vid);
562  }
563  }
564  }
565  }
566  ASSERTL1(EndOfDomain.size() == 2, "Did not find two ends of domain");
567 
568  for (auto &regIt : EndOfDomain)
569  {
570  if (GeometryToRegionsMap.count(regIt.first) != 0)
571  {
572  // Set up boundary condition up
573  auto iter = GeometryToRegionsMap.find(regIt.first);
574  ASSERTL1(iter != GeometryToRegionsMap.end(),
575  "Failied to find GeometryToRegionMap");
576 
577  int regionId = iter->second;
578  auto bregionsIter = bregions.find(regionId);
579  ASSERTL1(bregionsIter != bregions.end(),
580  "Failed to find boundary region");
581 
582  SpatialDomains::BoundaryRegionShPtr breg = bregionsIter->second;
583  returnval->AddBoundaryRegions(regionId, breg);
584 
585  auto bconditionsIter = bconditions.find(regionId);
586  ASSERTL1(bconditionsIter != bconditions.end(),
587  "Failed to find boundary collection");
588 
590  bconditionsIter->second;
591  returnval->AddBoundaryConditions(regionId, bcond);
592  }
593  }
594 
595  return returnval;
596 }
std::shared_ptr< BoundaryRegion > BoundaryRegionShPtr
Definition: Conditions.h:208
std::shared_ptr< BoundaryConditionMap > BoundaryConditionMapShPtr
Definition: Conditions.h:218

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

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 2704 of file DisContField.cpp.

2706 {
2707  ASSERTL0(mkey.GetMatrixType() == StdRegions::eHybridDGHelmBndLam,
2708  "Routine currently only tested for HybridDGHelmholtz");
2709 
2710  ASSERTL1(mkey.GetGlobalSysSolnType() != eDirectFullMatrix,
2711  "Full matrix global systems are not supported for HDG "
2712  "expansions");
2713 
2714  ASSERTL1(mkey.GetGlobalSysSolnType() == m_traceMap->GetGlobalSysSolnType(),
2715  "The local to global map is not set up for the requested "
2716  "solution type");
2717 
2718  GlobalLinSysSharedPtr glo_matrix;
2719  auto matrixIter = m_globalBndMat->find(mkey);
2720 
2721  if (matrixIter == m_globalBndMat->end())
2722  {
2723  glo_matrix = GenGlobalBndLinSys(mkey, m_traceMap);
2724  (*m_globalBndMat)[mkey] = glo_matrix;
2725  }
2726  else
2727  {
2728  glo_matrix = matrixIter->second;
2729  }
2730 
2731  return glo_matrix;
2732 }
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:2431
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().

◆ GetLocTraceToTraceMap()

void Nektar::MultiRegions::DisContField::GetLocTraceToTraceMap ( LocTraceToTraceMapSharedPtr LocTraceToTraceMap)
inline

Definition at line 121 of file DisContField.h.

123  {
124  LocTraceToTraceMap = m_locTraceToTraceMap;
125  }

References m_locTraceToTraceMap.

◆ GetNegatedFluxNormal()

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

Definition at line 2768 of file DisContField.cpp.

2769 {
2770  if (m_negatedFluxNormal.size() == 0)
2771  {
2772  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>>
2773  &elmtToTrace = m_traceMap->GetElmtToTrace();
2774 
2775  m_negatedFluxNormal.resize(2 * GetExpSize());
2776 
2777  for (int i = 0; i < GetExpSize(); ++i)
2778  {
2779 
2780  for (int v = 0; v < 2; ++v)
2781  {
2783  elmtToTrace[i][v]->as<LocalRegions::Expansion0D>();
2784 
2785  if (vertExp->GetLeftAdjacentElementExp()
2786  ->GetGeom()
2787  ->GetGlobalID() !=
2788  (*m_exp)[i]->GetGeom()->GetGlobalID())
2789  {
2790  m_negatedFluxNormal[2 * i + v] = true;
2791  }
2792  else
2793  {
2794  m_negatedFluxNormal[2 * i + v] = false;
2795  }
2796  }
2797  }
2798  }
2799 
2800  return m_negatedFluxNormal;
2801 }
std::vector< bool > m_negatedFluxNormal
Definition: DisContField.h:360
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

This routine determines if an element is to the "left" of the adjacent trace, which arises from the idea there is a local normal direction between two elements (i.e. on the trace) and one elements would then be the left.

This is typically required since we only define one normal direction along the trace which goes from the left adjacent element to the right adjacent element. It is also useful in DG discretisations to set up a local Riemann problem where we need to have the concept of a local normal direction which is unique between two elements.

There are two cases to be checked:

1) First is the trace on a boundary condition (perioidic or otherwise) or on a partition boundary. If a partition boundary then trace is always considered to be the left adjacent trace and the normal is pointing outward of the soltuion domain. We have to perform an additional case for a periodic boundary where wer chose the element with the lowest global id. If the trace is on a parallel partition we use a member of the traceMap where one side is chosen to contribute to a unique map and have a value which is not -1 so this is the identifier for the left adjacent side

2) If the element is a elemental boundary on one element the left adjacent element is defined by a link to the left element from the trace expansion and this is consistent with the defitiion of the normal which is determined by the largest id (in contrast to the periodic boundary definition !). This reversal of convention does not really matter providing the normals are defined consistently.

Definition at line 473 of file DisContField.cpp.

474 {
476  m_traceMap->GetElmtToTrace()[n][e];
477 
478  PeriodicMap periodicTraces = (m_expType == e1D) ? m_periodicVerts
479  : (m_expType == e2D) ? m_periodicEdges
480  : m_periodicFaces;
481 
482  bool fwd = true;
483  if (traceEl->GetLeftAdjacentElementTrace() == -1 ||
484  traceEl->GetRightAdjacentElementTrace() == -1)
485  {
486  // Boundary edge (1 connected element). Do nothing in
487  // serial.
488  auto it = m_boundaryTraces.find(traceEl->GetElmtId());
489 
490  // If the edge does not have a boundary condition set on
491  // it, then assume it is a partition edge or periodic.
492  if (it == m_boundaryTraces.end())
493  {
494  fwd = true;
495  }
496  }
497  else if (traceEl->GetLeftAdjacentElementTrace() != -1 &&
498  traceEl->GetRightAdjacentElementTrace() != -1)
499  {
500  // Non-boundary edge (2 connected elements).
501  fwd = (traceEl->GetLeftAdjacentElementExp().get()) == (*m_exp)[n].get();
502  }
503  else
504  {
505  ASSERTL2(false, "Unconnected trace element!");
506  }
507 
508  return fwd;
509 }
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed to...
Definition: ErrorUtil.hpp:272

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 4046 of file DisContField.cpp.

4048 {
4049  int i, e, ncoeff_edge;
4050  Array<OneD, const NekDouble> tmp_coeffs;
4051  Array<OneD, NekDouble> out_d(m_ncoeffs), out_tmp;
4052 
4053  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
4054  m_traceMap->GetElmtToTrace();
4055 
4056  StdRegions::Orientation edgedir;
4057 
4058  int cnt;
4059  int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs();
4060  Array<OneD, NekDouble> loc_lambda(LocBndCoeffs), edge_lambda;
4061 
4062  m_traceMap->GlobalToLocalBnd(m_trace->GetCoeffs(), loc_lambda);
4063 
4064  edge_lambda = loc_lambda;
4065 
4066  // Calculate Q using standard DG formulation.
4067  for (i = cnt = 0; i < GetExpSize(); ++i)
4068  {
4069  // Probably a better way of setting up lambda than this.
4070  // Note cannot use PutCoeffsInToElmts since lambda space
4071  // is mapped during the solve.
4072  int nEdges = (*m_exp)[i]->GetGeom()->GetNumEdges();
4073  Array<OneD, Array<OneD, NekDouble>> edgeCoeffs(nEdges);
4074 
4075  for (e = 0; e < nEdges; ++e)
4076  {
4077  edgedir = (*m_exp)[i]->GetTraceOrient(e);
4078  ncoeff_edge = elmtToTrace[i][e]->GetNcoeffs();
4079  edgeCoeffs[e] = Array<OneD, NekDouble>(ncoeff_edge);
4080  Vmath::Vcopy(ncoeff_edge, edge_lambda, 1, edgeCoeffs[e], 1);
4081  elmtToTrace[i][e]->SetCoeffsToOrientation(edgedir, edgeCoeffs[e],
4082  edgeCoeffs[e]);
4083  edge_lambda = edge_lambda + ncoeff_edge;
4084  }
4085 
4086  (*m_exp)[i]->DGDeriv(dir, tmp_coeffs = m_coeffs + m_coeff_offset[i],
4087  elmtToTrace[i], edgeCoeffs, out_tmp = out_d + cnt);
4088  cnt += (*m_exp)[i]->GetNcoeffs();
4089  }
4090 
4091  BwdTrans(out_d, m_phys);
4092  Vmath::Vsub(m_npoints, m_phys, 1, soln, 1, m_phys, 1);
4093  return L2(m_phys);
4094 }
int m_ncoeffs
The total number of local degrees of freedom. m_ncoeffs .
Definition: ExpList.h:1136
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:551
void BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function elementally evaluates the backward transformation of the global spectral/hp element exp...
Definition: ExpList.h:1849
Array< OneD, NekDouble > m_phys
The global expansion evaluated at the quadrature points.
Definition: ExpList.h:1175
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:419

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 2740 of file DisContField.cpp.

2741 {
2742  int i;
2743  bool returnval = true;
2744 
2745  for (i = 0; i < m_bndConditions.size(); ++i)
2746  {
2747  // check to see if boundary condition type is the same
2748  // and there are the same number of boundary
2749  // conditions in the boundary definition.
2750  if ((m_bndConditions[i]->GetBoundaryConditionType() !=
2751  In.m_bndConditions[i]->GetBoundaryConditionType()) ||
2752  (m_bndCondExpansions[i]->GetExpSize() !=
2753  In.m_bndCondExpansions[i]->GetExpSize()))
2754  {
2755  returnval = false;
2756  break;
2757  }
2758  }
2759 
2760  // Compare with all other processes. Return true only if all
2761  // processes report having the same boundary conditions.
2762  int vSame = (returnval ? 1 : 0);
2763  m_comm->GetRowComm()->AllReduce(vSame, LibUtilities::ReduceMin);
2764 
2765  return (vSame == 1);
2766 }

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 163 of file DisContField.cpp.

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

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_comm, 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 3305 of file DisContField.cpp.

3308 {
3309 
3310  ASSERTL0(m_expType != e1D, "This method is not setup or "
3311  "tested for 1D expansion");
3312 
3313  Array<OneD, NekDouble> Coeffs(m_trace->GetNcoeffs());
3314 
3315  m_trace->IProductWRTBase(Fwd, Coeffs);
3316  m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(0, Coeffs, outarray);
3317  m_trace->IProductWRTBase(Bwd, Coeffs);
3318  m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(1, Coeffs, outarray);
3319 }

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 3115 of file DisContField.cpp.

3117 {
3118  if (m_expType == e1D)
3119  {
3120  int n, offset, t_offset;
3121 
3122  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>>
3123  &elmtToTrace = m_traceMap->GetElmtToTrace();
3124 
3125  vector<bool> negatedFluxNormal = GetNegatedFluxNormal();
3126 
3127  for (n = 0; n < GetExpSize(); ++n)
3128  {
3129  // Number of coefficients on each element
3130  int e_ncoeffs = (*m_exp)[n]->GetNcoeffs();
3131 
3132  offset = GetCoeff_Offset(n);
3133 
3134  // Implementation for every points except Gauss points
3135  if ((*m_exp)[n]->GetBasis(0)->GetBasisType() !=
3137  {
3138  t_offset =
3139  GetTrace()->GetCoeff_Offset(elmtToTrace[n][0]->GetElmtId());
3140  if (negatedFluxNormal[2 * n])
3141  {
3142  outarray[offset] -= Fn[t_offset];
3143  }
3144  else
3145  {
3146  outarray[offset] += Fn[t_offset];
3147  }
3148 
3149  t_offset =
3150  GetTrace()->GetCoeff_Offset(elmtToTrace[n][1]->GetElmtId());
3151 
3152  if (negatedFluxNormal[2 * n + 1])
3153  {
3154  outarray[offset + (*m_exp)[n]->GetVertexMap(1)] -=
3155  Fn[t_offset];
3156  }
3157  else
3158  {
3159  outarray[offset + (*m_exp)[n]->GetVertexMap(1)] +=
3160  Fn[t_offset];
3161  }
3162  }
3163  else
3164  {
3165  int j;
3166  static DNekMatSharedPtr m_Ixm, m_Ixp;
3167  static int sav_ncoeffs = 0;
3168  if (!m_Ixm.get() || e_ncoeffs != sav_ncoeffs)
3169  {
3171  const LibUtilities::PointsKey BS_p(
3173  const LibUtilities::BasisKey BS_k(
3174  LibUtilities::eGauss_Lagrange, e_ncoeffs, BS_p);
3175 
3176  BASE = LibUtilities::BasisManager()[BS_k];
3177 
3178  Array<OneD, NekDouble> coords(1, 0.0);
3179 
3180  coords[0] = -1.0;
3181  m_Ixm = BASE->GetI(coords);
3182 
3183  coords[0] = 1.0;
3184  m_Ixp = BASE->GetI(coords);
3185 
3186  sav_ncoeffs = e_ncoeffs;
3187  }
3188 
3189  t_offset =
3190  GetTrace()->GetCoeff_Offset(elmtToTrace[n][0]->GetElmtId());
3191 
3192  if (negatedFluxNormal[2 * n])
3193  {
3194  for (j = 0; j < e_ncoeffs; j++)
3195  {
3196  outarray[offset + j] -=
3197  (m_Ixm->GetPtr())[j] * Fn[t_offset];
3198  }
3199  }
3200  else
3201  {
3202  for (j = 0; j < e_ncoeffs; j++)
3203  {
3204  outarray[offset + j] +=
3205  (m_Ixm->GetPtr())[j] * Fn[t_offset];
3206  }
3207  }
3208 
3209  t_offset =
3210  GetTrace()->GetCoeff_Offset(elmtToTrace[n][1]->GetElmtId());
3211 
3212  if (negatedFluxNormal[2 * n + 1])
3213  {
3214  for (j = 0; j < e_ncoeffs; j++)
3215  {
3216  outarray[offset + j] -=
3217  (m_Ixp->GetPtr())[j] * Fn[t_offset];
3218  }
3219  }
3220  else
3221  {
3222  for (j = 0; j < e_ncoeffs; j++)
3223  {
3224  outarray[offset + j] +=
3225  (m_Ixp->GetPtr())[j] * Fn[t_offset];
3226  }
3227  }
3228  }
3229  }
3230  }
3231  else // other dimensions
3232  {
3233  // Basis definition on each element
3234  LibUtilities::BasisSharedPtr basis = (*m_exp)[0]->GetBasis(0);
3235  if ((m_expType != e1D) &&
3236  (basis->GetBasisType() != LibUtilities::eGauss_Lagrange))
3237  {
3238  Array<OneD, NekDouble> Fcoeffs(m_trace->GetNcoeffs());
3239  m_trace->IProductWRTBase(Fn, Fcoeffs);
3240 
3241  m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(Fcoeffs,
3242  outarray);
3243  }
3244  else
3245  {
3246  int e, n, offset, t_offset;
3247  Array<OneD, NekDouble> e_outarray;
3248  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>>
3249  &elmtToTrace = m_traceMap->GetElmtToTrace();
3250 
3251  if (m_expType == e2D)
3252  {
3253  for (n = 0; n < GetExpSize(); ++n)
3254  {
3255  offset = GetCoeff_Offset(n);
3256  for (e = 0; e < (*m_exp)[n]->GetNtraces(); ++e)
3257  {
3258  t_offset = GetTrace()->GetPhys_Offset(
3259  elmtToTrace[n][e]->GetElmtId());
3260  (*m_exp)[n]->AddEdgeNormBoundaryInt(
3261  e, elmtToTrace[n][e], Fn + t_offset,
3262  e_outarray = outarray + offset);
3263  }
3264  }
3265  }
3266  else if (m_expType == e3D)
3267  {
3268  for (n = 0; n < GetExpSize(); ++n)
3269  {
3270  offset = GetCoeff_Offset(n);
3271  for (e = 0; e < (*m_exp)[n]->GetNtraces(); ++e)
3272  {
3273  t_offset = GetTrace()->GetPhys_Offset(
3274  elmtToTrace[n][e]->GetElmtId());
3275  (*m_exp)[n]->AddFaceNormBoundaryInt(
3276  e, elmtToTrace[n][e], Fn + t_offset,
3277  e_outarray = outarray + offset);
3278  }
3279  }
3280  }
3281  }
3282  }
3283 }
std::vector< bool > & GetNegatedFluxNormal(void)
std::shared_ptr< ExpList > & GetTrace()
Definition: ExpList.h:2319
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:2232
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:59
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:75

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 4324 of file DisContField.cpp.

4328 {
4329  Array<OneD, NekDouble> FCoeffs(m_trace->GetNcoeffs());
4330 
4331  m_trace->IProductWRTBase(FwdFlux, FCoeffs);
4332  m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(1, FCoeffs, outarray);
4333  m_trace->IProductWRTBase(BwdFlux, FCoeffs);
4334  m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(0, FCoeffs, outarray);
4335 }

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 3005 of file DisContField.cpp.

3008 {
3009  // Basis definition on each element
3010  LibUtilities::BasisSharedPtr basis = (*m_exp)[0]->GetBasis(0);
3011  if (basis->GetBasisType() != LibUtilities::eGauss_Lagrange)
3012  {
3013  Array<OneD, NekDouble> tracevals(
3014  m_locTraceToTraceMap->GetNLocTracePts(), 0.0);
3015 
3016  Array<OneD, NekDouble> invals =
3017  tracevals + m_locTraceToTraceMap->GetNFwdLocTracePts();
3018 
3019  m_locTraceToTraceMap->InterpTraceToLocTrace(1, Bwd, invals);
3020 
3021  m_locTraceToTraceMap->InterpTraceToLocTrace(0, Fwd, tracevals);
3022 
3023  m_locTraceToTraceMap->AddLocTracesToField(tracevals, field);
3024  }
3025  else
3026  {
3028  "v_AddTraceQuadPhysToField not coded for eGauss_Lagrange");
3029  }
3030 }

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 3492 of file DisContField.cpp.

3496 {
3497  int i;
3498  int npoints;
3499 
3500  MultiRegions::ExpListSharedPtr locExpList;
3501 
3502  for (i = 0; i < m_bndCondExpansions.size(); ++i)
3503  {
3504  if (time == 0.0 || m_bndConditions[i]->IsTimeDependent())
3505  {
3506  m_bndCondBndWeight[i] = 1.0;
3507  locExpList = m_bndCondExpansions[i];
3508 
3509  npoints = locExpList->GetNpoints();
3510  Array<OneD, NekDouble> x0(npoints, 0.0);
3511  Array<OneD, NekDouble> x1(npoints, 0.0);
3512  Array<OneD, NekDouble> x2(npoints, 0.0);
3513 
3514  locExpList->GetCoords(x0, x1, x2);
3515 
3516  if (x2_in != NekConstants::kNekUnsetDouble &&
3518  {
3519  Vmath::Fill(npoints, x2_in, x1, 1);
3520  Vmath::Fill(npoints, x3_in, x2, 1);
3521  }
3522  else if (x2_in != NekConstants::kNekUnsetDouble)
3523  {
3524  Vmath::Fill(npoints, x2_in, x2, 1);
3525  }
3526 
3527  // treat 1D expansions separately since we only
3528  // require an evaluation at a point rather than
3529  // any projections or inner products that are not
3530  // available in a PointExp
3531  if (m_expType == e1D)
3532  {
3533  if (m_bndConditions[i]->GetBoundaryConditionType() ==
3535  {
3536  m_bndCondExpansions[i]->SetCoeff(
3537  0, (std::static_pointer_cast<
3538  SpatialDomains ::DirichletBoundaryCondition>(
3539  m_bndConditions[i])
3540  ->m_dirichletCondition)
3541  .Evaluate(x0[0], x1[0], x2[0], time));
3542  m_bndCondExpansions[i]->SetPhys(
3543  0, m_bndCondExpansions[i]->GetCoeff(0));
3544  }
3545  else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3547  {
3548  m_bndCondExpansions[i]->SetCoeff(
3549  0, (std::static_pointer_cast<
3550  SpatialDomains ::NeumannBoundaryCondition>(
3551  m_bndConditions[i])
3552  ->m_neumannCondition)
3553  .Evaluate(x0[0], x1[0], x2[0], time));
3554  }
3555  else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3557  {
3558  m_bndCondExpansions[i]->SetCoeff(
3559  0, (std::static_pointer_cast<
3560  SpatialDomains ::RobinBoundaryCondition>(
3561  m_bndConditions[i])
3562  ->m_robinFunction)
3563  .Evaluate(x0[0], x1[0], x2[0], time));
3564  }
3565  else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3567  {
3568  continue;
3569  }
3570  else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3572  {
3573  }
3574  else
3575  {
3577  "This type of BC not implemented yet");
3578  }
3579  }
3580  else // 2D and 3D versions
3581  {
3582  if (m_bndConditions[i]->GetBoundaryConditionType() ==
3584  {
3586  std::static_pointer_cast<
3587  SpatialDomains::DirichletBoundaryCondition>(
3588  m_bndConditions[i]);
3589 
3590  Array<OneD, NekDouble> valuesFile(npoints, 1.0),
3591  valuesExp(npoints, 1.0);
3592 
3593  string filebcs = bcPtr->m_filename;
3594  string exprbcs = bcPtr->m_expr;
3595 
3596  if (filebcs != "")
3597  {
3598  ExtractFileBCs(filebcs, bcPtr->GetComm(), varName,
3599  locExpList);
3600  valuesFile = locExpList->GetPhys();
3601  }
3602 
3603  if (exprbcs != "")
3604  {
3605  LibUtilities::Equation condition =
3606  std::static_pointer_cast<
3607  SpatialDomains::DirichletBoundaryCondition>(
3608  m_bndConditions[i])
3609  ->m_dirichletCondition;
3610 
3611  condition.Evaluate(x0, x1, x2, time, valuesExp);
3612  }
3613 
3614  Vmath::Vmul(npoints, valuesExp, 1, valuesFile, 1,
3615  locExpList->UpdatePhys(), 1);
3616 
3617  locExpList->FwdTransBndConstrained(
3618  locExpList->GetPhys(), locExpList->UpdateCoeffs());
3619  }
3620  else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3622  {
3624  std::static_pointer_cast<
3625  SpatialDomains::NeumannBoundaryCondition>(
3626  m_bndConditions[i]);
3627  string filebcs = bcPtr->m_filename;
3628  if (filebcs != "")
3629  {
3630  ExtractFileBCs(filebcs, bcPtr->GetComm(), varName,
3631  locExpList);
3632  }
3633  else
3634  {
3635  LibUtilities::Equation condition =
3636  std::static_pointer_cast<
3637  SpatialDomains::NeumannBoundaryCondition>(
3638  m_bndConditions[i])
3639  ->m_neumannCondition;
3640  condition.Evaluate(x0, x1, x2, time,
3641  locExpList->UpdatePhys());
3642  }
3643 
3644  locExpList->IProductWRTBase(locExpList->GetPhys(),
3645  locExpList->UpdateCoeffs());
3646  }
3647  else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3649  {
3651  std::static_pointer_cast<
3652  SpatialDomains::RobinBoundaryCondition>(
3653  m_bndConditions[i]);
3654 
3655  string filebcs = bcPtr->m_filename;
3656 
3657  if (filebcs != "")
3658  {
3659  ExtractFileBCs(filebcs, bcPtr->GetComm(), varName,
3660  locExpList);
3661  }
3662  else
3663  {
3664  LibUtilities::Equation condition =
3665  std::static_pointer_cast<
3666  SpatialDomains::RobinBoundaryCondition>(
3667  m_bndConditions[i])
3668  ->m_robinFunction;
3669  condition.Evaluate(x0, x1, x2, time,
3670  locExpList->UpdatePhys());
3671  }
3672 
3673  locExpList->IProductWRTBase(locExpList->GetPhys(),
3674  locExpList->UpdateCoeffs());
3675  }
3676  else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3678  {
3679  continue;
3680  }
3681  else
3682  {
3684  "This type of BC not implemented yet");
3685  }
3686  }
3687  }
3688  }
3689 }
void ExtractFileBCs(const std::string &fileName, LibUtilities::CommSharedPtr comm, const std::string &varName, const std::shared_ptr< ExpList > locExpList)
Definition: ExpList.cpp:3304
NekDouble GetCoeff(int i)
Get the i th value (coefficient) of m_coeffs.
Definition: ExpList.h:2173
static const NekDouble kNekUnsetDouble
std::shared_ptr< DirichletBoundaryCondition > DirichletBCShPtr
Definition: Conditions.h:213
std::shared_ptr< NeumannBoundaryCondition > NeumannBCShPtr
Definition: Conditions.h:214
std::shared_ptr< RobinBoundaryCondition > RobinBCShPtr
Definition: Conditions.h:215
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:209
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 3032 of file DisContField.cpp.

3033 {
3034  ASSERTL1(m_physState == true, "local physical space is not true ");
3035  v_ExtractTracePhys(m_phys, outarray);
3036 }
virtual void v_ExtractTracePhys(Array< OneD, NekDouble > &outarray)
bool m_physState
The state of the array m_phys.
Definition: ExpList.h:1184

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 3052 of file DisContField.cpp.

3055 {
3056  LibUtilities::BasisSharedPtr basis = (*m_exp)[0]->GetBasis(0);
3057  if ((basis->GetBasisType() != LibUtilities::eGauss_Lagrange))
3058  {
3059  Vmath::Zero(outarray.size(), outarray, 1);
3060 
3061  Array<OneD, NekDouble> tracevals(
3062  m_locTraceToTraceMap->GetNFwdLocTracePts());
3063  m_locTraceToTraceMap->FwdLocTracesFromField(inarray, tracevals);
3064  m_locTraceToTraceMap->InterpLocTracesToTrace(0, tracevals, outarray);
3065  m_traceMap->GetAssemblyCommDG()->PerformExchange(outarray, outarray);
3066  }
3067  else
3068  {
3069 
3070  // Loop over elemente and collect forward expansion
3071  int nexp = GetExpSize();
3072  int n, p, offset, phys_offset;
3073  Array<OneD, NekDouble> t_tmp;
3074 
3075  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>>
3076  &elmtToTrace = m_traceMap->GetElmtToTrace();
3077 
3078  ASSERTL1(outarray.size() >= m_trace->GetNpoints(),
3079  "input array is of insufficient length");
3080 
3081  for (n = 0; n < nexp; ++n)
3082  {
3083  phys_offset = GetPhys_Offset(n);
3084 
3085  for (p = 0; p < (*m_exp)[n]->GetNtraces(); ++p)
3086  {
3087  offset =
3088  m_trace->GetPhys_Offset(elmtToTrace[n][p]->GetElmtId());
3089  (*m_exp)[n]->GetTracePhysVals(p, elmtToTrace[n][p],
3090  inarray + phys_offset,
3091  t_tmp = outarray + offset);
3092  }
3093  }
3094  }
3095 }
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:2240
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:492

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 2898 of file DisContField.cpp.

2901 {
2902  // Fill boundary conditions into missing elements
2903  if (PutFwdInBwdOnBCs) // just set Bwd value to be Fwd value on BCs
2904  {
2905  // Fill boundary conditions into missing elements
2906  for (int n = 0, cnt = 0; n < m_bndCondExpansions.size(); ++n)
2907  {
2908  if (m_bndConditions[n]->GetBoundaryConditionType() ==
2910  {
2911  auto ne = m_bndCondExpansions[n]->GetExpSize();
2912  for (int e = 0; e < ne; ++e)
2913  {
2914  auto npts =
2915  m_bndCondExpansions[n]->GetExp(e)->GetTotPoints();
2916  auto id2 = m_trace->GetPhys_Offset(
2917  m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
2918  Vmath::Vcopy(npts, &Fwd[id2], 1, &Bwd[id2], 1);
2919  }
2920 
2921  cnt += ne;
2922  }
2923  else if (m_bndConditions[n]->GetBoundaryConditionType() ==
2925  m_bndConditions[n]->GetBoundaryConditionType() ==
2927  {
2928  auto ne = m_bndCondExpansions[n]->GetExpSize();
2929  for (int e = 0; e < ne; ++e)
2930  {
2931  auto npts =
2932  m_bndCondExpansions[n]->GetExp(e)->GetTotPoints();
2933  auto id2 = m_trace->GetPhys_Offset(
2934  m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
2935  Vmath::Vcopy(npts, &Fwd[id2], 1, &Bwd[id2], 1);
2936  }
2937  cnt += ne;
2938  }
2939  else if (m_bndConditions[n]->GetBoundaryConditionType() !=
2941  {
2942  ASSERTL0(false, "Method not set up for this "
2943  "boundary condition.");
2944  }
2945  }
2946  }
2947  else
2948  {
2949  for (int n = 0, cnt = 0; n < m_bndCondExpansions.size(); ++n)
2950  {
2951  if (m_bndConditions[n]->GetBoundaryConditionType() ==
2953  {
2954  auto ne = m_bndCondExpansions[n]->GetExpSize();
2955  for (int e = 0; e < ne; ++e)
2956  {
2957  auto npts =
2958  m_bndCondExpansions[n]->GetExp(e)->GetTotPoints();
2959  auto id1 = m_bndCondExpansions[n]->GetPhys_Offset(e);
2960  auto id2 = m_trace->GetPhys_Offset(
2961  m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
2962  Vmath::Vcopy(npts,
2963  &(m_bndCondExpansions[n]->GetPhys())[id1], 1,
2964  &Bwd[id2], 1);
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 =
2977  m_bndCondExpansions[n]->GetExp(e)->GetTotPoints();
2978  auto id1 = m_bndCondExpansions[n]->GetPhys_Offset(e);
2979  ASSERTL0((m_bndCondExpansions[n]->GetPhys())[id1] == 0.0,
2980  "Method not set up for non-zero "
2981  "Neumann boundary condition");
2982  auto id2 = m_trace->GetPhys_Offset(
2983  m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
2984  Vmath::Vcopy(npts, &Fwd[id2], 1, &Bwd[id2], 1);
2985  }
2986 
2987  cnt += ne;
2988  }
2989  else if (m_bndConditions[n]->GetBoundaryConditionType() ==
2991  {
2992  // Do nothing
2993  }
2994  else if (m_bndConditions[n]->GetBoundaryConditionType() !=
2996  {
2998  "Method not set up for this boundary "
2999  "condition.");
3000  }
3001  }
3002  }
3003 }
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:2195

References ASSERTL0, Nektar::SpatialDomains::eDirichlet, Nektar::ErrorUtil::efatal, Nektar::SpatialDomains::eNeumann, Nektar::SpatialDomains::eNotDefined, Nektar::SpatialDomains::ePeriodic, Nektar::SpatialDomains::eRobin, Nektar::MultiRegions::ExpList::GetPhys(), 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 3694 of file DisContField.cpp.

3696 {
3697  int cnt;
3698  int npts;
3699  int e = 0;
3700 
3701  // Fill boundary conditions into missing elements
3702  int id2 = 0;
3703 
3704  for (int n = cnt = 0; n < m_bndCondExpansions.size(); ++n)
3705  {
3706 
3707  if (m_bndConditions[n]->GetBoundaryConditionType() ==
3709  {
3710  for (e = 0; e < m_bndCondExpansions[n]->GetExpSize(); ++e)
3711  {
3712  npts = m_bndCondExpansions[n]->GetExp(e)->GetTotPoints();
3713  id2 = m_trace->GetPhys_Offset(
3714  m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
3715  Vmath::Fill(npts, m_bndCondBndWeight[n], &weightave[id2], 1);
3716  Vmath::Fill(npts, 0.0, &weightjmp[id2], 1);
3717  }
3718 
3719  cnt += e;
3720  }
3721  else if (m_bndConditions[n]->GetBoundaryConditionType() ==
3723  m_bndConditions[n]->GetBoundaryConditionType() ==
3725  {
3726  for (e = 0; e < m_bndCondExpansions[n]->GetExpSize(); ++e)
3727  {
3728  npts = m_bndCondExpansions[n]->GetExp(e)->GetTotPoints();
3729  id2 = m_trace->GetPhys_Offset(
3730  m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
3731  Vmath::Fill(npts, m_bndCondBndWeight[n], &weightave[id2], 1);
3732  Vmath::Fill(npts, 0.0, &weightjmp[id2], 1);
3733  }
3734 
3735  cnt += e;
3736  }
3737  else if (m_bndConditions[n]->GetBoundaryConditionType() !=
3739  {
3741  "Method not set up for this boundary condition.");
3742  }
3743  }
3744 }

References Nektar::SpatialDomains::eDirichlet, Nektar::ErrorUtil::efatal, Nektar::SpatialDomains::eNeumann, Nektar::SpatialDomains::ePeriodic, Nektar::SpatialDomains::eRobin, Vmath::Fill(), 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 390 of file DisContField.h.

391 {
392  return m_bndCondBndWeight;
393 }

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 3903 of file DisContField.cpp.

3906 {
3907  int n, cnt, nq;
3908  int offsetOld, offsetNew;
3909  std::vector<unsigned int> eIDs;
3910 
3911  Array<OneD, int> ElmtID, TraceID;
3912  GetBoundaryToElmtMap(ElmtID, TraceID);
3913 
3914  // Skip other boundary regions
3915  for (cnt = n = 0; n < i; ++n)
3916  {
3917  cnt += m_bndCondExpansions[n]->GetExpSize();
3918  }
3919 
3920  // Populate eIDs with information from BoundaryToElmtMap
3921  for (n = 0; n < m_bndCondExpansions[i]->GetExpSize(); ++n)
3922  {
3923  eIDs.push_back(ElmtID[cnt + n]);
3924  }
3925 
3926  // Create expansion list
3928  *this, eIDs, DeclareCoeffPhysArrays, Collections::eNoCollection);
3929 
3930  // Copy phys and coeffs to new explist
3931  if (DeclareCoeffPhysArrays)
3932  {
3933  Array<OneD, NekDouble> tmp1, tmp2;
3934  for (n = 0; n < result->GetExpSize(); ++n)
3935  {
3936  nq = GetExp(ElmtID[cnt + n])->GetTotPoints();
3937  offsetOld = GetPhys_Offset(ElmtID[cnt + n]);
3938  offsetNew = result->GetPhys_Offset(n);
3939  Vmath::Vcopy(nq, tmp1 = GetPhys() + offsetOld, 1,
3940  tmp2 = result->UpdatePhys() + offsetNew, 1);
3941 
3942  nq = GetExp(ElmtID[cnt + n])->GetNcoeffs();
3943  offsetOld = GetCoeff_Offset(ElmtID[cnt + n]);
3944  offsetNew = result->GetCoeff_Offset(n);
3945  Vmath::Vcopy(nq, tmp1 = GetCoeffs() + offsetOld, 1,
3946  tmp2 = result->UpdateCoeffs() + offsetNew, 1);
3947  }
3948  }
3949 }
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:2100
const std::shared_ptr< LocalRegions::ExpansionVector > GetExp() const
This function returns the vector of elements in the expansion.
Definition: ExpList.h:2223

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::Collections::eNoCollection, 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 3748 of file DisContField.cpp.

3750 {
3751 
3752  if (m_BCtoElmMap.size() == 0)
3753  {
3754  switch (m_expType)
3755  {
3756  case e1D:
3757  {
3758  map<int, int> VertGID;
3759  int i, n, id;
3760  int bid, cnt, Vid;
3761  int nbcs = m_bndConditions.size();
3762 
3763  // make sure arrays are of sufficient length
3764  m_BCtoElmMap = Array<OneD, int>(nbcs, -1);
3765  m_BCtoTraceMap = Array<OneD, int>(nbcs);
3766 
3767  // setup map of all global ids along boundary
3768  cnt = 0;
3769  for (n = 0; n < m_bndCondExpansions.size(); ++n)
3770  {
3771  Vid = m_bndCondExpansions[n]
3772  ->GetExp(0)
3773  ->GetGeom()
3774  ->GetVertex(0)
3775  ->GetVid();
3776  VertGID[Vid] = cnt++;
3777  }
3778 
3779  // Loop over elements and find verts that match;
3781  for (cnt = n = 0; n < GetExpSize(); ++n)
3782  {
3783  exp = (*m_exp)[n];
3784  for (i = 0; i < exp->GetNverts(); ++i)
3785  {
3786  id = exp->GetGeom()->GetVid(i);
3787 
3788  if (VertGID.count(id) > 0)
3789  {
3790  bid = VertGID.find(id)->second;
3791  ASSERTL1(m_BCtoElmMap[bid] == -1,
3792  "Edge already set");
3793  m_BCtoElmMap[bid] = n;
3794  m_BCtoTraceMap[bid] = i;
3795  cnt++;
3796  }
3797  }
3798  }
3799  ASSERTL1(cnt == nbcs,
3800  "Failed to visit all boundary condtiions");
3801  }
3802  break;
3803  case e2D:
3804  {
3805  map<int, int> globalIdMap;
3806  int i, n;
3807  int cnt;
3808  int nbcs = 0;
3809 
3810  // Populate global ID map (takes global geometry
3811  // ID to local expansion list ID).
3812  for (i = 0; i < GetExpSize(); ++i)
3813  {
3814  globalIdMap[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
3815  }
3816 
3817  // Determine number of boundary condition expansions.
3818  for (i = 0; i < m_bndConditions.size(); ++i)
3819  {
3820  nbcs += m_bndCondExpansions[i]->GetExpSize();
3821  }
3822 
3823  // Initialize arrays
3824  m_BCtoElmMap = Array<OneD, int>(nbcs);
3825  m_BCtoTraceMap = Array<OneD, int>(nbcs);
3826 
3828  cnt = 0;
3829  for (n = 0; n < m_bndCondExpansions.size(); ++n)
3830  {
3831  for (i = 0; i < m_bndCondExpansions[n]->GetExpSize();
3832  ++i, ++cnt)
3833  {
3834  exp1d = m_bndCondExpansions[n]
3835  ->GetExp(i)
3836  ->as<LocalRegions::Expansion1D>();
3837 
3838  // Use edge to element map from MeshGraph.
3840  m_graph->GetElementsFromEdge(exp1d->GetGeom1D());
3841 
3842  m_BCtoElmMap[cnt] =
3843  globalIdMap[(*tmp)[0].first->GetGlobalID()];
3844  m_BCtoTraceMap[cnt] = (*tmp)[0].second;
3845  }
3846  }
3847  }
3848  break;
3849  case e3D:
3850  {
3851  map<int, int> globalIdMap;
3852  int i, n;
3853  int cnt;
3854  int nbcs = 0;
3855 
3856  // Populate global ID map (takes global geometry ID to local
3857  // expansion list ID).
3859  for (i = 0; i < GetExpSize(); ++i)
3860  {
3861  globalIdMap[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
3862  }
3863 
3864  // Determine number of boundary condition expansions.
3865  for (i = 0; i < m_bndConditions.size(); ++i)
3866  {
3867  nbcs += m_bndCondExpansions[i]->GetExpSize();
3868  }
3869 
3870  // Initialize arrays
3871  m_BCtoElmMap = Array<OneD, int>(nbcs);
3872  m_BCtoTraceMap = Array<OneD, int>(nbcs);
3873 
3875  for (cnt = n = 0; n < m_bndCondExpansions.size(); ++n)
3876  {
3877  for (i = 0; i < m_bndCondExpansions[n]->GetExpSize();
3878  ++i, ++cnt)
3879  {
3880  exp2d = m_bndCondExpansions[n]
3881  ->GetExp(i)
3882  ->as<LocalRegions::Expansion2D>();
3883 
3885  m_graph->GetElementsFromFace(exp2d->GetGeom2D());
3886  m_BCtoElmMap[cnt] =
3887  globalIdMap[tmp->at(0).first->GetGlobalID()];
3888  m_BCtoTraceMap[cnt] = tmp->at(0).second;
3889  }
3890  }
3891  }
3892  break;
3893  default:
3894  ASSERTL1(false, "Not setup for this expansion");
3895  break;
3896  }
3897  }
3898 
3899  ElmtID = m_BCtoElmMap;
3900  TraceID = m_BCtoTraceMap;
3901 }
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:50

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 384 of file DisContField.h.

386 {
387  v_GetFwdBwdTracePhys(m_phys, Fwd, Bwd);
388 }
virtual void v_GetFwdBwdTracePhys(Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd)
Definition: DisContField.h:384

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 2828 of file DisContField.cpp.

2834 {
2835  // Is this zeroing necessary?
2836  // Zero forward/backward vectors.
2837  Vmath::Zero(Fwd.size(), Fwd, 1);
2838  Vmath::Zero(Bwd.size(), Bwd, 1);
2839 
2840  // Basis definition on each element
2841  LibUtilities::BasisSharedPtr basis = (*m_exp)[0]->GetBasis(0);
2842  if ((basis->GetBasisType() != LibUtilities::eGauss_Lagrange))
2843  {
2844  // blocked routine
2845  Array<OneD, NekDouble> tracevals(
2846  m_locTraceToTraceMap->GetNLocTracePts());
2847 
2848  m_locTraceToTraceMap->LocTracesFromField(field, tracevals);
2849  m_locTraceToTraceMap->InterpLocTracesToTrace(0, tracevals, Fwd);
2850 
2851  Array<OneD, NekDouble> invals =
2852  tracevals + m_locTraceToTraceMap->GetNFwdLocTracePts();
2853  m_locTraceToTraceMap->InterpLocTracesToTrace(1, invals, Bwd);
2854  }
2855  else
2856  {
2857  // Loop over elements and collect forward expansion
2858  auto nexp = GetExpSize();
2859  Array<OneD, NekDouble> e_tmp;
2861 
2862  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>>
2863  &elmtToTrace = m_traceMap->GetElmtToTrace();
2864 
2865  for (int n = 0, cnt = 0; n < nexp; ++n)
2866  {
2867  exp = (*m_exp)[n];
2868  auto phys_offset = GetPhys_Offset(n);
2869 
2870  for (int e = 0; e < exp->GetNtraces(); ++e, ++cnt)
2871  {
2872  auto offset =
2873  m_trace->GetPhys_Offset(elmtToTrace[n][e]->GetElmtId());
2874 
2875  e_tmp =
2876  (m_leftAdjacentTraces[cnt]) ? Fwd + offset : Bwd + offset;
2877 
2878  exp->GetTracePhysVals(e, elmtToTrace[n][e], field + phys_offset,
2879  e_tmp);
2880  }
2881  }
2882  }
2883 
2885 
2886  if (FillBnd)
2887  {
2888  v_FillBwdWithBoundCond(Fwd, Bwd, PutFwdInBwdOnBCs);
2889  }
2890 
2891  if (DoExchange)
2892  {
2893  // Do parallel exchange for forwards/backwards spaces.
2894  m_traceMap->GetAssemblyCommDG()->PerformExchange(Fwd, Bwd);
2895  }
2896 }
virtual void v_PeriodicBwdCopy(const Array< OneD, const NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd)
Definition: DisContField.h:400
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_GetLeftAdjacentTraces()

std::vector< bool > & Nektar::MultiRegions::DisContField::v_GetLeftAdjacentTraces ( void  )
inlineprotectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 372 of file DisContField.h.

373 {
374  return m_leftAdjacentTraces;
375 }

References m_leftAdjacentTraces.

◆ 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 4306 of file DisContField.cpp.

4310 {
4311  if (locTraceFwd != NullNekDouble1DArray)
4312  {
4313  m_locTraceToTraceMap->InterpLocTracesToTraceTranspose(0, Fwd,
4314  locTraceFwd);
4315  }
4316 
4317  if (locTraceBwd != NullNekDouble1DArray)
4318  {
4319  m_locTraceToTraceMap->InterpLocTracesToTraceTranspose(1, Bwd,
4320  locTraceBwd);
4321  }
4322 }
static Array< OneD, NekDouble > NullNekDouble1DArray

References m_locTraceToTraceMap, and Nektar::NullNekDouble1DArray.

◆ v_GetLocTraceToTraceMap()

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

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 225 of file DisContField.h.

227  {
228  return m_locTraceToTraceMap;
229  }

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 350 of file DisContField.h.

353  {
354  periodicVerts = m_periodicVerts;
355  periodicEdges = m_periodicEdges;
356  periodicFaces = m_periodicFaces;
357  }

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 3979 of file DisContField.cpp.

3980 {
3981  int i, cnt;
3982  map<int, RobinBCInfoSharedPtr> returnval;
3983  Array<OneD, int> ElmtID, TraceID;
3984  GetBoundaryToElmtMap(ElmtID, TraceID);
3985 
3986  for (cnt = i = 0; i < m_bndCondExpansions.size(); ++i)
3987  {
3988  MultiRegions::ExpListSharedPtr locExpList;
3989 
3990  if (m_bndConditions[i]->GetBoundaryConditionType() ==
3992  {
3993  int e, elmtid;
3994  Array<OneD, NekDouble> Array_tmp;
3995 
3996  locExpList = m_bndCondExpansions[i];
3997 
3998  int npoints = locExpList->GetNpoints();
3999  Array<OneD, NekDouble> x0(npoints, 0.0);
4000  Array<OneD, NekDouble> x1(npoints, 0.0);
4001  Array<OneD, NekDouble> x2(npoints, 0.0);
4002  Array<OneD, NekDouble> coeffphys(npoints);
4003 
4004  locExpList->GetCoords(x0, x1, x2);
4005 
4006  LibUtilities::Equation coeffeqn =
4007  std::static_pointer_cast<
4008  SpatialDomains::RobinBoundaryCondition>(m_bndConditions[i])
4009  ->m_robinPrimitiveCoeff;
4010 
4011  // evalaute coefficient
4012  coeffeqn.Evaluate(x0, x1, x2, 0.0, coeffphys);
4013 
4014  for (e = 0; e < locExpList->GetExpSize(); ++e)
4015  {
4016  RobinBCInfoSharedPtr rInfo =
4018  TraceID[cnt + e],
4019  Array_tmp = coeffphys + locExpList->GetPhys_Offset(e));
4020 
4021  elmtid = ElmtID[cnt + e];
4022  // make link list if necessary
4023  if (returnval.count(elmtid) != 0)
4024  {
4025  rInfo->next = returnval.find(elmtid)->second;
4026  }
4027  returnval[elmtid] = rInfo;
4028  }
4029  }
4030  cnt += m_bndCondExpansions[i]->GetExpSize();
4031  }
4032 
4033  return returnval;
4034 }
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 210 of file DisContField.h.

211  {
213  {
214  SetUpDG();
215  }
216 
217  return m_trace;
218  }

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 220 of file DisContField.h.

221  {
222  return m_traceMap;
223  }

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 3321 of file DisContField.cpp.

3328 {
3329  boost::ignore_unused(varfactors, dirForcing);
3330  int i, n, cnt, nbndry;
3331  int nexp = GetExpSize();
3332 
3333  Array<OneD, NekDouble> f(m_ncoeffs);
3334  DNekVec F(m_ncoeffs, f, eWrapper);
3335  Array<OneD, NekDouble> e_f, e_l;
3336 
3337  //----------------------------------
3338  // Setup RHS Inner product if required
3339  //----------------------------------
3340  if (PhysSpaceForcing)
3341  {
3342  IProductWRTBase(inarray, f);
3343  Vmath::Neg(m_ncoeffs, f, 1);
3344  }
3345  else
3346  {
3347  Vmath::Smul(m_ncoeffs, -1.0, inarray, 1, f, 1);
3348  }
3349 
3350  //----------------------------------
3351  // Solve continuous Boundary System
3352  //----------------------------------
3353  int GloBndDofs = m_traceMap->GetNumGlobalBndCoeffs();
3354  int NumDirBCs = m_traceMap->GetNumLocalDirBndCoeffs();
3355  int e_ncoeffs;
3356 
3357  GlobalMatrixKey HDGLamToUKey(StdRegions::eHybridDGLamToU,
3358  NullAssemblyMapSharedPtr, factors, varcoeff);
3359  const DNekScalBlkMatSharedPtr &HDGLamToU = GetBlockMatrix(HDGLamToUKey);
3360 
3361  // Retrieve number of local trace space coefficients N_{\lambda},
3362  // and set up local elemental trace solution \lambda^e.
3363  int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs();
3364  Array<OneD, NekDouble> bndrhs(LocBndCoeffs, 0.0);
3365  Array<OneD, NekDouble> loclambda(LocBndCoeffs, 0.0);
3366  DNekVec LocLambda(LocBndCoeffs, loclambda, eWrapper);
3367 
3368  //----------------------------------
3369  // Evaluate Trace Forcing
3370  // Kirby et al, 2010, P23, Step 5.
3371  //----------------------------------
3372  // Determing <u_lam,f> terms using HDGLamToU matrix
3373  for (cnt = n = 0; n < nexp; ++n)
3374  {
3375  nbndry = (*m_exp)[n]->NumDGBndryCoeffs();
3376 
3377  e_ncoeffs = (*m_exp)[n]->GetNcoeffs();
3378  e_f = f + m_coeff_offset[n];
3379  e_l = bndrhs + cnt;
3380 
3381  // use outarray as tmp space
3382  DNekVec Floc(nbndry, e_l, eWrapper);
3383  DNekVec ElmtFce(e_ncoeffs, e_f, eWrapper);
3384  Floc = Transpose(*(HDGLamToU->GetBlock(n, n))) * ElmtFce;
3385 
3386  cnt += nbndry;
3387  }
3388 
3389  Array<OneD, const int> bndCondMap =
3390  m_traceMap->GetBndCondCoeffsToLocalTraceMap();
3391  Array<OneD, const NekDouble> Sign = m_traceMap->GetLocalToGlobalBndSign();
3392 
3393  // Copy Dirichlet boundary conditions and weak forcing
3394  // into trace space
3395  int locid;
3396  cnt = 0;
3397  for (i = 0; i < m_bndCondExpansions.size(); ++i)
3398  {
3399  Array<OneD, const NekDouble> bndcoeffs =
3400  m_bndCondExpansions[i]->GetCoeffs();
3401 
3402  if (m_bndConditions[i]->GetBoundaryConditionType() ==
3404  {
3405  for (int j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); ++j)
3406  {
3407  locid = bndCondMap[cnt + j];
3408  loclambda[locid] = Sign[locid] * bndcoeffs[j];
3409  }
3410  }
3411  else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3413  m_bndConditions[i]->GetBoundaryConditionType() ==
3415  {
3416  // Add weak boundary condition to trace forcing
3417  for (int j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); ++j)
3418  {
3419  locid = bndCondMap[cnt + j];
3420  bndrhs[locid] += Sign[locid] * bndcoeffs[j];
3421  }
3422  }
3423 
3424  cnt += (m_bndCondExpansions[i])->GetNcoeffs();
3425  }
3426 
3427  //----------------------------------
3428  // Solve trace problem: \Lambda = K^{-1} F
3429  // K is the HybridDGHelmBndLam matrix.
3430  //----------------------------------
3431  if (GloBndDofs - NumDirBCs > 0)
3432  {
3433  GlobalLinSysKey key(StdRegions::eHybridDGHelmBndLam, m_traceMap,
3434  factors, varcoeff);
3436  LinSys->Solve(bndrhs, loclambda, m_traceMap);
3437 
3438  // For consistency with previous version put global
3439  // solution into m_trace->m_coeffs
3440  m_traceMap->LocalToGlobal(loclambda, m_trace->UpdateCoeffs());
3441  }
3442 
3443  //----------------------------------
3444  // Internal element solves
3445  //----------------------------------
3446  GlobalMatrixKey invHDGhelmkey(StdRegions::eInvHybridDGHelmholtz,
3447  NullAssemblyMapSharedPtr, factors, varcoeff);
3448 
3449  const DNekScalBlkMatSharedPtr &InvHDGHelm = GetBlockMatrix(invHDGhelmkey);
3450  DNekVec out(m_ncoeffs, outarray, eWrapper);
3451  Vmath::Zero(m_ncoeffs, outarray, 1);
3452 
3453  // out = u_f + u_lam = (*InvHDGHelm)*f + (LamtoU)*Lam
3454  out = (*InvHDGHelm) * F + (*HDGLamToU) * LocLambda;
3455 }
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)
This function calculates the inner product of a function with respect to all local expansion modes .
Definition: ExpList.h:1787
int GetNcoeffs(void) const
Returns the total number of local degrees of freedom .
Definition: ExpList.h:1641
const DNekScalBlkMatSharedPtr & GetBlockMatrix(const GlobalMatrixKey &gkey)
Definition: ExpList.cpp:2010
static AssemblyMapSharedPtr NullAssemblyMapSharedPtr
Definition: AssemblyMap.h:53
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:79
NekVector< NekDouble > DNekVec
Definition: NekTypeDefs.hpp:48
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:518
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:248

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 400 of file DisContField.h.

402 {
403  for (int n = 0; n < m_periodicFwdCopy.size(); ++n)
404  {
405  Bwd[m_periodicBwdCopy[n]] = Fwd[m_periodicFwdCopy[n]];
406  }
407 }

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 3954 of file DisContField.cpp.

3955 {
3956  ExpList::v_Reset();
3957 
3958  // Reset boundary condition expansions.
3959  for (int n = 0; n < m_bndCondExpansions.size(); ++n)
3960  {
3961  m_bndCondExpansions[n]->Reset();
3962  }
3963 }
virtual void v_Reset()
Reset geometry information, metrics, matrix managers and geometry information.
Definition: ExpList.cpp:2670

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 395 of file DisContField.h.

396 {
397  m_bndCondBndWeight[index] = value;
398 }

References m_bndCondBndWeight.

◆ v_UpdateBndCondExpansion()

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

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 287 of file DisContField.h.

288  {
289  return m_bndCondExpansions[i];
290  }

References m_bndCondExpansions.

◆ v_UpdateBndConditions()

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

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 293 of file DisContField.h.

294  {
295  return m_bndConditions;
296  }

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 163 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 152 of file DisContField.h.

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

◆ m_leftAdjacentTraces

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

◆ m_locTraceToTraceMap

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

◆ m_negatedFluxNormal

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

Definition at line 360 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 130 of file DisContField.h.

◆ m_periodicBwdCopy

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

Definition at line 186 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 173 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 178 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 185 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 168 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