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...
 
 ~DisContField () override
 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)
 
void GetFwdBwdTracePhys (const Array< OneD, const NekDouble > &field, Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd, const Array< OneD, const SpatialDomains::BoundaryConditionShPtr > &bndCond, const Array< OneD, const ExpListSharedPtr > &BndCondExp)
 This method extracts the "forward" and "backward" trace data from the array field and puts the data into output vectors Fwd and Bwd. More...
 
- Public Member Functions inherited from Nektar::MultiRegions::ExpList
 ExpList (const ExpansionType Type=eNoType)
 The default constructor using a type. More...
 
 ExpList (const ExpList &in, const bool DeclareCoeffPhysArrays=true)
 The copy constructor. More...
 
 ExpList (const ExpList &in, const std::vector< unsigned int > &eIDs, const bool DeclareCoeffPhysArrays=true, const Collections::ImplementationType ImpType=Collections::eNoImpType)
 Constructor copying only elements defined in eIds. More...
 
 ExpList (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &graph, const bool DeclareCoeffPhysArrays=true, const std::string &var="DefaultVar", const Collections::ImplementationType ImpType=Collections::eNoImpType)
 Generate an ExpList from a meshgraph graph and session file. More...
 
 ExpList (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::ExpansionInfoMap &expansions, const bool DeclareCoeffPhysArrays=true, const Collections::ImplementationType ImpType=Collections::eNoImpType)
 Sets up a list of local expansions based on an expansion Map. More...
 
 ExpList (const SpatialDomains::PointGeomSharedPtr &geom)
 Specialised constructors for 0D Expansions Wrapper around LocalRegion::PointExp - used in PrePacing.cpp. More...
 
 ExpList (const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, const ExpListSharedPtr > &bndConstraint, const Array< OneD, const SpatialDomains ::BoundaryConditionShPtr > &bndCond, const LocalRegions::ExpansionVector &locexp, const SpatialDomains::MeshGraphSharedPtr &graph, const 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, 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 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 ResetMatrices ()
 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)
 
std::shared_ptr< InterfaceMapDG > & GetInterfaceMap (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 > &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...
 
ExpListSharedPtrv_GetTrace () override
 
AssemblyMapDGSharedPtrv_GetTraceMap (void) override
 
InterfaceMapDGSharedPtrv_GetInterfaceMap (void) override
 
const LocTraceToTraceMapSharedPtrv_GetLocTraceToTraceMap (void) const override
 
std::vector< bool > & v_GetLeftAdjacentTraces (void) override
 
void v_AddTraceIntegral (const Array< OneD, const NekDouble > &Fn, Array< OneD, NekDouble > &outarray) override
 Add trace contributions into elemental coefficient spaces. More...
 
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...
 
void v_AddTraceQuadPhysToField (const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &field) override
 
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...
 
void v_ExtractTracePhys (Array< OneD, NekDouble > &outarray) override
 
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)
 
std::map< int, RobinBCInfoSharedPtrv_GetRobinBCInfo () override
 
const Array< OneD, const MultiRegions::ExpListSharedPtr > & v_GetBndCondExpansions () override
 
const Array< OneD, const SpatialDomains::BoundaryConditionShPtr > & v_GetBndConditions () override
 
MultiRegions::ExpListSharedPtrv_UpdateBndCondExpansion (int i) override
 
Array< OneD, SpatialDomains::BoundaryConditionShPtr > & v_UpdateBndConditions () override
 
void v_GetBoundaryToElmtMap (Array< OneD, int > &ElmtID, Array< OneD, int > &TraceID) override
 
void v_GetBndElmtExpansion (int i, std::shared_ptr< ExpList > &result, const bool DeclareCoeffPhysArrays) override
 
void v_Reset () override
 Reset this field, so that geometry information can be updated. More...
 
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...
 
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...
 
void v_PeriodicBwdCopy (const Array< OneD, const NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd) override
 
void v_FillBwdWithBwdWeight (Array< OneD, NekDouble > &weightave, Array< OneD, NekDouble > &weightjmp) override
 Fill the weight with m_bndCondBndWeight. More...
 
void v_GetFwdBwdTracePhys (Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd) override
 
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
 
void v_FillBwdWithBoundCond (const Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd, bool PutFwdInBwdOnBCs) override
 
void FillBwdWithBoundCond (const Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd, const Array< OneD, const SpatialDomains::BoundaryConditionShPtr > &bndConditions, const Array< OneD, const ExpListSharedPtr > &BndCondExpansions, bool PutFwdInBwdOnBCs)
 
const Array< OneD, const NekDouble > & v_GetBndCondBwdWeight () override
 
void v_SetBndCondBwdWeight (const int index, const NekDouble value) override
 
void v_GetPeriodicEntities (PeriodicMap &periodicVerts, PeriodicMap &periodicEdges, PeriodicMap &periodicFaces) override
 Obtain a copy of the periodic edges and vertices for this field. More...
 
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 const Array< OneD, const std::shared_ptr< ExpList > > & v_GetBndCondExpansions (void)
 
virtual const Array< OneD, const NekDouble > & v_GetBndCondBwdWeight ()
 
virtual void v_SetBndCondBwdWeight (const int index, const NekDouble value)
 
virtual std::shared_ptr< ExpList > & v_UpdateBndCondExpansion (int i)
 
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 std::shared_ptr< ExpList > & v_GetTrace ()
 
virtual std::shared_ptr< AssemblyMapDG > & v_GetTraceMap ()
 
virtual std::shared_ptr< InterfaceMapDG > & v_GetInterfaceMap ()
 
virtual const Array< OneD, const int > & v_GetTraceBndMap ()
 
virtual const std::shared_ptr< LocTraceToTraceMap > & v_GetLocTraceToTraceMap (void) const
 
virtual std::vector< bool > & v_GetLeftAdjacentTraces (void)
 
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 > &Fn, Array< OneD, NekDouble > &outarray)
 
virtual void v_AddFwdBwdTraceIntegral (const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &outarray)
 
virtual void v_GetFwdBwdTracePhys (Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd)
 
virtual void v_GetFwdBwdTracePhys (const Array< OneD, const NekDouble > &field, Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd, bool FillBnd=true, bool PutFwdInBwdOnBCs=false, bool DoExchange=true)
 
virtual void v_FillBwdWithBoundCond (const Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd, bool PutFwdInBwdOnBCs)
 
virtual void v_AddTraceQuadPhysToField (const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &field)
 
virtual void v_AddTraceQuadPhysToOffDiag (const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &field)
 
virtual void v_GetLocTraceFromTracePts (const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &locTraceFwd, Array< OneD, NekDouble > &locTraceBwd)
 
virtual void v_FillBwdWithBwdWeight (Array< OneD, NekDouble > &weightave, Array< OneD, NekDouble > &weightjmp)
 
virtual void v_PeriodicBwdCopy (const Array< OneD, const NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd)
 
virtual const std::vector< bool > & v_GetLeftAdjacentFaces (void) const
 
virtual void v_ExtractTracePhys (Array< OneD, NekDouble > &outarray)
 
virtual void v_ExtractTracePhys (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
virtual void v_MultiplyByInvMassMatrix (const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
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)
 
virtual GlobalLinSysKey v_LinearAdvectionDiffusionReactionSolve (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)
 
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_Reset ()
 Reset geometry information, metrics, matrix managers and geometry information. More...
 
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_IProductWRTDerivBase (const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
 
virtual void v_IProductWRTDerivBase (const Array< OneD, const Array< OneD, 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_Curl (Array< OneD, Array< OneD, NekDouble > > &Vel, Array< OneD, Array< OneD, NekDouble > > &Q)
 
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_GetBoundaryToElmtMap (Array< OneD, int > &ElmtID, Array< OneD, int > &EdgeID)
 
virtual void v_GetBndElmtExpansion (int i, std::shared_ptr< ExpList > &result, const bool DeclareCoeffPhysArrays)
 
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 const Array< OneD, const SpatialDomains::BoundaryConditionShPtr > & v_GetBndConditions ()
 
virtual Array< OneD, SpatialDomains::BoundaryConditionShPtr > & v_UpdateBndConditions ()
 
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)
 
virtual std::map< int, RobinBCInfoSharedPtrv_GetRobinBCInfo (void)
 
virtual void v_GetPeriodicEntities (PeriodicMap &periodicVerts, PeriodicMap &periodicEdges, PeriodicMap &periodicFaces)
 
virtual LibUtilities::BasisSharedPtr v_GetHomogeneousBasis (void)
 
virtual void v_SetHomo1DSpecVanVisc (Array< OneD, NekDouble > visc)
 
virtual std::shared_ptr< ExpList > & v_GetPlane (int n)
 
virtual void v_AddTraceIntegralToOffDiag (const Array< OneD, const NekDouble > &FwdFlux, const Array< OneD, const NekDouble > &BwdFlux, Array< OneD, NekDouble > &outarray)
 

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

Constructor & Destructor Documentation

◆ DisContField() [1/6]

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

Default constructor.

Constructs an empty expansion list with no boundary conditions.

Definition at line 62 of file DisContField.cpp.

65{
66}
Array< OneD, SpatialDomains::BoundaryConditionShPtr > m_bndConditions
An array which contains the information about the boundary condition structure definition on the diff...
Definition: DisContField.h:143
ExpListSharedPtr m_trace
Trace space storage for points between elements.
Definition: DisContField.h:167
Array< OneD, MultiRegions::ExpListSharedPtr > m_bndCondExpansions
An object which contains the discretised boundary conditions.
Definition: DisContField.h:156
ExpList(const ExpansionType Type=eNoType)
The default constructor using a type.
Definition: ExpList.cpp:91
static ExpListSharedPtr NullExpListSharedPtr
Definition: ExpList.h:1499

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

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

615 : ExpList(pSession, domain, graph1D, true, variable, SetToOneSpaceDimension,
616 comm, ImpType)
617{
618 if (variable.compare("DefaultVar") != 0)
619 {
621 GetDomainBCs(domain, Allbcs, variable);
622
624 EvaluateBoundaryConditions(0.0, variable);
626 FindPeriodicTraces(*DomBCs, variable);
627 }
628
629 SetUpDG(variable);
630}
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:1060
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 636 of file DisContField.cpp.

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

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

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

784 : ExpList(In)
785{
786}

◆ ~DisContField()

Nektar::MultiRegions::DisContField::~DisContField ( )
override

Destructor.

Definition at line 791 of file DisContField.cpp.

792{
793}

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

4247{
4248 int i, cnt, e, ncoeff_trace;
4249 Array<OneD, NekDouble> force, out_tmp, qrhs, qrhs1;
4250 Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
4251 m_traceMap->GetElmtToTrace();
4252
4254
4255 int nq_elmt, nm_elmt;
4256 int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs();
4257 Array<OneD, NekDouble> loc_lambda(LocBndCoeffs), trace_lambda;
4258 Array<OneD, NekDouble> tmp_coeffs;
4259 m_traceMap->GlobalToLocalBnd(m_trace->GetCoeffs(), loc_lambda);
4260
4261 trace_lambda = loc_lambda;
4262
4263 int dim = (m_expType == e2D) ? 2 : 3;
4264
4265 int num_points[] = {0, 0, 0};
4266 int num_modes[] = {0, 0, 0};
4267
4268 // Calculate Q using standard DG formulation.
4269 for (i = cnt = 0; i < GetExpSize(); ++i)
4270 {
4271 nq_elmt = (*m_exp)[i]->GetTotPoints();
4272 nm_elmt = (*m_exp)[i]->GetNcoeffs();
4273 qrhs = Array<OneD, NekDouble>(nq_elmt);
4274 qrhs1 = Array<OneD, NekDouble>(nq_elmt);
4275 force = Array<OneD, NekDouble>(2 * nm_elmt);
4276 out_tmp = force + nm_elmt;
4278
4279 for (int j = 0; j < dim; ++j)
4280 {
4281 num_points[j] = (*m_exp)[i]->GetBasis(j)->GetNumPoints();
4282 num_modes[j] = (*m_exp)[i]->GetBasis(j)->GetNumModes();
4283 }
4284
4285 // Probably a better way of setting up lambda than this. Note
4286 // cannot use PutCoeffsInToElmts since lambda space is mapped
4287 // during the solve.
4288 int nTraces = (*m_exp)[i]->GetNtraces();
4289 Array<OneD, Array<OneD, NekDouble>> traceCoeffs(nTraces);
4290
4291 for (e = 0; e < (*m_exp)[i]->GetNtraces(); ++e)
4292 {
4293 edgedir = (*m_exp)[i]->GetTraceOrient(e);
4294 ncoeff_trace = elmtToTrace[i][e]->GetNcoeffs();
4295 traceCoeffs[e] = Array<OneD, NekDouble>(ncoeff_trace);
4296 Vmath::Vcopy(ncoeff_trace, trace_lambda, 1, traceCoeffs[e], 1);
4297 if (dim == 2)
4298 {
4299 elmtToTrace[i][e]->SetCoeffsToOrientation(
4300 edgedir, traceCoeffs[e], traceCoeffs[e]);
4301 }
4302 else
4303 {
4304 (*m_exp)[i]
4305 ->as<LocalRegions::Expansion3D>()
4306 ->SetFaceToGeomOrientation(e, traceCoeffs[e]);
4307 }
4308 trace_lambda = trace_lambda + ncoeff_trace;
4309 }
4310
4311 // creating orthogonal expansion (checking if we have quads or
4312 // triangles)
4313 LibUtilities::ShapeType shape = (*m_exp)[i]->DetShapeType();
4314 switch (shape)
4315 {
4317 {
4318 const LibUtilities::PointsKey PkeyQ1(
4320 const LibUtilities::PointsKey PkeyQ2(
4322 LibUtilities::BasisKey BkeyQ1(LibUtilities::eOrtho_A,
4323 num_modes[0], PkeyQ1);
4324 LibUtilities::BasisKey BkeyQ2(LibUtilities::eOrtho_A,
4325 num_modes[1], PkeyQ2);
4327 std::dynamic_pointer_cast<SpatialDomains::QuadGeom>(
4328 (*m_exp)[i]->GetGeom());
4330 BkeyQ1, BkeyQ2, qGeom);
4331 }
4332 break;
4334 {
4335 const LibUtilities::PointsKey PkeyT1(
4337 const LibUtilities::PointsKey PkeyT2(
4338 num_points[1], LibUtilities::eGaussRadauMAlpha1Beta0);
4339 LibUtilities::BasisKey BkeyT1(LibUtilities::eOrtho_A,
4340 num_modes[0], PkeyT1);
4341 LibUtilities::BasisKey BkeyT2(LibUtilities::eOrtho_B,
4342 num_modes[1], PkeyT2);
4344 std::dynamic_pointer_cast<SpatialDomains::TriGeom>(
4345 (*m_exp)[i]->GetGeom());
4347 BkeyT1, BkeyT2, tGeom);
4348 }
4349 break;
4351 {
4352 const LibUtilities::PointsKey PkeyH1(
4354 const LibUtilities::PointsKey PkeyH2(
4356 const LibUtilities::PointsKey PkeyH3(
4358 LibUtilities::BasisKey BkeyH1(LibUtilities::eOrtho_A,
4359 num_modes[0], PkeyH1);
4360 LibUtilities::BasisKey BkeyH2(LibUtilities::eOrtho_A,
4361 num_modes[1], PkeyH2);
4362 LibUtilities::BasisKey BkeyH3(LibUtilities::eOrtho_A,
4363 num_modes[2], PkeyH3);
4365 std::dynamic_pointer_cast<SpatialDomains::HexGeom>(
4366 (*m_exp)[i]->GetGeom());
4368 BkeyH1, BkeyH2, BkeyH3, hGeom);
4369 }
4370 break;
4372 {
4373 const LibUtilities::PointsKey PkeyT1(
4375 const LibUtilities::PointsKey PkeyT2(
4376 num_points[1], LibUtilities::eGaussRadauMAlpha1Beta0);
4377 const LibUtilities::PointsKey PkeyT3(
4378 num_points[2], LibUtilities::eGaussRadauMAlpha2Beta0);
4379 LibUtilities::BasisKey BkeyT1(LibUtilities::eOrtho_A,
4380 num_modes[0], PkeyT1);
4381 LibUtilities::BasisKey BkeyT2(LibUtilities::eOrtho_B,
4382 num_modes[1], PkeyT2);
4383 LibUtilities::BasisKey BkeyT3(LibUtilities::eOrtho_C,
4384 num_modes[2], PkeyT3);
4386 std::dynamic_pointer_cast<SpatialDomains::TetGeom>(
4387 (*m_exp)[i]->GetGeom());
4389 BkeyT1, BkeyT2, BkeyT3, tGeom);
4390 }
4391 break;
4392 default:
4394 "Wrong shape type, HDG postprocessing is not "
4395 "implemented");
4396 };
4397
4398 // DGDeriv
4399 // (d/dx w, d/dx q_0)
4400 (*m_exp)[i]->DGDeriv(0, tmp_coeffs = coeffs + m_coeff_offset[i],
4401 elmtToTrace[i], traceCoeffs, out_tmp);
4402 (*m_exp)[i]->BwdTrans(out_tmp, qrhs);
4403 ppExp->IProductWRTDerivBase(0, qrhs, force);
4404
4405 // + (d/dy w, d/dy q_1)
4406 (*m_exp)[i]->DGDeriv(1, tmp_coeffs = coeffs + m_coeff_offset[i],
4407 elmtToTrace[i], traceCoeffs, out_tmp);
4408
4409 (*m_exp)[i]->BwdTrans(out_tmp, qrhs);
4410 ppExp->IProductWRTDerivBase(1, qrhs, out_tmp);
4411
4412 Vmath::Vadd(nm_elmt, force, 1, out_tmp, 1, force, 1);
4413
4414 // determine force[0] = (1,u)
4415 (*m_exp)[i]->BwdTrans(tmp_coeffs = coeffs + m_coeff_offset[i], qrhs);
4416 force[0] = (*m_exp)[i]->Integral(qrhs);
4417
4418 // multiply by inverse Laplacian matrix
4419 // get matrix inverse
4420 LocalRegions::MatrixKey lapkey(StdRegions::eInvLaplacianWithUnityMean,
4421 ppExp->DetShapeType(), *ppExp);
4422 DNekScalMatSharedPtr lapsys = ppExp->GetLocMatrix(lapkey);
4423
4424 NekVector<NekDouble> in(nm_elmt, force, eWrapper);
4425 NekVector<NekDouble> out(nm_elmt);
4426
4427 out = (*lapsys) * in;
4428
4429 // Transforming back to modified basis
4430 Array<OneD, NekDouble> work(nq_elmt);
4431 ppExp->BwdTrans(out.GetPtr(), work);
4432 (*m_exp)[i]->FwdTrans(work, tmp_coeffs = outarray + m_coeff_offset[i]);
4433 }
4434}
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition: ErrorUtil.hpp:202
int GetExpSize(void)
This function returns the number of elements in the expansion.
Definition: ExpList.h:2038
Array< OneD, int > m_coeff_offset
Offset of elemental data into the array m_coeffs.
Definition: ExpList.h:1123
ExpansionType m_expType
Expansion type.
Definition: ExpList.h:1051
@ eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:51
@ eOrtho_A
Principle Orthogonal Functions .
Definition: BasisType.h:42
@ eOrtho_C
Principle Orthogonal Functions .
Definition: BasisType.h:46
@ eOrtho_B
Principle Orthogonal Functions .
Definition: BasisType.h:44
std::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:66
std::shared_ptr< QuadGeom > QuadGeomSharedPtr
Definition: HexGeom.h:45
std::shared_ptr< HexGeom > HexGeomSharedPtr
Definition: HexGeom.h:84
std::shared_ptr< TetGeom > TetGeomSharedPtr
Definition: TetGeom.h:82
std::shared_ptr< TriGeom > TriGeomSharedPtr
Definition: TriGeom.h:56
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.hpp:180
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825

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

◆ FillBwdWithBoundCond()

void Nektar::MultiRegions::DisContField::FillBwdWithBoundCond ( const Array< OneD, NekDouble > &  Fwd,
Array< OneD, NekDouble > &  Bwd,
const Array< OneD, const SpatialDomains::BoundaryConditionShPtr > &  bndConditions,
const Array< OneD, const ExpListSharedPtr > &  bndCondExpansions,
bool  PutFwdInBwdOnBCs 
)
protected

function allowing different BCs to be passed to routine

Definition at line 3029 of file DisContField.cpp.

3035{
3036
3037 // Fill boundary conditions into missing elements
3038 if (PutFwdInBwdOnBCs) // just set Bwd value to be Fwd value on BCs
3039 {
3040 // Fill boundary conditions into missing elements
3041 for (int n = 0, cnt = 0; n < bndCondExpansions.size(); ++n)
3042 {
3043 if (bndConditions[n]->GetBoundaryConditionType() ==
3045 {
3046 auto ne = bndCondExpansions[n]->GetExpSize();
3047 for (int e = 0; e < ne; ++e)
3048 {
3049 auto npts = bndCondExpansions[n]->GetExp(e)->GetTotPoints();
3050 auto id2 = m_trace->GetPhys_Offset(
3051 m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
3052 Vmath::Vcopy(npts, &Fwd[id2], 1, &Bwd[id2], 1);
3053 }
3054
3055 cnt += ne;
3056 }
3057 else if (bndConditions[n]->GetBoundaryConditionType() ==
3059 bndConditions[n]->GetBoundaryConditionType() ==
3061 {
3062 auto ne = bndCondExpansions[n]->GetExpSize();
3063 for (int e = 0; e < ne; ++e)
3064 {
3065 auto npts = bndCondExpansions[n]->GetExp(e)->GetTotPoints();
3066 auto id2 = m_trace->GetPhys_Offset(
3067 m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
3068 Vmath::Vcopy(npts, &Fwd[id2], 1, &Bwd[id2], 1);
3069 }
3070 cnt += ne;
3071 }
3072 else if (bndConditions[n]->GetBoundaryConditionType() !=
3074 {
3075 ASSERTL0(false, "Method not set up for this "
3076 "boundary condition.");
3077 }
3078 }
3079 }
3080 else
3081 {
3082 for (int n = 0, cnt = 0; n < bndCondExpansions.size(); ++n)
3083 {
3084 if (bndConditions[n]->GetBoundaryConditionType() ==
3086 {
3087 auto ne = bndCondExpansions[n]->GetExpSize();
3088 for (int e = 0; e < ne; ++e)
3089 {
3090 auto npts = bndCondExpansions[n]->GetExp(e)->GetTotPoints();
3091 auto id1 = bndCondExpansions[n]->GetPhys_Offset(e);
3092 auto id2 = m_trace->GetPhys_Offset(
3093 m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
3094
3095 Vmath::Vcopy(npts, &(bndCondExpansions[n]->GetPhys())[id1],
3096 1, &Bwd[id2], 1);
3097 }
3098 cnt += ne;
3099 }
3100 else if (bndConditions[n]->GetBoundaryConditionType() ==
3102 bndConditions[n]->GetBoundaryConditionType() ==
3104 {
3105 auto ne = m_bndCondExpansions[n]->GetExpSize();
3106 for (int e = 0; e < ne; ++e)
3107 {
3108 auto npts = bndCondExpansions[n]->GetExp(e)->GetTotPoints();
3109 auto id1 = bndCondExpansions[n]->GetPhys_Offset(e);
3110
3111 ASSERTL0((bndCondExpansions[n]->GetPhys())[id1] == 0.0,
3112 "Method not set up for non-zero "
3113 "Neumann boundary condition");
3114 auto id2 = m_trace->GetPhys_Offset(
3115 m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
3116
3117 Vmath::Vcopy(npts, &Fwd[id2], 1, &Bwd[id2], 1);
3118 }
3119
3120 cnt += ne;
3121 }
3122 else if (bndConditions[n]->GetBoundaryConditionType() ==
3124 {
3125 // Do nothing
3126 }
3127 else if (bndConditions[n]->GetBoundaryConditionType() !=
3129 {
3131 "Method not set up for this boundary "
3132 "condition.");
3133 }
3134 }
3135 }
3136}
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
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:2030

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_trace, m_traceMap, NEKERROR, and Vmath::Vcopy().

Referenced by GetFwdBwdTracePhys(), and v_FillBwdWithBoundCond().

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

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

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(), Vmath::Vsum(), Nektar::UnitTests::w(), and Nektar::UnitTests::z().

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

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

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

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

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

Referenced by DisContField().

◆ GetFwdBwdTracePhys()

void Nektar::MultiRegions::DisContField::GetFwdBwdTracePhys ( const Array< OneD, const NekDouble > &  field,
Array< OneD, NekDouble > &  Fwd,
Array< OneD, NekDouble > &  Bwd,
const Array< OneD, const SpatialDomains::BoundaryConditionShPtr > &  BndCond,
const Array< OneD, const ExpListSharedPtr > &  BndCondExp 
)

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

Definition at line 2912 of file DisContField.cpp.

2917{
2918 ASSERTL1((*m_exp)[0]->GetBasis(0)->GetBasisType() !=
2920 "Routine not set up for Gauss Lagrange points distribution");
2921
2922 // blocked routine
2923 Array<OneD, NekDouble> tracevals(m_locTraceToTraceMap->GetNLocTracePts());
2924
2925 m_locTraceToTraceMap->LocTracesFromField(field, tracevals);
2926 m_locTraceToTraceMap->InterpLocTracesToTrace(0, tracevals, Fwd);
2927
2928 Array<OneD, NekDouble> invals =
2929 tracevals + m_locTraceToTraceMap->GetNFwdLocTracePts();
2930 m_locTraceToTraceMap->InterpLocTracesToTrace(1, invals, Bwd);
2931
2933
2934 FillBwdWithBoundCond(Fwd, Bwd, BndCond, BndCondExp, false);
2935
2936 // Do parallel exchange for forwards/backwards spaces.
2937 m_traceMap->GetAssemblyCommDG()->PerformExchange(Fwd, Bwd);
2938
2939 // Do exchange of interface traces (local and parallel)
2940 // We may have to split this out into separate local and
2941 // parallel for IP method???
2942 m_interfaceMap->ExchangeTrace(Fwd, Bwd);
2943}
void v_PeriodicBwdCopy(const Array< OneD, const NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd) override
void FillBwdWithBoundCond(const Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd, const Array< OneD, const SpatialDomains::BoundaryConditionShPtr > &bndConditions, const Array< OneD, const ExpListSharedPtr > &BndCondExpansions, bool PutFwdInBwdOnBCs)
@ eGauss_Lagrange
Lagrange Polynomials using the Gauss points.
Definition: BasisType.h:57

References ASSERTL1, Nektar::LibUtilities::eGauss_Lagrange, FilterPython_Function::field, FillBwdWithBoundCond(), Nektar::MultiRegions::ExpList::m_exp, m_interfaceMap, m_locTraceToTraceMap, m_traceMap, and v_PeriodicBwdCopy().

◆ GetGlobalBndLinSys()

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

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

Definition at line 2704 of file DisContField.cpp.

2706{
2707 ASSERTL0(mkey.GetMatrixType() == StdRegions::eHybridDGHelmBndLam,
2708 "Routine currently only tested for HybridDGHelmholtz");
2709
2710 ASSERTL1(mkey.GetGlobalSysSolnType() != eDirectFullMatrix,
2711 "Full matrix global systems are not supported for HDG "
2712 "expansions");
2713
2714 ASSERTL1(mkey.GetGlobalSysSolnType() == m_traceMap->GetGlobalSysSolnType(),
2715 "The local to global map is not set up for the requested "
2716 "solution type");
2717
2718 GlobalLinSysSharedPtr glo_matrix;
2719 auto matrixIter = m_globalBndMat->find(mkey);
2720
2721 if (matrixIter == m_globalBndMat->end())
2722 {
2723 glo_matrix = GenGlobalBndLinSys(mkey, m_traceMap);
2724 (*m_globalBndMat)[mkey] = glo_matrix;
2725 }
2726 else
2727 {
2728 glo_matrix = matrixIter->second;
2729 }
2730
2731 return glo_matrix;
2732}
std::shared_ptr< GlobalLinSys > GenGlobalBndLinSys(const GlobalLinSysKey &mkey, const AssemblyMapSharedPtr &locToGloMap)
Generate a GlobalLinSys from information provided by the key "mkey" and the mapping provided in LocTo...
Definition: ExpList.cpp:3093
std::shared_ptr< GlobalLinSys > GlobalLinSysSharedPtr
Pointer to a GlobalLinSys object.
Definition: GlobalLinSys.h:51

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

125 {
126 LocTraceToTraceMap = m_locTraceToTraceMap;
127 }

References m_locTraceToTraceMap.

◆ GetNegatedFluxNormal()

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

Definition at line 2768 of file DisContField.cpp.

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

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

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

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

4176{
4177 int i, e, ncoeff_edge;
4178 Array<OneD, const NekDouble> tmp_coeffs;
4179 Array<OneD, NekDouble> out_d(m_ncoeffs), out_tmp;
4180
4181 Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
4182 m_traceMap->GetElmtToTrace();
4183
4185
4186 int cnt;
4187 int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs();
4188 Array<OneD, NekDouble> loc_lambda(LocBndCoeffs), edge_lambda;
4189 m_traceMap->GlobalToLocalBnd(m_trace->GetCoeffs(), loc_lambda);
4190
4191 edge_lambda = loc_lambda;
4192
4193 // Calculate Q using standard DG formulation.
4194 for (i = cnt = 0; i < GetExpSize(); ++i)
4195 {
4196 // Probably a better way of setting up lambda than this.
4197 // Note cannot use PutCoeffsInToElmts since lambda space
4198 // is mapped during the solve.
4199 int nEdges = (*m_exp)[i]->GetGeom()->GetNumEdges();
4200 Array<OneD, Array<OneD, NekDouble>> edgeCoeffs(nEdges);
4201
4202 for (e = 0; e < nEdges; ++e)
4203 {
4204 edgedir = (*m_exp)[i]->GetTraceOrient(e);
4205 ncoeff_edge = elmtToTrace[i][e]->GetNcoeffs();
4206 edgeCoeffs[e] = Array<OneD, NekDouble>(ncoeff_edge);
4207 Vmath::Vcopy(ncoeff_edge, edge_lambda, 1, edgeCoeffs[e], 1);
4208 elmtToTrace[i][e]->SetCoeffsToOrientation(edgedir, edgeCoeffs[e],
4209 edgeCoeffs[e]);
4210 edge_lambda = edge_lambda + ncoeff_edge;
4211 }
4212
4213 (*m_exp)[i]->DGDeriv(dir, tmp_coeffs = coeffs + m_coeff_offset[i],
4214 elmtToTrace[i], edgeCoeffs, out_tmp = out_d + cnt);
4215 cnt += (*m_exp)[i]->GetNcoeffs();
4216 }
4217
4218 Array<OneD, NekDouble> phys(m_npoints);
4219 BwdTrans(out_d, phys);
4220 Vmath::Vsub(m_npoints, phys, 1, soln, 1, phys, 1);
4221 return L2(phys);
4222}
int m_ncoeffs
The total number of local degrees of freedom. m_ncoeffs .
Definition: ExpList.h:1063
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:514
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:1717
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.hpp:220

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

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

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

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

◆ SetUpDG()

void Nektar::MultiRegions::DisContField::SetUpDG ( const std::string  variable = "DefaultVar",
const Collections::ImplementationType  ImpType = Collections::eNoImpType 
)
protected

Set up all DG member variables and maps.

Definition at line 161 of file DisContField.cpp.

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

3424{
3425 ASSERTL0(m_expType != e1D, "This method is not setup or "
3426 "tested for 1D expansion");
3427 Array<OneD, NekDouble> Coeffs(m_trace->GetNcoeffs());
3428 m_trace->IProductWRTBase(Fwd, Coeffs);
3429 m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(0, Coeffs, outarray);
3430 m_trace->IProductWRTBase(Bwd, Coeffs);
3431 m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(1, Coeffs, outarray);
3432}

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

3252{
3253 if (m_expType == e1D)
3254 {
3255 int n, offset, t_offset;
3256 Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>>
3257 &elmtToTrace = m_traceMap->GetElmtToTrace();
3258 vector<bool> negatedFluxNormal = GetNegatedFluxNormal();
3259 for (n = 0; n < GetExpSize(); ++n)
3260 {
3261 // Number of coefficients on each element
3262 int e_ncoeffs = (*m_exp)[n]->GetNcoeffs();
3263 offset = GetCoeff_Offset(n);
3264 // Implementation for every points except Gauss points
3265 if ((*m_exp)[n]->GetBasis(0)->GetBasisType() !=
3267 {
3268 t_offset =
3269 GetTrace()->GetCoeff_Offset(elmtToTrace[n][0]->GetElmtId());
3270 if (negatedFluxNormal[2 * n])
3271 {
3272 outarray[offset] -= Fn[t_offset];
3273 }
3274 else
3275 {
3276 outarray[offset] += Fn[t_offset];
3277 }
3278 t_offset =
3279 GetTrace()->GetCoeff_Offset(elmtToTrace[n][1]->GetElmtId());
3280 if (negatedFluxNormal[2 * n + 1])
3281 {
3282 outarray[offset + (*m_exp)[n]->GetVertexMap(1)] -=
3283 Fn[t_offset];
3284 }
3285 else
3286 {
3287 outarray[offset + (*m_exp)[n]->GetVertexMap(1)] +=
3288 Fn[t_offset];
3289 }
3290 }
3291 else
3292 {
3293 int j;
3294 static DNekMatSharedPtr m_Ixm, m_Ixp;
3295 static int sav_ncoeffs = 0;
3296 if (!m_Ixm.get() || e_ncoeffs != sav_ncoeffs)
3297 {
3299 const LibUtilities::PointsKey BS_p(
3301 const LibUtilities::BasisKey BS_k(
3302 LibUtilities::eGauss_Lagrange, e_ncoeffs, BS_p);
3303 BASE = LibUtilities::BasisManager()[BS_k];
3304 Array<OneD, NekDouble> coords(1, 0.0);
3305 coords[0] = -1.0;
3306 m_Ixm = BASE->GetI(coords);
3307 coords[0] = 1.0;
3308 m_Ixp = BASE->GetI(coords);
3309 sav_ncoeffs = e_ncoeffs;
3310 }
3311 t_offset =
3312 GetTrace()->GetCoeff_Offset(elmtToTrace[n][0]->GetElmtId());
3313 if (negatedFluxNormal[2 * n])
3314 {
3315 for (j = 0; j < e_ncoeffs; j++)
3316 {
3317 outarray[offset + j] -=
3318 (m_Ixm->GetPtr())[j] * Fn[t_offset];
3319 }
3320 }
3321 else
3322 {
3323 for (j = 0; j < e_ncoeffs; j++)
3324 {
3325 outarray[offset + j] +=
3326 (m_Ixm->GetPtr())[j] * Fn[t_offset];
3327 }
3328 }
3329 t_offset =
3330 GetTrace()->GetCoeff_Offset(elmtToTrace[n][1]->GetElmtId());
3331 if (negatedFluxNormal[2 * n + 1])
3332 {
3333 for (j = 0; j < e_ncoeffs; j++)
3334 {
3335 outarray[offset + j] -=
3336 (m_Ixp->GetPtr())[j] * Fn[t_offset];
3337 }
3338 }
3339 else
3340 {
3341 for (j = 0; j < e_ncoeffs; j++)
3342 {
3343 outarray[offset + j] +=
3344 (m_Ixp->GetPtr())[j] * Fn[t_offset];
3345 }
3346 }
3347 }
3348 }
3349 }
3350 else // other dimensions
3351 {
3352 // Basis definition on each element
3353 LibUtilities::BasisSharedPtr basis = (*m_exp)[0]->GetBasis(0);
3354 if ((m_expType != e1D) &&
3355 (basis->GetBasisType() != LibUtilities::eGauss_Lagrange))
3356 {
3357 Array<OneD, NekDouble> Fcoeffs(m_trace->GetNcoeffs());
3358 m_trace->IProductWRTBase(Fn, Fcoeffs);
3359 m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(Fcoeffs,
3360 outarray);
3361 }
3362 else
3363 {
3364 int e, n, offset, t_offset;
3365 Array<OneD, NekDouble> e_outarray;
3366 Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>>
3367 &elmtToTrace = m_traceMap->GetElmtToTrace();
3368 if (m_expType == e2D)
3369 {
3370 for (n = 0; n < GetExpSize(); ++n)
3371 {
3372 offset = GetCoeff_Offset(n);
3373 for (e = 0; e < (*m_exp)[n]->GetNtraces(); ++e)
3374 {
3375 t_offset = GetTrace()->GetPhys_Offset(
3376 elmtToTrace[n][e]->GetElmtId());
3377 (*m_exp)[n]->AddEdgeNormBoundaryInt(
3378 e, elmtToTrace[n][e], Fn + t_offset,
3379 e_outarray = outarray + offset);
3380 }
3381 }
3382 }
3383 else if (m_expType == e3D)
3384 {
3385 for (n = 0; n < GetExpSize(); ++n)
3386 {
3387 offset = GetCoeff_Offset(n);
3388 for (e = 0; e < (*m_exp)[n]->GetNtraces(); ++e)
3389 {
3390 t_offset = GetTrace()->GetPhys_Offset(
3391 elmtToTrace[n][e]->GetElmtId());
3392 (*m_exp)[n]->AddFaceNormBoundaryInt(
3393 e, elmtToTrace[n][e], Fn + t_offset,
3394 e_outarray = outarray + offset);
3395 }
3396 }
3397 }
3398 }
3399 }
3400}
std::vector< bool > & GetNegatedFluxNormal(void)
std::shared_ptr< ExpList > & GetTrace()
Definition: ExpList.h:2153
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:2078
BasisManagerT & BasisManager(void)
std::shared_ptr< Basis > BasisSharedPtr
@ eGaussGaussLegendre
1D Gauss-Gauss-Legendre quadrature points
Definition: PointsType.h:46
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 4454 of file DisContField.cpp.

4458{
4459 Array<OneD, NekDouble> FCoeffs(m_trace->GetNcoeffs());
4460
4461 m_trace->IProductWRTBase(FwdFlux, FCoeffs);
4462 m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(1, FCoeffs, outarray);
4463 m_trace->IProductWRTBase(BwdFlux, FCoeffs);
4464 m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(0, FCoeffs, outarray);
4465}

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

3151{
3152 // Basis definition on each element
3153 LibUtilities::BasisSharedPtr basis = (*m_exp)[0]->GetBasis(0);
3154 if (basis->GetBasisType() != LibUtilities::eGauss_Lagrange)
3155 {
3156 Array<OneD, NekDouble> tracevals(
3157 m_locTraceToTraceMap->GetNLocTracePts(), 0.0);
3158
3159 Array<OneD, NekDouble> invals =
3160 tracevals + m_locTraceToTraceMap->GetNFwdLocTracePts();
3161
3162 m_locTraceToTraceMap->InterpTraceToLocTrace(1, Bwd, invals);
3163
3164 m_locTraceToTraceMap->InterpTraceToLocTrace(0, Fwd, tracevals);
3165
3166 m_locTraceToTraceMap->AddLocTracesToField(tracevals, field);
3167 }
3168 else
3169 {
3171 "v_AddTraceQuadPhysToField not coded for eGauss_Lagrange");
3172 }
3173}

References Nektar::ErrorUtil::efatal, Nektar::LibUtilities::eGauss_Lagrange, FilterPython_Function::field, 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 3602 of file DisContField.cpp.

3606{
3607 int i;
3608 int npoints;
3609
3611
3612 for (i = 0; i < m_bndCondExpansions.size(); ++i)
3613 {
3614 if (time == 0.0 || m_bndConditions[i]->IsTimeDependent())
3615 {
3616 m_bndCondBndWeight[i] = 1.0;
3617 locExpList = m_bndCondExpansions[i];
3618
3619 npoints = locExpList->GetNpoints();
3620 Array<OneD, NekDouble> x0(npoints, 0.0);
3621 Array<OneD, NekDouble> x1(npoints, 0.0);
3622 Array<OneD, NekDouble> x2(npoints, 0.0);
3623
3624 locExpList->GetCoords(x0, x1, x2);
3625
3626 if (x2_in != NekConstants::kNekUnsetDouble &&
3628 {
3629 Vmath::Fill(npoints, x2_in, x1, 1);
3630 Vmath::Fill(npoints, x3_in, x2, 1);
3631 }
3632 else if (x2_in != NekConstants::kNekUnsetDouble)
3633 {
3634 Vmath::Fill(npoints, x2_in, x2, 1);
3635 }
3636
3637 // treat 1D expansions separately since we only
3638 // require an evaluation at a point rather than
3639 // any projections or inner products that are not
3640 // available in a PointExp
3641 if (m_expType == e1D)
3642 {
3643 if (m_bndConditions[i]->GetBoundaryConditionType() ==
3645 {
3646
3647 m_bndCondExpansions[i]->SetCoeff(
3648 0, (std::static_pointer_cast<
3649 SpatialDomains::DirichletBoundaryCondition>(
3650 m_bndConditions[i])
3651 ->m_dirichletCondition)
3652 .Evaluate(x0[0], x1[0], x2[0], time));
3653 m_bndCondExpansions[i]->SetPhys(
3654 0, m_bndCondExpansions[i]->GetCoeff(0));
3655 }
3656 else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3658 {
3659 m_bndCondExpansions[i]->SetCoeff(
3660 0, (std::static_pointer_cast<
3661 SpatialDomains::NeumannBoundaryCondition>(
3662 m_bndConditions[i])
3663 ->m_neumannCondition)
3664 .Evaluate(x0[0], x1[0], x2[0], time));
3665 }
3666 else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3668 {
3669 m_bndCondExpansions[i]->SetCoeff(
3670 0, (std::static_pointer_cast<
3671 SpatialDomains::RobinBoundaryCondition>(
3672 m_bndConditions[i])
3673 ->m_robinFunction)
3674 .Evaluate(x0[0], x1[0], x2[0], time));
3675 }
3676 else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3678 {
3679 continue;
3680 }
3681 else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3683 {
3684 }
3685 else
3686 {
3688 "This type of BC not implemented yet");
3689 }
3690 }
3691 else // 2D and 3D versions
3692 {
3693 if (m_bndConditions[i]->GetBoundaryConditionType() ==
3695 {
3697 std::static_pointer_cast<
3698 SpatialDomains::DirichletBoundaryCondition>(
3699 m_bndConditions[i]);
3700
3701 Array<OneD, NekDouble> valuesFile(npoints, 1.0),
3702 valuesExp(npoints, 1.0);
3703
3704 string bcfilename = bcPtr->m_filename;
3705 string exprbcs = bcPtr->m_expr;
3706
3707 if (bcfilename != "")
3708 {
3709 locExpList->ExtractCoeffsFromFile(
3710 bcfilename, bcPtr->GetComm(), varName,
3711 locExpList->UpdateCoeffs());
3712 locExpList->BwdTrans(locExpList->GetCoeffs(),
3713 locExpList->UpdatePhys());
3714 valuesFile = locExpList->GetPhys();
3715 }
3716
3717 if (exprbcs != "")
3718 {
3719 LibUtilities::Equation condition =
3720 std::static_pointer_cast<
3721 SpatialDomains::DirichletBoundaryCondition>(
3722 m_bndConditions[i])
3723 ->m_dirichletCondition;
3724
3725 condition.Evaluate(x0, x1, x2, time, valuesExp);
3726 }
3727
3728 Vmath::Vmul(npoints, valuesExp, 1, valuesFile, 1,
3729 locExpList->UpdatePhys(), 1);
3730 locExpList->FwdTransBndConstrained(
3731 locExpList->GetPhys(), locExpList->UpdateCoeffs());
3732 }
3733 else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3735 {
3737 std::static_pointer_cast<
3738 SpatialDomains::NeumannBoundaryCondition>(
3739 m_bndConditions[i]);
3740 string bcfilename = bcPtr->m_filename;
3741
3742 if (bcfilename != "")
3743 {
3744 locExpList->ExtractCoeffsFromFile(
3745 bcfilename, bcPtr->GetComm(), varName,
3746 locExpList->UpdateCoeffs());
3747 locExpList->BwdTrans(locExpList->GetCoeffs(),
3748 locExpList->UpdatePhys());
3749 }
3750 else
3751 {
3752 LibUtilities::Equation condition =
3753 std::static_pointer_cast<
3754 SpatialDomains::NeumannBoundaryCondition>(
3755 m_bndConditions[i])
3756 ->m_neumannCondition;
3757 condition.Evaluate(x0, x1, x2, time,
3758 locExpList->UpdatePhys());
3759 }
3760
3761 locExpList->IProductWRTBase(locExpList->GetPhys(),
3762 locExpList->UpdateCoeffs());
3763 }
3764 else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3766 {
3768 std::static_pointer_cast<
3769 SpatialDomains::RobinBoundaryCondition>(
3770 m_bndConditions[i]);
3771
3772 string bcfilename = bcPtr->m_filename;
3773
3774 if (bcfilename != "")
3775 {
3776 locExpList->ExtractCoeffsFromFile(
3777 bcfilename, bcPtr->GetComm(), varName,
3778 locExpList->UpdateCoeffs());
3779 locExpList->BwdTrans(locExpList->GetCoeffs(),
3780 locExpList->UpdatePhys());
3781 }
3782 else
3783 {
3784 LibUtilities::Equation condition =
3785 std::static_pointer_cast<
3786 SpatialDomains::RobinBoundaryCondition>(
3787 m_bndConditions[i])
3788 ->m_robinFunction;
3789 condition.Evaluate(x0, x1, x2, time,
3790 locExpList->UpdatePhys());
3791 }
3792
3793 locExpList->IProductWRTBase(locExpList->GetPhys(),
3794 locExpList->UpdateCoeffs());
3795 }
3796 else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3798 {
3799 continue;
3800 }
3801 else
3802 {
3804 "This type of BC not implemented yet");
3805 }
3806 }
3807 }
3808 }
3809}
NekDouble GetCoeff(int i)
Get the i th value (coefficient) of m_coeffs.
Definition: ExpList.h:2010
static const NekDouble kNekUnsetDouble
std::shared_ptr< DirichletBoundaryCondition > DirichletBCShPtr
Definition: Conditions.h:214
std::shared_ptr< NeumannBoundaryCondition > NeumannBCShPtr
Definition: Conditions.h:215
std::shared_ptr< RobinBoundaryCondition > RobinBCShPtr
Definition: Conditions.h:216
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.hpp:72
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.hpp:54

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

3176{
3177 ASSERTL1(m_physState == true, "local physical space is not true ");
3178 v_ExtractTracePhys(m_phys, outarray);
3179}
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:1107
Array< OneD, NekDouble > m_phys
The global expansion evaluated at the quadrature points.
Definition: ExpList.h:1099

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

3197{
3198 LibUtilities::BasisSharedPtr basis = (*m_exp)[0]->GetBasis(0);
3199 if ((basis->GetBasisType() != LibUtilities::eGauss_Lagrange))
3200 {
3201 Vmath::Zero(outarray.size(), outarray, 1);
3202 Array<OneD, NekDouble> tracevals(
3203 m_locTraceToTraceMap->GetNFwdLocTracePts());
3204 m_locTraceToTraceMap->FwdLocTracesFromField(inarray, tracevals);
3205 m_locTraceToTraceMap->InterpLocTracesToTrace(0, tracevals, outarray);
3206 m_traceMap->GetAssemblyCommDG()->PerformExchange(outarray, outarray);
3207 }
3208 else
3209 {
3210 // Loop over elemente and collect forward expansion
3211 int nexp = GetExpSize();
3212 int n, p, offset, phys_offset;
3213 Array<OneD, NekDouble> t_tmp;
3214 Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>>
3215 &elmtToTrace = m_traceMap->GetElmtToTrace();
3216 ASSERTL1(outarray.size() >= m_trace->GetNpoints(),
3217 "input array is of insufficient length");
3218 for (n = 0; n < nexp; ++n)
3219 {
3220 phys_offset = GetPhys_Offset(n);
3221 for (p = 0; p < (*m_exp)[n]->GetNtraces(); ++p)
3222 {
3223 offset =
3224 m_trace->GetPhys_Offset(elmtToTrace[n][p]->GetElmtId());
3225 (*m_exp)[n]->GetTracePhysVals(p, elmtToTrace[n][p],
3226 inarray + phys_offset,
3227 t_tmp = outarray + offset);
3228 }
3229 }
3230 }
3231}
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:2085
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.hpp:273

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

3022{
3024 PutFwdInBwdOnBCs);
3025}

References FillBwdWithBoundCond(), m_bndCondExpansions, and m_bndConditions.

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

3816{
3817 int cnt;
3818 int npts;
3819 int e = 0;
3820
3821 // Fill boundary conditions into missing elements
3822 int id2 = 0;
3823
3824 for (int n = cnt = 0; n < m_bndCondExpansions.size(); ++n)
3825 {
3826
3827 if (m_bndConditions[n]->GetBoundaryConditionType() ==
3829 {
3830 for (e = 0; e < m_bndCondExpansions[n]->GetExpSize(); ++e)
3831 {
3832 npts = m_bndCondExpansions[n]->GetExp(e)->GetTotPoints();
3833 id2 = m_trace->GetPhys_Offset(
3834 m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
3835 Vmath::Fill(npts, m_bndCondBndWeight[n], &weightave[id2], 1);
3836 Vmath::Fill(npts, 0.0, &weightjmp[id2], 1);
3837 }
3838
3839 cnt += e;
3840 }
3841 else if (m_bndConditions[n]->GetBoundaryConditionType() ==
3843 m_bndConditions[n]->GetBoundaryConditionType() ==
3845 {
3846 for (e = 0; e < m_bndCondExpansions[n]->GetExpSize(); ++e)
3847 {
3848 npts = m_bndCondExpansions[n]->GetExp(e)->GetTotPoints();
3849 id2 = m_trace->GetPhys_Offset(
3850 m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
3851 Vmath::Fill(npts, m_bndCondBndWeight[n], &weightave[id2], 1);
3852 Vmath::Fill(npts, 0.0, &weightjmp[id2], 1);
3853 }
3854
3855 cnt += e;
3856 }
3857 else if (m_bndConditions[n]->GetBoundaryConditionType() !=
3859 {
3861 "Method not set up for this boundary condition.");
3862 }
3863 }
3864}

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

3139{
3140 return m_bndCondBndWeight;
3141}

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

2839{
2840 return m_bndCondExpansions;
2841}

References m_bndCondExpansions.

◆ v_GetBndConditions()

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

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 2843 of file DisContField.cpp.

2845{
2846 return m_bndConditions;
2847}

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

4034{
4035 int n, cnt;
4036 std::vector<unsigned int> eIDs;
4037
4038 Array<OneD, int> ElmtID, TraceID;
4039 GetBoundaryToElmtMap(ElmtID, TraceID);
4040
4041 // Skip other boundary regions
4042 for (cnt = n = 0; n < i; ++n)
4043 {
4044 cnt += m_bndCondExpansions[n]->GetExpSize();
4045 }
4046
4047 // Populate eIDs with information from BoundaryToElmtMap
4048 for (n = 0; n < m_bndCondExpansions[i]->GetExpSize(); ++n)
4049 {
4050 eIDs.push_back(ElmtID[cnt + n]);
4051 }
4052
4053 // Create expansion list
4055 *this, eIDs, DeclareCoeffPhysArrays, Collections::eNoCollection);
4056 // Copy phys and coeffs to new explist
4057 if (DeclareCoeffPhysArrays)
4058 {
4059 int nq;
4060 int offsetOld, offsetNew;
4061 Array<OneD, NekDouble> tmp1, tmp2;
4062 for (n = 0; n < result->GetExpSize(); ++n)
4063 {
4064 nq = GetExp(ElmtID[cnt + n])->GetTotPoints();
4065 offsetOld = GetPhys_Offset(ElmtID[cnt + n]);
4066 offsetNew = result->GetPhys_Offset(n);
4067 Vmath::Vcopy(nq, tmp1 = GetPhys() + offsetOld, 1,
4068 tmp2 = result->UpdatePhys() + offsetNew, 1);
4069 nq = GetExp(ElmtID[cnt + n])->GetNcoeffs();
4070 offsetOld = GetCoeff_Offset(ElmtID[cnt + n]);
4071 offsetNew = result->GetCoeff_Offset(n);
4072 Vmath::Vcopy(nq, tmp1 = GetCoeffs() + offsetOld, 1,
4073 tmp2 = result->UpdateCoeffs() + offsetNew, 1);
4074 }
4075 }
4076}
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:1944
const std::shared_ptr< LocalRegions::ExpansionVector > GetExp() const
This function returns the vector of elements in the expansion.
Definition: ExpList.h:2070

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

3870{
3871
3872 if (m_BCtoElmMap.size() == 0)
3873 {
3874 switch (m_expType)
3875 {
3876 case e1D:
3877 {
3878 map<int, int> VertGID;
3879 int i, n, id;
3880 int bid, cnt, Vid;
3881 int nbcs = 0;
3882 // Determine number of boundary condition expansions.
3883 for (i = 0; i < m_bndConditions.size(); ++i)
3884 {
3885 nbcs += m_bndCondExpansions[i]->GetExpSize();
3886 }
3887
3888 // make sure arrays are of sufficient length
3889 m_BCtoElmMap = Array<OneD, int>(nbcs, -1);
3890 m_BCtoTraceMap = Array<OneD, int>(nbcs);
3891
3892 // setup map of all global ids along boundary
3893 cnt = 0;
3894 for (n = 0; n < m_bndCondExpansions.size(); ++n)
3895 {
3896 for (i = 0; i < m_bndCondExpansions[n]->GetExpSize(); ++i)
3897 {
3898 Vid = m_bndCondExpansions[n]
3899 ->GetExp(i)
3900 ->GetGeom()
3901 ->GetVertex(0)
3902 ->GetVid();
3903 VertGID[Vid] = cnt++;
3904 }
3905 }
3906
3907 // Loop over elements and find verts that match;
3909 for (cnt = n = 0; n < GetExpSize(); ++n)
3910 {
3911 exp = (*m_exp)[n];
3912 for (i = 0; i < exp->GetNverts(); ++i)
3913 {
3914 id = exp->GetGeom()->GetVid(i);
3915
3916 if (VertGID.count(id) > 0)
3917 {
3918 bid = VertGID.find(id)->second;
3919 ASSERTL1(m_BCtoElmMap[bid] == -1,
3920 "Edge already set");
3921 m_BCtoElmMap[bid] = n;
3922 m_BCtoTraceMap[bid] = i;
3923 cnt++;
3924 }
3925 }
3926 }
3927 ASSERTL1(cnt == nbcs,
3928 "Failed to visit all boundary condtiions");
3929 }
3930 break;
3931 case e2D:
3932 {
3933 map<int, int> globalIdMap;
3934 int i, n;
3935 int cnt;
3936 int nbcs = 0;
3937
3938 // Populate global ID map (takes global geometry
3939 // ID to local expansion list ID).
3940 for (i = 0; i < GetExpSize(); ++i)
3941 {
3942 globalIdMap[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
3943 }
3944
3945 // Determine number of boundary condition expansions.
3946 for (i = 0; i < m_bndConditions.size(); ++i)
3947 {
3948 nbcs += m_bndCondExpansions[i]->GetExpSize();
3949 }
3950
3951 // Initialize arrays
3952 m_BCtoElmMap = Array<OneD, int>(nbcs);
3953 m_BCtoTraceMap = Array<OneD, int>(nbcs);
3954
3956 cnt = 0;
3957 for (n = 0; n < m_bndCondExpansions.size(); ++n)
3958 {
3959 for (i = 0; i < m_bndCondExpansions[n]->GetExpSize();
3960 ++i, ++cnt)
3961 {
3962 exp1d = m_bndCondExpansions[n]
3963 ->GetExp(i)
3964 ->as<LocalRegions::Expansion1D>();
3965
3966 // Use edge to element map from MeshGraph.
3968 m_graph->GetElementsFromEdge(exp1d->GetGeom1D());
3969
3970 m_BCtoElmMap[cnt] =
3971 globalIdMap[(*tmp)[0].first->GetGlobalID()];
3972 m_BCtoTraceMap[cnt] = (*tmp)[0].second;
3973 }
3974 }
3975 }
3976 break;
3977 case e3D:
3978 {
3979 map<int, int> globalIdMap;
3980 int i, n;
3981 int cnt;
3982 int nbcs = 0;
3983
3984 // Populate global ID map (takes global geometry ID to local
3985 // expansion list ID).
3987 for (i = 0; i < GetExpSize(); ++i)
3988 {
3989 globalIdMap[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
3990 }
3991
3992 // Determine number of boundary condition expansions.
3993 for (i = 0; i < m_bndConditions.size(); ++i)
3994 {
3995 nbcs += m_bndCondExpansions[i]->GetExpSize();
3996 }
3997
3998 // Initialize arrays
3999 m_BCtoElmMap = Array<OneD, int>(nbcs);
4000 m_BCtoTraceMap = Array<OneD, int>(nbcs);
4001
4003 for (cnt = n = 0; n < m_bndCondExpansions.size(); ++n)
4004 {
4005 for (i = 0; i < m_bndCondExpansions[n]->GetExpSize();
4006 ++i, ++cnt)
4007 {
4008 exp2d = m_bndCondExpansions[n]
4009 ->GetExp(i)
4010 ->as<LocalRegions::Expansion2D>();
4011
4013 m_graph->GetElementsFromFace(exp2d->GetGeom2D());
4014 m_BCtoElmMap[cnt] =
4015 globalIdMap[tmp->at(0).first->GetGlobalID()];
4016 m_BCtoTraceMap[cnt] = tmp->at(0).second;
4017 }
4018 }
4019 }
4020 break;
4021 default:
4022 ASSERTL1(false, "Not setup for this expansion");
4023 break;
4024 }
4025 }
4026
4027 ElmtID = m_BCtoElmMap;
4028 TraceID = m_BCtoTraceMap;
4029}
Array< OneD, int > m_BCtoTraceMap
Definition: DisContField.h:59
std::shared_ptr< Expansion2D > Expansion2DSharedPtr
Definition: Expansion1D.h:46
std::shared_ptr< Expansion1D > Expansion1DSharedPtr
Definition: Expansion1D.h:50
std::shared_ptr< Expansion3D > Expansion3DSharedPtr
Definition: Expansion2D.h:47

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

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

References Nektar::MultiRegions::ExpList::m_phys, and v_GetFwdBwdTracePhys().

Referenced by v_GetFwdBwdTracePhys().

◆ 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

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 2945 of file DisContField.cpp.

2950{
2951 // Is this zeroing necessary?
2952 // Zero forward/backward vectors.
2953 Vmath::Zero(Fwd.size(), Fwd, 1);
2954 Vmath::Zero(Bwd.size(), Bwd, 1);
2955
2956 // Basis definition on each element
2957 LibUtilities::BasisSharedPtr basis = (*m_exp)[0]->GetBasis(0);
2958 if ((basis->GetBasisType() != LibUtilities::eGauss_Lagrange))
2959 {
2960 // blocked routine
2961 Array<OneD, NekDouble> tracevals(
2962 m_locTraceToTraceMap->GetNLocTracePts());
2963
2964 m_locTraceToTraceMap->LocTracesFromField(field, tracevals);
2965 m_locTraceToTraceMap->InterpLocTracesToTrace(0, tracevals, Fwd);
2966
2967 Array<OneD, NekDouble> invals =
2968 tracevals + m_locTraceToTraceMap->GetNFwdLocTracePts();
2969 m_locTraceToTraceMap->InterpLocTracesToTrace(1, invals, Bwd);
2970 }
2971 else
2972 {
2973 // Loop over elements and collect forward expansion
2974 auto nexp = GetExpSize();
2975 Array<OneD, NekDouble> e_tmp;
2977
2978 Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>>
2979 &elmtToTrace = m_traceMap->GetElmtToTrace();
2980
2981 for (int n = 0, cnt = 0; n < nexp; ++n)
2982 {
2983 exp = (*m_exp)[n];
2984 auto phys_offset = GetPhys_Offset(n);
2985
2986 for (int e = 0; e < exp->GetNtraces(); ++e, ++cnt)
2987 {
2988 auto offset =
2989 m_trace->GetPhys_Offset(elmtToTrace[n][e]->GetElmtId());
2990
2991 e_tmp =
2992 (m_leftAdjacentTraces[cnt]) ? Fwd + offset : Bwd + offset;
2993
2994 exp->GetTracePhysVals(e, elmtToTrace[n][e], field + phys_offset,
2995 e_tmp);
2996 }
2997 }
2998 }
2999
3001
3002 if (FillBnd)
3003 {
3004 v_FillBwdWithBoundCond(Fwd, Bwd, PutFwdInBwdOnBCs);
3005 }
3006
3007 if (DoExchange)
3008 {
3009 // Do parallel exchange for forwards/backwards spaces.
3010 m_traceMap->GetAssemblyCommDG()->PerformExchange(Fwd, Bwd);
3011
3012 // Do exchange of interface traces (local and parallel)
3013 // We may have to split this out into separate local and
3014 // parallel for IP method???
3015 m_interfaceMap->ExchangeTrace(Fwd, Bwd);
3016 }
3017}
void v_FillBwdWithBoundCond(const Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd, bool PutFwdInBwdOnBCs) override

References Nektar::LibUtilities::eGauss_Lagrange, FilterPython_Function::field, 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_GetInterfaceMap()

InterfaceMapDGSharedPtr & Nektar::MultiRegions::DisContField::v_GetInterfaceMap ( void  )
overrideprotectedvirtual

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 2818 of file DisContField.cpp.

2819{
2820 return m_interfaceMap;
2821}

References m_interfaceMap.

◆ v_GetLeftAdjacentTraces()

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

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 2832 of file DisContField.cpp.

2833{
2834 return m_leftAdjacentTraces;
2835}

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

4440{
4441 if (locTraceFwd != NullNekDouble1DArray)
4442 {
4443 m_locTraceToTraceMap->InterpLocTracesToTraceTranspose(0, Fwd,
4444 locTraceFwd);
4445 }
4446
4447 if (locTraceBwd != NullNekDouble1DArray)
4448 {
4449 m_locTraceToTraceMap->InterpLocTracesToTraceTranspose(1, Bwd,
4450 locTraceBwd);
4451 }
4452}
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 2823 of file DisContField.cpp.

2825{
2826 return m_locTraceToTraceMap;
2827}

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

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

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

4107{
4108 int i, cnt;
4109 map<int, RobinBCInfoSharedPtr> returnval;
4110 Array<OneD, int> ElmtID, TraceID;
4111 GetBoundaryToElmtMap(ElmtID, TraceID);
4112
4113 for (cnt = i = 0; i < m_bndCondExpansions.size(); ++i)
4114 {
4116
4117 if (m_bndConditions[i]->GetBoundaryConditionType() ==
4119 {
4120 int e, elmtid;
4121 Array<OneD, NekDouble> Array_tmp;
4122
4123 locExpList = m_bndCondExpansions[i];
4124
4125 int npoints = locExpList->GetNpoints();
4126 Array<OneD, NekDouble> x0(npoints, 0.0);
4127 Array<OneD, NekDouble> x1(npoints, 0.0);
4128 Array<OneD, NekDouble> x2(npoints, 0.0);
4129 Array<OneD, NekDouble> coeffphys(npoints);
4130
4131 locExpList->GetCoords(x0, x1, x2);
4132
4133 LibUtilities::Equation coeffeqn =
4134 std::static_pointer_cast<
4135 SpatialDomains::RobinBoundaryCondition>(m_bndConditions[i])
4136 ->m_robinPrimitiveCoeff;
4137
4138 // evalaute coefficient
4139 coeffeqn.Evaluate(x0, x1, x2, 0.0, coeffphys);
4140
4141 for (e = 0; e < locExpList->GetExpSize(); ++e)
4142 {
4143 RobinBCInfoSharedPtr rInfo =
4145 TraceID[cnt + e],
4146 Array_tmp = coeffphys + locExpList->GetPhys_Offset(e));
4147
4148 elmtid = ElmtID[cnt + e];
4149 // make link list if necessary
4150 if (returnval.count(elmtid) != 0)
4151 {
4152 rInfo->next = returnval.find(elmtid)->second;
4153 }
4154 returnval[elmtid] = rInfo;
4155 }
4156 }
4157 cnt += m_bndCondExpansions[i]->GetExpSize();
4158 }
4159
4160 return returnval;
4161}
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 2803 of file DisContField.cpp.

2804{
2806 {
2807 SetUpDG();
2808 }
2809
2810 return m_trace;
2811}

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

2814{
2815 return m_traceMap;
2816}

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.

Definition at line 3434 of file DisContField.cpp.

3441{
3442 int i, n, cnt, nbndry;
3443 int nexp = GetExpSize();
3444 Array<OneD, NekDouble> f(m_ncoeffs);
3445 DNekVec F(m_ncoeffs, f, eWrapper);
3446 Array<OneD, NekDouble> e_f, e_l;
3447 //----------------------------------
3448 // Setup RHS Inner product if required
3449 //----------------------------------
3450 if (PhysSpaceForcing)
3451 {
3452 IProductWRTBase(inarray, f);
3453 Vmath::Neg(m_ncoeffs, f, 1);
3454 }
3455 else
3456 {
3457 Vmath::Smul(m_ncoeffs, -1.0, inarray, 1, f, 1);
3458 }
3459 //----------------------------------
3460 // Solve continuous Boundary System
3461 //----------------------------------
3462 int GloBndDofs = m_traceMap->GetNumGlobalBndCoeffs();
3463 int NumDirBCs = m_traceMap->GetNumLocalDirBndCoeffs();
3464 int e_ncoeffs;
3465 GlobalMatrixKey HDGLamToUKey(StdRegions::eHybridDGLamToU,
3467 const DNekScalBlkMatSharedPtr &HDGLamToU = GetBlockMatrix(HDGLamToUKey);
3468 // Retrieve number of local trace space coefficients N_{\lambda},
3469 // and set up local elemental trace solution \lambda^e.
3470 int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs();
3471 Array<OneD, NekDouble> bndrhs(LocBndCoeffs, 0.0);
3472 Array<OneD, NekDouble> loclambda(LocBndCoeffs, 0.0);
3473 DNekVec LocLambda(LocBndCoeffs, loclambda, eWrapper);
3474 //----------------------------------
3475 // Evaluate Trace Forcing
3476 // Kirby et al, 2010, P23, Step 5.
3477 //----------------------------------
3478 // Determing <u_lam,f> terms using HDGLamToU matrix
3479 for (cnt = n = 0; n < nexp; ++n)
3480 {
3481 nbndry = (*m_exp)[n]->NumDGBndryCoeffs();
3482 e_ncoeffs = (*m_exp)[n]->GetNcoeffs();
3483 e_f = f + m_coeff_offset[n];
3484 e_l = bndrhs + cnt;
3485 // use outarray as tmp space
3486 DNekVec Floc(nbndry, e_l, eWrapper);
3487 DNekVec ElmtFce(e_ncoeffs, e_f, eWrapper);
3488 Floc = Transpose(*(HDGLamToU->GetBlock(n, n))) * ElmtFce;
3489 cnt += nbndry;
3490 }
3491 Array<OneD, const int> bndCondMap =
3492 m_traceMap->GetBndCondCoeffsToLocalTraceMap();
3493 Array<OneD, const NekDouble> Sign = m_traceMap->GetLocalToGlobalBndSign();
3494 // Copy Dirichlet boundary conditions and weak forcing
3495 // into trace space
3496 int locid;
3497 cnt = 0;
3498 for (i = 0; i < m_bndCondExpansions.size(); ++i)
3499 {
3500 const Array<OneD, const NekDouble> bndcoeffs =
3501 m_bndCondExpansions[i]->GetCoeffs();
3502
3503 if (m_bndConditions[i]->GetBoundaryConditionType() ==
3505 {
3506 for (int j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); ++j)
3507 {
3508 locid = bndCondMap[cnt + j];
3509 loclambda[locid] = Sign[locid] * bndcoeffs[j];
3510 }
3511 }
3512 else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3514 m_bndConditions[i]->GetBoundaryConditionType() ==
3516 {
3517 // Add weak boundary condition to trace forcing
3518 for (int j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); ++j)
3519 {
3520 locid = bndCondMap[cnt + j];
3521 bndrhs[locid] += Sign[locid] * bndcoeffs[j];
3522 }
3523 }
3524 else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3526 {
3527 // skip increment of cnt if ePeriodic
3528 // because bndCondMap does not include ePeriodic
3529 continue;
3530 }
3531 cnt += (m_bndCondExpansions[i])->GetNcoeffs();
3532 }
3533
3534 //----------------------------------
3535 // Solve trace problem: \Lambda = K^{-1} F
3536 // K is the HybridDGHelmBndLam matrix.
3537 //----------------------------------
3538 if (GloBndDofs - NumDirBCs > 0)
3539 {
3540 GlobalLinSysKey key(StdRegions::eHybridDGHelmBndLam, m_traceMap,
3541 factors, varcoeff);
3543 LinSys->Solve(bndrhs, loclambda, m_traceMap);
3544
3545 // For consistency with previous version put global
3546 // solution into m_trace->m_coeffs
3547 m_traceMap->LocalToGlobal(loclambda, m_trace->UpdateCoeffs());
3548 }
3549
3550 //----------------------------------
3551 // Internal element solves
3552 //----------------------------------
3553 GlobalMatrixKey invHDGhelmkey(StdRegions::eInvHybridDGHelmholtz,
3555
3556 const DNekScalBlkMatSharedPtr &InvHDGHelm = GetBlockMatrix(invHDGhelmkey);
3557 DNekVec out(m_ncoeffs, outarray, eWrapper);
3558 Vmath::Zero(m_ncoeffs, outarray, 1);
3559
3560 // out = u_f + u_lam = (*InvHDGHelm)*f + (LamtoU)*Lam
3561 out = (*InvHDGHelm) * F + (*HDGLamToU) * LocLambda;
3562
3563 // Return empty GlobalLinSysKey
3564 return NullGlobalLinSysKey;
3565}
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:1650
int GetNcoeffs(void) const
Returns the total number of local degrees of freedom .
Definition: ExpList.h:1512
const DNekScalBlkMatSharedPtr & GetBlockMatrix(const GlobalMatrixKey &gkey)
Definition: ExpList.cpp:2656
static AssemblyMapSharedPtr NullAssemblyMapSharedPtr
Definition: AssemblyMap.h:51
static GlobalLinSysKey NullGlobalLinSysKey(StdRegions::eNoMatrixType)
StdRegions::ConstFactorMap factors
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:79
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
NekVector< NekDouble > DNekVec
Definition: NekTypeDefs.hpp:48
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.hpp:292
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.hpp:100

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::VarcoeffHashingTest::factors, 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 2860 of file DisContField.cpp.

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

References m_periodicBwdCopy, and m_periodicFwdCopy.

Referenced by GetFwdBwdTracePhys(), and 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 4081 of file DisContField.cpp.

4082{
4084
4085 // Reset boundary condition expansions.
4086 for (int n = 0; n < m_bndCondExpansions.size(); ++n)
4087 {
4088 m_bndCondExpansions[n]->Reset();
4089 }
4090}
virtual void v_Reset()
Reset geometry information, metrics, matrix managers and geometry information.
Definition: ExpList.cpp:3339

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

3144{
3145 m_bndCondBndWeight[index] = value;
3146}

References m_bndCondBndWeight.

◆ v_UpdateBndCondExpansion()

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

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 2849 of file DisContField.cpp.

2850{
2851 return m_bndCondExpansions[i];
2852}

References m_bndCondExpansions.

◆ v_UpdateBndConditions()

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

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 2854 of file DisContField.cpp.

2856{
2857 return m_bndConditions;
2858}

References m_bndConditions.

Member Data Documentation

◆ m_BCtoElmMap

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

Definition at line 58 of file DisContField.h.

Referenced by v_GetBoundaryToElmtMap().

◆ m_BCtoTraceMap

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

Definition at line 59 of file DisContField.h.

Referenced by v_GetBoundaryToElmtMap().

◆ m_bndCondBndWeight

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

◆ m_bndCondExpansions

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

◆ m_bndConditions

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

◆ m_boundaryTraces

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

A set storing the global IDs of any boundary Verts.

Definition at line 175 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 164 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 161 of file DisContField.h.

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

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

◆ m_periodicBwdCopy

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

Definition at line 198 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 185 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 190 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 197 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 180 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