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 > &coeffs, 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 (const Array< OneD, const NekDouble > &coeffs, 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...
 
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 (int i, NekDouble val)
 Set the i th value of m_phys to value val. 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...
 
GlobalLinSysKey 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...
 
GlobalLinSysKey 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 int npts, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool Shuff=true, bool UnShuff=true)
 
void HomogeneousBwdTrans (const int npts, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool Shuff=true, bool UnShuff=true)
 
void DealiasedProd (const int num_dofs, const Array< OneD, NekDouble > &inarray1, const Array< OneD, NekDouble > &inarray2, Array< OneD, NekDouble > &outarray)
 
void DealiasedDotProd (const int num_dofs, 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 (const Array< OneD, NekDouble > coeffs)
 Fill Bnd Condition expansion from the values stored in expansion. More...
 
void FillBndCondFromField (const int nreg, const Array< OneD, NekDouble > coeffs)
 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_\infty\) error of the global 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 ()
 Calculates the \(H^1\) error of the global spectral/hp element approximation. More...
 
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...
 
size_t 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::ExpansionSharedPtrGetExpFromGeomId (int n)
 This function returns (a shared pointer to) the local elemental expansion of the \(n^{\mathrm{th}}\) element given a global geometry ID. 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 local contiguous list of coeffs correspoinding to element n. More...
 
int GetPhys_Offset (int n) const
 Get the start offset position for a local contiguous list of quadrature points in a full array 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)
 
virtual 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, std::unordered_map< int, int > zIdToPlane=std::unordered_map< int, int >())
 Extract the data in fielddata into the coeffs. More...
 
void ExtractCoeffsFromFile (const std::string &fileName, LibUtilities::CommSharedPtr comm, const std::string &varName, Array< OneD, NekDouble > &coeffs)
 
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)
 
int GetPoolCount (std::string)
 
void UnsetGlobalLinSys (GlobalLinSysKey, bool)
 
LibUtilities::NekManager< GlobalLinSysKey, GlobalLinSys > & GetGlobalLinSysManager (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)
 
const std::unordered_map< int, int > & GetElmtToExpId (void)
 This function returns the map of index inside m_exp to geom id. More...
 
int GetElmtToExpId (int elmtId)
 This function returns the index inside m_exp for a given geom id. More...
 

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...
 
void SetUpDG (const std::string="DefaultVar", const Collections::ImplementationType ImpType=Collections::eNoImpType)
 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 ExpListSharedPtrv_GetTrace () override
 
virtual AssemblyMapDGSharedPtrv_GetTraceMap (void) override
 
virtual const LocTraceToTraceMapSharedPtrv_GetLocTraceToTraceMap (void) const override
 
virtual std::vector< bool > & v_GetLeftAdjacentTraces (void) override
 
virtual void v_AddTraceIntegral (const Array< OneD, const NekDouble > &Fn, Array< OneD, NekDouble > &outarray) override
 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) override
 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) override
 
virtual void v_ExtractTracePhys (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
 This method extracts the trace (verts in 1D) from the field inarray and puts the values in outarray. More...
 
virtual void v_ExtractTracePhys (Array< OneD, NekDouble > &outarray) override
 
virtual void v_GetLocTraceFromTracePts (const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &locTraceFwd, Array< OneD, NekDouble > &locTraceBwd) override
 
void GenerateFieldBnd1D (SpatialDomains::BoundaryConditions &bcs, const std::string variable)
 
virtual std::map< int, RobinBCInfoSharedPtrv_GetRobinBCInfo () override
 
virtual const Array< OneD, const MultiRegions::ExpListSharedPtr > & v_GetBndCondExpansions () override
 
virtual const Array< OneD, const SpatialDomains::BoundaryConditionShPtr > & v_GetBndConditions () override
 
virtual MultiRegions::ExpListSharedPtrv_UpdateBndCondExpansion (int i) override
 
virtual Array< OneD, SpatialDomains::BoundaryConditionShPtr > & v_UpdateBndConditions () override
 
virtual void v_GetBoundaryToElmtMap (Array< OneD, int > &ElmtID, Array< OneD, int > &TraceID) override
 
virtual void v_GetBndElmtExpansion (int i, std::shared_ptr< ExpList > &result, const bool DeclareCoeffPhysArrays) override
 
virtual void v_Reset () override
 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) override
 Evaluate all boundary conditions at a given time.. More...
 
virtual GlobalLinSysKey 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) override
 Solve the Helmholtz equation. More...
 
virtual void v_PeriodicBwdCopy (const Array< OneD, const NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd) override
 
virtual void v_FillBwdWithBwdWeight (Array< OneD, NekDouble > &weightave, Array< OneD, NekDouble > &weightjmp) override
 Fill the weight with m_bndCondBndWeight. More...
 
virtual void v_GetFwdBwdTracePhys (Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd) override
 
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) override
 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) override
 
virtual const Array< OneD, const NekDouble > & v_GetBndCondBwdWeight () override
 
virtual void v_SetBndCondBwdWeight (const int index, const NekDouble value) override
 
virtual void v_GetPeriodicEntities (PeriodicMap &periodicVerts, PeriodicMap &periodicEdges, PeriodicMap &periodicFaces) override
 Obtain a copy of the periodic edges and vertices for this field. More...
 
virtual void v_AddTraceIntegralToOffDiag (const Array< OneD, const NekDouble > &FwdFlux, const Array< OneD, const NekDouble > &BwdFlux, Array< OneD, NekDouble > &outarray) override
 
- 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 size_t 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 GlobalLinSysKey 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 (const Array< OneD, NekDouble > coeffs)
 
virtual void v_FillBndCondFromField (const int nreg, const Array< OneD, NekDouble > coeffs)
 
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=NullNekDouble1DArray, 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 int npts, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool Shuff=true, bool UnShuff=true)
 
virtual void v_HomogeneousBwdTrans (const int npts, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool Shuff=true, bool UnShuff=true)
 
virtual void v_DealiasedProd (const int num_dofs, const Array< OneD, NekDouble > &inarray1, const Array< OneD, NekDouble > &inarray2, Array< OneD, NekDouble > &outarray)
 
virtual void v_DealiasedDotProd (const int num_dofs, 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, std::unordered_map< int, int > zIdToPlane)
 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_WriteVtkPieceData (std::ostream &outfile, int expansion, std::string var)
 
virtual void v_WriteVtkPieceHeader (std::ostream &outfile, int expansion, int istrip)
 
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)
 
virtual int v_GetPoolCount (std::string)
 
virtual void v_UnsetGlobalLinSys (GlobalLinSysKey, bool)
 
virtual LibUtilities::NekManager< GlobalLinSysKey, GlobalLinSys > & v_GetGlobalLinSysManager (void)
 
void ExtractFileBCs (const std::string &fileName, LibUtilities::CommSharedPtr comm, const std::string &varName, const std::shared_ptr< ExpList > locExpList)
 
virtual LibUtilities::BasisSharedPtr v_GetHomogeneousBasis (void)
 
virtual void v_SetHomo1DSpecVanVisc (Array< OneD, NekDouble > visc)
 
virtual std::shared_ptr< ExpList > & v_GetPlane (int n)
 

Protected Attributes

size_t m_numDirBndCondExpansions
 The number of boundary segments on which Dirichlet boundary conditions are imposed. More...
 
Array< OneD, SpatialDomains::BoundaryConditionShPtrm_bndConditions
 An array which contains the information about the boundary condition structure definition on the different boundary regions. More...
 
Array< OneD, MultiRegions::ExpListSharedPtrm_bndCondExpansions
 An object which contains the discretised boundary conditions. More...
 
Array< OneD, NekDoublem_bndCondBndWeight
 
InterfaceMapDGSharedPtr m_interfaceMap
 Interfaces mapping for trace space. 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)
 

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 57 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 structure definition on the diff...
Definition: DisContField.h:137
ExpListSharedPtr m_trace
Trace space storage for points between elements.
Definition: DisContField.h:161
Array< OneD, MultiRegions::ExpListSharedPtr > m_bndCondExpansions
An object which contains the discretised boundary conditions.
Definition: DisContField.h:150
ExpList(const ExpansionType Type=eNoType)
The default constructor using a type.
Definition: ExpList.cpp:102
static ExpListSharedPtr NullExpListSharedPtr
Definition: ExpList.h:1479

◆ 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  int i, cnt;
116  Array<OneD, int> ElmtID, TraceID;
117  GetBoundaryToElmtMap(ElmtID, TraceID);
118 
119  for (cnt = i = 0; i < m_bndCondExpansions.size(); ++i)
120  {
122  locExpList = m_bndCondExpansions[i];
123 
124  for (int e = 0; e < locExpList->GetExpSize(); ++e)
125  {
126  (*m_exp)[ElmtID[cnt + e]]->SetTraceExp(TraceID[cnt + e],
127  locExpList->GetExp(e));
128  locExpList->GetExp(e)->SetAdjacentElementExp(
129  TraceID[cnt + e], (*m_exp)[ElmtID[cnt + e]]);
130  }
131  cnt += m_bndCondExpansions[i]->GetExpSize();
132  }
133 
134  if (SetUpJustDG)
135  {
136  SetUpDG(variable, ImpType);
137  }
138  else
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 SetUpDG(const std::string="DefaultVar", const Collections::ImplementationType ImpType=Collections::eNoImpType)
Set up all DG member variables and maps.
void FindPeriodicTraces(const SpatialDomains::BoundaryConditions &bcs, const std::string variable)
Generate a associative map of periodic vertices in a mesh.
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:2222
std::shared_ptr< LocalRegions::ExpansionVector > m_exp
The list of local expansions.
Definition: ExpList.h:1111
void EvaluateBoundaryConditions(const NekDouble time=0.0, const std::string varName="", const NekDouble=NekConstants::kNekUnsetDouble, const NekDouble=NekConstants::kNekUnsetDouble)
Definition: ExpList.h:2211
void ApplyGeomInfo()
Apply geometry information to each expansion.
Definition: ExpList.cpp:2661
LibUtilities::SessionReaderSharedPtr m_session
Session.
Definition: ExpList.h:1051
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 610 of file DisContField.cpp.

618  : ExpList(pSession, domain, graph1D, true, variable, SetToOneSpaceDimension,
619  comm, ImpType)
620 {
621  if (variable.compare("DefaultVar") != 0)
622  {
624  GetDomainBCs(domain, Allbcs, variable);
625 
626  GenerateBoundaryConditionExpansion(m_graph, *DomBCs, variable);
627  EvaluateBoundaryConditions(0.0, variable);
628  ApplyGeomInfo();
629  FindPeriodicTraces(*DomBCs, variable);
630  }
631 
632  SetUpDG(variable);
633 }
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:1053
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 639 of file DisContField.cpp.

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

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

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

787  : ExpList(In)
788 {
789 }

◆ ~DisContField()

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

Destructor.

Definition at line 794 of file DisContField.cpp.

795 {
796 }

Member Function Documentation

◆ EvaluateHDGPostProcessing()

void Nektar::MultiRegions::DisContField::EvaluateHDGPostProcessing ( const Array< OneD, const NekDouble > &  coeffs,
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 4204 of file DisContField.cpp.

4207 {
4208  int i, cnt, e, ncoeff_trace;
4209  Array<OneD, NekDouble> force, out_tmp, qrhs, qrhs1;
4210  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
4211  m_traceMap->GetElmtToTrace();
4212 
4213  StdRegions::Orientation edgedir;
4214 
4215  int nq_elmt, nm_elmt;
4216  int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs();
4217  Array<OneD, NekDouble> loc_lambda(LocBndCoeffs), trace_lambda;
4218  Array<OneD, NekDouble> tmp_coeffs;
4219  m_traceMap->GlobalToLocalBnd(m_trace->GetCoeffs(), loc_lambda);
4220 
4221  trace_lambda = loc_lambda;
4222 
4223  int dim = (m_expType == e2D) ? 2 : 3;
4224 
4225  int num_points[] = {0, 0, 0};
4226  int num_modes[] = {0, 0, 0};
4227 
4228  // Calculate Q using standard DG formulation.
4229  for (i = cnt = 0; i < GetExpSize(); ++i)
4230  {
4231  nq_elmt = (*m_exp)[i]->GetTotPoints();
4232  nm_elmt = (*m_exp)[i]->GetNcoeffs();
4233  qrhs = Array<OneD, NekDouble>(nq_elmt);
4234  qrhs1 = Array<OneD, NekDouble>(nq_elmt);
4235  force = Array<OneD, NekDouble>(2 * nm_elmt);
4236  out_tmp = force + nm_elmt;
4238 
4239  for (int j = 0; j < dim; ++j)
4240  {
4241  num_points[j] = (*m_exp)[i]->GetBasis(j)->GetNumPoints();
4242  num_modes[j] = (*m_exp)[i]->GetBasis(j)->GetNumModes();
4243  }
4244 
4245  // Probably a better way of setting up lambda than this. Note
4246  // cannot use PutCoeffsInToElmts since lambda space is mapped
4247  // during the solve.
4248  int nTraces = (*m_exp)[i]->GetNtraces();
4249  Array<OneD, Array<OneD, NekDouble>> traceCoeffs(nTraces);
4250 
4251  for (e = 0; e < (*m_exp)[i]->GetNtraces(); ++e)
4252  {
4253  edgedir = (*m_exp)[i]->GetTraceOrient(e);
4254  ncoeff_trace = elmtToTrace[i][e]->GetNcoeffs();
4255  traceCoeffs[e] = Array<OneD, NekDouble>(ncoeff_trace);
4256  Vmath::Vcopy(ncoeff_trace, trace_lambda, 1, traceCoeffs[e], 1);
4257  if (dim == 2)
4258  {
4259  elmtToTrace[i][e]->SetCoeffsToOrientation(
4260  edgedir, traceCoeffs[e], traceCoeffs[e]);
4261  }
4262  else
4263  {
4264  (*m_exp)[i]
4265  ->as<LocalRegions::Expansion3D>()
4266  ->SetFaceToGeomOrientation(e, traceCoeffs[e]);
4267  }
4268  trace_lambda = trace_lambda + ncoeff_trace;
4269  }
4270 
4271  // creating orthogonal expansion (checking if we have quads or
4272  // triangles)
4273  LibUtilities::ShapeType shape = (*m_exp)[i]->DetShapeType();
4274  switch (shape)
4275  {
4277  {
4278  const LibUtilities::PointsKey PkeyQ1(
4279  num_points[0], LibUtilities::eGaussLobattoLegendre);
4280  const LibUtilities::PointsKey PkeyQ2(
4281  num_points[1], LibUtilities::eGaussLobattoLegendre);
4282  LibUtilities::BasisKey BkeyQ1(LibUtilities::eOrtho_A,
4283  num_modes[0], PkeyQ1);
4284  LibUtilities::BasisKey BkeyQ2(LibUtilities::eOrtho_A,
4285  num_modes[1], PkeyQ2);
4287  std::dynamic_pointer_cast<SpatialDomains::QuadGeom>(
4288  (*m_exp)[i]->GetGeom());
4290  BkeyQ1, BkeyQ2, qGeom);
4291  }
4292  break;
4294  {
4295  const LibUtilities::PointsKey PkeyT1(
4296  num_points[0], LibUtilities::eGaussLobattoLegendre);
4297  const LibUtilities::PointsKey PkeyT2(
4298  num_points[1], LibUtilities::eGaussRadauMAlpha1Beta0);
4299  LibUtilities::BasisKey BkeyT1(LibUtilities::eOrtho_A,
4300  num_modes[0], PkeyT1);
4301  LibUtilities::BasisKey BkeyT2(LibUtilities::eOrtho_B,
4302  num_modes[1], PkeyT2);
4304  std::dynamic_pointer_cast<SpatialDomains::TriGeom>(
4305  (*m_exp)[i]->GetGeom());
4307  BkeyT1, BkeyT2, tGeom);
4308  }
4309  break;
4311  {
4312  const LibUtilities::PointsKey PkeyH1(
4313  num_points[0], LibUtilities::eGaussLobattoLegendre);
4314  const LibUtilities::PointsKey PkeyH2(
4315  num_points[1], LibUtilities::eGaussLobattoLegendre);
4316  const LibUtilities::PointsKey PkeyH3(
4317  num_points[2], LibUtilities::eGaussLobattoLegendre);
4318  LibUtilities::BasisKey BkeyH1(LibUtilities::eOrtho_A,
4319  num_modes[0], PkeyH1);
4320  LibUtilities::BasisKey BkeyH2(LibUtilities::eOrtho_A,
4321  num_modes[1], PkeyH2);
4322  LibUtilities::BasisKey BkeyH3(LibUtilities::eOrtho_A,
4323  num_modes[2], PkeyH3);
4325  std::dynamic_pointer_cast<SpatialDomains::HexGeom>(
4326  (*m_exp)[i]->GetGeom());
4328  BkeyH1, BkeyH2, BkeyH3, hGeom);
4329  }
4330  break;
4332  {
4333  const LibUtilities::PointsKey PkeyT1(
4334  num_points[0], LibUtilities::eGaussLobattoLegendre);
4335  const LibUtilities::PointsKey PkeyT2(
4336  num_points[1], LibUtilities::eGaussRadauMAlpha1Beta0);
4337  const LibUtilities::PointsKey PkeyT3(
4338  num_points[2], LibUtilities::eGaussRadauMAlpha2Beta0);
4339  LibUtilities::BasisKey BkeyT1(LibUtilities::eOrtho_A,
4340  num_modes[0], PkeyT1);
4341  LibUtilities::BasisKey BkeyT2(LibUtilities::eOrtho_B,
4342  num_modes[1], PkeyT2);
4343  LibUtilities::BasisKey BkeyT3(LibUtilities::eOrtho_C,
4344  num_modes[2], PkeyT3);
4346  std::dynamic_pointer_cast<SpatialDomains::TetGeom>(
4347  (*m_exp)[i]->GetGeom());
4349  BkeyT1, BkeyT2, BkeyT3, tGeom);
4350  }
4351  break;
4352  default:
4354  "Wrong shape type, HDG postprocessing is not "
4355  "implemented");
4356  };
4357 
4358  // DGDeriv
4359  // (d/dx w, d/dx q_0)
4360  (*m_exp)[i]->DGDeriv(0, tmp_coeffs = coeffs + m_coeff_offset[i],
4361  elmtToTrace[i], traceCoeffs, out_tmp);
4362  (*m_exp)[i]->BwdTrans(out_tmp, qrhs);
4363  ppExp->IProductWRTDerivBase(0, qrhs, force);
4364 
4365  // + (d/dy w, d/dy q_1)
4366  (*m_exp)[i]->DGDeriv(1, tmp_coeffs = coeffs + m_coeff_offset[i],
4367  elmtToTrace[i], traceCoeffs, out_tmp);
4368 
4369  (*m_exp)[i]->BwdTrans(out_tmp, qrhs);
4370  ppExp->IProductWRTDerivBase(1, qrhs, out_tmp);
4371 
4372  Vmath::Vadd(nm_elmt, force, 1, out_tmp, 1, force, 1);
4373 
4374  // determine force[0] = (1,u)
4375  (*m_exp)[i]->BwdTrans(tmp_coeffs = coeffs + m_coeff_offset[i], qrhs);
4376  force[0] = (*m_exp)[i]->Integral(qrhs);
4377 
4378  // multiply by inverse Laplacian matrix
4379  // get matrix inverse
4380  LocalRegions::MatrixKey lapkey(StdRegions::eInvLaplacianWithUnityMean,
4381  ppExp->DetShapeType(), *ppExp);
4382  DNekScalMatSharedPtr lapsys = ppExp->GetLocMatrix(lapkey);
4383 
4384  NekVector<NekDouble> in(nm_elmt, force, eWrapper);
4385  NekVector<NekDouble> out(nm_elmt);
4386 
4387  out = (*lapsys) * in;
4388 
4389  // Transforming back to modified basis
4390  Array<OneD, NekDouble> work(nq_elmt);
4391  ppExp->BwdTrans(out.GetPtr(), work);
4392  (*m_exp)[i]->FwdTrans(work, tmp_coeffs = outarray + m_coeff_offset[i]);
4393  }
4394 }
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition: ErrorUtil.hpp:209
int GetExpSize(void)
This function returns the number of elements in the expansion.
Definition: ExpList.h:1996
Array< OneD, int > m_coeff_offset
Offset of elemental data into the array m_coeffs.
Definition: ExpList.h:1120
ExpansionType m_expType
Exapnsion type.
Definition: ExpList.h:1044
@ 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< HexGeom > HexGeomSharedPtr
Definition: HexGeom.h:87
std::shared_ptr< TetGeom > TetGeomSharedPtr
Definition: TetGeom.h:85
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_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 868 of file DisContField.cpp.

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

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

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

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

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

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

126  {
127  LocTraceToTraceMap = m_locTraceToTraceMap;
128  }

References m_locTraceToTraceMap.

◆ GetNegatedFluxNormal()

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

Definition at line 2774 of file DisContField.cpp.

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

479 {
481  m_traceMap->GetElmtToTrace()[n][e];
482 
483  PeriodicMap periodicTraces = (m_expType == e1D) ? m_periodicVerts
484  : (m_expType == e2D) ? m_periodicEdges
485  : m_periodicFaces;
486 
487  bool fwd = true;
488  if (traceEl->GetLeftAdjacentElementTrace() == -1 ||
489  traceEl->GetRightAdjacentElementTrace() == -1)
490  {
491  // Boundary edge (1 connected element). Do nothing in
492  // serial.
493  auto it = m_boundaryTraces.find(traceEl->GetElmtId());
494 
495  // If the edge does not have a boundary condition set on
496  // it, then assume it is a partition edge or periodic.
497  if (it == m_boundaryTraces.end())
498  {
499  fwd = true;
500  }
501  }
502  else if (traceEl->GetLeftAdjacentElementTrace() != -1 &&
503  traceEl->GetRightAdjacentElementTrace() != -1)
504  {
505  // Non-boundary edge (2 connected elements).
506  fwd = (traceEl->GetLeftAdjacentElementExp().get()) == (*m_exp)[n].get();
507  }
508  else
509  {
510  ASSERTL2(false, "Unconnected trace element!");
511  }
512 
513  return fwd;
514 }
#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 > &  coeffs,
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 4133 of file DisContField.cpp.

4136 {
4137  int i, e, ncoeff_edge;
4138  Array<OneD, const NekDouble> tmp_coeffs;
4139  Array<OneD, NekDouble> out_d(m_ncoeffs), out_tmp;
4140 
4141  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
4142  m_traceMap->GetElmtToTrace();
4143 
4144  StdRegions::Orientation edgedir;
4145 
4146  int cnt;
4147  int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs();
4148  Array<OneD, NekDouble> loc_lambda(LocBndCoeffs), edge_lambda;
4149  m_traceMap->GlobalToLocalBnd(m_trace->GetCoeffs(), loc_lambda);
4150 
4151  edge_lambda = loc_lambda;
4152 
4153  // Calculate Q using standard DG formulation.
4154  for (i = cnt = 0; i < GetExpSize(); ++i)
4155  {
4156  // Probably a better way of setting up lambda than this.
4157  // Note cannot use PutCoeffsInToElmts since lambda space
4158  // is mapped during the solve.
4159  int nEdges = (*m_exp)[i]->GetGeom()->GetNumEdges();
4160  Array<OneD, Array<OneD, NekDouble>> edgeCoeffs(nEdges);
4161 
4162  for (e = 0; e < nEdges; ++e)
4163  {
4164  edgedir = (*m_exp)[i]->GetTraceOrient(e);
4165  ncoeff_edge = elmtToTrace[i][e]->GetNcoeffs();
4166  edgeCoeffs[e] = Array<OneD, NekDouble>(ncoeff_edge);
4167  Vmath::Vcopy(ncoeff_edge, edge_lambda, 1, edgeCoeffs[e], 1);
4168  elmtToTrace[i][e]->SetCoeffsToOrientation(edgedir, edgeCoeffs[e],
4169  edgeCoeffs[e]);
4170  edge_lambda = edge_lambda + ncoeff_edge;
4171  }
4172 
4173  (*m_exp)[i]->DGDeriv(dir, tmp_coeffs = coeffs + m_coeff_offset[i],
4174  elmtToTrace[i], edgeCoeffs, out_tmp = out_d + cnt);
4175  cnt += (*m_exp)[i]->GetNcoeffs();
4176  }
4177 
4178  Array<OneD, NekDouble> phys(m_npoints);
4179  BwdTrans(out_d, phys);
4180  Vmath::Vsub(m_npoints, phys, 1, soln, 1, phys, 1);
4181  return L2(phys);
4182 }
int m_ncoeffs
The total number of local degrees of freedom. m_ncoeffs .
Definition: ExpList.h:1056
NekDouble L2(const Array< OneD, const NekDouble > &inarray, const Array< OneD, const NekDouble > &soln=NullNekDouble1DArray)
This function calculates the error of the global This function calculates the error with respect to...
Definition: ExpList.h:506
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:1677
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_ncoeffs, Nektar::MultiRegions::ExpList::m_npoints, 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 2746 of file DisContField.cpp.

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

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",
const Collections::ImplementationType  ImpType = Collections::eNoImpType 
)
protected

Set up all DG member variables and maps.

Definition at line 163 of file DisContField.cpp.

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

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

3384 {
3385  ASSERTL0(m_expType != e1D, "This method is not setup or "
3386  "tested for 1D expansion");
3387  Array<OneD, NekDouble> Coeffs(m_trace->GetNcoeffs());
3388  m_trace->IProductWRTBase(Fwd, Coeffs);
3389  m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(0, Coeffs, outarray);
3390  m_trace->IProductWRTBase(Bwd, Coeffs);
3391  m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(1, Coeffs, outarray);
3392 }

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 
)
overrideprotectedvirtual

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

3212 {
3213  if (m_expType == e1D)
3214  {
3215  int n, offset, t_offset;
3216  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>>
3217  &elmtToTrace = m_traceMap->GetElmtToTrace();
3218  vector<bool> negatedFluxNormal = GetNegatedFluxNormal();
3219  for (n = 0; n < GetExpSize(); ++n)
3220  {
3221  // Number of coefficients on each element
3222  int e_ncoeffs = (*m_exp)[n]->GetNcoeffs();
3223  offset = GetCoeff_Offset(n);
3224  // Implementation for every points except Gauss points
3225  if ((*m_exp)[n]->GetBasis(0)->GetBasisType() !=
3227  {
3228  t_offset =
3229  GetTrace()->GetCoeff_Offset(elmtToTrace[n][0]->GetElmtId());
3230  if (negatedFluxNormal[2 * n])
3231  {
3232  outarray[offset] -= Fn[t_offset];
3233  }
3234  else
3235  {
3236  outarray[offset] += Fn[t_offset];
3237  }
3238  t_offset =
3239  GetTrace()->GetCoeff_Offset(elmtToTrace[n][1]->GetElmtId());
3240  if (negatedFluxNormal[2 * n + 1])
3241  {
3242  outarray[offset + (*m_exp)[n]->GetVertexMap(1)] -=
3243  Fn[t_offset];
3244  }
3245  else
3246  {
3247  outarray[offset + (*m_exp)[n]->GetVertexMap(1)] +=
3248  Fn[t_offset];
3249  }
3250  }
3251  else
3252  {
3253  int j;
3254  static DNekMatSharedPtr m_Ixm, m_Ixp;
3255  static int sav_ncoeffs = 0;
3256  if (!m_Ixm.get() || e_ncoeffs != sav_ncoeffs)
3257  {
3259  const LibUtilities::PointsKey BS_p(
3261  const LibUtilities::BasisKey BS_k(
3262  LibUtilities::eGauss_Lagrange, e_ncoeffs, BS_p);
3263  BASE = LibUtilities::BasisManager()[BS_k];
3264  Array<OneD, NekDouble> coords(1, 0.0);
3265  coords[0] = -1.0;
3266  m_Ixm = BASE->GetI(coords);
3267  coords[0] = 1.0;
3268  m_Ixp = BASE->GetI(coords);
3269  sav_ncoeffs = e_ncoeffs;
3270  }
3271  t_offset =
3272  GetTrace()->GetCoeff_Offset(elmtToTrace[n][0]->GetElmtId());
3273  if (negatedFluxNormal[2 * n])
3274  {
3275  for (j = 0; j < e_ncoeffs; j++)
3276  {
3277  outarray[offset + j] -=
3278  (m_Ixm->GetPtr())[j] * Fn[t_offset];
3279  }
3280  }
3281  else
3282  {
3283  for (j = 0; j < e_ncoeffs; j++)
3284  {
3285  outarray[offset + j] +=
3286  (m_Ixm->GetPtr())[j] * Fn[t_offset];
3287  }
3288  }
3289  t_offset =
3290  GetTrace()->GetCoeff_Offset(elmtToTrace[n][1]->GetElmtId());
3291  if (negatedFluxNormal[2 * n + 1])
3292  {
3293  for (j = 0; j < e_ncoeffs; j++)
3294  {
3295  outarray[offset + j] -=
3296  (m_Ixp->GetPtr())[j] * Fn[t_offset];
3297  }
3298  }
3299  else
3300  {
3301  for (j = 0; j < e_ncoeffs; j++)
3302  {
3303  outarray[offset + j] +=
3304  (m_Ixp->GetPtr())[j] * Fn[t_offset];
3305  }
3306  }
3307  }
3308  }
3309  }
3310  else // other dimensions
3311  {
3312  // Basis definition on each element
3313  LibUtilities::BasisSharedPtr basis = (*m_exp)[0]->GetBasis(0);
3314  if ((m_expType != e1D) &&
3315  (basis->GetBasisType() != LibUtilities::eGauss_Lagrange))
3316  {
3317  Array<OneD, NekDouble> Fcoeffs(m_trace->GetNcoeffs());
3318  m_trace->IProductWRTBase(Fn, Fcoeffs);
3319  m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(Fcoeffs,
3320  outarray);
3321  }
3322  else
3323  {
3324  int e, n, offset, t_offset;
3325  Array<OneD, NekDouble> e_outarray;
3326  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>>
3327  &elmtToTrace = m_traceMap->GetElmtToTrace();
3328  if (m_expType == e2D)
3329  {
3330  for (n = 0; n < GetExpSize(); ++n)
3331  {
3332  offset = GetCoeff_Offset(n);
3333  for (e = 0; e < (*m_exp)[n]->GetNtraces(); ++e)
3334  {
3335  t_offset = GetTrace()->GetPhys_Offset(
3336  elmtToTrace[n][e]->GetElmtId());
3337  (*m_exp)[n]->AddEdgeNormBoundaryInt(
3338  e, elmtToTrace[n][e], Fn + t_offset,
3339  e_outarray = outarray + offset);
3340  }
3341  }
3342  }
3343  else if (m_expType == e3D)
3344  {
3345  for (n = 0; n < GetExpSize(); ++n)
3346  {
3347  offset = GetCoeff_Offset(n);
3348  for (e = 0; e < (*m_exp)[n]->GetNtraces(); ++e)
3349  {
3350  t_offset = GetTrace()->GetPhys_Offset(
3351  elmtToTrace[n][e]->GetElmtId());
3352  (*m_exp)[n]->AddFaceNormBoundaryInt(
3353  e, elmtToTrace[n][e], Fn + t_offset,
3354  e_outarray = outarray + offset);
3355  }
3356  }
3357  }
3358  }
3359  }
3360 }
std::vector< bool > & GetNegatedFluxNormal(void)
std::shared_ptr< ExpList > & GetTrace()
Definition: ExpList.h:2112
int GetCoeff_Offset(int n) const
Get the start offset position for a local contiguous list of coeffs correspoinding to element n.
Definition: ExpList.h:2037
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 
)
overrideprotectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 4414 of file DisContField.cpp.

4418 {
4419  Array<OneD, NekDouble> FCoeffs(m_trace->GetNcoeffs());
4420 
4421  m_trace->IProductWRTBase(FwdFlux, FCoeffs);
4422  m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(1, FCoeffs, outarray);
4423  m_trace->IProductWRTBase(BwdFlux, FCoeffs);
4424  m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(0, FCoeffs, outarray);
4425 }

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 
)
overrideprotectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 3108 of file DisContField.cpp.

3111 {
3112  // Basis definition on each element
3113  LibUtilities::BasisSharedPtr basis = (*m_exp)[0]->GetBasis(0);
3114  if (basis->GetBasisType() != LibUtilities::eGauss_Lagrange)
3115  {
3116  Array<OneD, NekDouble> tracevals(
3117  m_locTraceToTraceMap->GetNLocTracePts(), 0.0);
3118 
3119  Array<OneD, NekDouble> invals =
3120  tracevals + m_locTraceToTraceMap->GetNFwdLocTracePts();
3121 
3122  m_locTraceToTraceMap->InterpTraceToLocTrace(1, Bwd, invals);
3123 
3124  m_locTraceToTraceMap->InterpTraceToLocTrace(0, Fwd, tracevals);
3125 
3126  m_locTraceToTraceMap->AddLocTracesToField(tracevals, field);
3127  }
3128  else
3129  {
3131  "v_AddTraceQuadPhysToField not coded for eGauss_Lagrange");
3132  }
3133 }

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 
)
overrideprotectedvirtual

Evaluate all boundary conditions at a given time..

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 3562 of file DisContField.cpp.

3566 {
3567  int i;
3568  int npoints;
3569 
3570  MultiRegions::ExpListSharedPtr locExpList;
3571 
3572  for (i = 0; i < m_bndCondExpansions.size(); ++i)
3573  {
3574  if (time == 0.0 || m_bndConditions[i]->IsTimeDependent())
3575  {
3576  m_bndCondBndWeight[i] = 1.0;
3577  locExpList = m_bndCondExpansions[i];
3578 
3579  npoints = locExpList->GetNpoints();
3580  Array<OneD, NekDouble> x0(npoints, 0.0);
3581  Array<OneD, NekDouble> x1(npoints, 0.0);
3582  Array<OneD, NekDouble> x2(npoints, 0.0);
3583 
3584  locExpList->GetCoords(x0, x1, x2);
3585 
3586  if (x2_in != NekConstants::kNekUnsetDouble &&
3588  {
3589  Vmath::Fill(npoints, x2_in, x1, 1);
3590  Vmath::Fill(npoints, x3_in, x2, 1);
3591  }
3592  else if (x2_in != NekConstants::kNekUnsetDouble)
3593  {
3594  Vmath::Fill(npoints, x2_in, x2, 1);
3595  }
3596 
3597  // treat 1D expansions separately since we only
3598  // require an evaluation at a point rather than
3599  // any projections or inner products that are not
3600  // available in a PointExp
3601  if (m_expType == e1D)
3602  {
3603  if (m_bndConditions[i]->GetBoundaryConditionType() ==
3605  {
3606 
3607  m_bndCondExpansions[i]->SetCoeff(
3608  0, (std::static_pointer_cast<
3609  SpatialDomains::DirichletBoundaryCondition>(
3610  m_bndConditions[i])
3611  ->m_dirichletCondition)
3612  .Evaluate(x0[0], x1[0], x2[0], time));
3613  m_bndCondExpansions[i]->SetPhys(
3614  0, m_bndCondExpansions[i]->GetCoeff(0));
3615  }
3616  else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3618  {
3619  m_bndCondExpansions[i]->SetCoeff(
3620  0, (std::static_pointer_cast<
3621  SpatialDomains::NeumannBoundaryCondition>(
3622  m_bndConditions[i])
3623  ->m_neumannCondition)
3624  .Evaluate(x0[0], x1[0], x2[0], time));
3625  }
3626  else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3628  {
3629  m_bndCondExpansions[i]->SetCoeff(
3630  0, (std::static_pointer_cast<
3631  SpatialDomains::RobinBoundaryCondition>(
3632  m_bndConditions[i])
3633  ->m_robinFunction)
3634  .Evaluate(x0[0], x1[0], x2[0], time));
3635  }
3636  else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3638  {
3639  continue;
3640  }
3641  else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3643  {
3644  }
3645  else
3646  {
3648  "This type of BC not implemented yet");
3649  }
3650  }
3651  else // 2D and 3D versions
3652  {
3653  if (m_bndConditions[i]->GetBoundaryConditionType() ==
3655  {
3657  std::static_pointer_cast<
3658  SpatialDomains::DirichletBoundaryCondition>(
3659  m_bndConditions[i]);
3660 
3661  Array<OneD, NekDouble> valuesFile(npoints, 1.0),
3662  valuesExp(npoints, 1.0);
3663 
3664  string bcfilename = bcPtr->m_filename;
3665  string exprbcs = bcPtr->m_expr;
3666 
3667  if (bcfilename != "")
3668  {
3669  locExpList->ExtractCoeffsFromFile(
3670  bcfilename, bcPtr->GetComm(), varName,
3671  locExpList->UpdateCoeffs());
3672  locExpList->BwdTrans(locExpList->GetCoeffs(),
3673  locExpList->UpdatePhys());
3674  valuesFile = locExpList->GetPhys();
3675  }
3676 
3677  if (exprbcs != "")
3678  {
3679  LibUtilities::Equation condition =
3680  std::static_pointer_cast<
3681  SpatialDomains::DirichletBoundaryCondition>(
3682  m_bndConditions[i])
3683  ->m_dirichletCondition;
3684 
3685  condition.Evaluate(x0, x1, x2, time, valuesExp);
3686  }
3687 
3688  Vmath::Vmul(npoints, valuesExp, 1, valuesFile, 1,
3689  locExpList->UpdatePhys(), 1);
3690  locExpList->FwdTransBndConstrained(
3691  locExpList->GetPhys(), locExpList->UpdateCoeffs());
3692  }
3693  else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3695  {
3697  std::static_pointer_cast<
3698  SpatialDomains::NeumannBoundaryCondition>(
3699  m_bndConditions[i]);
3700  string bcfilename = bcPtr->m_filename;
3701 
3702  if (bcfilename != "")
3703  {
3704  locExpList->ExtractCoeffsFromFile(
3705  bcfilename, bcPtr->GetComm(), varName,
3706  locExpList->UpdateCoeffs());
3707  locExpList->BwdTrans(locExpList->GetCoeffs(),
3708  locExpList->UpdatePhys());
3709  }
3710  else
3711  {
3712  LibUtilities::Equation condition =
3713  std::static_pointer_cast<
3714  SpatialDomains::NeumannBoundaryCondition>(
3715  m_bndConditions[i])
3716  ->m_neumannCondition;
3717  condition.Evaluate(x0, x1, x2, time,
3718  locExpList->UpdatePhys());
3719  }
3720 
3721  locExpList->IProductWRTBase(locExpList->GetPhys(),
3722  locExpList->UpdateCoeffs());
3723  }
3724  else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3726  {
3728  std::static_pointer_cast<
3729  SpatialDomains::RobinBoundaryCondition>(
3730  m_bndConditions[i]);
3731 
3732  string bcfilename = bcPtr->m_filename;
3733 
3734  if (bcfilename != "")
3735  {
3736  locExpList->ExtractCoeffsFromFile(
3737  bcfilename, bcPtr->GetComm(), varName,
3738  locExpList->UpdateCoeffs());
3739  locExpList->BwdTrans(locExpList->GetCoeffs(),
3740  locExpList->UpdatePhys());
3741  }
3742  else
3743  {
3744  LibUtilities::Equation condition =
3745  std::static_pointer_cast<
3746  SpatialDomains::RobinBoundaryCondition>(
3747  m_bndConditions[i])
3748  ->m_robinFunction;
3749  condition.Evaluate(x0, x1, x2, time,
3750  locExpList->UpdatePhys());
3751  }
3752 
3753  locExpList->IProductWRTBase(locExpList->GetPhys(),
3754  locExpList->UpdateCoeffs());
3755  }
3756  else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3758  {
3759  continue;
3760  }
3761  else
3762  {
3764  "This type of BC not implemented yet");
3765  }
3766  }
3767  }
3768  }
3769 }
NekDouble GetCoeff(int i)
Get the i th value (coefficient) of m_coeffs.
Definition: ExpList.h:1968
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(), 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)
overrideprotectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 3135 of file DisContField.cpp.

3136 {
3137  ASSERTL1(m_physState == true, "local physical space is not true ");
3138  v_ExtractTracePhys(m_phys, outarray);
3139 }
virtual void v_ExtractTracePhys(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
This method extracts the trace (verts in 1D) from the field inarray and puts the values in outarray.
bool m_physState
The state of the array m_phys.
Definition: ExpList.h:1100
Array< OneD, NekDouble > m_phys
The global expansion evaluated at the quadrature points.
Definition: ExpList.h:1092

References ASSERTL1, Nektar::MultiRegions::ExpList::m_phys, Nektar::MultiRegions::ExpList::m_physState, and v_ExtractTracePhys().

◆ v_ExtractTracePhys() [2/2]

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

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

3157 {
3158  LibUtilities::BasisSharedPtr basis = (*m_exp)[0]->GetBasis(0);
3159  if ((basis->GetBasisType() != LibUtilities::eGauss_Lagrange))
3160  {
3161  Vmath::Zero(outarray.size(), outarray, 1);
3162  Array<OneD, NekDouble> tracevals(
3163  m_locTraceToTraceMap->GetNFwdLocTracePts());
3164  m_locTraceToTraceMap->FwdLocTracesFromField(inarray, tracevals);
3165  m_locTraceToTraceMap->InterpLocTracesToTrace(0, tracevals, outarray);
3166  m_traceMap->GetAssemblyCommDG()->PerformExchange(outarray, outarray);
3167  }
3168  else
3169  {
3170  // Loop over elemente and collect forward expansion
3171  int nexp = GetExpSize();
3172  int n, p, offset, phys_offset;
3173  Array<OneD, NekDouble> t_tmp;
3174  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>>
3175  &elmtToTrace = m_traceMap->GetElmtToTrace();
3176  ASSERTL1(outarray.size() >= m_trace->GetNpoints(),
3177  "input array is of insufficient length");
3178  for (n = 0; n < nexp; ++n)
3179  {
3180  phys_offset = GetPhys_Offset(n);
3181  for (p = 0; p < (*m_exp)[n]->GetNtraces(); ++p)
3182  {
3183  offset =
3184  m_trace->GetPhys_Offset(elmtToTrace[n][p]->GetElmtId());
3185  (*m_exp)[n]->GetTracePhysVals(p, elmtToTrace[n][p],
3186  inarray + phys_offset,
3187  t_tmp = outarray + offset);
3188  }
3189  }
3190  }
3191 }
int GetPhys_Offset(int n) const
Get the start offset position for a local contiguous list of quadrature points in a full array corres...
Definition: ExpList.h:2044
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().

Referenced by v_ExtractTracePhys().

◆ v_FillBwdWithBoundCond()

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

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 2988 of file DisContField.cpp.

2991 {
2992  // Fill boundary conditions into missing elements
2993  if (PutFwdInBwdOnBCs) // just set Bwd value to be Fwd value on BCs
2994  {
2995  // Fill boundary conditions into missing elements
2996  for (int n = 0, cnt = 0; n < m_bndCondExpansions.size(); ++n)
2997  {
2998  if (m_bndConditions[n]->GetBoundaryConditionType() ==
3000  {
3001  auto ne = m_bndCondExpansions[n]->GetExpSize();
3002  for (int e = 0; e < ne; ++e)
3003  {
3004  auto npts =
3005  m_bndCondExpansions[n]->GetExp(e)->GetTotPoints();
3006  auto id2 = m_trace->GetPhys_Offset(
3007  m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
3008  Vmath::Vcopy(npts, &Fwd[id2], 1, &Bwd[id2], 1);
3009  }
3010 
3011  cnt += ne;
3012  }
3013  else if (m_bndConditions[n]->GetBoundaryConditionType() ==
3015  m_bndConditions[n]->GetBoundaryConditionType() ==
3017  {
3018  auto ne = m_bndCondExpansions[n]->GetExpSize();
3019  for (int e = 0; e < ne; ++e)
3020  {
3021  auto npts =
3022  m_bndCondExpansions[n]->GetExp(e)->GetTotPoints();
3023  auto id2 = m_trace->GetPhys_Offset(
3024  m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
3025  Vmath::Vcopy(npts, &Fwd[id2], 1, &Bwd[id2], 1);
3026  }
3027  cnt += ne;
3028  }
3029  else if (m_bndConditions[n]->GetBoundaryConditionType() !=
3031  {
3032  ASSERTL0(false, "Method not set up for this "
3033  "boundary condition.");
3034  }
3035  }
3036  }
3037  else
3038  {
3039  for (int n = 0, cnt = 0; n < m_bndCondExpansions.size(); ++n)
3040  {
3041  if (m_bndConditions[n]->GetBoundaryConditionType() ==
3043  {
3044  auto ne = m_bndCondExpansions[n]->GetExpSize();
3045  for (int e = 0; e < ne; ++e)
3046  {
3047  auto npts =
3048  m_bndCondExpansions[n]->GetExp(e)->GetTotPoints();
3049  auto id1 = m_bndCondExpansions[n]->GetPhys_Offset(e);
3050  auto id2 = m_trace->GetPhys_Offset(
3051  m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
3052 
3053  Vmath::Vcopy(npts,
3054  &(m_bndCondExpansions[n]->GetPhys())[id1], 1,
3055  &Bwd[id2], 1);
3056  }
3057  cnt += ne;
3058  }
3059  else if (m_bndConditions[n]->GetBoundaryConditionType() ==
3061  m_bndConditions[n]->GetBoundaryConditionType() ==
3063  {
3064  auto ne = m_bndCondExpansions[n]->GetExpSize();
3065  for (int e = 0; e < ne; ++e)
3066  {
3067  auto npts =
3068  m_bndCondExpansions[n]->GetExp(e)->GetTotPoints();
3069  auto id1 = m_bndCondExpansions[n]->GetPhys_Offset(e);
3070 
3071  ASSERTL0((m_bndCondExpansions[n]->GetPhys())[id1] == 0.0,
3072  "Method not set up for non-zero "
3073  "Neumann boundary condition");
3074  auto id2 = m_trace->GetPhys_Offset(
3075  m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
3076 
3077  Vmath::Vcopy(npts, &Fwd[id2], 1, &Bwd[id2], 1);
3078  }
3079 
3080  cnt += ne;
3081  }
3082  else if (m_bndConditions[n]->GetBoundaryConditionType() ==
3084  {
3085  // Do nothing
3086  }
3087  else if (m_bndConditions[n]->GetBoundaryConditionType() !=
3089  {
3091  "Method not set up for this boundary "
3092  "condition.");
3093  }
3094  }
3095  }
3096 }
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:1988

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 
)
overrideprotectedvirtual

Fill the weight with m_bndCondBndWeight.

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 3774 of file DisContField.cpp.

3776 {
3777  int cnt;
3778  int npts;
3779  int e = 0;
3780 
3781  // Fill boundary conditions into missing elements
3782  int id2 = 0;
3783 
3784  for (int n = cnt = 0; n < m_bndCondExpansions.size(); ++n)
3785  {
3786 
3787  if (m_bndConditions[n]->GetBoundaryConditionType() ==
3789  {
3790  for (e = 0; e < m_bndCondExpansions[n]->GetExpSize(); ++e)
3791  {
3792  npts = m_bndCondExpansions[n]->GetExp(e)->GetTotPoints();
3793  id2 = m_trace->GetPhys_Offset(
3794  m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
3795  Vmath::Fill(npts, m_bndCondBndWeight[n], &weightave[id2], 1);
3796  Vmath::Fill(npts, 0.0, &weightjmp[id2], 1);
3797  }
3798 
3799  cnt += e;
3800  }
3801  else if (m_bndConditions[n]->GetBoundaryConditionType() ==
3803  m_bndConditions[n]->GetBoundaryConditionType() ==
3805  {
3806  for (e = 0; e < m_bndCondExpansions[n]->GetExpSize(); ++e)
3807  {
3808  npts = m_bndCondExpansions[n]->GetExp(e)->GetTotPoints();
3809  id2 = m_trace->GetPhys_Offset(
3810  m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
3811  Vmath::Fill(npts, m_bndCondBndWeight[n], &weightave[id2], 1);
3812  Vmath::Fill(npts, 0.0, &weightjmp[id2], 1);
3813  }
3814 
3815  cnt += e;
3816  }
3817  else if (m_bndConditions[n]->GetBoundaryConditionType() !=
3819  {
3821  "Method not set up for this boundary condition.");
3822  }
3823  }
3824 }

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 ( )
overrideprotectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 3098 of file DisContField.cpp.

3099 {
3100  return m_bndCondBndWeight;
3101 }

References m_bndCondBndWeight.

◆ v_GetBndCondExpansions()

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

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 2839 of file DisContField.cpp.

2840 {
2841  return m_bndCondExpansions;
2842 }

References m_bndCondExpansions.

◆ v_GetBndConditions()

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

Reimplemented from Nektar::MultiRegions::ExpList.

Reimplemented in Nektar::MultiRegions::ContField.

Definition at line 2845 of file DisContField.cpp.

2846 {
2847  return m_bndConditions;
2848 }

References m_bndConditions.

◆ v_GetBndElmtExpansion()

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

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 3991 of file DisContField.cpp.

3994 {
3995  int n, cnt;
3996  std::vector<unsigned int> eIDs;
3997 
3998  Array<OneD, int> ElmtID, TraceID;
3999  GetBoundaryToElmtMap(ElmtID, TraceID);
4000 
4001  // Skip other boundary regions
4002  for (cnt = n = 0; n < i; ++n)
4003  {
4004  cnt += m_bndCondExpansions[n]->GetExpSize();
4005  }
4006 
4007  // Populate eIDs with information from BoundaryToElmtMap
4008  for (n = 0; n < m_bndCondExpansions[i]->GetExpSize(); ++n)
4009  {
4010  eIDs.push_back(ElmtID[cnt + n]);
4011  }
4012 
4013  // Create expansion list
4015  *this, eIDs, DeclareCoeffPhysArrays, Collections::eNoCollection);
4016  // Copy phys and coeffs to new explist
4017  if (DeclareCoeffPhysArrays)
4018  {
4019  int nq;
4020  int offsetOld, offsetNew;
4021  Array<OneD, NekDouble> tmp1, tmp2;
4022  for (n = 0; n < result->GetExpSize(); ++n)
4023  {
4024  nq = GetExp(ElmtID[cnt + n])->GetTotPoints();
4025  offsetOld = GetPhys_Offset(ElmtID[cnt + n]);
4026  offsetNew = result->GetPhys_Offset(n);
4027  Vmath::Vcopy(nq, tmp1 = GetPhys() + offsetOld, 1,
4028  tmp2 = result->UpdatePhys() + offsetNew, 1);
4029  nq = GetExp(ElmtID[cnt + n])->GetNcoeffs();
4030  offsetOld = GetCoeff_Offset(ElmtID[cnt + n]);
4031  offsetNew = result->GetCoeff_Offset(n);
4032  Vmath::Vcopy(nq, tmp1 = GetCoeffs() + offsetOld, 1,
4033  tmp2 = result->UpdateCoeffs() + offsetNew, 1);
4034  }
4035  }
4036 }
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:1902
const std::shared_ptr< LocalRegions::ExpansionVector > GetExp() const
This function returns the vector of elements in the expansion.
Definition: ExpList.h:2029

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 
)
overrideprotectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 3828 of file DisContField.cpp.

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

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 2870 of file DisContField.cpp.

2872 {
2873  v_GetFwdBwdTracePhys(m_phys, Fwd, Bwd);
2874 }
virtual void v_GetFwdBwdTracePhys(Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd) override

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 
)
overrideprotectedvirtual

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

2919 {
2920  // Is this zeroing necessary?
2921  // Zero forward/backward vectors.
2922  Vmath::Zero(Fwd.size(), Fwd, 1);
2923  Vmath::Zero(Bwd.size(), Bwd, 1);
2924 
2925  // Basis definition on each element
2926  LibUtilities::BasisSharedPtr basis = (*m_exp)[0]->GetBasis(0);
2927  if ((basis->GetBasisType() != LibUtilities::eGauss_Lagrange))
2928  {
2929  // blocked routine
2930  Array<OneD, NekDouble> tracevals(
2931  m_locTraceToTraceMap->GetNLocTracePts());
2932 
2933  m_locTraceToTraceMap->LocTracesFromField(field, tracevals);
2934  m_locTraceToTraceMap->InterpLocTracesToTrace(0, tracevals, Fwd);
2935 
2936  Array<OneD, NekDouble> invals =
2937  tracevals + m_locTraceToTraceMap->GetNFwdLocTracePts();
2938  m_locTraceToTraceMap->InterpLocTracesToTrace(1, invals, Bwd);
2939  }
2940  else
2941  {
2942  // Loop over elements and collect forward expansion
2943  auto nexp = GetExpSize();
2944  Array<OneD, NekDouble> e_tmp;
2946 
2947  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>>
2948  &elmtToTrace = m_traceMap->GetElmtToTrace();
2949 
2950  for (int n = 0, cnt = 0; n < nexp; ++n)
2951  {
2952  exp = (*m_exp)[n];
2953  auto phys_offset = GetPhys_Offset(n);
2954 
2955  for (int e = 0; e < exp->GetNtraces(); ++e, ++cnt)
2956  {
2957  auto offset =
2958  m_trace->GetPhys_Offset(elmtToTrace[n][e]->GetElmtId());
2959 
2960  e_tmp =
2961  (m_leftAdjacentTraces[cnt]) ? Fwd + offset : Bwd + offset;
2962 
2963  exp->GetTracePhysVals(e, elmtToTrace[n][e], field + phys_offset,
2964  e_tmp);
2965  }
2966  }
2967  }
2968 
2970 
2971  if (FillBnd)
2972  {
2973  v_FillBwdWithBoundCond(Fwd, Bwd, PutFwdInBwdOnBCs);
2974  }
2975 
2976  if (DoExchange)
2977  {
2978  // Do parallel exchange for forwards/backwards spaces.
2979  m_traceMap->GetAssemblyCommDG()->PerformExchange(Fwd, Bwd);
2980 
2981  // Do exchange of interface traces (local and parallel)
2982  // We may have to split this out into separate local and
2983  // parallel for IP method???
2984  m_interfaceMap->ExchangeTrace(Fwd, Bwd);
2985  }
2986 }
virtual void v_PeriodicBwdCopy(const Array< OneD, const NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd) override
virtual void v_FillBwdWithBoundCond(const Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd, bool PutFwdInBwdOnBCs) override

References Nektar::LibUtilities::eGauss_Lagrange, Nektar::MultiRegions::ExpList::GetExpSize(), Nektar::MultiRegions::ExpList::GetPhys_Offset(), m_interfaceMap, 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  )
overrideprotectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 2833 of file DisContField.cpp.

2834 {
2835  return m_leftAdjacentTraces;
2836 }

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 
)
overrideprotectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 4396 of file DisContField.cpp.

4400 {
4401  if (locTraceFwd != NullNekDouble1DArray)
4402  {
4403  m_locTraceToTraceMap->InterpLocTracesToTraceTranspose(0, Fwd,
4404  locTraceFwd);
4405  }
4406 
4407  if (locTraceBwd != NullNekDouble1DArray)
4408  {
4409  m_locTraceToTraceMap->InterpLocTracesToTraceTranspose(1, Bwd,
4410  locTraceBwd);
4411  }
4412 }
static Array< OneD, NekDouble > NullNekDouble1DArray

References m_locTraceToTraceMap, and Nektar::NullNekDouble1DArray.

◆ v_GetLocTraceToTraceMap()

const LocTraceToTraceMapSharedPtr & Nektar::MultiRegions::DisContField::v_GetLocTraceToTraceMap ( void  ) const
overrideprotectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 2824 of file DisContField.cpp.

2826 {
2827  return m_locTraceToTraceMap;
2828 }

References m_locTraceToTraceMap.

◆ v_GetPeriodicEntities()

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

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

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 2880 of file DisContField.cpp.

2883 {
2884  periodicVerts = m_periodicVerts;
2885  periodicEdges = m_periodicEdges;
2886  periodicFaces = m_periodicFaces;
2887 }

References m_periodicEdges, m_periodicFaces, and m_periodicVerts.

◆ v_GetRobinBCInfo()

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

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

4067 {
4068  int i, cnt;
4069  map<int, RobinBCInfoSharedPtr> returnval;
4070  Array<OneD, int> ElmtID, TraceID;
4071  GetBoundaryToElmtMap(ElmtID, TraceID);
4072 
4073  for (cnt = i = 0; i < m_bndCondExpansions.size(); ++i)
4074  {
4075  MultiRegions::ExpListSharedPtr locExpList;
4076 
4077  if (m_bndConditions[i]->GetBoundaryConditionType() ==
4079  {
4080  int e, elmtid;
4081  Array<OneD, NekDouble> Array_tmp;
4082 
4083  locExpList = m_bndCondExpansions[i];
4084 
4085  int npoints = locExpList->GetNpoints();
4086  Array<OneD, NekDouble> x0(npoints, 0.0);
4087  Array<OneD, NekDouble> x1(npoints, 0.0);
4088  Array<OneD, NekDouble> x2(npoints, 0.0);
4089  Array<OneD, NekDouble> coeffphys(npoints);
4090 
4091  locExpList->GetCoords(x0, x1, x2);
4092 
4093  LibUtilities::Equation coeffeqn =
4094  std::static_pointer_cast<
4095  SpatialDomains::RobinBoundaryCondition>(m_bndConditions[i])
4096  ->m_robinPrimitiveCoeff;
4097 
4098  // evalaute coefficient
4099  coeffeqn.Evaluate(x0, x1, x2, 0.0, coeffphys);
4100 
4101  for (e = 0; e < locExpList->GetExpSize(); ++e)
4102  {
4103  RobinBCInfoSharedPtr rInfo =
4105  TraceID[cnt + e],
4106  Array_tmp = coeffphys + locExpList->GetPhys_Offset(e));
4107 
4108  elmtid = ElmtID[cnt + e];
4109  // make link list if necessary
4110  if (returnval.count(elmtid) != 0)
4111  {
4112  rInfo->next = returnval.find(elmtid)->second;
4113  }
4114  returnval[elmtid] = rInfo;
4115  }
4116  }
4117  cnt += m_bndCondExpansions[i]->GetExpSize();
4118  }
4119 
4120  return returnval;
4121 }
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()

ExpListSharedPtr & Nektar::MultiRegions::DisContField::v_GetTrace ( )
overrideprotectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 2809 of file DisContField.cpp.

2810 {
2812  {
2813  SetUpDG();
2814  }
2815 
2816  return m_trace;
2817 }

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

◆ v_GetTraceMap()

AssemblyMapDGSharedPtr & Nektar::MultiRegions::DisContField::v_GetTraceMap ( void  )
overrideprotectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 2819 of file DisContField.cpp.

2820 {
2821  return m_traceMap;
2822 }

References m_traceMap.

◆ v_HelmSolve()

GlobalLinSysKey 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 
)
overrideprotectedvirtual

Solve the Helmholtz equation.

Reimplemented from Nektar::MultiRegions::ExpList.

Reimplemented in Nektar::MultiRegions::ContField.

Definition at line 3394 of file DisContField.cpp.

3400 {
3401  boost::ignore_unused(varfactors, dirForcing);
3402  int i, n, cnt, nbndry;
3403  int nexp = GetExpSize();
3404  Array<OneD, NekDouble> f(m_ncoeffs);
3405  DNekVec F(m_ncoeffs, f, eWrapper);
3406  Array<OneD, NekDouble> e_f, e_l;
3407  //----------------------------------
3408  // Setup RHS Inner product if required
3409  //----------------------------------
3410  if (PhysSpaceForcing)
3411  {
3412  IProductWRTBase(inarray, f);
3413  Vmath::Neg(m_ncoeffs, f, 1);
3414  }
3415  else
3416  {
3417  Vmath::Smul(m_ncoeffs, -1.0, inarray, 1, f, 1);
3418  }
3419  //----------------------------------
3420  // Solve continuous Boundary System
3421  //----------------------------------
3422  int GloBndDofs = m_traceMap->GetNumGlobalBndCoeffs();
3423  int NumDirBCs = m_traceMap->GetNumLocalDirBndCoeffs();
3424  int e_ncoeffs;
3425  GlobalMatrixKey HDGLamToUKey(StdRegions::eHybridDGLamToU,
3426  NullAssemblyMapSharedPtr, factors, varcoeff);
3427  const DNekScalBlkMatSharedPtr &HDGLamToU = GetBlockMatrix(HDGLamToUKey);
3428  // Retrieve number of local trace space coefficients N_{\lambda},
3429  // and set up local elemental trace solution \lambda^e.
3430  int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs();
3431  Array<OneD, NekDouble> bndrhs(LocBndCoeffs, 0.0);
3432  Array<OneD, NekDouble> loclambda(LocBndCoeffs, 0.0);
3433  DNekVec LocLambda(LocBndCoeffs, loclambda, eWrapper);
3434  //----------------------------------
3435  // Evaluate Trace Forcing
3436  // Kirby et al, 2010, P23, Step 5.
3437  //----------------------------------
3438  // Determing <u_lam,f> terms using HDGLamToU matrix
3439  for (cnt = n = 0; n < nexp; ++n)
3440  {
3441  nbndry = (*m_exp)[n]->NumDGBndryCoeffs();
3442  e_ncoeffs = (*m_exp)[n]->GetNcoeffs();
3443  e_f = f + m_coeff_offset[n];
3444  e_l = bndrhs + cnt;
3445  // use outarray as tmp space
3446  DNekVec Floc(nbndry, e_l, eWrapper);
3447  DNekVec ElmtFce(e_ncoeffs, e_f, eWrapper);
3448  Floc = Transpose(*(HDGLamToU->GetBlock(n, n))) * ElmtFce;
3449  cnt += nbndry;
3450  }
3451  Array<OneD, const int> bndCondMap =
3452  m_traceMap->GetBndCondCoeffsToLocalTraceMap();
3453  Array<OneD, const NekDouble> Sign = m_traceMap->GetLocalToGlobalBndSign();
3454  // Copy Dirichlet boundary conditions and weak forcing
3455  // into trace space
3456  int locid;
3457  cnt = 0;
3458  for (i = 0; i < m_bndCondExpansions.size(); ++i)
3459  {
3460  const Array<OneD, const NekDouble> bndcoeffs =
3461  m_bndCondExpansions[i]->GetCoeffs();
3462 
3463  if (m_bndConditions[i]->GetBoundaryConditionType() ==
3465  {
3466  for (int j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); ++j)
3467  {
3468  locid = bndCondMap[cnt + j];
3469  loclambda[locid] = Sign[locid] * bndcoeffs[j];
3470  }
3471  }
3472  else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3474  m_bndConditions[i]->GetBoundaryConditionType() ==
3476  {
3477  // Add weak boundary condition to trace forcing
3478  for (int j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); ++j)
3479  {
3480  locid = bndCondMap[cnt + j];
3481  bndrhs[locid] += Sign[locid] * bndcoeffs[j];
3482  }
3483  }
3484  else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3486  {
3487  // skip increment of cnt if ePeriodic
3488  // because bndCondMap does not include ePeriodic
3489  continue;
3490  }
3491  cnt += (m_bndCondExpansions[i])->GetNcoeffs();
3492  }
3493 
3494  //----------------------------------
3495  // Solve trace problem: \Lambda = K^{-1} F
3496  // K is the HybridDGHelmBndLam matrix.
3497  //----------------------------------
3498  if (GloBndDofs - NumDirBCs > 0)
3499  {
3500  GlobalLinSysKey key(StdRegions::eHybridDGHelmBndLam, m_traceMap,
3501  factors, varcoeff);
3503  LinSys->Solve(bndrhs, loclambda, m_traceMap);
3504 
3505  // For consistency with previous version put global
3506  // solution into m_trace->m_coeffs
3507  m_traceMap->LocalToGlobal(loclambda, m_trace->UpdateCoeffs());
3508  }
3509 
3510  //----------------------------------
3511  // Internal element solves
3512  //----------------------------------
3513  GlobalMatrixKey invHDGhelmkey(StdRegions::eInvHybridDGHelmholtz,
3514  NullAssemblyMapSharedPtr, factors, varcoeff);
3515 
3516  const DNekScalBlkMatSharedPtr &InvHDGHelm = GetBlockMatrix(invHDGhelmkey);
3517  DNekVec out(m_ncoeffs, outarray, eWrapper);
3518  Vmath::Zero(m_ncoeffs, outarray, 1);
3519 
3520  // out = u_f + u_lam = (*InvHDGHelm)*f + (LamtoU)*Lam
3521  out = (*InvHDGHelm) * F + (*HDGLamToU) * LocLambda;
3522 
3523  // Return empty GlobalLinSysKey
3524  return NullGlobalLinSysKey;
3525 }
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:1630
int GetNcoeffs(void) const
Returns the total number of local degrees of freedom .
Definition: ExpList.h:1492
const DNekScalBlkMatSharedPtr & GetBlockMatrix(const GlobalMatrixKey &gkey)
Definition: ExpList.cpp:2016
static AssemblyMapSharedPtr NullAssemblyMapSharedPtr
Definition: AssemblyMap.h:53
static GlobalLinSysKey NullGlobalLinSysKey(StdRegions::eNoMatrixType)
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::ePeriodic, 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, Nektar::MultiRegions::NullGlobalLinSysKey(), 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 
)
overrideprotectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 2861 of file DisContField.cpp.

2863 {
2864  for (int n = 0; n < m_periodicFwdCopy.size(); ++n)
2865  {
2866  Bwd[m_periodicBwdCopy[n]] = Fwd[m_periodicFwdCopy[n]];
2867  }
2868 }

References m_periodicBwdCopy, and m_periodicFwdCopy.

Referenced by v_GetFwdBwdTracePhys().

◆ v_Reset()

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

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

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 4041 of file DisContField.cpp.

4042 {
4043  ExpList::v_Reset();
4044 
4045  // Reset boundary condition expansions.
4046  for (int n = 0; n < m_bndCondExpansions.size(); ++n)
4047  {
4048  m_bndCondExpansions[n]->Reset();
4049  }
4050 }
virtual void v_Reset()
Reset geometry information, metrics, matrix managers and geometry information.
Definition: ExpList.cpp:2673

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

◆ v_SetBndCondBwdWeight()

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

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 3103 of file DisContField.cpp.

3104 {
3105  m_bndCondBndWeight[index] = value;
3106 }

References m_bndCondBndWeight.

◆ v_UpdateBndCondExpansion()

MultiRegions::ExpListSharedPtr & Nektar::MultiRegions::DisContField::v_UpdateBndCondExpansion ( int  i)
overrideprotectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 2850 of file DisContField.cpp.

2851 {
2852  return m_bndCondExpansions[i];
2853 }

References m_bndCondExpansions.

◆ v_UpdateBndConditions()

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

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 2856 of file DisContField.cpp.

2857 {
2858  return m_bndConditions;
2859 }

References m_bndConditions.

Member Data Documentation

◆ m_BCtoElmMap

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

Definition at line 60 of file DisContField.h.

Referenced by v_GetBoundaryToElmtMap().

◆ m_BCtoTraceMap

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

Definition at line 61 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 169 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 158 of file DisContField.h.

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

◆ m_interfaceMap

InterfaceMapDGSharedPtr Nektar::MultiRegions::DisContField::m_interfaceMap
protected

Interfaces mapping for trace space.

Definition at line 155 of file DisContField.h.

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

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

Referenced by GetNegatedFluxNormal().

◆ m_numDirBndCondExpansions

size_t Nektar::MultiRegions::DisContField::m_numDirBndCondExpansions
protected

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

Definition at line 133 of file DisContField.h.

◆ m_periodicBwdCopy

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

Definition at line 192 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 179 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 184 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 191 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 174 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