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:3323
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 4246 of file DisContField.cpp.

4249{
4250 int i, cnt, e, ncoeff_trace;
4251 Array<OneD, NekDouble> force, out_tmp, qrhs, qrhs1;
4252 Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
4253 m_traceMap->GetElmtToTrace();
4254
4256
4257 int nq_elmt, nm_elmt;
4258 int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs();
4259 Array<OneD, NekDouble> loc_lambda(LocBndCoeffs), trace_lambda;
4260 Array<OneD, NekDouble> tmp_coeffs;
4261 m_traceMap->GlobalToLocalBnd(m_trace->GetCoeffs(), loc_lambda);
4262
4263 trace_lambda = loc_lambda;
4264
4265 int dim = (m_expType == e2D) ? 2 : 3;
4266
4267 int num_points[] = {0, 0, 0};
4268 int num_modes[] = {0, 0, 0};
4269
4270 // Calculate Q using standard DG formulation.
4271 for (i = cnt = 0; i < GetExpSize(); ++i)
4272 {
4273 nq_elmt = (*m_exp)[i]->GetTotPoints();
4274 nm_elmt = (*m_exp)[i]->GetNcoeffs();
4275 qrhs = Array<OneD, NekDouble>(nq_elmt);
4276 qrhs1 = Array<OneD, NekDouble>(nq_elmt);
4277 force = Array<OneD, NekDouble>(2 * nm_elmt);
4278 out_tmp = force + nm_elmt;
4280
4281 for (int j = 0; j < dim; ++j)
4282 {
4283 num_points[j] = (*m_exp)[i]->GetBasis(j)->GetNumPoints();
4284 num_modes[j] = (*m_exp)[i]->GetBasis(j)->GetNumModes();
4285 }
4286
4287 // Probably a better way of setting up lambda than this. Note
4288 // cannot use PutCoeffsInToElmts since lambda space is mapped
4289 // during the solve.
4290 int nTraces = (*m_exp)[i]->GetNtraces();
4291 Array<OneD, Array<OneD, NekDouble>> traceCoeffs(nTraces);
4292
4293 for (e = 0; e < (*m_exp)[i]->GetNtraces(); ++e)
4294 {
4295 edgedir = (*m_exp)[i]->GetTraceOrient(e);
4296 ncoeff_trace = elmtToTrace[i][e]->GetNcoeffs();
4297 traceCoeffs[e] = Array<OneD, NekDouble>(ncoeff_trace);
4298 Vmath::Vcopy(ncoeff_trace, trace_lambda, 1, traceCoeffs[e], 1);
4299 if (dim == 2)
4300 {
4301 elmtToTrace[i][e]->SetCoeffsToOrientation(
4302 edgedir, traceCoeffs[e], traceCoeffs[e]);
4303 }
4304 else
4305 {
4306 (*m_exp)[i]
4307 ->as<LocalRegions::Expansion3D>()
4308 ->SetFaceToGeomOrientation(e, traceCoeffs[e]);
4309 }
4310 trace_lambda = trace_lambda + ncoeff_trace;
4311 }
4312
4313 // creating orthogonal expansion (checking if we have quads or
4314 // triangles)
4315 LibUtilities::ShapeType shape = (*m_exp)[i]->DetShapeType();
4316 switch (shape)
4317 {
4319 {
4320 const LibUtilities::PointsKey PkeyQ1(
4322 const LibUtilities::PointsKey PkeyQ2(
4324 LibUtilities::BasisKey BkeyQ1(LibUtilities::eOrtho_A,
4325 num_modes[0], PkeyQ1);
4326 LibUtilities::BasisKey BkeyQ2(LibUtilities::eOrtho_A,
4327 num_modes[1], PkeyQ2);
4329 std::dynamic_pointer_cast<SpatialDomains::QuadGeom>(
4330 (*m_exp)[i]->GetGeom());
4332 BkeyQ1, BkeyQ2, qGeom);
4333 }
4334 break;
4336 {
4337 const LibUtilities::PointsKey PkeyT1(
4339 const LibUtilities::PointsKey PkeyT2(
4340 num_points[1], LibUtilities::eGaussRadauMAlpha1Beta0);
4341 LibUtilities::BasisKey BkeyT1(LibUtilities::eOrtho_A,
4342 num_modes[0], PkeyT1);
4343 LibUtilities::BasisKey BkeyT2(LibUtilities::eOrtho_B,
4344 num_modes[1], PkeyT2);
4346 std::dynamic_pointer_cast<SpatialDomains::TriGeom>(
4347 (*m_exp)[i]->GetGeom());
4349 BkeyT1, BkeyT2, tGeom);
4350 }
4351 break;
4353 {
4354 const LibUtilities::PointsKey PkeyH1(
4356 const LibUtilities::PointsKey PkeyH2(
4358 const LibUtilities::PointsKey PkeyH3(
4360 LibUtilities::BasisKey BkeyH1(LibUtilities::eOrtho_A,
4361 num_modes[0], PkeyH1);
4362 LibUtilities::BasisKey BkeyH2(LibUtilities::eOrtho_A,
4363 num_modes[1], PkeyH2);
4364 LibUtilities::BasisKey BkeyH3(LibUtilities::eOrtho_A,
4365 num_modes[2], PkeyH3);
4367 std::dynamic_pointer_cast<SpatialDomains::HexGeom>(
4368 (*m_exp)[i]->GetGeom());
4370 BkeyH1, BkeyH2, BkeyH3, hGeom);
4371 }
4372 break;
4374 {
4375 const LibUtilities::PointsKey PkeyT1(
4377 const LibUtilities::PointsKey PkeyT2(
4378 num_points[1], LibUtilities::eGaussRadauMAlpha1Beta0);
4379 const LibUtilities::PointsKey PkeyT3(
4380 num_points[2], LibUtilities::eGaussRadauMAlpha2Beta0);
4381 LibUtilities::BasisKey BkeyT1(LibUtilities::eOrtho_A,
4382 num_modes[0], PkeyT1);
4383 LibUtilities::BasisKey BkeyT2(LibUtilities::eOrtho_B,
4384 num_modes[1], PkeyT2);
4385 LibUtilities::BasisKey BkeyT3(LibUtilities::eOrtho_C,
4386 num_modes[2], PkeyT3);
4388 std::dynamic_pointer_cast<SpatialDomains::TetGeom>(
4389 (*m_exp)[i]->GetGeom());
4391 BkeyT1, BkeyT2, BkeyT3, tGeom);
4392 }
4393 break;
4394 default:
4396 "Wrong shape type, HDG postprocessing is not "
4397 "implemented");
4398 };
4399
4400 // DGDeriv
4401 // (d/dx w, d/dx q_0)
4402 (*m_exp)[i]->DGDeriv(0, tmp_coeffs = coeffs + m_coeff_offset[i],
4403 elmtToTrace[i], traceCoeffs, out_tmp);
4404 (*m_exp)[i]->BwdTrans(out_tmp, qrhs);
4405 ppExp->IProductWRTDerivBase(0, qrhs, force);
4406
4407 // + (d/dy w, d/dy q_1)
4408 (*m_exp)[i]->DGDeriv(1, tmp_coeffs = coeffs + m_coeff_offset[i],
4409 elmtToTrace[i], traceCoeffs, out_tmp);
4410
4411 (*m_exp)[i]->BwdTrans(out_tmp, qrhs);
4412 ppExp->IProductWRTDerivBase(1, qrhs, out_tmp);
4413
4414 Vmath::Vadd(nm_elmt, force, 1, out_tmp, 1, force, 1);
4415
4416 // determine force[0] = (1,u)
4417 (*m_exp)[i]->BwdTrans(tmp_coeffs = coeffs + m_coeff_offset[i], qrhs);
4418 force[0] = (*m_exp)[i]->Integral(qrhs);
4419
4420 // multiply by inverse Laplacian matrix
4421 // get matrix inverse
4422 LocalRegions::MatrixKey lapkey(StdRegions::eInvLaplacianWithUnityMean,
4423 ppExp->DetShapeType(), *ppExp);
4424 DNekScalMatSharedPtr lapsys = ppExp->GetLocMatrix(lapkey);
4425
4426 NekVector<NekDouble> in(nm_elmt, force, eWrapper);
4427 NekVector<NekDouble> out(nm_elmt);
4428
4429 out = (*lapsys) * in;
4430
4431 // Transforming back to modified basis
4432 Array<OneD, NekDouble> work(nq_elmt);
4433 ppExp->BwdTrans(out.GetPtr(), work);
4434 (*m_exp)[i]->FwdTrans(work, tmp_coeffs = outarray + m_coeff_offset[i]);
4435 }
4436}
#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 3031 of file DisContField.cpp.

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

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, ASSERTL1, Nektar::LibUtilities::Interpreter::DefineFunction(), Nektar::SpatialDomains::PointGeom::dist(), Nektar::MultiRegions::e1D, Nektar::MultiRegions::e2D, Nektar::MultiRegions::e3D, Nektar::StdRegions::eBackwards, Nektar::StdRegions::eDir1BwdDir1_Dir2BwdDir2, Nektar::StdRegions::eDir1BwdDir1_Dir2FwdDir2, Nektar::StdRegions::eDir1BwdDir2_Dir2BwdDir1, Nektar::StdRegions::eDir1BwdDir2_Dir2FwdDir1, Nektar::StdRegions::eDir1FwdDir1_Dir2BwdDir2, Nektar::StdRegions::eDir1FwdDir1_Dir2FwdDir2, Nektar::StdRegions::eDir1FwdDir2_Dir2BwdDir1, Nektar::StdRegions::eDir1FwdDir2_Dir2FwdDir1, Nektar::ErrorUtil::efatal, Nektar::StdRegions::eForwards, Nektar::StdRegions::eNoOrientation, Nektar::SpatialDomains::ePeriodic, Nektar::LibUtilities::Interpreter::Evaluate(), Nektar::MultiRegions::ExpList::GetBoundaryCondition(), Nektar::SpatialDomains::BoundaryConditions::GetBoundaryConditions(), Nektar::SpatialDomains::BoundaryConditions::GetBoundaryRegions(), Nektar::SpatialDomains::QuadGeom::GetFaceOrientation(), Nektar::SpatialDomains::TriGeom::GetFaceOrientation(), Nektar::MultiRegions::RotPeriodicInfo::m_angle, Nektar::MultiRegions::ExpList::m_comm, Nektar::MultiRegions::RotPeriodicInfo::m_dir, Nektar::MultiRegions::ExpList::m_expType, Nektar::MultiRegions::ExpList::m_graph, m_periodicEdges, m_periodicFaces, m_periodicVerts, Nektar::MultiRegions::RotPeriodicInfo::m_tol, NEKERROR, CellMLToNektar.cellml_metadata::p, Nektar::LibUtilities::ReduceMax, Nektar::LibUtilities::ReduceSum, Nektar::SpatialDomains::PointGeom::Rotate(), sign, tinysimd::sqrt(), 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:208
std::shared_ptr< BoundaryConditionMap > BoundaryConditionMapShPtr
Definition: Conditions.h:218

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

Referenced by DisContField().

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

2919{
2920 ASSERTL1((*m_exp)[0]->GetBasis(0)->GetBasisType() !=
2922 "Routine not set up for Gauss Lagrange points distribution");
2923
2924 // blocked routine
2925 Array<OneD, NekDouble> tracevals(m_locTraceToTraceMap->GetNLocTracePts());
2926
2927 m_locTraceToTraceMap->LocTracesFromField(field, tracevals);
2928 m_locTraceToTraceMap->InterpLocTracesToTrace(0, tracevals, Fwd);
2929
2930 Array<OneD, NekDouble> invals =
2931 tracevals + m_locTraceToTraceMap->GetNFwdLocTracePts();
2932 m_locTraceToTraceMap->InterpLocTracesToTrace(1, invals, Bwd);
2933
2935
2936 FillBwdWithBoundCond(Fwd, Bwd, BndCond, BndCondExp, false);
2937
2938 // Do parallel exchange for forwards/backwards spaces.
2939 m_traceMap->GetAssemblyCommDG()->PerformExchange(Fwd, Bwd);
2940
2941 // Do exchange of interface traces (local and parallel)
2942 // We may have to split this out into separate local and
2943 // parallel for IP method???
2944 m_interfaceMap->ExchangeTrace(Fwd, Bwd);
2945}
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 2706 of file DisContField.cpp.

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

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

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

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

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

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

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

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

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

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

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

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

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

3178{
3179 ASSERTL1(m_physState == true, "local physical space is not true ");
3180 v_ExtractTracePhys(m_phys, outarray);
3181}
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 3196 of file DisContField.cpp.

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

3024{
3026 PutFwdInBwdOnBCs);
3027}

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

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

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

3141{
3142 return m_bndCondBndWeight;
3143}

References m_bndCondBndWeight.

◆ v_GetBndCondExpansions()

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

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 2839 of file DisContField.cpp.

2841{
2842 return m_bndCondExpansions;
2843}

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

2847{
2848 return m_bndConditions;
2849}

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

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

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

2873{
2874 v_GetFwdBwdTracePhys(m_phys, Fwd, Bwd);
2875}
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 2947 of file DisContField.cpp.

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

2821{
2822 return m_interfaceMap;
2823}

References m_interfaceMap.

◆ v_GetLeftAdjacentTraces()

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

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 2834 of file DisContField.cpp.

2835{
2836 return m_leftAdjacentTraces;
2837}

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

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

2827{
2828 return m_locTraceToTraceMap;
2829}

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

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

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

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

2806{
2808 {
2809 SetUpDG();
2810 }
2811
2812 return m_trace;
2813}

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

2816{
2817 return m_traceMap;
2818}

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

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

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

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

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

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

3146{
3147 m_bndCondBndWeight[index] = value;
3148}

References m_bndCondBndWeight.

◆ v_UpdateBndCondExpansion()

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

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 2851 of file DisContField.cpp.

2852{
2853 return m_bndCondExpansions[i];
2854}

References m_bndCondExpansions.

◆ v_UpdateBndConditions()

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

Reimplemented from Nektar::MultiRegions::ExpList.

Definition at line 2856 of file DisContField.cpp.

2858{
2859 return m_bndConditions;
2860}

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